#  Impulse response estimation via flexible local projections
\-
Latest draft available at this link††thanks: We are thankful to Vito Polito
for helpful comments and suggestions.

Haroon Mumtaz Queen Mary, University of London, School of Economics and
Finance, Mile End road, London E1 4LJ, UK. e-mail<EMAIL_ADDRESS>Michele
Piffer King’s Business School, King’s College London, Bush House, 30 Aldwych,
London WC2B 4BG, UK. e-mail<EMAIL_ADDRESS>
###### Abstract

This paper introduces a flexible local projection that generalises the model
by Jordà (2005) to a non-parametric setting using Bayesian Additive Regression
Trees. Monte Carlo experiments show that our BART-LP model is able to capture
non-linearities in the impulse responses. Our first application shows that the
fiscal multiplier is stronger in recession than expansion only in response to
contractionary fiscal shocks, but not in response to expansionary fiscal
shocks. We then show that financial shocks generate effects on the economy
that increase more than proportionately in the size of the shock when the
shock is negative, but not when the shock is positive.

JEL classification: C14, C11, C32, E52.
Keywords: Non-linear models, non-parametric techniques, identification.

##  1 Introduction

Estimation of impulse responses (IRFs) via local projections (LP) by Jordà
(2005) has become increasingly common in applied Macroeconometric analysis. A
key feature of the local projection estimator is that it estimates IRFs of
variable $y_{t}$ to an innovation to variable $x_{t}$ directly via linear
regressions of the form $y_{t+h}=\beta_{h}x_{t}+d_{h}w_{t}+u_{t+h}$, where
$w_{t}$ denotes control variables. Given their flexibility, considerable
attention has been given to investigate the properties of the LP estimator,
see, for instance, Stock & Watson (2018) and Plagborg-Møller & Wolf (2021).

In their most popular specification, LP estimators impose a linearity between
$y_{t+h}$ and $(x_{t},w_{t})$. This limitation implies that linear LPs cannot
be used to study non-linear effects of the shocks of interest, for instance
non-linearities on the sign or size of the shock, or on the economic
conditions when the shock occurs. Some extensions of the linear LP estimator
have been proposed, but they all rely on the functional form introduced to
model the non-linearity. Jordà (2005) proposes the use of quadratic and cubic
terms. Auerbach & Gorodnichenko (2013a, b) and Ramey & Zubairy (2018) use a
smooth transition function and a threshold function, respectively. Ruisi
(2019) and Lusompa (2021) use a time varying extension of LP based on
parametric state-space models, while Inoue et al. (2022) provides a more
general framework for modelling structural shifts.

In this paper we propose a flexible non-linear extension of the LP estimator
that does not require assumptions on the functional form of the LP regression
equation. We propose a non-parametric LP estimator that uses the Bayesian
Additive Regression Trees (BART) model to approximate the unknown function
$m_{h}(z_{t})$ in the more general equation $y_{t+h}=m_{h}(z_{t})+u_{t+h}$,
with $z_{t}=(x_{t},w_{t})$. Introduced by Chipman et al. (2010), BART uses
regression trees as its building block. Regression trees split the space of
explanatory variables $z_{t}$ into sub-groups based on rules of the form
$z_{t,j}<C$, where $j=1,2,..,K$. The function $m_{h}(z_{t})$ is approximated
as a sum of a large number of small trees. Chipman et al. (2010) show that
BART is able to approximate highly non-linear functions accurately.

We first illustrate how BART techniques can be applied in a non-linear LP
estimator, and we refer to this new methodology as BART-LP. We show that BART-
LP can handle autocorrelation in the error terms, a problem already discussed
in the literature of linear LP estimators. We then document the performance of
BART-LP using Monte Carlo analysis. We build our simulations on three models.
First, we use the same SVAR-GARCH model employed in the simulation by Jordà
(2005), where a structural shock generates non-linear effects that affect the
variance of the shock. Second, we use a recursive Threshold VAR model in which
shocks generate stronger effects in one of the two regimes of the model.
Third, we use a sign-dependent moving average model in which the true monetary
policy shock generates different effects depending on the sign of the shock.
In all cases, the BART-LP is capable of recovering the true impulse responses,
while the linear LP typically estimates a weighted average of the true
underlying non-linear impulse responses. We focus the discussion on the IRF
estimation, which is separate from identification of the structural shocks. As
also in Jordà (2005), we do not investigate the topic of identification in our
framework.

We then apply the new methodology to shed some light on two separate ongoing
debates in the literature. We first apply our methodology to US fiscal
spending shocks. While Auerbach & Gorodnichenko (2013b) has argued that US
fiscal multipliers are stronger in recession than in boom, their evidence has
been called into question by Ramey & Zubairy (2018). We use our BART-LP
procedure to show that the answer to this question depends on the sign of the
shock. We show that the multiplier in response to a positive shock does not
change significantly depending on whether the shock hits the economy in an
expansion or in a recession. However, in response to a negative shock the
multiplier is stronger in a recession. Compared to the models used by Auerbach
& Gorodnichenko (2013b) and Ramey & Zubairy (2018), our methodology can detect
both non-linearities over the state of the economy and non-linearities over
the sign and size of the shock. Our results hence help reconcile the
difference found by Auerbach & Gorodnichenko (2013b) and Ramey & Zubairy
(2018), who used models that are more suitable to detect non-linearities over
the state of the economy, but not over the sign of the shocks. Last, we
revisit the analysis by Forni et al. (2021) on financial shocks. We confirm
their result that negative financial shocks generate detrimental effects on
output that increase more than proportionally in the size of the shock. These
effects cannot otherwise be detected by a linear model, which is doomed to
underestimate the effect of strong negative financial shocks.

This paper relates to the literature that studies how BART techniques can be
used in Macroeconometrics. Huber & Rossini (2021) introduce a VAR model where
the dynamics of the endogenous variables are modelled using BART. The authors
model the impact of uncertainty shocks using their proposed model. Huber et
al. (2020) extend the BART-VAR to a mixed frequency setting and evaluate the
forecasting performance of the model. Clark et al. (2021) show that
multivariate BART regression models perform well in terms of tail forecasting.
To the best of our knowledge, our paper is the first one to use BART in an LP
framework. The paper is also part of a broad literature that studies the
advantages of IRF estimation using LP estimators, relative to constructing
IRFs on vector autoregressive models. Several contributions document the
performance of LP estimators, including Gonçalves et al. (2022), Kilian & Kim
(2011), Alloza et al. (2019), Breitung et al. (2019) Herbst & Johannsen
(2021), and Bruns & Lütkepohl (2022). While LP estimators are usually proposed
in a frequentist setting, we follow Miranda-Agrippino & Ricco (2021) and take
a Bayesian approach to LP, yet in a non-linear framework.

The paper is organised as follows. Section 2 presents the empirical model.
Section 3 reports the results from the simulation exercise. Section 4 shows
the application to fiscal and financial shocks. Section 5 concludes.

##  2 Flexible local projections

In this section we outline the methodology, which we refer to as the BART-LP
model, or flexible local projections. We discuss how BART-LP approximates the
unknown conditional expectation function of local projection models, discuss
the prior, and outline the posterior sampler. We then discuss how to compute
generalized impulse responses to structural shocks within our framework.

###  2.1 The BART approximation

We work with the equation

$y_{t+h}=m_{h}\big{(}x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}\big{)}+\epsilon_{t+h}^{(h)},$
(1)

where $y_{t+h}$ denotes the scalar variable of interest, $h=0,1,...,H$ is the
impulse response horizon, and $\epsilon_{t+h}^{(h)}$ satisfies
$E\big{(}\epsilon_{t+h}^{(h)}|x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}\big{)}=0$.
We aim to study how $y_{t+h}$ responds to a change in the scalar variable
$x_{t}$. The vector $\boldsymbol{z}_{t}$ contains observable control
variables, possibly including lagged values of $y_{t}$ and contemporaneous
and/or lagged values of other variables. The vector
$\boldsymbol{w}_{t+h}^{(h)}$ contains additional control variables in the form
of estimated residuals, as explained in Section 2.3. While
$(x_{t},\boldsymbol{z}_{t})$ are the same for every regression model $h$,
$\boldsymbol{w}_{t+h}^{(h)}$ can potentially change. The function $m_{h}(.)$
captures the true unknown conditional expectation function. The residual
$\epsilon^{(h)}_{t+h}$ is assumed to be normally distributed with variance
$\sigma_{t+h}^{2~{}(h)}$. As noted in Jordà (2005), $\epsilon^{(h)}_{t+h}$ is
serially correlated for $h\geq 1$.

It is common in the literature to assume a functional form for
$m_{h}\big{(}x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}\big{)}$. The
most popular applications of LPs use a function of the type

$\displaystyle y_{t+h}$
$\displaystyle=~{}~{}~{}~{}~{}~{}~{}g(q_{t})\big{[}\alpha^{(h)}_{0}x_{t}+\boldsymbol{\alpha}_{1}^{\prime(h)}\boldsymbol{z}_{t}+\boldsymbol{\alpha}_{2}^{\prime(h)}\boldsymbol{w}_{t+h}^{(h)}\big{]}+$
(2)
$\displaystyle+\big{(}1-g(q_{t})\big{)}\big{[}\beta^{(h)}_{0}x_{t}+\boldsymbol{\beta}_{1}^{\prime(h)}\boldsymbol{z}_{t}+\boldsymbol{\beta}_{2}^{\prime(h)}\boldsymbol{w}_{t+h}^{(h)}\big{]}+\epsilon_{t+h}^{(h)}.$
(3)

The special case of a linear model sets $g(q_{t})=1$, $\forall~{}t$, and
estimates the impulse response to a shock to $x_{t}$ using the estimates for
$\\{\alpha^{(h)}_{0}\\}_{h=0}^{H}$ (Jordà, 2005). Non-linear applications
usually specify a transition variable $q_{t}$, assume a specific non-linear
functional form for $g(.)$, and compute non-linear impulse responses as a
function of the estimates for
$\\{\alpha^{(h)}_{0},\beta^{(h)}_{0}\\}_{h=0}^{H}$. For example, the smooth
transition formulation by Auerbach & Gorodnichenko (2013a, b) sets $g(.)$
equal to the logistic function, while the threshold formulation by Ramey &
Zubairy (2018) and Alpanda et al. (2021) sets $g(.)$ equal to the indicator
function.

Our paper differs from the existing literature by approximating the unknown
conditional expectation function
$m_{h}\big{(}x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}\big{)}$ using
Bayesian Additive Regression Trees (BART). It is assumed that

$m_{h}(x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)})\approx
f_{h}(x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}|\Gamma^{(h)},\boldsymbol{\mu}^{(h)})=\sum\limits_{j=1}^{J}f_{h,j}(x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}|\Gamma_{j}^{(h)},\boldsymbol{\mu}_{j}^{(h)}),$
(4)

where
$f_{h,j}(x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}|\Gamma^{(h)}_{j},\boldsymbol{\mu}_{j}^{(h)})$
denotes a single regression tree $j$ at horizon $h$ and
$f_{h}(x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}|\Gamma^{(h)},\boldsymbol{\mu}^{(h)})$
denotes the sum of $J$ regression trees. The parameters of the regression
trees are the tree structures
$\Gamma^{(h)}=[\Gamma^{(h)}_{1},..,\Gamma^{(h)}_{J}]$ and the terminal nodes
(or leaves)
$\boldsymbol{\mu}^{(h)}=(\boldsymbol{\mu}_{1}^{(h)^{\prime}},..,\boldsymbol{\mu}_{J}^{(h)^{\prime}})^{\prime}$,
with $\boldsymbol{\mu}_{j}^{(h)}$ of dimension $B_{j}\times 1$ and $B_{j}$ the
number of terminal nodes of tree $j$. As an illustration, each regression tree
divides the space of each explanatory variable by using binary splitting
rules. Denoting $X_{i}$ as the i-th entry of the vector
$\big{(}x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}\big{)}$, these
rules are defined as:

$\displaystyle X_{i}$ $\displaystyle\leq c,$ (5) $\displaystyle X_{i}$
$\displaystyle>c,$

with $c$ the threshold value. Observations are assigned according to these
splitting rules, and the terminal nodes return the fitted value conditional on
the split. The fitted value of the dependent variable, based on a single
regression tree, is then given by

$f_{h,j}(x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}|\Gamma_{j}^{(h)},\boldsymbol{\mu}_{j}^{(h)})=\sum_{b=1}^{B_{j}}I(X_{i})\mu_{j,b}^{(h)},$
(6)

where $I(.)$ denotes an indicator function that equals 1 if $X_{i}$ belongs to
the set defined by the splitting rule implicit in $\Gamma_{j}^{(h)}$. Note
that the complexity of each tree is determined by $B_{j}$, the number of
terminal nodes. We refer the reader to the Online Appendix for an illustrative
example of a regression tree, and to Hill et al. (2020) for a comprehensive
review.

The model in equation (4) approximates
$m_{h}\big{(}x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}\big{)}$ using
a sum $J$ trees. Each tree in the sum is restricted to be small a priori to
avoid overfitting, and thus explains a small proportion of $y_{t+h}$ and is a
‘weak learner’. Chipman et al. (2010) show that a low value of $J$ reduces
predictive accuracy. As $J$ increases, predictive performance initially
improves, but this improvement tapers off, eventually. In practice, studies
such as Huber et al. (2020) note that the difference in predictive accuracy is
negligible for $J>150$ and complex functions can be easily approximated using
$J=200$ or $250$.

The BART approximation of the relationship between $y_{t+h}$ and $x_{t}$ has
implications for the properties of the impulse responses. As the regression
trees split the space of the covariates via rules of the type shown in
equation (5), the estimated predictions on the right-hand side of equation
(15) are dependent on their history. Similarly, the shock $d$ to variable
$x_{t}$ can lead to predictions that proportionally differ if the size and
sign of the shock leads to the covariate space where the relationship between
$y_{t+h}$ and $x_{t}$ is substantially different from the ‘average’ impact. We
now discuss estimation, and then return to a detailed discussion of non-linear
impulse responses in the BART-LP model in Section 2.4.

###  2.2 Estimation

The model in equation (1) can be estimated using the MCMC algorithm described
in Chipman et al. (2010), which we summarize here for completeness.

####  2.2.1 Priors

The prior distributions proposed by Chipman et al. (2010) play a crucial role,
as they are devised to reduce the possibility of overfitting. The joint prior
for the parameters of the $J$ trees of the model at each horizon $h$ is
factored as follows:

$p\big{(}(\Gamma_{1}^{(h)},\boldsymbol{\mu}_{1}^{(h)}),(\Gamma_{2}^{(h)},\boldsymbol{\mu}_{2}^{(h)}),...,(\Gamma_{J}^{(h)},\boldsymbol{\mu}_{J}^{(h)})\big{)}=\prod_{j=1}^{J}p(\boldsymbol{\mu}_{j}^{(h)}|\Gamma_{j}^{(h)})p(\Gamma^{(h)}_{j}),$
(7)

where
$p(\boldsymbol{\mu}^{(h)}_{j}|\Gamma^{(h)}_{j})=\prod_{b=1}^{B_{j}}p(\mu_{b,j}^{(h)}|\Gamma_{j}^{(h)})$.

The prior for the tree structure $\Gamma_{j}^{(h)}$ depends on the probability
that the node at depth $d=0,1,2,..$ is not a terminal node. This prior
probability is given by $\alpha(1+d)^{-\beta}$ where $\alpha\in(0,1)$ and
$\beta>0$. Higher values of $\beta$ and smaller values of $\alpha$ reduce this
probability and impose a stronger belief that the tree has a simple (i.e.
shorter) structure. We follow the recommendation by Chipman et al. (2010) and
set $\alpha=0.95$ and $\beta=2$. The prior for the threshold value $c$ implies
that this parameter is assumed to be uniform over the range of the values
taken by the variables. In the default setting, the choice of splitting
variable is also assumed to be uniform across the regressors.

To define $p(\boldsymbol{\mu}^{(h)}_{j}|\Gamma^{(h)}_{j})$ Chipman et al.
(2010) first transform the dependent variable so that it lies between $-0.5$
and $0.5$. As a consequence,
$m_{h}\big{(}x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}\big{)}$ is
also expected to lie between these values. The prior
$p(\boldsymbol{\mu}^{(h)}_{j}|\Gamma^{(h)}_{j})$ is assumed to be normal
$N(0,S)$. The variance $S$ is set as $\frac{1}{2\kappa(J^{0.5})}$, with
$\kappa$ set to 2, the value recommended by Chipman et al. (2010). Under this
default prior, there is a 95% probability that the conditional mean of the
dependent variable lies between $-0.5$ and $0.5$.

A conjugate inverse $\chi^{2}$ prior is used for the variance
$\sigma_{t+h}^{2~{}(h)}$. The hyperparameters of the prior distribution are
set by using an estimate $\hat{\sigma}_{t+s}^{2~{}(h)}$ of the variance
obtained from a linear regression. If the true model is non-linear
$\hat{\sigma}_{t+s}^{2~{}(h)}$ will be biased upwards. Under the default
prior, the hyperparameters are chosen so that
$\Pr(\sigma_{t+s}^{(h)}<\hat{\sigma}_{t+s}^{(h)})=0.9$.

The total number of trees $J$ is fixed.

####  2.2.2 MCMC algorithm

The MCMC algorithm devised by Chipman et al. (2010) samples from the
conditional posterior distributions of $\sigma_{t+h}^{2~{}(h)}$ and the
parameters of the regression trees in each iteration.111Intuitive descriptions
of this MCMC algorithm can be found in Clark et al. (2021) and Hill et al.
(2020). Each iteration of the algorithm samples from the following conditional
posteriors:

  1. 1.

conditional on the trees, the error variance can be easily drawn from the
inverse Gamma distribution;

  2. 2.

the conditional posterior distribution of the tree structure is not known in
closed form and a Metropolis-Hastings algorithm is used. Define $R^{(h)}_{j}$
as the residual:

$R^{(h)}_{j}=y_{t+h}-\sum_{i\neq
j}f(x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}|\Gamma^{(h)}_{j},\boldsymbol{\mu}^{(h)}_{j}).$
(8)

The $j-th$ tree is proposed using the density
$q(\Gamma_{j}^{new},\Gamma_{j}^{old})$. Chipman et al. (2010) use a proposal
density that incorporates 4 moves: (i) splitting the node into two new nodes
(grow), (ii) transforming adjacent nodes to terminal node (prune), (iii)
changing the decision rule of an interior node (change), (iv) swapping a
decision rule between a node that is above and the node before it (swap). The
probabilities associated with these moves are fixed at $0.25,0.25,0.4$ and
$0.1$ respectively. The proposed tree structure $\Gamma_{j}^{new}$ is accepted
with probability

$\alpha=\frac{q(\Gamma_{j}^{new},\Gamma_{j}^{old})p(R_{j}^{(h)}|\Gamma_{j}^{new},\sigma_{t+h}^{2~{}h})p(\Gamma_{j}^{new})}{q(\Gamma_{j}^{old},\Gamma_{j}^{new})p(R_{j}^{(h)}|\Gamma_{j}^{old},\sigma_{t+h}^{2})p(\Gamma_{j}^{old})},$
(9)

where $p(R_{j}|\Gamma_{j},\sigma_{t+h}^{2~{}h})$ is the conditional likelihood
and $p(\Gamma_{j})$ denotes the prior. This step is repeated for $j=1,2,..,J$
trees;

  3. 3.

the conditional posterior distribution of the terminal node parameters is
Gaussian with the parameters known in closed form. Therefore, the draw of
$\boldsymbol{\mu}^{(h)}_{j}$ for $j=1,2,..,J$ can be carried out in a
straightforward manner;

  4. 4.

given a draw of the model parameters conditioning on
$(x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)})$, the predicted value
can be computed as

$E\left(y_{t+h}|x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}\right)=\sum_{j=1}^{J}f_{h}(x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}|\Gamma^{(h)}_{j},\boldsymbol{\mu}^{(h)}_{j})=\sum_{j=1}^{J}\sum_{b=1}^{B_{j}}I(x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)})\mu_{j,b}^{(h)},$
(10)

with $I(x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)})$ an indicator
function equal to 1 if $(x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)})$
belongs to the splitting rule implied by
$[\Gamma_{1}^{(h)},...,\Gamma_{J}^{(h)}]$.

###  2.3 Autocorrelation

The residual term in LP models is known to be autocorrelated, a feature that
must be taken into account in the estimation. In the case of linear local
projections, it has been shown that the residual at horizon $h$ follows a
$MA(h-1)$ process, see for instance Lusompa (2021). Lusompa (2021) suggests a
GLS procedure whereby the autocorrelation is eliminated by including leads of
the LP residuals from horizon $h=0$ in the conditioning set.222Lusompa (2021)
suggests an efficient strategy that transforms the dependent variable of the
LP regressions and does not require one to explicitly include the horizon $0$
residuals as regressors.

The non-parametric setting considered in this paper encompasses non-linear
models. For the purpose of illustration, consider a simple non-parametric
AR(1) model

$y_{t+1}=v(y_{t};A_{1})+e_{t+1}.$ (11)

Iterating the process forward 3 periods as an example gives

$\displaystyle y_{t+2}$
$\displaystyle=v\big{(}v(y_{t};A_{1})+e_{t+1};A_{1}\big{)}+e_{t+2},$
$\displaystyle y_{t+3}$
$\displaystyle=v\big{(}v(v(y_{t};A_{1})+e_{t+1};A_{1}\big{)}+e_{t+2};A_{1})+e_{t+3}.$

It is useful to compare this with a BART-LP for this horizon:

$y_{t+3}=f_{3}(y_{t}|\Gamma_{3},\boldsymbol{\mu}_{3})+\epsilon_{t+3}^{(3)}.$
(12)

The function $f_{3}(y_{t}|\Gamma_{3},\boldsymbol{\mu}_{3})$ approximates the
non-linear relationship between $y_{t}$ and its lead, but does not account for
the dependence between the dependent variable and lagged shocks. Thus in this
setting, the residual $\epsilon_{t+3}^{(3)}$ is a non-linear function of
$e_{t+1}$ and $e_{t+2}$, and has a non-linear autocorrelation structure.

In general, the Voltera expansion of any non-linear time-series shows its
complex dependence on past shocks:

$y_{t}=\sum\limits_{i=0}^{\infty}\phi^{i}e_{t-i}+\sum\limits_{i=0}^{\infty}\sum\limits_{i=j}^{\infty}\zeta_{ij}e_{t-i}e_{t-j}+\sum\limits_{i=0}^{\infty}\sum\limits_{i=0}^{\infty}\sum\limits_{k=0}^{\infty}\phi_{ij}e_{t-i}e_{t-j}e_{t-k}+...$
(13)

To account for this autocorrelation, we propose to include an estimate of the
shocks,
$\boldsymbol{w}_{t+h}^{(h)}=(\hat{e}_{t+1},\hat{e}_{t+2,}...,\hat{e}_{t+h-1})$,
as additional covariates in the $h-$period BART-LP. Following Lusompa (2021)
we construct $\boldsymbol{w}_{t+h}^{(h)}$ at every horizon $h\geq 1$ using the
residuals of the period $h=0$ flexible local projection
$y_{t}=f_{0}(x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t}^{(0)}|\Gamma^{(0)},\boldsymbol{\mu}^{(0)})+\epsilon_{t}^{(0)}$,
with $\boldsymbol{w}_{t}^{(0)}=\boldsymbol{0}$. Then, the flexible local
projection for period $h$ is specified as

$y_{t+h}=f_{h}(x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}|\Gamma^{(h)},\boldsymbol{\mu}^{(h)})+\epsilon^{(h)}_{t+h}.$
(14)

The BART approximation of the true non-linear function
$m_{h}(x_{t},\boldsymbol{z}_{t},\hat{w}_{t+h}^{(h)})$ proxies the non-linear
dependence of $y_{t+h}$ on $e_{t+1},e_{t+2,}...,e_{t+h-1}$ and ameliorates the
autocorrelation in $\epsilon_{t+h}^{(h)}$.

###  2.4 Generalized structural impulse responses

The computation of impulse responses to structural shocks typically relies on
two conceptually different pillars: an estimation procedure for the impulse
responses, and an identification scheme for the structural shock of interest.
Our paper only aims to advance the literature on the former, and does not
investigate the topic of identification in a non-parametric framework.

Define $\boldsymbol{y}_{t}$ a $k\times 1$ vector of variables of interest,
with $x_{t}$ one of the entries of $\boldsymbol{y}_{t}$. Define $\phi_{i,h}$
the impulse response of variable $y_{i,t}$ to a structural shock to variable
$x_{t}$ of size $d$, with $i=1,..,k$. Following Koop et al. (1996),
$\phi_{i,h}$ is given by

$\phi_{i,h}=E(y_{i,t+h}|x_{t}=d;\boldsymbol{n}_{t})-E(y_{i,t+h}|x_{t},\boldsymbol{n}_{t}).$
(15)

We compute generalized impulse responses numerically using the following
algorithm:

  1. 1.

estimate the model for $h=0$,

$y_{i,t}=f_{0}\big{(}x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t}^{(0)}|\Gamma^{(0)},\boldsymbol{\mu}^{(0)}\big{)}+\epsilon_{i,t}^{(0)},$
(16)

with $\boldsymbol{w}_{t}^{(0)}=\boldsymbol{0}$, and store $D$ vectors of
dimension $T\times 1$ containing the estimated residuals associated with $D$
posterior draws, $\\{\hat{\epsilon}_{i,t}^{(0),d}\\}_{t=1}^{T}$, $d=1,..,D$;

  2. 2.

for the generic draw $d$ of the residuals
$\\{\hat{\epsilon}_{i,t}^{(0),d}\\}_{t=1}^{T}$ and for the generic horizon
$h$, estimate the model

$y_{i,t+h}=f_{h}\big{(}x_{t},\boldsymbol{z}_{t},\boldsymbol{w}_{t+h}^{(h)}|\Gamma^{(h)},\boldsymbol{\mu}^{(h)}\big{)}+\epsilon_{i,t+h}^{(h)},$
(17)

with

$\boldsymbol{w}_{t+h}^{(h)}=(\hat{\epsilon}_{i,t}^{(0),d},\hat{\epsilon}_{i,t-1}^{(0),d},...,\hat{\epsilon}_{i,t-h+1}^{(0),d});$
(18)

  3. 3.

compute the predicted values $(\hat{y}_{i,h}^{0,d},~{}\hat{y}_{i,h}^{1,d})$
associated with one posterior draw from model (17) conditioning on

$\displaystyle\big{(}x_{t}=\bar{x},~{}~{}~{}~{}~{}~{}$
$\displaystyle\boldsymbol{z}_{t}=\bar{\boldsymbol{z}},~{}\boldsymbol{w}_{t+h}^{(h)}=\bar{\boldsymbol{w}}^{(h)}\big{)},$
(19) $\displaystyle\big{(}x_{t}=\bar{x}+d,~{}$
$\displaystyle\boldsymbol{z}_{t}=\bar{\boldsymbol{z}},~{}\boldsymbol{w}_{t+h}^{(h)}=\bar{\boldsymbol{w}}^{(h)}\big{)},$
(20)

respectively, with
$\big{(}\bar{x},\bar{\boldsymbol{z}},\bar{\boldsymbol{w}}^{(h)}\big{)}$
defined below. Compute
$\psi_{i,h}^{d}=\hat{y}_{i,h}^{1,d}-\hat{y}_{i,h}^{0,d}$ and store

$\big{(}\hat{y}_{i,h}^{0,d},~{}~{}\hat{y}_{i,h}^{1,d},~{}~{}\psi_{i,h}^{d}\big{)};$
(21)

  4. 4.

repeat steps 2-3 for $h=0,..,H$, store
$\\{\big{(}\hat{y}_{i,h}^{0,d},~{}~{}\hat{y}_{i,h}^{1,d},~{}~{}\psi_{i,h}^{d}\big{)}\\}_{h=0}^{H}$
given the same posterior draw $d$ for
$\\{\hat{\epsilon}_{i,t}^{(0),d}\\}_{t=1}^{T}$;

  5. 5.

repeat steps 2-4 for $d=1,..,D$;

  6. 6.

compute the average across posterior draws

$\psi_{i,h}=\frac{1}{D}\sum_{d=1}^{D}\hat{y}_{i,h}^{1,d}-\frac{1}{D}\sum_{d=1}^{D}\hat{y}_{i,h}^{0,d}.$
(22)

The exact implementation of the above procedure depends on how the structural
shocks are identified. Different options are in principle available. One
option is to follow Plagborg-Møller & Wolf (2021) and Barnichon & Brownlees
(2019) and replicate a recursive identification scheme by using appropriate
control variables. For example, if the aim is to estimate the response of GDP
to an interest rate shock that is restricted to have a zero contemporaneous
impact on GDP and CPI in a trivariate model, one can set $x_{t}$ equal to the
policy interest rate and add contemporaneous GDP and CPI into
$\boldsymbol{z}_{t}$, which will also include $L$ lags of the variables.
Alternatively, one can set $x_{t}$ equal to either a proxy for the structural
shock of interest, or equal to the true realizations of the shocks, if
available, and set
$\boldsymbol{z}_{t}=(\boldsymbol{y}_{t-1}^{\prime},..,\boldsymbol{y}_{t-L}^{\prime})^{\prime}$.
Either way, the generalized impulse responses are computed as
$\phi_{i,h}=\psi_{i,h}$, with $\psi_{i,h}$ from equation (22).

An alternative approach is to follow Jordà (2005) more closely and use a
separate impulse vector. One can estimate a preliminary SVAR on
$\boldsymbol{y}_{t}$ and use a preferred identification approach to estimate
the impulse vector $\boldsymbol{d}$ of variables $\boldsymbol{y}_{t}$ to a
shock to variable $x_{t}$ of size $d$. The above algorithm can then be run by
setting $x_{t}=0$ and
$\boldsymbol{z}_{t}=(\boldsymbol{y}_{t-1}^{\prime},..,\boldsymbol{y}_{t-L}^{\prime})^{\prime}$.
Step 1 of the algorithm is still required to generate
$\\{\hat{\epsilon}_{i,t}^{(0),d}\\}_{t=1}^{T}$, $d=1,..,D$. Steps 2 to 4 are
then run for $h=1,2,..,H$, replacing equations (19)-(20) with

$\displaystyle\big{(}x_{t}=0,~{}$
$\displaystyle\boldsymbol{z}_{t}=\bar{\boldsymbol{z}},~{}~{}~{}~{}~{}~{}\boldsymbol{w}_{t+h}^{(h)}=\bar{\boldsymbol{w}}^{(h)}\big{)},$
(23) $\displaystyle\big{(}x_{t}=0,~{}$
$\displaystyle\boldsymbol{z}_{t}=\bar{\boldsymbol{z}}+\boldsymbol{d},~{}\boldsymbol{w}_{t+h}^{(h)}=\bar{\boldsymbol{w}}^{(h)}\big{)}.$
(24)

The generalized impulse responses are then computed as $\phi_{i,h}=d_{i}$ for
$h=0$ with $d_{i}$ the $i-th$ entry of $\boldsymbol{d}$, and as
$\phi_{i,h}=\psi_{i,h}$ from (22) for $h=1,..,H$.

Last, to be operational, the above procedure requires specifying the
conditioning values
$\big{(}\bar{x},\bar{\boldsymbol{z}},\bar{\boldsymbol{w}}^{(h)}\big{)}$ as
well as the shock $d$. The former can be set differently depending on whether
the shock is simulated to hit the economy at any point in time or on a subset
of periods. Our baseline specification sets
$\big{(}\bar{x},\bar{\boldsymbol{z}},\bar{\boldsymbol{w}}^{(h)}\big{)}$ equal
to the sample average, which is computed within each regression model $h$.
$d$, instead, can be set in accordance to the sign and size of the intended
shock.

##  3 Monte Carlo simulation

We use Monte Carlo simulations to assess if BART-LP can recover the true non-
linear patters of a data generating process. We use three different data
generating processes: a SVAR-GARCH model, a Threshold VAR model, and a sign-
dependent Moving Average model.

###  3.1 SVAR-GARCH model

The first model we use for simulations is the model from Section III.B in
Jordà (2005),

$\displaystyle\begin{pmatrix}y_{1t}\\\ y_{2t}\\\ y_{3t}\end{pmatrix}$
$\displaystyle=A\begin{pmatrix}y_{1t-1}\\\ y_{2t-1}\\\
y_{3t-1}\end{pmatrix}+Bh_{t}+\begin{pmatrix}\sqrt{h_{t}}\epsilon_{1t}\\\
\epsilon_{2t}\\\ \epsilon_{3t}\end{pmatrix},$ (25)
$\displaystyle(\epsilon_{1t},\epsilon_{2t},\epsilon_{3t})^{\prime}$
$\displaystyle=\boldsymbol{\epsilon}_{t}\sim N(\boldsymbol{0},I_{3}),$ (26)
$\displaystyle h_{t}$
$\displaystyle=0.5+0.5h_{t-1}+0.3\sqrt{h_{t}}\epsilon_{1t},$ (27)
$\displaystyle B=\begin{pmatrix}-1.75\\\ -1.5\\\ 1.75\end{pmatrix}~{}~{}~{}$
$\displaystyle~{}~{}~{}~{}~{}A=\begin{pmatrix}0.5&-0.25&0.25\\\
0.75&0.25&0.25\\\ -0.25&-0.25&0.75\end{pmatrix}.$ (28)

In this model, each shock $i$ affects only variable $i$ contemporaneously,
while all shocks affect all variables after one period. Contrary to the second
and the third shock, the first shock features time-varying variance.333We code
the simulation exercise following the exact code available in the replication
files from Jordà (2005).

Figure 1: Monte Carlo simulation - SVAR-GARCH

We use a simulation exercise to study how well the BART-LP methodology
recovers the true impulse responses associated with a shock to $y_{1t}$ of
size $1$. We first use generalized impulse responses to simulate the true
impulse responses associated with $\epsilon_{1t}=1$. We then generate
simulated data for 300 periods and discard the first 100, as in Jordà (2005).
The remaining 200 observations are used to estimate the pointwise linear LP
impulse response, as well as the impulse response computed from the BART-LP
method, setting $d=1$. We replicate the exercise 100 times and store the
generated estimates. The estimated model uses 2 lags, and sets $x_{t}$ equal
to the true realizations of $\epsilon_{1t}$. We set $J=250$.

Figure 1 reports the results of the analysis. The squared black line shows the
true generalized impulse responses. The shock increases $y_{1t}$ on impact by
1, and generates no contemporaneous response in $(y_{2t},y_{3t})$. It
subsequently generates an oscillating pattern in $(y_{1t},y_{2t})$ and a hump-
shaped response of $y_{3t}$. The dashed line shows the pointwise mean impulse
response from the linear LP model, with the mean computed on the pointwise
estimates along the 100 iterations. The continuous line and the shaded area
report the pointwise median and 95% band computed over the 100 median
responses from the BART-LP method, where within each iteration, the median
response from BART-LP is computed over 2,000 posterior draws. We see that
BART-LP does remarkably well in replicating the response of all variables,
capturing both the timing and the magnitude of the response correctly. By
contrast, the response estimated with the linear model estimated an attenuated
effect, underestimating the effect of the shock in the short horizon.

###  3.2 Threshold VAR model

The second model we use for simulations is

$\displaystyle\boldsymbol{y}_{t}$
$\displaystyle=[\Pi_{1}\boldsymbol{y}_{t-1}+B_{1}\boldsymbol{\epsilon}_{t}]\cdot\text{I}(y_{3,t-1}\leq
0)~{}+$ (29a)
$\displaystyle~{}+[\Pi_{2}\boldsymbol{y}_{t-1}+B_{2}\boldsymbol{\epsilon}_{t}]\cdot\text{I}(y_{3,t-1}>0),$
(29b) $\displaystyle\boldsymbol{\epsilon}_{t}$ $\displaystyle\sim
N(\boldsymbol{0},I),$ (29c)

with $\boldsymbol{y}_{t}=(y_{1t},y_{2t},y_{3t})^{\prime}$,
$\boldsymbol{\epsilon}_{t}=(\epsilon_{1t},\epsilon_{2t},\epsilon_{3t})^{\prime}$
and

$\displaystyle\Pi_{1}$ $\displaystyle=\begin{pmatrix}0.25&0.25&-0.25\\\
-0.25&0.25&-0.25\\\
0.25&0.25&0.15\end{pmatrix},~{}~{}~{}~{}~{}~{}B_{1}=\begin{pmatrix}0.10&0&0\\\
-0.20&0.15&0\\\ 0.10&-0.10&1\end{pmatrix},$ (30) $\displaystyle\Pi_{2}$
$\displaystyle=\begin{pmatrix}0.50&1.25&-1.75\\\ -0.25&0.50&-1.25\\\
0.25&0.25&0.15\end{pmatrix},~{}~{}~{}~{}~{}~{}B_{2}=\begin{pmatrix}0.10&0&0\\\
-0.20&0.15&0\\\ 0.10&-0.10&0.40\end{pmatrix}.$ (31)

The model is a recursive Threshold Vector Autoregressive model that jumps
across two regimes depending on the endogenous evolution of the third
variable. See, for instance, Castelnuovo & Pellegrino (2018).

We first compute the true generalized impulse responses to a positive one-
standard-deviation shock to a $y_{3t}$ in regime 1 by setting the initial
condition of the impulse response to $\boldsymbol{y}_{0}=\boldsymbol{0}$. We
then study how well the linear-LP and the BART-LP estimators recover the true
impulse responses. We generate a dataset of 300 observations, discard the
first 100 and use the remaining 200 to estimate the impulse responses. We
compute the point estimates of the linear LP, and compute the pointwise median
over 2,000 posterior draws from the BART-LP estimates. We then repeat the
exercise over 100 iterations. The estimated models use 4 lags and set
$(x_{t},\boldsymbol{z}_{t})$ to replicate the recursive ordering within LP
models, as discussed in Section 2.4. We set $J=250$.

Figure 2: Monte Carlo simulation - TVAR

The shock hits the system when the system is in regime 1. Under linearity, the
model would stay in regime 1, with an impulse response uniquely pinned down by
$(\Pi_{1},B_{1})$. Instead, the non-linearity of the model implies an
endogenous evolution across regimes, in accordance with the endogenous
response of $y_{3t}$. Figure 2 shows the results of the exercise. The pink
continuous line shows the true linear impulse response associated with regime
1. The squared black line captures the true generalized impulse responses. The
two lines differ. The dashed line shows the estimated response from the
linear-LP. The figure shows that this response lies between the true
generalized and the true linear impulse responses. By contrast, the impulse
response estimated via BART-LP better estimates the true generalized impulse
response. The model correctly captures the evolution of all three variables.
It detects that the endogenous evolution across regimes makes the response of
the first two variables more pronounced compared to a linear model that
remains in regime 1.

###  3.3 Sign-dependent Moving Average model

The third model we use for simulation is

$\displaystyle\boldsymbol{y}_{t}$
$\displaystyle=\sum_{l=0}^{20}\boldsymbol{\beta}_{gdp,~{}l}\cdot\epsilon_{gdp,t-l}~{}+~{}\sum_{l=0}^{20}\boldsymbol{\beta}_{\pi,l}\cdot\epsilon_{\pi,t-l}~{}+~{}$
(32)
$\displaystyle~{}~{}+\sum_{l=0}^{20}\Big{[}\boldsymbol{\beta}_{ff,l}^{+}\cdot\text{I}(\epsilon_{ff,t-l}\geq
0)~{}+~{}\boldsymbol{\beta}_{ff,l}^{-}\cdot\text{I}(\epsilon_{ff,t-l}<0)\Big{]}\cdot\epsilon_{ff,t-l},$
(33) $\displaystyle\boldsymbol{\epsilon}_{t}$ $\displaystyle\sim
N(\boldsymbol{0},I),$ (34)

with $\boldsymbol{y}_{t}=\big{(}gdp_{t},~{}\pi_{t},~{}ff_{t}\big{)}^{\prime}$
a vector containing real GDP, inflation, and the federal funds rate, and
$\boldsymbol{\epsilon}=\big{(}\epsilon_{gdp,t},~{}\epsilon_{\pi,t},~{}\epsilon_{t,ff}\big{)}^{\prime}$
a vector of structural shocks. The model is a moving average process of order
20 driven by three structural shocks, a GDP shock, an inflation shock, and a
monetary policy shock. The monetary policy shock affects the variables
differently at each horizon $t+h$, depending on whether the monetary policy
shock at time $t$ is positive or negative. The true impulse responses are
captured by $\\{\boldsymbol{\beta}_{gdp,~{}l}\\}_{l=0}^{20}$ for the GDP
shock, by $\\{\boldsymbol{\beta}_{\pi,~{}l}\\}_{l=0}^{20}$ for the inflation
shock, and by
$\\{\boldsymbol{\beta}_{ff,~{}l}^{+},~{}\boldsymbol{\beta}_{ff,~{}l}^{-}\\}_{l=0}^{20}$
for the positive and negative monetary policy shock.

We calibrate the model following an approach similar to Barnichon & Brownlees
(2019). We first estimate a recursive VAR model on US real GDP, inflation and
the federal funds rate to estimate linear impulse responses to a GDP shock, an
inflation shock and a monetary policy shock. We then set
$\\{\boldsymbol{\beta}_{gdp,l},~{}\boldsymbol{\beta}_{\pi,l},~{}\boldsymbol{\beta}_{ff,l}^{-}\\}_{l=0}^{20}$
equal to the estimated impulse responses to three shock: the GDP shock, the
inflation shock, and the monetary policy shock. Last, we set
$\boldsymbol{\beta}_{ff,l}^{+}=\boldsymbol{\beta}_{ff,l}^{-}$ for every $l$,
with two exceptions: (1) at horizons $l=2,3$ the first entry of
$\boldsymbol{\beta}_{ff,l}^{+}$ equals 3 times the first entry
$\boldsymbol{\beta}_{ff,l}^{-}$, and (2) at horizons $l=7,8,..,20$ the second
entry of $\boldsymbol{\beta}_{ff,l}^{+}$ equal 3 times the second entry of
$\boldsymbol{\beta}_{ff,l}^{-}$. This implies a data generating process in
which positive monetary policy shocks affect output and inflation 3 times more
than negative shocks in the short term (for output) and in the medium term
(for inflation) horizons.

Figure 3: Monte Carlo simulation - Sign-dependent Moving Average

We design the simulation exercise as follows. We generate an artificial
dataset of size $T=400$, discard the first 100 observations and use the
remaining 300 for estimation. We estimate the BART-LP model by controlling for
the true shocks driving the data, together with up to 2 lags of the shocks and
the endogenous variables. Last, we replicate the analysis 100 times and store
the pointwise responses from the linear and the BART LPs, computing separately
the generalized impulse responses associated with $d\pm 1$. We set $J=250$.

Figure 3 shows the impulse responses. The $(+)$ and $(-)$ red lines show the
true impulse responses associated with a positive and a negative monetary
policy shock, respectively. The responses to a negative shock are reported
with flipped sign to improve the comparison. By construction, the true
response of the federal funds rate does not change in the sign of the shock,
while the true response of GDP and inflation is 3 times stronger in response
to a positive shock in periods in which the estimated linear responses imply
strong responses of GDP and inflation. The black continuous line shows the
pointwise median response associated with the linear LP. Note that it sits
approximately halfway through the true positive and negative responses. The
blue and grey shaded areas report the 68% pointwise posterior bands associated
with the BART-LP model following a positive and a negative shock,
respectively. The BART-LP methodology estimates the true response of GDP
correctly both qualitatively and quantitatively. It correctly detects that
positive shocks have stronger short run effects on GDP than negative shocks,
and also estimates the true values of the responses accurately. As for
inflation, BART-LP correctly estimates the differential response
qualitatively, detecting the stronger effect of positive shocks. While the
response to a positive shock is somewhat underestimated, the response to a
negative shock sits approximately in the middle of the estimated bands. All in
the all, the model correctly detects that positive and negative shocks do not
have the same effect.

Table 1: Monte Carlo Probability that the effect is stronger for a positive shock | GDP | Inflation | Fed funds
---|---|---|---
$h$ | $T=100$ | $T=200$ | $T=300$ | $T=100$ | $T=200$ | $T=300$ | $T=100$ | $T=200$ | $T=300$
0 | 49 | 48 | 54 | 54 | 52 | 59 | 58 | 47 | 56
2 | 100 | 99 | 100 | 59 | 60 | 54 | 56 | 44 | 46
3 | 81 | 86 | 91 | 50 | 53 | 47 | 54 | 39 | 46
6 | 54 | 43 | 38 | 55 | 49 | 49 | 47 | 44 | 49
14 | 49 | 51 | 45 | 73 | 93 | 89 | 43 | 52 | 53
15 | 50 | 54 | 40 | 68 | 81 | 82 | 49 | 48 | 49

We further assess the ability of the model to detect differences in positive
and negative monetary policy shocks by replicating the analysis over datasets
of different lengths. We set the estimation sample period equal to
$T=100,200,300$ observations, keeping the initial discarded observations to
100. For each sample size, we replicate the analysis 100 times and compute the
percentage of Monte Carlo iterations in which the effect is stronger after a
positive rather than a negative shock in absolute value. Table 1 reports the
results for few illustrative horizons, indicating in bold the horizons and
variables for which the data generating process features a non-linearity in
the response to positive and negative shocks. On impact, the data generating
process implies no non-linearity between positive and negative shocks. Indeed,
the estimated impulse responses are well distributed across iterations, with
around 50% of the iterations detecting a stronger effect of positive shocks
and the remaining 50% detecting the opposite. At horizons 2-3, close to all
iterations detect that the effects of a monetary shock are stronger on GDP,
consistent with the data generating process, while correctly not finding
evidence of non-linearities in the response of inflation and the policy rate.
6 horizons from the shocks the model correctly detects no non-linearity in any
variable. At horizons 14-15 it correctly detects a stronger effect on
inflation associated with positive compared to negative shocks, and no non-
linearity in the remaining variables.

##  4 Empirical analysis

In this section, we apply the proposed model to two recent issues that have
featured prominently in the empirical literature on non-linear macroeconomic
dynamics.

###  4.1 Fiscal shocks during recessions and expansions

In a seminal contribution Auerbach & Gorodnichenko (2013a) use a smooth
transition VAR model to show that the response of output to government
spending shocks is larger during recessions. However, this evidence was
disputed by Ramey & Zubairy (2018) who use a 120 year sample of quarterly data
to estimate the response to military spending news shocks and show that there
is no systematic difference between the spending multiplier in normal periods
and those characterised by slack. Both papers use a non-linear LP model, but
postulate a functional form for the transition, using either a logistic
function or an indicator function.

We build on section VI-B of Ramey & Zubairy (2018) and re-visit the analysis
of state-dependent fiscal multipliers using our flexible BART-LP model from
equation (1). We set $y_{t+h}$ equal to $h$-period ahead of (1) military
spending news $(news_{t})$, (2) real per-capita government spending $(g_{t})$,
and (3) real per-capita GDP $(y_{t})$.444The data set collated by Ramey &
Zubairy (2018) is quarterly and runs from 1889 Q3 to 2015 Q3. The variables
used in our analysis are downloaded from the website of the Journal of
Political Economy. Government spending and GDP are transformed via the
procedure described in Gordon & Krenn (2010).  We first estimate a linear SVAR
model in these three variables and estimated the impulse vector
$\boldsymbol{d}$ associated with the first shock. We then compute impulse
responses as explained in Section 2.4. We include 4 lags of the variables into
the model. We set the number of trees to 250 and use 2,000 posterior draws,
with a burn-in of 1,000 draws.

We study non-linear effects of government spending shocks as follows. We
define periods of recessions as periods where $y_{t}$ is less than its 20th
percentile, while expansions denote periods where $y_{t}$ is above the 80th
percentile. We then compute generalized impulse responses by generalizing over
randomly drawn values of the conditioning values. Compared to Auerbach &
Gorodnichenko (2013a) and Ramey & Zubairy (2018), our framework also allows
for the computation of impulse responses that potentially differ for the size
of the shock. For this reason, we compute the generalized impulse responses to
either a contractionary or an expansionary fiscal spending shock, both studied
either in a recession or an expansion.

Figure 4: Response to positive (top panel) and negative (bottom panel)
military news shocks during recessions and expansions  Figure 5: Cumulated
multiplier

The top panel of Figure 4 presents the response to an expansionary spending
shock normalised to increase $g$ by 0.1 units on impact. The response of $g$
and $y$ is estimated to be more persistent during recessions. However, the
response of these two variables moves closely together in the two regimes. As
a consequence, the spending multiplier (defined as the ratio of the cumulated
response of $y$ to the cumulated response of $g$) is similar across regimes.
This can be seen from the left panel of Figure 5, which displays the estimated
multipliers and shows that evidence for a systematically larger multiplier
during recessions is weak. However, it is interesting to note that
contractionary spending shocks produce different dynamics. The bottom panel of
Figure 4 shows that the response of $y$ to these shocks during expansions is
less persistent. As a consequence the multiplier for spending cuts during
these periods is systematically smaller than the estimate during recessions.
This distinction helps reconcile the different results found by Auerbach &
Gorodnichenko (2013a) and Ramey & Zubairy (2018) by documenting that state-
specific non-linearities can be different depending on the sign if the shock.

###  4.2 Non-linear impact of Financial shocks

In a recent contribution Forni et al. (2021) use a non-linear vector moving
average model to show that large negative financial shocks have a
proportionally larger impact on the US economy. We revisit this question using
the flexible local projection proposed in this paper. We use the model from
equation (1) setting $y_{t+h}$ equal to (1) growth of industrial production
$(IP_{t})$, (2) CPI inflation $(CPI_{t})$, (3) unemployment rate $(U_{t})$,
(4) Excess bond premium $(EBP_{t})$ of Gilchrist & Zakrajšek (2012), (5) stock
returns $(STOCK_{t})$ and (6) federal funds rate $(FFR_{t})$. The data is
monthly and the sample runs from 1973M1 to 2022 M2.555The excess bond premium
is downloaded from the federal reserve website. Stock returns are calculated
using the Standards and Poor total return index obtained from Global Financial
database. The remaining variables are taken from the FRED database.

We follow the identification approach by Forni et al. (2021) and identify the
monetary policy shock as the shock that affects fast moving variables but not
slow moving variables. We estimate a preliminary linear VAR model using the
above ordering of the variables and estimate the impact vector associated with
the entry of the excess bond premium in a recursive identification. We then
add 6 lags of all variables. We simulate positive and negative shocks,
considering either a shock that moves the excess bond premium by $\pm$50 basis
points, or by $\pm$300 basis points.

Figure 6: Response to contractionary and expansionary financial shocks that
change the excess bond premium by 50 basis points. The response to
expansionary shocks has been multiplied by -1 for the purpose of comparison.
Figure 7: Response to contractionary and expansionary financial shocks that
change the excess bond premium by 300 basis points. The response to
expansionary shocks has been multiplied by -1 for the purpose of comparison

The results of the analysis are shown in Figure 6 and Figure 7, which report
the impulse responses associated with the small or the large exogenous
variation of the excess bond premium, respectively. The response to
expansionary shocks has been multiplied by -1, to improve the comparison.
Figure 6 shows that when the magnitude of the financial shock is relatively
small, systematic evidence for differences in responses across the sign of the
shock is largely absent. However, for large shocks there is clear evidence of
sign non-linearity. Figure 7 shows that shocks that increase $EBP_{t}$ by 300
basis points are associated with large declines in output, inflation, stock
returns and the interest rate while the unemployment rate rises. The responses
to expansionary shocks of the same size are substantially smaller in
magnitude. These results broadly confirm the findings reported by Forni et al.
(2021).

##  5 Conclusions

Local projections are widely used in Macroeconometrics, as they provide a
flexible tool to estimate impulse responses to structural shocks of interest.
However, the most popular linear specification of local projections introduces
the assumption of a linear relationship among variables within each horizon
$h$ considered. This paper introduces a flexible local projection that
generalises the model of Jordà (2005) to a non-parametric setting by using
Bayesian Additive Regression Trees (BART). Using Monte Carlo experiments, we
show that the model is able to capture impulse response non-linearities driven
by state-dependence, sign-dependence or size-dependence.

We apply our methodology to US fiscal and financial shocks. We show that while
it is true that the fiscal multiplier is stronger in recession, as advocated
by Auerbach & Gorodnichenko (2013a), this holds true only in response to a
contractionary fiscal shock. In response to an expansionary shock, we confirm
the result by Ramey & Zubairy (2018), namely that the fiscal multiplier does
not change significantly between recession and expansion. We then confirm the
results by Forni et al. (2021) that financial shocks have non-linear effects
on the economy. A financial shock that increases the cost of borrowing
generates contractionary effects that increase more than proportionately in
the size of the shock. This suggests that strong negative financial shocks
generate stronger effects that would be otherwise predicted using linear
models.

## References

  * (1)
  * Alloza et al. (2019) Alloza, M., Gonzalo, J. & Sanz, C. (2019), ‘Dynamic effects of persistent shocks’.
  * Alpanda et al. (2021) Alpanda, S., Granziera, E. & Zubairy, S. (2021), ‘State dependence of monetary policy across business, credit and interest rate cycles’, European Economic Review 140, 103936.
  * Auerbach & Gorodnichenko (2013a) Auerbach, A. J. & Gorodnichenko, Y. (2013a), Fiscal multipliers in recession and expansion, in A. Alesina & F. Giavazzi, eds, ‘Fiscal policy after the financial crisis’, University of Chicago Press, pp. 63–98.
  * Auerbach & Gorodnichenko (2013b) Auerbach, A. J. & Gorodnichenko, Y. (2013b), ‘Output spillovers from fiscal policy’, American Economic Review 103(3), 141–46.
  * Barnichon & Brownlees (2019) Barnichon, R. & Brownlees, C. (2019), ‘Impulse Response Estimation by Smooth Local Projections’, The Review of Economics and Statistics 101(3), 522–530.
  * Breitung et al. (2019) Breitung, J., Brüggemann, R. et al. (2019), ‘Projection estimators for structural impulse responses’, University of Konstanz Department of Economics Working Paper Series 5.
  * Bruns & Lütkepohl (2022) Bruns, M. & Lütkepohl, H. (2022), ‘Comparison of local projection estimators for proxy vector autoregressions’, Journal of Economic Dynamics and Control 134, 104277.
  * Castelnuovo & Pellegrino (2018) Castelnuovo, E. & Pellegrino, G. (2018), ‘Uncertainty-dependent effects of monetary policy shocks: A New-Keynesian interpretation’, Journal of Economic Dynamics and Control 93, 277–296.
  * Chipman et al. (2010) Chipman, H. A., George, E. I. & McCulloch, R. E. (2010), ‘BART: Bayesian additive regression trees’, The Annals of Applied Statistics 4(1), 266–298.
  * Clark et al. (2021) Clark, T. E., Huber, F., Koop, G., Marcellino, M. & Pfarrhofer, M. (2021), ‘Tail forecasting with multivariate Bayesian additive regression trees’.
  * Forni et al. (2021) Forni, M., Gambetti, L., Maffei-Faccioli, N. & Sala, L. (2021), ‘Nonlinear transmission of financial shocks: Some new evidence’.
  * Gilchrist & Zakrajšek (2012) Gilchrist, S. & Zakrajšek, E. (2012), ‘Credit spreads and business cycle fluctuations’, American Economic Review 102(4), 1692–1720.
  * Gonçalves et al. (2022) Gonçalves, S., Herrera, A. M., Kilian, L. & Pesavento, E. (2022), ‘When do state-dependent local projections work?’.
  * Gordon & Krenn (2010) Gordon, R. & Krenn, R. (2010), ‘The end of the Great Depression 1939-41: Policy contributions and fiscal multipliers’, (16380).
  * Herbst & Johannsen (2021) Herbst, E. & Johannsen, B. K. (2021), ‘Bias in local projections’.
  * Hill et al. (2020) Hill, J., Linero, A. & Murray, J. (2020), ‘Bayesian additive regression trees: A review and look forward’, Annual Review of Statistics and Its Application 7(1), 251–278.
  * Huber et al. (2020) Huber, F., Koop, G., Onorante, L., Pfarrhofer, M. & Schreiner, J. (2020), ‘Nowcasting in a pandemic using non-parametric mixed frequency VARs’, Journal of Econometrics .
  * Huber & Rossini (2021) Huber, F. & Rossini, L. (2021), ‘Inference in Bayesian additive vector autoregressive tree models’.
  * Inoue et al. (2022) Inoue, A., Rossi, B. & Wang, Y. (2022), ‘Local Projections in Unstable Environments’.
  * Jordà (2005) Jordà, Ò. (2005), ‘Estimation and inference of impulse responses by local projections’, American Economic Review 95(1), 161–182.
  * Kilian & Kim (2011) Kilian, L. & Kim, Y. J. (2011), ‘How reliable are local projection estimators of impulse responses?’, Review of Economics and Statistics 93(4), 1460–1466.
  * Koop et al. (1996) Koop, G., Pesaran, M. H. & Potter, S. M. (1996), ‘Impulse response analysis in nonlinear multivariate models’, Journal of Econometrics 74(1), 119–147.
  * Lusompa (2021) Lusompa, A. (2021), ‘Local Projections, Autocorrelation, and Efficiency’, (RWP 21-01).
  * Miranda-Agrippino & Ricco (2021) Miranda-Agrippino, S. & Ricco, G. (2021), ‘Bayesian Local Projections’.
  * Plagborg-Møller & Wolf (2021) Plagborg-Møller, M. & Wolf, C. K. (2021), ‘Local Projections and VARs Estimate the Same Impulse Responses’, Econometrica 89(2), 955–980.
  * Ramey & Zubairy (2018) Ramey, V. A. & Zubairy, S. (2018), ‘Government spending multipliers in good times and in bad: Evidence from US historical data’, Journal of Political Economy 126(2), 850–901.
  * Ruisi (2019) Ruisi, G. (2019), ‘Time-Varying Local Projections’, (891).
  * Stock & Watson (2018) Stock, J. H. & Watson, M. W. (2018), ‘Identification and estimation of dynamic causal effects in macroeconomics using external instruments’, The Economic Journal 128(610), 917–948.


# Bitext Mining for Low-Resource Languages via Contrastive Learning

Weiting Tan  Philipp Koehn
Center for Language and Speech Processing
Computer Science Department
Johns Hopkins University
{wtan12<EMAIL_ADDRESS>
###### Abstract

Mining high-quality bitexts for low-resource languages is challenging. This
paper shows that sentence representation of language models fine-tuned with
multiple negatives ranking loss, a contrastive objective, helps retrieve clean
bitexts. Experiments show that parallel data mined from our approach
substantially outperform the previous state-of-the-art method on low resource
languages Khmer and Pashto.

##  1 Introduction

Modern neural machine translation (NMT) system’s success largely depends on
the amount of high quality parallel training data.
ParaCrawl111https://paracrawl.eu/, One of the popular projects to mine
bitexts, crawls webpages and retrieve sentence pairs for various languages. In
this paper, we improve the quality of bitexts mined from ParaCrawl for two of
the low resource languages (data size smaller than 10 million), Khmer (km) and
Pashto (ps). ParaCrawl mines corpus with a pipeline of four major steps:

  1. 1.

Website Crawling: crawl websites and collect contents of web pages.

  2. 2.

Document Alignment: from collected web pages, find contents that align with
each other in different languages. Since web page has blocks of contents, this
step aligns content on document level

  3. 3.

Sentence Alignment: from aligned documents, retrieve aligned sentences by
finding matched sentences pairs of two languages.

  4. 4.

Sentence Filtering: from aligned sentences, filter out noisy sentence pairs
and use the rest as clean parallel data for downstream tasks such as training
NMT systems.

---
Figure 1: Heatmap of cosine similarity for EN-PS sentences computed by LASER
(left) and our fine-tuned Embedding (right). X-axis is Pashto sentence (left
to right: About Us, Departments, Services, Regulation). Matching pairs have
higher similarity score by our embedding. LASER embedding also mistakenly give
"Archive" high score with "Regulation".

While each step of the pipeline can be improved, we focus on sentence
alignment and sentence filtering steps, both of which could benefit from an
improved sentence scoring function. In this paper, we apply contrastive
learning Chen et al. (2020); Henderson et al. (2017) to fine-tune a sentence
transformer model and use it to align and filter sentences for Pashto and
Khmer. Our contrastively fine-tuned sentence transformer achieves better
results than previous state-of-the-art sentence representation LASER Artetxe
and Schwenk (2019b) as well as other top-performance filtering systems
Açarçiçek et al. (2020); Lu et al. (2020) submitted to WMT 2020 Corpus
Filtering Task Koehn et al. (2020). We release a toolkit222Code available at:
https://github.com/steventan0110/align-filter for replicating our experiments
and it can be applied to other language pairs with document-aligned data.

---
Figure 2: Fine-tune Sentence Transformer (BERT) with Multiple Negatives
Ranking Loss

##  2 Related Work

###  2.1 Sentence Alignment

Earlier sentence aligners use heuristics such as sentence length and word
frequency to find matching sentence pairs in two documents. One such tool,
Hunalign Varga et al. (2007), is still used today to align sentences from
document pairs in the ParaCrawl project. More recent aligners such as
Bleualign Sennrich and Volk (2010) use translation system to get both
documents into the same language and find matching sentence pairs. Another
recent aligner, Vecalign333https://github.com/thompsonb/vecalign Thompson and
Koehn (2019) uses time series warping Salvador and Chan (2007) to improve the
algorithm run-time for alignment and can be used with any cross-lingual
sentence embedding such as LASER to compute a similarity score.

###  2.2 Sentence Pair Filtering

Different filtering methods have been proposed in the past few years, for
instance in the context of the shared tasks organized by WMT Buck and Koehn
(2016); Koehn et al. (2018, 2019, 2020) and BUCC Zweigenbaum et al. (2017,
2018). The cosine similarity score computed by LASER is widely used to filter
sentences because (1) pre-trained LASER embedding models are publicly
available and (2) computing cosine similarity based on embeddings is fast.
Another high performance method is dual conditional cross-entropy Junczys-
Dowmunt (2018), which uses translation system to find maximal symmetric
agreement among sentence pairs. With recent advance on pre-trained language
models such as BERT Devlin et al. (2018), Roberta Liu et al. (2019) and XLM-
Roberta Conneau et al. (2019), proxy learning (sentence filtering as a binary
classification task) Açarçiçek et al. (2020) also performs very well.

##  3 Methodology

###  3.1 Fine-tune Sentence Transformer

Açarçiçek et al. (2020) proposed proxy learning for the sentence filtering
task. They used a pre-trained language model with a binary classification head
to detect high quality sentence pairs. Inspired by proxy learning, sentence
transformers Reimers and Gurevych (2019), and contrastive learning Chen et al.
(2020), we fine-tuned sentence transformers following figure 2 to learn a
sentence embedding for source and target languages. To fine-tune sentence
transformers (we use sentence-BERT, or SBERT), we construct positive and
negative pairs to train models with Multiple Negative Ranking Loss (MNR;
Henderson et al., 2017). Given N aligned sentence pairs
$\\{(s_{1},t_{1}),\cdots(s_{N},t_{N})\\}$, each aligned pair is a positive
sample. To construct negative samples, for any given source sentence $s_{i}$,
we use window size W to take neighbors of $t_{i}$ to form negative samples
$\\{(s_{i},t_{i-w}),(s_{i},t_{i-w+1}),\cdots,(s_{i},t_{i+w})\\}$. We also
randomly sample R sentences from target side to form negative samples with
$s_{i}$. Following MNR Loss, the training objective is computed as

$\mathcal{J}(s,t,\theta)=-\frac{1}{K}\sum_{i=1}^{K}\left[d(s_{i},t_{i})-\log\sum_{j=1}^{2W+R}e^{d(s_{i},t_{j})}\right]$
(1)

Where K is the batch size, $(s_{i},t_{i})$ is the aligned source-target
sentence pair (positive sample) and $(s_{i},t_{j})$ is the negative sample.
The distance or similarity score is measured by cosine similarity:

$d(s_{i},t_{i})=\cos(r_{s_{i}},r_{t_{i}})$ (2)

where $r_{s_{i}}$ is the high dimensional sentence representation of $s_{i}$
encoded by the pre-trained language model $\theta$. By minimizing
$\mathcal{J}(s,t,\theta)$, the model learns to maximize the difference between
similarity scores of positive and negative pairs. Therefore, models fine-tuned
with MNR not only recognize similar sentences but can also discard noisy
sentences. The advantage of our fine-tuned models against other contrastively
trained systems like Açarçiçek et al. (2020) is that our representation can be
quickly computed for millions of sentences and then used for alignment or
filtering tasks.

###  3.2 Sentence Alignment

The task of sentence alignment is to find matching sentence pairs in each
aligned document pair
$\\{D^{\text{src}},D^{\text{tgt}}\\}=\\{\\{s_{1},s_{2}\cdots,s_{m},\\},\\{t_{1},t_{2}\cdots,s_{n},\\}\\}$.
We hope to retrieve k sentence pairs
$\\{(s_{i_{1}},t_{j_{1}}),\cdots,(s_{i_{k}},t_{i_{k}})\\}$, where each index
$i_{k}(j_{k})$ correspond to a set of indexes in
$D^{\text{src}}(D^{\text{tgt}})$. For example, $i_{1}=(1,2,3),j_{1}=(1)$
stands for aligning $\\{s_{1},s_{2},s_{3}\\}$ from source to $\\{t_{1}\\}$
from target. We use Vecalign as the alignment algorithm because it is designed
to work with any high dimensional sentence embedding and it uses approximate
dynamic programming algorithm that run in $O(NM)$ time ($N$ and $M$ are the
number of sentences in source and target document). Therefore we use Vecalign
to quickly align sentences in document-aligned corpus by feeding in LASER or
our fine-tuned SBERT embedding. For details of Vecalign, we direct readers to
the original paper 444https://aclanthology.org/D19-1136.pdf.

###  3.3 Sentence Filtering

We replicated the filtering system from HUAWEI Açarçiçek et al. (2020) which
rank $1^{st}$ on corpus filtering task for Pashto and $2^{nd}$ for Khmer in
WMT 2020. HUAWEI’s system directly fine-tune language models with a binary
classification head555We follow their practice to fine-tune XLM-Roberta model
so we can filter the corpus by ranking scores predicted by the model. We also
experimented with sentence representation (LASER and our fine-tuned SBERT) to
filter corpus. Since we need to compute a similarity score based on two high
dimensional vectors, we resort to margin score Artetxe and Schwenk (2019a), a
similarity function that is shown to alleviate the "hubness" problem
Radovanović et al. (2010); Lazaridou et al. (2015). For each given sentence
pair $(x,y)$ and the encoded representation $(r_{x},r_{y})$, the score is
computed as

$\begin{split}&d(r_{x},r_{y})=\text{margin}(\cos(r_{x},r_{y}),\\\
&\sum_{z\in\text{NN}_{k}(x)}\frac{\cos(r_{x},r_{z})}{2k}+\sum_{z\in\text{NN}_{k}(y)}\frac{\cos(r_{y},r_{z})}{2k})\end{split}$
(3)

where $\text{NN}_{k}(x)$ is the k-nearest neighbor of x in the corpus666We use
Faiss Johnson et al. (2019) to compute the score and use its default value
k=4. In practice, we use ratio as the margin function, namely
$\text{margin}(a,b)=\frac{a}{b}$.

##  4 Experiments and Results

###  4.1 Mined Datasets Description

We use the evaluation setup of the WMT 2020 shared task on parallel corpus
filtering777https://www.statmt.org/wmt20/parallel-corpus-filtering.html Koehn
et al. (2020). Alignment and filtering methods are evaluated by training MT
systems on the resulting parallel corpora and assessing their quality with
BLEU Papineni et al. (2002). We start with the sentence-aligned corpus and the
document-aligned corpus provided by WMT. We denote the sentence-aligned corpus
as HUNALIGN since it is aligned by
Hunalign888http://mokk.bme.hu/en/resources/hunalign/ tool. We use Vecalign
with LASER embeddings and our fine-tuned SBERT embeddings on the released
document-aligned corpus, producing two more versions of sentence-aligned
corpora: LASER-ALIGN and SBERT-ALIGN.

For each of the three versions of parallel corpus (LASER-ALIGN, SBERT-ALIGN,
and HUNALIGN), we de-noise it with three filtering methods: LASER-FILTER,
SBERT-FILTER, HUAWEI-FILTER. Both LASER-FILTER and SBERT-FILTER rank sentences
with margin score function, with the only difference being the sentence
representation. HUAWEI-FILTER is our replication of the filtering system from
HUAWEI as described in section 3.3.

---
Figure 3: Alignment method comparison among sentence-aligned corpus HUNALIGN,
LASER-ALIGN, and SBERT-ALIGN. For each corpus, three filtering methods (LASER,
SBERT, HUAWEI) are experimented and we plot the highest score out of three.
The highest score of WMT2020 Corpus Filtering task is also shown here in red
line for comparison (We direct readers to original papers submitted to WMT2020
for detailed description of their filtering methods).

###  4.2 Results and Analysis

To evaluate performance of different methods, we rely on the BLEU score
Papineni et al. (2002) of the neural machine translation model trained
following the flores baseline
setting999https://github.com/facebookresearch/flores. Complete experimental
results (Table 1 and 2) and detailed description of preprocessing and fine-
tuning steps are included in appendix. In this section, Figures 3 and 4 are
used to help visualize results. Figure 3 shows the best BLEU score achieved by
each sentence-aligned corpus. For each corpus, we experimented with three
filtering techniques and only the highest score of the three is plotted. For
both languages (ps and km), the best score is from SBERT-ALIGN (10.43 for
Pashto and 11.83 for Khmer), which is about 1 BLEU point boost compared to the
wining system in WMT20. The substantial difference between semantic-
representation-based methods (LASER/SBERT-align) and heuristic-based method
(Hunalign) can be easily found in the plot. Between LASER-ALIGN and SBERT-
ALIGN, the advantage of the latter seems small but here we only plot the
highest score out of three filtering methods. In fact, most time HUAWEI-FILTER
works the best (and it is the most computational-intensive filtering method of
the three). Therefore, we plot figure 4 to better compare LASER and fine-tuned
SBERT. We uses LASER or SBERT for both sentence alignment and filtering steps.
SBERT-based alignment and filtering works much better than the LASER-based
method (about +2 BLEU), demonstrating SBERT’s effectiveness as alignment &
filtering technique.

Combining results above, we show that SBERT is a better representation of low
resource languages, a better quality-scoring mechanism for sentence alignment
and filtering. We believe there is still room for further improvement, since
SBERT is only fine-tuned on target language and English but WMT-released
document-aligned corpus are noisy, with boilerplate and sentences in the wrong
language. A natural next step is to fine-tune SBERT with our proposed
technique on much larger amount of data, covering more languages.

---
Figure 4: LASER and fine-tuned SBERT Comparison. For both languages, using
SBERT (red and green lines) to align and filter sentences results in >1 BLEU
score than LASER (orange and blue lines).

##  5 Conclusion

We empirically show that SBERT, fine-tuned with Multiple Negative Ranking
Loss, is a good sentence representation of low resource languages. Using our
fine-tuned SBERT as a sentence-aligner (with Vecalign) or filter (with margin-
based score) produces better training data for downstream neural machine
translation models.

## References

  * Açarçiçek et al. (2020) Haluk Açarçiçek, Talha Çolakoğlu, Pınar Ece Aktan Hatipoğlu, Chong Hsuan Huang, and Wei Peng. 2020.  Filtering noisy parallel corpus using transformers with proxy task learning.  In _Proceedings of the Fifth Conference on Machine Translation_ , pages 940–946, Online. Association for Computational Linguistics.
  * Artetxe and Schwenk (2019a) Mikel Artetxe and Holger Schwenk. 2019a.  Margin-based parallel corpus mining with multilingual sentence embeddings.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_. Association for Computational Linguistics.
  * Artetxe and Schwenk (2019b) Mikel Artetxe and Holger Schwenk. 2019b.  Massively multilingual sentence embeddings for zero-shot cross-lingual transfer and beyond.  _Transactions of the Association for Computational Linguistics_ , 7:597–610.
  * Buck and Koehn (2016) Christian Buck and Philipp Koehn. 2016.  Findings of the WMT 2016 bilingual document alignment shared task.  In _Proceedings of the First Conference on Machine Translation: Volume 2, Shared Task Papers_ , pages 554–563, Berlin, Germany. Association for Computational Linguistics.
  * Chen et al. (2020) Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton. 2020.  A simple framework for contrastive learning of visual representations.
  * Conneau et al. (2019) Alexis Conneau, Kartikay Khandelwal, Naman Goyal, Vishrav Chaudhary, Guillaume Wenzek, Francisco Guzmán, Edouard Grave, Myle Ott, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Unsupervised cross-lingual representation learning at scale.
  * Devlin et al. (2018) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2018.  Bert: Pre-training of deep bidirectional transformers for language understanding.
  * Henderson et al. (2017) Matthew Henderson, Rami Al-Rfou, Brian Strope, Yun-hsuan Sung, Laszlo Lukacs, Ruiqi Guo, Sanjiv Kumar, Balint Miklos, and Ray Kurzweil. 2017.  Efficient natural language response suggestion for smart reply.
  * Johnson et al. (2019) Jeff Johnson, Matthijs Douze, and Hervé Jégou. 2019.  Billion-scale similarity search with GPUs.  _IEEE Transactions on Big Data_ , 7(3):535–547.
  * Junczys-Dowmunt (2018) Marcin Junczys-Dowmunt. 2018.  Dual conditional cross-entropy filtering of noisy parallel corpora.
  * Koehn et al. (2020) Philipp Koehn, Vishrav Chaudhary, Ahmed El-Kishky, Naman Goyal, Peng-Jen Chen, and Francisco Guzmán. 2020.  Findings of the WMT 2020 shared task on parallel corpus filtering and alignment.  In _Proceedings of the Fifth Conference on Machine Translation_ , pages 726–742, Online. Association for Computational Linguistics.
  * Koehn et al. (2019) Philipp Koehn, Francisco Guzmán, Vishrav Chaudhary, and Juan Pino. 2019.  Findings of the WMT 2019 shared task on parallel corpus filtering for low-resource conditions.  In _Proceedings of the Fourth Conference on Machine Translation (Volume 3: Shared Task Papers, Day 2)_ , pages 54–72, Florence, Italy. Association for Computational Linguistics.
  * Koehn et al. (2018) Philipp Koehn, Huda Khayrallah, Kenneth Heafield, and Mikel L. Forcada. 2018.  Findings of the WMT 2018 shared task on parallel corpus filtering.  In _Proceedings of the Third Conference on Machine Translation: Shared Task Papers_ , pages 726–739, Belgium, Brussels. Association for Computational Linguistics.
  * Lazaridou et al. (2015) Angeliki Lazaridou, Georgiana Dinu, and Marco Baroni. 2015.  Hubness and pollution: Delving into cross-space mapping for zero-shot learning.  In _Proceedings of the 53rd Annual Meeting of the Association for Computational Linguistics and the 7th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , pages 270–280, Beijing, China. Association for Computational Linguistics.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized bert pretraining approach.
  * Lu et al. (2020) Jun Lu, Xin Ge, Yangbin Shi, and Yuqi Zhang. 2020.  Alibaba submission to the WMT20 parallel corpus filtering task.  In _Proceedings of the Fifth Conference on Machine Translation_ , pages 979–984, Online. Association for Computational Linguistics.
  * Papineni et al. (2002) Kishore Papineni, Salim Roukos, Todd Ward, and Wei-Jing Zhu. 2002.  Bleu: a method for automatic evaluation of machine translation.  In _Proceedings of the 40th Annual Meeting of the Association for Computational Linguistics_ , pages 311–318, Philadelphia, Pennsylvania, USA. Association for Computational Linguistics.
  * Radovanović et al. (2010) Miloš Radovanović, Alexandros Nanopoulos, and Mirjana Ivanović. 2010.  Hubs in space: Popular nearest neighbors in high-dimensional data.  _Journal of Machine Learning Research_ , 11(86):2487–2531.
  * Reimers and Gurevych (2019) Nils Reimers and Iryna Gurevych. 2019.  Sentence-BERT: Sentence embeddings using Siamese BERT-networks.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , pages 3982–3992, Hong Kong, China. Association for Computational Linguistics.
  * Salvador and Chan (2007) Stan Salvador and Philip Chan. 2007.  Toward accurate dynamic time warping in linear time and space.  _Intell. Data Anal._ , 11(5):561–580.
  * Sennrich and Volk (2010) Rico Sennrich and Martin Volk. 2010.  MT-based sentence alignment for OCR-generated parallel texts.  In _Proceedings of the 9th Conference of the Association for Machine Translation in the Americas: Research Papers_ , Denver, Colorado, USA. Association for Machine Translation in the Americas.
  * Thompson and Koehn (2019) Brian Thompson and Philipp Koehn. 2019.  Vecalign: Improved sentence alignment in linear time and space.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , pages 1342–1348, Hong Kong, China. Association for Computational Linguistics.
  * Varga et al. (2007) Dániel Varga, Péter Halácsy, András Kornai, Nagy Viktor, Nagy Laszlo, Németh László, and Tron Viktor. 2007.  Parallel corpora for medium density languages.
  * Zweigenbaum et al. (2017) Pierre Zweigenbaum, Serge Sharoff, and Reinhard Rapp. 2017.  Overview of the second BUCC shared task: Spotting parallel sentences in comparable corpora.  In _Proceedings of the 10th Workshop on Building and Using Comparable Corpora_ , pages 60–67, Vancouver, Canada. Association for Computational Linguistics.
  * Zweigenbaum et al. (2018) Pierre Zweigenbaum, Serge Sharoff, and Reinhard Rapp. 2018.  Overview of the Third BUCC Shared Task: Spotting Parallel Sentences in Comparable Corpora.  In _Workshop on Building and Using Comparable Corpora_ , Miyazaki, Japan.

##  6 Appendix

###  6.1 Evaluation Dataset Size

Our experiments’ results are shown in Table 1 and 2, where three types of
alignment method and three types are filtering methods are experimented (in
total 9 combinations). For each of the 9 possible alignment-filtering method,
4 versions are created based on how much data is sub-sampled. We use the
threshold 2, 3, 5, 7 million ($\\#$tokens on English side) following the
practice from WMT 2020 Corpus Filtering Task. Note that for Khmer, we see that
the BLEU scores are still going up for some alignment-filtering methods (for
example, for SBERT-ALIGN HUAWEI-FILTER, BLEU goes up from 11.13 to 11.83).
Therefore we also experimented with sub-sampling 9-million datasets and
verified that BLEU score did not increase anymore.

###  6.2 Preprocessing

For sentence alignment step, we did not employ any pre-processing techniques
because most documents contain noisy sentences and removing those sentences
would make it harder to align sentences. After retrieving sentence-aligned
datasets (LASER-ALIGN, SBERT-ALIGN, HUNALIGN), we pre-process the datasets
before sentence filtering step. First, we de-duplicate the datasets, which
filters out about 90% data (since most aligned sentences are duplicate pairs).
Second, we remove the sentence pair that has over 90% overlap between source
and target side sentences. Lastly, We also use fasttext language
id101010https://fasttext.cc/docs/en/language-identification.html to check
every aligned sentence pair and remove it if its English side is not predicted
as en. Note that this is a very lenient filter given the noisy sentence-
aligned dataset we retrieved. In fact, language id filtering plays an
important role for LASER-FILTER for km-en task. The BLEU score under LASER-
FILTER is significantly worse than the other two filtering methods, especially
when sub-sample size is small (2 or 3 million tokens). This is because LASER
would select many sentence pairs that is not Khmer as top-scoring pairs. When
filtering out sentences based on language id for both English and Khmer,
LASER-FILTER can achieve better performances (though still worse than our
SBERT-FILTER and HUAWEI-FILTER results), similar to the scores reported in WMT
2020 Corpus Filtering Task.

###  6.3 Fine-tune Sentence-BERT

To build SBERT-ALIGN corpus, we fine-tune the SBERT model as described in
section 3.1 and figure 2. To fine-tune SBERT, we need a parallel corpus to
sample positive and negative pairs from. We experimented with both HUNALIGN
and LASER-ALIGN corpus. It is unsurprising that LASER-ALIGN works better
because it has more correctly aligned sentences. Thus, we fine-tuned SBERT
based on the LASER-ALIGN corpus and then use it to align sentences from
document-aligned data, producing the sentence-aligned corpus SBERT-ALIGN.

Alignment Type | LASER Filter | SBERT Filter | HUAWEI Filter
---|---|---|---
2M | 3M | 5M | 7M | 2M | 3M | 5M | 7M | 2M | 3M | 5M | 7M
HUNALIGN | 6.19 | 7.53 | 7.71 | 7.86 | 7.96 | 8.73 | 9.05 | 8.44 | 7.48 | 8.70 | 9.26 | 8.33
LASER-ALIGN | 6.04 | 7.73 | 7.03 | 7.99 | 8.41 | 9.17 | 9.82 | 8.76 | 9.35 | 10.18 | 10.38 | 8.70
SBERT-ALIGN | 6.13 | 7.72 | 8.66 | 8.44 | 8.13 | 9.08 | 10.07 | 9.00 | 9.19 | 10.02 | 10.43 | 9.61
Table 1: BLEU Score of Neural Machine Translation model trained on ps-en datasets that are mined with different sentence-alignment and sentence-filtering methods. Each row corresponds to one sentence-aligner and each column corresponds to one sentence-filter. Under each filter methods, four sub-columns indicate the number of millions of English tokens sub-sampled from the corpus following the practice from WMT20 Shared Task on Corpus Filtering Açarçiçek et al. (2020). The bold numbers are the highest scores for each alignment type, which are used to plot figure 3 Alignment Type | LASER Filter | SBERT Filter | HUAWEI Filter
---|---|---|---
2M | 3M | 5M | 7M | 2M | 3M | 5M | 7M | 2M | 3M | 5M | 7M
HUNALIGN | 5.60 | 6.70 | 7.59 | 8.37 | 7.94 | 8.27 | 8.33 | 8.02 | 8.77 | 9.25 | 10.38 | 10.03
LASER-ALIGN | 5.29 | 6.91 | 9.05 | 9.33 | 9.03 | 9.44 | 10.13 | 10.79 | 8.96 | 9.86 | 11.11 | 11.68
SBERT-ALIGN | 4.96 | 7.02 | 8.65 | 9.14 | 9.22 | 9.45 | 10.18 | 10.12 | 8.85 | 10.15 | 11.13 | 11.83
Table 2: BLEU Score of Neural Machine Translation model trained on km-en
datasets that are mined with different sentence-alignment and sentence-
filtering methods. The bold numbers are the highest scores for each alignment
type, which are used to plot figure 3