Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
91,684 characters · 18 sections · 97 citation commands
Local Projections vs. VARs: Lessons From Thousands of DGPs
}
Keywords: external instrument, impulse response function, local projection, proxy variable, structural vector autoregression. JEL codes: C32, C36.
Since Jorda2005 introduced the popular local projection (LP) impulse response estimator, there has been a debate about its benefits and drawbacks relative to Vector Autoregression (VAR) estimation Sims1980. Recently, Plagborg2020 proved that these two methods in fact estimate precisely the same impulse responses asymptotically, provided that the lag length used for estimation tends to infinity. This result holds regardless of identification scheme and regardless of the underlying data generating process (DGP). Nevertheless, the question of which estimator to choose in finite samples remains open. It is also an urgent question, since researchers have remarked that LPs and VARs can give conflicting results when applied to central economic questions such as the effects of monetary or fiscal stimulus Ramey2016,Nakamura2018.
Whereas the LP estimator utilizes the sample autocovariances flexibly by directly projecting an outcome at the future horizon $h$ on current covariates, a VAR($p$) estimator instead extrapolates longer-run impulse responses from the first $p$ sample autocovariances. Hence, though the estimates from the two methods agree approximately at horizons $h \leq p$, they can disagree substantially at intermediate and long horizons.\footnote{See Plagborg2020 for a formal result.} Intuitively, the extrapolation employed by VARs should yield a lower variance but potentially a higher bias than for LPs, perfectly analogous to the trade-off between direct and iterated reduced-form forecasts Schorfheide2005,Kilian2017.\footnote{The trade-off is also conceptually similar to the relationship between polynomial series estimators and kernel estimators in cross-sectional nonparametric regression.} How much more should one care about bias than variance to optimally choose the LP estimator over the VAR estimator in realistic sample sizes? And how does the trade-off depend on the DGP? Unfortunately, these questions are challenging to answer analytically, due to the dynamic and nonlinear nature of the time series estimators, as well as the breadth of DGPs encountered in applied practice.
In this paper we illuminate the bias-variance trade-off in impulse response estimation through a comprehensive simulation study, applying LP and VAR methods to thousands of empirically relevant DGPs. Our goal is to identify which estimators perform well on average across many DGPs and thus may serve as practical default procedures. Rather than insisting on the usual binary distinction between “local projections” and “VARs”, we furthermore consider an entire menu of related estimation approaches that employ bias correction, shrinkage, or model-averaging. We find that the usual least-squares LP estimator tends to have lower bias than the least-squares VAR estimator, as expected, but also that this bias reduction comes at the cost of substantially higher variance. Out of all the procedures we analyze, bias-corrected LP is the most attractive estimator if and only if the researcher overwhelmingly prioritizes bias. If, however, the researcher also cares about precision (as in the conventional mean squared error criterion), then VAR methods are the most attractive; in particular, Bayesian VARs perform well at short and long horizons, while it is difficult to beat the least-squares VAR estimator at intermediate and long horizons.
Our simulation study considers an extensive array of DGPs, obtained by drawing specifications at random from a large-scale, empirically calibrated dynamic factor model (DFM).\footnote{Our overall approach is inspired by Lazarus2018, who are instead interested in the question of how to select among different long-run variance estimators.} We fit the DFM to the data set of Stock2016, which contains a large number of quarterly U.S. macroeconomic time series spanning a wide variety of variable categories. As emphasized by Stock2016, such DFMs can accurately capture the joint co-movements of conventional macroeconomic data, and so our simulation results will be informative about the universe of standard U.S. time series. This estimated DFM exhibits realistic and complex dynamics in the short and long run, including cointegrating relationships among the latent factors. From the encompassing 207-variable DFM we then draw 6,000 random subsets of five variables (subject to constraints that emulate applied practice); all results reported below are essentially unchanged if instead we limit attention to 17 of the most commonly used macro series out of the 207. The randomly drawn subsets of time series constitute the set of DGPs that we consider for our simulation study. As the calibrated DFM is known to us, we can compute the true impulse responses, and therefore also estimator biases and mean squared errors. Importantly, none of these many DGPs can be exactly represented as a finite-order VAR model, yielding a non-trivial bias-variance trade-off between the LP and VAR estimators. Moreover, our DGPs exhibit substantial heterogeneity in how well they can be approximated by VAR models, in persistence and shape of impulse response functions, and in the invertibility of the structural shocks, consistent with the heterogeneity faced by applied researchers. While our results inevitably depend on the specification of the encompassing model, we believe that an estimation method that works well across our multitude of empirically calibrated DGPs has substantial promise as a default procedure.
We study the ability of several variants of LP and VAR methods to accurately estimate impulse response functions. Consistent with the majority of applied work, the estimators are applied to data in levels, rather than transforming to stationarity prior to estimation. Since VARs with very large lag lengths are asymptotically equivalent to LPs Plagborg2020,Xu2023, we focus on VAR estimators with moderate lag length choices, as conventionally found in the literature. In addition to the popular least-squares LP and VAR estimators, we further enrich the bias-variance possibility frontier by considering: (i) small-sample bias correction of the VAR coefficients Pope1990,Kilian1998 and LP impulse response estimates Herbst2021; (ii) penalized LP Barnichon2019, which smooths out impulse response functions; (iii) Bayesian VAR estimation, with priors selected as in Giannone2015; and (iv) model averaging of univariate and multivariate VAR models of various lag lengths Hansen2016. For each estimation method, we consider three oft-used structural identification schemes: observed shocks, instrumental variables (IVs)/proxies, and recursive identification. For IV identification, we further distinguish between internal IV methods Ramey2011,Plagborg2020 and external IV methods Stock2008,Stock2012bpea,Mertens2013. We then evaluate the performance of these estimators through the lens of loss functions with varying weights on bias and variance.
Applying the estimation methods to simulated data from the thousands of DGPs, a clear and unavoidable bias-variance trade-off emerges. We highlight four main lessons:
Our findings provide a novel perspective on recent work emphasizing the potential dangers of VAR model mis-specification Ramey2016,Nakamura2018. We consider DGPs that do not admit finite-order VAR representations, so VAR methods indeed suffer from larger bias, as cautioned there. Reducing that bias via direct projection, however, tends to incur a steep cost in terms of increased sampling variance at intermediate and long horizons. Researchers who prefer to employ LP estimators should therefore be prepared to pay that price, and furthermore should apply the Herbst2021 bias correction procedure when their data is persistent, as is usually the case.
\paragraph{Literature.} Our simulation study is inspired by the seminal work of Marcellino2006 on direct and iterated multi-step forecasts, though we focus instead on structural impulse responses. While simulation studies in the forecasting literature often analyze low-dimensional specifications, we consider multi-variable systems, consistent with standard practice in the applied structural macroeconometrics literature. The structural perspective also requires us to contend with issues such as the variety of different popular shock identification schemes, normalization of impulse responses, and the special role of external instrumental variables.
Our large-scale model set-up differs from prior simulation studies of LP and VAR methods, which have considered at most a handful of DGPs. Examples here include Jorda2005, Meier2005, Kilian2011, Brugnolini2018, Choi2019, Austin2020, and Bruns2021. These papers either obtain their DGPs from stylized, low-dimensional VARMA models, calibrated DSGE models, and/or a few empirically calibrated VAR models. Our encompassing DGP is instead designed to closely mimic applied practice: we consider a non-stationary DFM with rich common and idiosyncratic dynamics that accurately captures key properties of the kinds of aggregate time series typically used in standard macroeconometric analyses. Our analysis also differs in the following respects: we consider shrinkage estimation procedures as competitors to the least-squares estimators; we study several popular structural identification schemes; and we examine how our conclusions vary with the impulse response horizon and the researcher's loss function. All these features are essential to the above-mentioned main lessons that we draw from our results.
Even though the simulation results are at the heart of our analysis, we start off by illustrating the bias-variance trade-off through an analytical example that builds on Schorfheide2005. That paper develops a general theory of the asymptotic bias and variance of direct and iterated (reduced-form) forecasts under local mis-specification. While these theoretical results are valuable for analytically distilling the forces at work, they do not by themselves resolve the bias-variance trade-off faced by practitioners, as this trade-off invariably depends in a complicated fashion on many features of the DGP.
Finally, we stress that our paper focuses solely on point estimation, as opposed to inference or hypothesis testing. See Inoue2020, MontielOlea2020, and Xu2023 for theoretical as well as simulation results on VAR and LP confidence interval procedures. Moreover, we focus exclusively on impulse response estimands, rather than variance decompositions or historical decompositions.
\paragraph{Outline.} (ref) illustrates the bias-variance trade-off for LP and VAR estimators using a simple analytical example. (ref) describes the empirically calibrated dynamic factor model that we use to generate our many DGPs. (ref) defines the menu of LP- and VAR-based estimation procedures. (ref) contains our main simulation results and robustness checks. (ref) summarizes the lessons for applied researchers and then offers guidance for future research. The appendix contains implementation details. A supplemental appendix with proofs and further simulation results as well as a Matlab code suite are available online.\footnote{\url{https://github.com/dake-li/lp_var_simul} }
This section motivates our simulation study with an analytical discussion of the bias-variance trade-off between LP and VAR impulse response estimators. (ref) analyzes these estimators in the context of a simple toy model that cleanly illustrates the trade-off, and (ref) connects this analytical discussion to the rest of the paper.
Plagborg2020 show that the impulse response estimands of VAR and LP estimators with $p$ lags generally differ at horizons $h > p$: the VAR extrapolates from the first $p$ sample autocovariances, while LP exploits all autocovariances out to horizon $h+p$. This observation suggests the presence of a bias-variance trade-off whenever the true DGP is not a finite-order VAR, perfectly analogous to the choice between “direct” and “iterated” predictions in multi-step forecasting Marcellino2006. We here formalize this basic intuition by extending the arguments of Schorfheide2005 to structural impulse response estimation in a simple, albeit non-stationary DGP.
\paragraph{Model.} Consider a simple sequence of drifting DGPs for the scalar time series $y_t$:
where $\varepsilon_t \equiv (\varepsilon_{1,t},\varepsilon_{2,t})'$ is an i.i.d. white noise process with $\operatorname*{Var}(\varepsilon_t) = \operatorname*{diag}(1,\sigma_2^2)$, and $y_0=0$. We assume that the researcher observes $w_t \equiv (\varepsilon_{1,t},y_t)'$, i.e., she observes the shock $\varepsilon_{1,t}$ but not $\varepsilon_{2,t}$. The above DGP drifts towards a unit-root VAR(1) process in $w_t$ at rate $T^{-1/2}$, where $T$ is the sample size. We show below that this ensures a non-trivial bias-variance trade-off in the limit $T\to\infty$. The DGP captures the notion that finite-order autoregressive models are often a good---but not exact---approximation to the true underlying DGP. The degree of autoregressive mis-specification is governed by the parameter $\alpha$.\footnote{To interpret its units, consider a distributed lag regression of $\Delta y_t$ on $\varepsilon_{1,t}$, $\varepsilon_{1,t-1}$, and $\varepsilon_{1,t-2}$. Then it is standard to show that the t-statistic for significance of the second lag converges in distribution to $N(\alpha/\sigma_2,1)$.}
We are interested in the impulse responses of $y_t$ with respect to a unit impulse in $\varepsilon_{1,t}$. The true impulse response function implied by the model (ref) equals $\theta_{h,T} \equiv 1+\tau \mathbbm{1}(h \geq 1) + \alpha T^{-1/2}\mathbbm{1}(h \geq 2)$ at horizon $h$. This impulse response function reflects---in stark fashion---the common empirical finding that signal-to-noise ratios are especially low at longer horizons, here $h \geq 2$, in the sense that the increment $\theta_{2,T}-\theta_{1,T}$ is of the same asymptotic order as the standard errors of the LP and VAR estimators, as shown formally below.
\paragraph{Estimators.} For now, we consider two estimators of $\theta_{h,T}$.
\paragraph{Trade-off.} Along the stated asymptote, the researcher faces a clear bias-variance trade-off between the LP and VAR impulse response estimators:
At horizons $h \in \lbrace 0,1\rbrace$, there is no bias-variance trade-off: on impact, the two estimators are numerically equivalent; at $h=1$, both are asymptotically unbiased with identical asymptotic variance, consistent with Plagborg2020. Intuitively, the equivalence at $h=1$ reflects the fact that the VAR(1) estimator does not extrapolate, instead reporting the direct projection of $y_{t+1}$ on $w_t$, exactly as LP does Plagborg2020.
At horizons $h \geq 2$ (i.e., exceeding the lag length used for estimation), the bias-variance trade-off is non-trivial. Specifically, the asymptotic biases satisfy $|\operatorname{aBias}_{\text{VAR},h}|=|\alpha| > 0 =|\operatorname{aBias}_{\text{LP},h}|$ whenever $\alpha \neq 0$, while the asymptotic variances satisfy $\operatorname{aVar}_{\text{LP},h}-\operatorname{aVar}_{\text{VAR},h} = 1+\sigma_2^2+(h-2)[(1+\tau)^2+\sigma_2^2]>0$. Intuitively, LP directly projects $y_{t+h}$ on the shock $\varepsilon_{1,t}$, which is uncorrelated with any lagged controls, so the asymptotic bias is always zero. In contrast, the VAR(1) estimator extrapolates the response at horizon $h$: the model's structure implies that the precisely estimated autocovariances at lag 1 suffice to compute impulse responses at longer horizons. Though this tight parametric extrapolation yields a low variance relative to LP, it incurs a bias due to dynamic mis-specification when $\alpha \neq 0$. In the simple DGP (ref), the asymptotic bias of the VAR estimator could be eliminated by simply increasing the lag length to 2 or higher, but in practice it may be difficult to determine the appropriate lag length, as we demonstrate below in (ref). The fact that both the asymptotic bias and variance of the VAR impulse response estimator are constant at horizons $h \geq 2$ is a special feature of the stylized DGP (ref).\footnote{In this DGP, the VAR coefficients on lagged $y_t$ are estimated super-consistently due to the unit root, so to first order, estimation uncertainty arises only from the coefficients on lagged $\varepsilon_{1,t}$.} Nevertheless, our simulation study below will demonstrate the robustness of the qualitative predictions that (i) the bias of VAR is high at intermediate and long horizons relative to LP, and (ii) the difference between the variance of LP and that of VAR tends to increase as a function of the horizon.\footnote{We derived similar analytical results for a stationary DGP in a previous working paper version of this article Li2022.}
How does the optimal choice of estimator depend on the researcher's preferences concerning bias and variance? To evaluate the performance of a given estimator $\hat{\theta}_h$ of $\theta_{h,T}$, we will throughout this paper consider loss functions of the form\footnote{The objective function (ref) is not a loss function in the usual decision theoretic sense (which would call it a risk function when $\omega=\frac{1}{2}$). We proceed with the non-standard terminology for ease of exposition.}
For $\omega = \frac{1}{2}$, this is proportional to the mean squared error (MSE). For $\omega > \frac{1}{2}$, the researcher is more concerned about (squared) bias than variance, and for $\omega=1$ the researcher exclusively cares about bias. Substituting the asymptotic bias and variance expressions in (ref) into the above loss function, we find that LP is preferred over VAR (asymptotically) if and only if the researcher prioritizes bias sufficiently heavily at the expense of variance, namely when $\omega \geq \omega_h^* \equiv 1-\alpha^2/(\alpha^2 + \operatorname{aVar}_{\text{LP},h}-\operatorname{aVar}_{\text{VAR},h}) \in (0,1)$ (focusing here on the interesting case $h \geq 2$ and $\alpha \neq 0$). We remark that---even in the very simple DGP (ref)---the indifference weight $\omega_h^*$ depends sensitively on all the model parameters and the horizon $h$.
Because analytical bias-variance calculations will invariably end up depending in complicated ways on a multitude of parameters, we will in the rest of this paper use simulations to explore the nature of the bias-variance trade-off across a rich and empirically relevant set of DGPs. In the language of (ref), these DGPs will inform us about empirically plausible degrees of mis-specification $\alpha$, impulse response function shapes $\tau$, and relative shock importances $\sigma_2^2$, and therefore about the practically relevant bias weight $\omega_h^*$ necessary to justify the use of one linear projection technique over another one. Moreover, we will also consider several variants of the standard least-squares LP and VAR estimators, thus allowing us to further trace out the bias-variance possibility frontier.
This section presents our DGPs. We define the empirically calibrated encompassing model in (ref), from which we draw thousands of DGPs with corresponding structural impulse response estimands, as described in (ref). We discuss implementation details in (ref), and provide summary statistics for the DGPs in (ref). Various modifications to this baseline set of DGPs are considered later in (ref).
We construct our simulation DGPs from an encompassing model that is known to accurately capture the time series properties of many U.S. macroeconomic time series: a dynamic factor model (DFM) fitted to the well-known Stock2016 data set. Because we seek to follow applied practice in using data in levels rather than first differences, we employ a non-stationary variant of the DFM estimated by Stock2016.
The DFM postulates that a large-dimensional $n_X \times 1$ vector $X_t$ of observed macroeconomic time series is driven by a low-dimensional $n_f \times 1$ vector $f_t$ of latent factors, as well as an $n_X \times 1$ vector $v_t$ of idiosyncratic components. The latent factors are assumed to follow a non-stationary Vector Error Correction Model (VECM) with VAR($p_f$) representation
where $\varepsilon_t=(\varepsilon_{1,t},\dots,\varepsilon_{n_f,t})'$ is an $n_f \times 1$ vector of aggregate shocks, which are i.i.d. and mutually uncorrelated, with $\operatorname*{Var}(\varepsilon_t) = I_{n_f}$. The $n_f \times n_f$ matrix $H$ determines the impact impulse responses of the factors with respect to the aggregate shocks. The observed macroeconomic series $X_t$ are given by
where the idiosyncratic component $v_{i,t}$ for macro observable $X_{i,t}$ follows the potentially non-stationary AR($p_v$) process
with $\xi_{i,t}$ i.i.d. across $t$ and $i$. We assume that all shocks and innovations are jointly normal and homoskedastic. We will next in (ref) describe how we construct our many lower-dimensional DGPs from this encompassing large-scale DFM; (ref) then follows up with implementation details, including in particular a discussion of how the parameters of this non-stationary DFM are calibrated to the Stock2016 data set.
We use the encompassing model (ref)--(ref) to build thousands of lower-dimensional DGPs for our simulation study. Specifically, for each DGP, we draw a random subset of $n_{\bar{w}}$ variables $\bar{w}_t$ from the large vector $X_t$, i.e., $\bar{w}_t \subset X_t$. The variables $\bar{w}_t$ follow the time series process implied by the encompassing model (ref)--(ref). In particular, $\bar{w}_t$ is driven by some combination of aggregate structural shocks $\varepsilon_t$ and idiosyncratic components $v_t$. We draw thousands of such random combinations of variables, thus yielding thousands of lower-dimensional DGPs. The details of how we select the variable combinations are postponed until (ref).
For each DGP drawn in this way, we consider three types of structural impulse response estimands, chosen to mimic as closely as possible popular schemes for identifying the effects of policy shocks in applied macroeconometrics Ramey2016,Stock2016. In the following, $y_t \in \bar{w}_t$ denotes a response variable of interest in the DGP, $i_t \in \bar{w}_t$ is a policy variable used to normalize the scale of the shock (if applicable), $z_t$ is an external instrument (if applicable), and $w_t$ denotes the vector of all observed time series in the DGP.
This section first discusses how we estimate the DFM and then specifies the particular DGPs and structural impulse responses that we consider in the simulation study.
\paragraph{DFM parameters.} We parametrize the DFM (ref)--(ref) by estimating the model on the Stock2016 data set. Recall that we model the variables in levels rather than first differences, unlike Stock2016. We provide a brief overview of our approach here, with details in (ref).
We begin with the vector of observables $X_t$. As in Stock2016, that vector contains quarterly observations on 207 time series for 1959Q1--2014Q4, mostly consisting of real activity variables, price measures, interest rates, asset and wealth variables, and productivity series.\footnote{Table 1 and the Data Appendix of Stock2016 list all variables and their categories.} Each series is seasonally adjusted as in Stock2016. However, unlike those authors we do not transform the series to stationarity. Instead, variables that they transform to (non-log or log) first differences, we now keep in (non-log or log) levels; and variables that they transform to log second differences (which are mostly price indices), we only transform to log first differences. All estimation procedures mentioned below control for series-specific and common deterministic linear time trends.
Our estimation approach is intended to allow for rich long- and short-run dynamics, opening the door for meaningful mis-specification of short-lag VARs. Following Bai2004 and Barigozzi2021, we estimate the non-stationary DFM by extracting factors from differenced (and subsequently de-meaned) data, cumulating these factors, and then fitting a VECM to the cumulated factors; the aforementioned papers show that this estimation strategy consistently estimates the true VECM parameters under weak conditions that allow for cointegration. We set the number of factors $n_f$ equal to 6 as in Stock2016. The VECM is estimated by quasi-maximum-likelihood without restricting the adjustment coefficients or cointegrating relations.\footnote{The reason we fit a VECM rather than an unrestricted VAR in levels to the factors is that the VAR estimator may underestimate persistence in finite samples, as is well known. While bias correction procedures exist, they may not always work well in practice. The VECM approach instead errs on the side of overstating the role of permanent shocks, consistent with our goal of allowing for rich long-run dynamics.} The cointegration rank of the VECM is selected by the Johansen1995 maximum eigenvalue test, which indicates that the latent factors are driven by four common stochastic trends. As in Stock2016, we fit AR($p_v$) processes by OLS to each idiosyncratic residual after removing the estimated factors, separately for each $i$. We use lag lengths $p_f=p_v = 4$ for both the factor process and the idiosyncratic component processes, which is at the upper end of what is preferred by the Akaike Information Criterion, consistent with our goal of allowing rich dynamics. The above-mentioned estimation procedure pins down all parameters of the DFM except for the structural impact response matrix $H$; we discuss below how we construct that matrix.
While the estimated encompassing DFM assumes a cointegrated VECM for the latent factors, the lower-dimensional DGPs that we subsequently extract from the DFM will not satisfy exact finite-order VECM or VAR processes, and will not be exactly cointegrated. This follows from the presence of the idiosyncratic components and the mismatch in dimensions between the latent factors and the observable series (as specified below). Indeed, we show below that most of the lower-dimensional DGPs feature a combination of exact unit roots (imposed in the factor VECM), several roots near unity (owing partly to the factor process and partly to the idiosyncratic component processes), as well as smaller roots that induce transitory dynamics. Our DGPs are therefore consistent with the common empirical finding that there is often substantial ambiguity about the appropriate VAR lag lengths, the exact magnitude of roots, and the presence or absence of cointegrating relationships.
\paragraph{DGP and estimand selection.} To provide a comprehensive picture of the bias-variance trade-off, we select thousands of different sets of observables $\bar{w}_t \subset X_t$. We consider two protocols for selecting these observables---one aimed at mimicking monetary policy shock applications, and one aimed at fiscal policy shock applications. Specifically, for each type of policy shock, we randomly draw 3,000 configurations of $n_{\bar{w}} = 5$ macroeconomic observables $\bar{w}_t$. Thus, we end up with a total of 6,000 DGPs. For the monetary policy DGPs we restrict $\bar{w}_t$ to always contain the federal funds rate, while for the fiscal policy DGPs we restrict $\bar{w}_t$ to contain federal government spending. These two series are chosen as the policy variables $i_t$ for the IV and recursive estimands. The remaining four variables in $\bar{w}_t$ are then selected uniformly at random from $X_t$, except we impose that at least one variable should be a measure of real activity, and at least one other variable a measure of prices.\footnote{Real activity series correspond to categories 1--3 in the classification in Table 1 of Stock2016, while price series correspond to category 6.} The impulse response variable $y_t$ is selected uniformly at random from the four series (other than $i_t$).
For each of the DGPs, we implement the structural impulse response estimands as follows:
Consistent with the experience of applied researchers, our DGPs exhibit substantial heterogeneity along several dimensions. (ref) displays the distribution of various population parameters across our 6,000 DGPs. The table focuses on impulse responses with respect to directly observed monetary policy and government spending shocks, though results for recursively defined shocks are similar, as shown in (ref).
First of all, the DGPs feature varying degrees of persistence. All DGPs have unit roots by construction; nevertheless, the DGPs differ in how heavily they load on the various non-stationary and stationary linear combinations of the latent factors. (ref) reports the ratio of the traces of the long-run variance matrix and variance matrix applied to differenced data, $\operatorname*{trace}(\mathit{LRV}(\Delta \bar{w}_t))/\operatorname*{trace}(\operatorname*{Var}(\Delta \bar{w}_t))$. This measure varies widely across the DGPs, with median equal to 1.02 (as when all series are simple random walks), and the 90th percentile equal to 3.54 (consistent with strong positive autocorrelation of the first differences). We will consider an alternative set of moderately persistent, stationary DGPs in one of our main robustness checks in (ref).
Second, the DGPs are heterogeneous in terms of how well they can be approximated by a low-order VAR. (ref) reports the ratio $\sum_{\ell=5}^{1000} \|A_\ell^w\|/\sum_{\ell=1}^{1000} \|A_\ell^w\|$, which measures the relative magnitude of the coefficient matrices $\lbrace A_\ell^w \rbrace_{\ell}$ in the VAR($\infty$) representation for $\lbrace \bar{w}_t \rbrace$ at or after lag 5 (with $\|\cdot\|$ here denoting the Frobenius matrix norm). The 10th and 90th percentiles equal 0.14 and 0.37, respectively. Hence, the analysis in (ref) suggests that the bias of low-order VAR procedures will vary substantially across the various DGPs that we consider in our simulations.
Third, for the IV specifications, we note that our DGPs differ in terms of shock invertibility and IV strength. The degree of invertibility is defined as the R-squared in a population projection of the shock of interest on current and lagged macro observables $\lbrace \bar{w}_{t-\ell} \rbrace_{\ell=0}^\infty$.\footnote{Projections on infinite collections of lagged variables are defined as the limit when the lag length tends to $\infty$, using a diffuse initialization of the Kalman filter.} The bias of some SVAR-based external instrument procedures depends on how far below 1 this measure is, as discussed further in (ref). The table shows that 90% of the DGPs have degrees of invertibility below 49%, i.e., substantial non-invertibility.\footnote{Leeper2013 argue that adding forward-looking variables to a VAR ameliorates the invertibility problem. However, if we restrict attention to the 1,457 DGPs that contain at least one time series in the “Asset Price & Sentiment” category (see (ref)), the 90th percentile of the degree of invertibility increases only marginally to 51%.} This is not surprising: the DFM (ref)--(ref) features a realistic amount of idiosyncratic noise $v_t$, making it challenging to accurately back out the aggregate shock $\varepsilon_{1,t}$ from a small number of time series $\bar{w}_t$. The strength of the IV is by construction borderline weak to moderate, as the population first stage F-statistic (from a regression of the policy variable $i_t$ on the IV $z_t$, controlling for lagged data) is calibrated to vary between approximately 10 and 30, given sample size $T=200$.
Finally, the true impulse response estimands exhibit a wide variety of shapes. (ref) shows that the impulse response functions peak at very different horizons and are typically not simple monotonically decaying or even hump-shaped functions: the median number of interior local extrema of the impulse response functions is 2 (a monotonic function would have 0; a hump-shaped function would have 1). Many impulse response functions change sign at some horizon, as evidenced by the average response (across horizons) typically being much smaller than the maximal response. Finally, the smoothness of the impulse response functions varies substantially: the R-squared value in a regression of the impulse responses $\lbrace \theta_h \rbrace_{h=0}^{20}$ on a quadratic polynomial $b_0+b_1 \times h + b_2 \times h^2$ has 10th and 90th percentiles given by 0.46 and 0.98, respectively. For further illustration, (ref) displays the true values of six impulse response functions, providing a representative picture of the heterogeneity. The figure illustrates that, while some impulse response functions approximately return to 0 at long horizons, many do not, and some have the largest response even beyond horizon $h = 20$.
We now give a brief overview of the different VAR- and LP-based estimation methods that we consider in the simulation study.\footnote{To visualize the various estimation methods, (ref) plots the estimated impulse response functions in a few data sets simulated from a single DGP.} Though all these methods aim at estimating the same population impulse responses defined in (ref), they differ in terms of their bias-variance properties, and in terms of their robustness to non-invertibility. Further implementation details are relegated to (ref). All estimators include an intercept.
\paragraph{Local projection approaches.} The basic idea behind local projections, as proposed by Jorda2005, is to estimate the impulse responses separately at each horizon by a direct regression of the future outcome on current covariates. We consider three such approaches:
\paragraph{VAR approaches.} Like local projections, a VAR with lag length $p$ flexibly estimates the impulse responses out to horizon $p$; however, the VAR extrapolates the responses at longer horizons $h>p$ using only the sample autocovariances out to lag $p$. As suggested by the analysis in (ref), this tends to generate impulse response estimates with lower variance but higher bias than LP estimates at intermediate and long horizons. We consider four such VAR-based approaches:
Observed shock identification is carried out by simply ordering the shock first in the recursive VAR. Recursive identification is implemented as usual in the VAR literature. We consider two different approaches to IV estimation:
We implement the “internal instruments” approach using all four types of VAR estimation techniques described earlier. For brevity, we only consider the least-squares version of the “external instrument” SVAR-IV estimator.
\paragraph{Lag length selection.} As a baseline, the LP and VAR estimators use $p=4$ lags for estimation (except of course VAR model averaging, which uses many different lag lengths). In our DGPs, the Akaike Information Criterion almost always selects very short lag lengths $\hat{p}_{AIC}$, as we discuss further in (ref) below. Thus, for all intents and purposes, our results may be interpreted as having been generated by the lag length selection rule $p = \max\lbrace \hat{p}_{AIC}, 4\rbrace$. Our reading of applied practice is that researchers typically include at least 4 lags in quarterly data. Results for $p=8$ are discussed in (ref).
This section presents our simulation results. We summarize the results through four lessons, presented in separate subsections. The first three lessons focus on observed shock identification. The fourth lesson is concerned with IV identification. We show in (ref) that these conclusions are qualitatively robust to several alterations of our baseline simulation specification (including less persistent DGPs and recursive identification). Finally, in (ref), we justify our focus on the average performance of estimators across DGPs, by arguing that there is limited scope for selecting among estimators in a data-dependent way.
Throughout this section we present results for our 6,000 monetary and fiscal policy shock DGPs considered jointly rather than separately. For each DGP, we simulate time series of length $T=200$ quarters and approximate the population bias and variance of the estimators by averaging across 5,000 Monte Carlo simulations. The main results (excluding robustness checks) take about one week to produce in Matlab on a research computing cluster with 300 parallel cores.
Our first takeaway is that researchers invariably face a bias-variance trade-off: because most of our DGPs are not well approximated by finite-order VAR models, least-squares LPs tend to have lower bias, while least-squares VAR estimators tend to have lower variance, consistent with the simple analytical example provided in (ref). Strictly speaking, these statements are only exactly true for the bias-corrected versions of the estimators Herbst2021,Pope1990,Kilian1998, as the high persistence of our DGPs imparts a sizable finite-sample bias in the estimators at intermediate and long horizons, particularly for LPs. This bias correction, however, is not a free lunch, as it increases variance.
(ref) depict the bias-variance trade-off at various horizons. These figures show the median (across our 6,000 DGPs) of the absolute bias $|E(\hat{\theta}_h-\theta_h)|$ or the standard deviation $\sqrt{\operatorname*{Var}(\hat{\theta}_h)}$, respectively, as a function of the horizon. The different lines correspond to different estimators $\hat{\theta}_h$, with least-squares LP and VAR being the thick lines. Before taking the median, we cancel out the units of the response variables by dividing the bias and standard deviation by $\sqrt{\frac{1}{21}\sum_{h=0}^{20}\theta_h^2}$, i.e., the root mean squared value of the true impulse response function out to horizon 20. Note that the scale of the vertical axis differs between the bias and standard deviation plots.
The figures show that least-squares LP and VAR estimators have similar bias and variance at horizons $h \leq p = 4$, but not at longer horizons $h > p$. The median biases then generally increase with the horizon, with the bias of VAR exceeding that of LP, except at long horizons.\footnote{Kilian2011 find in simulations that LP does not have lower bias than VAR estimators, but they consider a different variant of LP that uses an auxiliary VAR to identify the structural shocks.} While the median standard deviation of LP is increasing in the horizon, that of VAR instead displays a hump-shaped pattern. At long horizons, the median standard deviation of LP is about double that of VAR. These observations are broadly consistent with the asymptotic results in (ref), Schorfheide2005, and Plagborg2020.
Our results also show that the bias correction procedure of Herbst2021 is critical to achieving uniformly low bias for the LP approach. Though the asymptotic bias of LP is zero when the shock is observed, as discussed in (ref), the high persistence of our DGPs implies that the small-sample bias of least-squares LP is non-negligible at intermediate and long horizons, especially the latter. The bias-corrected version of LP proposed by Herbst2021 (thin line with small dots in the figures) eliminates about a third of the bias at all horizons. In comparison with the LP case, bias correction is not as critical for VAR estimation, though the bias-corrected VAR estimator (dashed line) does have a somewhat lower bias than the least-squares VAR estimator at long horizons. After bias correction, LP has lower (median) bias than VAR at all horizons, as predicted by asymptotic theory. We further show in (ref) below that such bias correction is not nearly as important in less persistent, stationary DGPs.
Bias correction is not a free lunch, however, as it is associated with a substantial increase in variance. (ref) shows that the bias-corrected LP and VAR estimators have uniformly higher median standard deviation than the uncorrected estimators. In fact, bias-corrected LP has not only the uniformly lowest median bias among the methods we consider, it also has the uniformly highest standard deviation. (ref) show head-to-head comparisons of the least-squares and bias-corrected estimators for the LP and VAR cases, respectively. The figures show the fraction of DGPs for which the least-squares estimator achieves a lower loss (ref) than the bias-corrected estimator, as a function of the horizon $h$ and the weight $\omega$ attached to squared bias in the loss function; to interpret the figures, recall that $\omega=0.5$ corresponds to MSE loss, while $\omega=1$ corresponds to an exclusive focus on bias at the expense of variance. The darker the plot, the more often is the least-squares estimator preferred over the bias-corrected one. Evidently, one has to attach a very high weight $\omega$ to bias in the loss function to prefer the bias-corrected estimator in more than 60% of DGPs; furthermore, a researcher with MSE loss would usually prefer the uncorrected estimators.
Our second takeaway is that bias-corrected LP is the single best estimator in our choice set if and only if the researcher's loss function overwhelmingly prioritizes bias. In contrast, uncorrected LP is never the best option if the goal is to minimize average loss across our DGPs. Under MSE loss, penalized LP typically outperforms the other LP procedures as it has substantially lower variance, though at the expense of a moderate increase in bias.
(ref) shows the optimal estimation method as a function of the horizon $h$ and the bias weight $\omega$. The colors and patterns indicate the estimation method that minimizes the average loss (ref) across DGPs, after normalizing the loss to cancel out units as in (ref). In this subsection we focus on the top part of (ref), i.e., where the weight $\omega$ on bias in the loss function is high. Bias-corrected LP emerges as the best estimator at most horizons in this case, as is to be expected given its excellent bias properties in (ref); nevertheless, the figure shows that the optimality of bias-corrected LP is predicated on $\omega$ exceeding roughly 0.9, or even higher at some horizons, corresponding to an overwhelming focus on minimizing bias rather than variance. In contrast, uncorrected least-squares LP is essentially dominated: it has greater bias than bias-corrected LP (notably at longer horizons), yet materially higher variance than least-squares VAR or other shrinkage methods, and so no part of (ref) is orange with diagonal lines. We discuss the rest of (ref) in the next subsection.
(ref) compare bias-corrected LP to bias-corrected VAR and to penalized LP, respectively. The former figure shows that bias-corrected LP is only preferred to bias-corrected VAR in at least 60% of DGPs when $\omega \geq 0.9$. In the latter figure, we see that the smoothing of impulse responses across horizons that the penalized LP estimator performs is usually attractive whenever $\omega \leq 0.9$, except at very short and very long horizons. By “betting on smoothness”, penalized LP achieves a substantial variance reduction relative to the un-penalized LP procedures, at the expense of a moderate increase in bias, see (ref). In fact, there is a region of (ref) with intermediate horizons and moderately high weight on bias where penalized LP (green with diagonal cross-hatching) is the single best estimator. These findings underscore our conclusion that, across the majority of the DGPs, the use of bias-corrected LP can only be justified by committing to a nearly exclusive focus on minimizing bias, with little regard for precision.
Our third takeaway is that VAR estimators are attractive to researchers who place at least moderate weight on variance in their loss function. But the choice of VAR method depends on the horizon: Bayesian VARs perform well at short horizons, least-squares VARs at intermediate horizons, and at long horizons the two are comparable. VAR model averaging, on the other hand, performs poorly regardless of bias-variance preferences.
Returning to (ref), we see that for bias weights $\omega$ below 0.9, the optimal estimation method is almost always either least-squares VAR (purple areas) or BVAR (solid-dotted blue). The key attractive property of BVAR is that it has the lowest (median) standard deviation at all horizons among the methods we consider, as seen in (ref), though it also has high bias relative to least-squares VAR at intermediate horizons, as shown in (ref). The relatively high bias at intermediate horizons is possibly due to the fact that its prior specification, which is conventional in the literature, is motivated by one-step-ahead and long-run forecasting properties, as opposed to medium-run properties.\footnote{Moreover, the Giannone2015 approach of choosing the prior hyper-parameters to maximize the marginal likelihood implicitly targets one-step-ahead forecasts (see Equation 5 in their paper).}
(ref) shows that the head-to-head performance of least-squares VAR vs. Bayesian VAR depends on the horizon. At short horizons $h \leq 4$, BVAR is preferred in the majority of DGPs, and indeed it is the overall best estimator for most loss functions that place non-trivial weight on variance (see (ref)). However, at intermediate horizons $h \in [5,12]$, least-squares VAR is preferred over BVAR in the clear majority of DGPs for most loss functions, and the former estimator is the overall preferred method for loss functions with $\omega \leq 0.8$. At long horizons $h \geq 13$, the two VAR methods are comparable and outperform all other methods, unless the weight on bias in the loss function is high.
Finally, we remark that bias-corrected VAR and VAR model averaging are rarely, if ever, optimal. Bias-corrected VAR (yellow with horizontal lines in (ref)) can be rationalized at short horizons if the concern for bias is high, but the difference compared to least-squares VAR is small at these horizons, as discussed in (ref). VAR model averaging performs poorly regardless of loss function and horizon, as it has substantial bias as well as a high standard deviation relative to other VAR-based estimators (see (ref)). Closer inspection reveals that the high standard deviation is a consequence of a very fat-tailed sampling distribution, with a non-negligible probability of erratic estimates.\footnote{We use Hansen2016's (Hansen2016) code off the shelf. It would be interesting to investigate whether the procedure could be modified to avoid erratic estimates, perhaps by regularizing the averaging weights.}
Our last takeaway is concerned with IV/proxy identification. Among the invertibility-robust “internal instruments” estimators, the bias-variance trade-off is very similar to that already discussed above for the case of an observed shock. The alternative “external instruments” SVAR-IV procedure, however, contributes starkly to the trade-off: it can be severely biased due to its lack of robustness to non-invertibility, but at the same time it also has substantially lower dispersion than the “internal instruments” procedures.
Since first and second moments of IV estimators may not exist theoretically Sawa1972, we in this subsection report median bias (i.e., in each DGP, the median of the estimation error) instead of (mean) bias, and the interquartile range instead of the standard deviation.\footnote{For completeness, (mean) bias and standard deviation are reported in (ref).} We refer to the latter as “dispersion.”
(ref) show the median bias and interquartile range of the various IV estimators. If we ignore the dotted line representing SVAR-IV, these figures are qualitatively similar to those presented in (ref). However, SVAR-IV stands out by exhibiting especially high median bias and especially low interquartile range at all horizons. This is consistent with the existing theoretical work referenced in (ref): unlike the “internal instruments” procedures, SVAR-IV is asymptotically biased when the shock is not invertible, and we saw in (ref) that the degree of invertibility is generally low in our DGPs.\footnote{Consistent with theory, we furthermore find that the median bias of SVAR-IV is particularly large relative to other estimation methods in the subset of DGPs with the smallest degree of invertibility. See (ref).} On the other hand, the SVAR-IV procedure has fewer parameters to estimate (as it excludes the IV $z_t$ from the reduced-form VAR regression), causing a reduction in dispersion relative to the other procedures. Though we view the high median bias of SVAR-IV across our DGPs as worrying, its low dispersion is intriguing and may in some cases trump the bias concerns.
This section argues that our main conclusions in (ref) are robust to several alterations of our baseline simulation specification. We pay particular attention to an exercise that replaces our non-stationary encompassing DFM with a stationary version. Various other robustness checks are listed subsequently, with details relegated to (ref).
\paragraph{Stationary DGPs.} While the majority of applied papers estimate VARs and LPs with a mix of non-stationary and stationary variables in levels Ramey2016, in some cases researchers transform all their data to stationarity prior to the analysis. To cover such applications, we have repeated our analysis using the stationary estimated DFM of Stock2016 as our encompassing model. We construct impulse response estimands as before and compare the performance of the same estimation methods, except that the BVAR estimator uses a prior that shrinks towards white noise rather than random walks. Details on the implementation and results are presented in (ref).
Our headline qualitative conclusions go through in the stationary DGPs. We observe the same bias-variance trade-off as in our main analysis, with LPs achieving lower bias than VARs at the cost of elevated variance. As a result, except for researchers that exclusively prioritize bias, least-squares VARs or some kind of shrinkage---in the form of Bayesian VARs or penalized LPs---are preferred. The two most notable differences from our baseline analysis are that (i) penalized LP outperforms BVAR for MSE loss at very short horizons, and (ii) due to the moderate persistence of the stationary DGPs, bias correction has less bite, and uncorrected least-squares LP has near-zero bias at all horizons.
\paragraph{Other robustness checks.} The following modifications to our baseline simulation specification all leave our main conclusions qualitatively unchanged.
It is natural to ask whether, instead of selecting estimators based on average performance across DGPs, the choice of estimator can be guided by the data at hand in each given DGP. We now show that this appears to be difficult, as conventional model selection or evaluation criteria are unable to detect even substantial mis-specification of the VAR(4) model in the vast majority of our DGPs. These findings are consistent with the previously documented poor performance of the VAR model averaging estimator. For simplicity, we focus here on observed shock identification.
First, the Akaike Information Criterion tends to select very short lag lengths $\hat{p}_{AIC}$ in our DGPs, as already mentioned earlier. The 90th percentile of $\hat{p}_{AIC}$ (across simulations) does not exceed 2 in any of our 6,000 DGPs, and it in fact equals 2 in only 68.3% of those DGPs. This frequently used model selection tool therefore essentially never indicates that the VAR(4) specification is mis-specified.
Second, the Lagrange Multiplier test of residual serial correlation has low power in most of our DGPs. We carry out this test by regressing the sample VAR residuals on their first lags, controlling for four lags of the observed variables, and employing the likelihood ratio test defined in Johansen1995. Using a 10% significance level for the test, only around 8% of the DGPs exhibit a rejection probability above 25%, and none of the DGPs have a rejection probability above 50%. Hence, this conventional specification test of the VAR(4) model is under-powered, despite the fact that many of our DGPs are in fact not well approximated by a VAR(4) model in population, as shown in (ref).
It is of course possible that other model selection criteria or specification tests will work better. However, at a minimum, the performance of the VAR model averaging estimator discussed in (ref) and the evidence presented in this subsection together suggest that it is not straightforward to develop effective data-dependent estimator selection rules for use on conventional macroeconomic time series data.
We conducted a large-scale simulation study of the performance of LP and VAR structural impulse response estimators, as well as several variants of these methods. We drew the following four main conclusions.
These conclusions inevitably depend on the choice of encompassing model and the specific implementation of the impulse response estimators. Our paper first and foremost has aimed to bring the bias-variance trade-off in impulse response estimation to the attention of applied researchers. Our particular quantification of this trade-off has sought to capture the wide range of applied settings faced by macroeconomists, by fitting a dynamic factor model with rich short-run and long-run dynamics to the well-known Stock2016 data set. Our online code repository (see (ref)) facilitates experimentation with alternative encompassing models or estimation procedures.
Our findings point to several potential areas for future research. First, we conjecture that the bias-variance trade-off may differ quantitatively in panel data settings, to the extent that the availability of a large cross section reduces the sampling variance of the estimators for a given time dimension, thus potentially making LP relatively more attractive than in the pure time series case. Second, our analysis has focused on the average performance of estimators across DGPs because we find that conventional model selection or evaluation tools are unable to detect substantial mis-specification of low-order VARs in our simulations; nevertheless, we view data-dependent estimator selection as an area ripe for further investigation. Third, it may be worth investigating whether the performance of the Bayesian VAR procedure at intermediate horizons can be improved by developing alternative prior specifications that are specifically aimed at structural impulse response estimation rather than forecasting, unlike the priors used in much of the literature. Fourth, for the case of IV/proxy identification, an interesting question is whether it is possible to develop alternative invertibility-robust estimation procedures that capture some of the variance improvement enjoyed by the non-robust SVAR-IV estimator. Fifth, we leave exploration of other structural shock identification schemes---such as sign restrictions, long-run restrictions, and non-recursive short-run restrictions---to future work. Sixth, while our simulations were calibrated to quarterly data, it would be illuminating to see whether our conclusions apply also to monthly calibrations.