The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
78,020 characters
Local Projection Inference is Simpler and More Robust Than You Think
\title{\texorpdfstring{\vspace{-2em}}{}Local Projection Inference is Simpler \texorpdfstring{\\}{} and More Robust Than You Think\thanks{Email: {\tt [email removed]}, {\tt [email removed]}. We are grateful for comments from two anonymous referees, Jushan Bai, Otavio Bartalotti, Guillaume Chevillon, Max Dovi, Marco Del Negro, Domenico Giannone, Nikolay Gospodinov, Michael Jansson, \`Oscar Jord\`a, Lutz Kilian, Michal Koles\'{a}r, Simon Lee, Sophocles Mavroeidis, Konrad Menzel, Ulrich M\"{u}ller, Serena Ng, Elena Pesavento, Mark Watson, Christian Wolf, Tao Zha, and numerous seminar participants. We would like to especially thank Atsushi Inoue and Anna Mikusheva for a very helpful discussion of our paper. Montiel Olea would like to thank Qifan Han and Giovanni Topa for excellent research assistance. Plagborg-M{\o}ller acknowledges that this material is based upon work supported by the NSF under Grant {\#}1851665. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF.}}
\author{Jos\'{e} Luis Montiel Olea \\ Columbia University \and Mikkel Plagborg-M{\o}ller \\ Princeton University}
\date{First version: March 17, 2020 \texorpdfstring{\\[0.5ex]}{}This version: December 4, 2020}
\maketitle
\begin{abstract}
Applied macroeconomists often compute confidence intervals for impulse responses using local projections, i.e., direct linear regressions of future outcomes on current covariates. This paper proves that local projection inference robustly handles two issues that commonly arise in applications: highly persistent data and the estimation of impulse responses at long horizons. We consider local projections that control for lags of the variables in the regression. We show that lag-augmented local projections with normal critical values are asymptotically valid uniformly over (i) both stationary and non-stationary data, and also over (ii) a wide range of response horizons. Moreover, lag augmentation obviates the need to correct standard errors for serial correlation in the regression residuals. Hence, local projection inference is arguably both simpler than previously thought and more robust than standard autoregressive inference, whose validity is known to depend sensitively on the persistence of the data and on the length of the horizon.
\end{abstract}
\noindent \emph{Keywords:} impulse response, local projection, long horizon, uniform inference.
\section{Introduction}
Impulse response functions are key objects of interest in empirical macroeconomic analysis. It is increasingly popular to estimate these parameters using the method of \emph{local projections} \citep{Jorda2005}: simple linear regressions of a future outcome on current covariates \citep{Ramey2016,Angrist2018,Nakamura2018,Stock2018}. Since local projection estimators are regression coefficients, they have a simple and intuitive interpretation. Moreover, inference can be carried out using textbook standard error formulae, adjusting for serial correlation in the (multi-step forecast) regression residuals.
Despite its popularity, there exist no theoretical results justifying the use of local projection \emph{inference} over autoregressive procedures. From an identification and estimation standpoint, \citet{Kilian2017} and \citet{PlagborgMoller2019} argue that neither local projections nor Vector Autoregressions (VARs) dominate the other in terms of mean squared error in finite samples, and in population the two methods are equivalent. However, from an inference perspective, the only available guidance on the relative performance of local projections comes in the form of a small number of simulation studies, which by necessity cannot cover the entire range of empirically relevant data generating processes.
In this paper we show that---in addition to its intuitive appeal---frequentist local projection inference is robust to two common features of macroeconomic applications: highly persistent data and the estimation of impulse responses at long horizons. Key to our result is that we consider \emph{lag-augmented} local projections, which use lags of the variables in the regression as controls. Formally, we prove that standard confidence intervals based on such lag-augmented local projections have correct asymptotic coverage \emph{uniformly} over the persistence in the data generating process and over a wide range of horizons.\footnote{We focus on marginal inference on individual impulse responses, not \emph{simultaneous} inference on a vector of several response horizons \citep{IK2016,MontielOlea2019}.} This means that confidence intervals remain valid even if the data exhibits unit roots, and even at horizons $h$ that are allowed to grow with the sample size $T$, e.g., $h=h_T \propto T^{\eta}$, $\eta \in [0,1)$. In fact, when persistence is not an issue, and the data is known to be stationary, local projection inference is also valid at \emph{long} horizons; i.e., horizons that are a non-negligible fraction of the sample size ($h_T \propto T$).
Lag-augmenting local projections not only robustifies inference, it also simplifies the computation of standard errors by obviating the adjustment for serial correlation in the residuals. It is common practice in the local projections literature to compute Heteroskedasticity and Autocorrelation Consistent/Robust (HAC/HAR) standard errors \citep{Jorda2005,Ramey2016,Kilian2017,Stock2018}. Instead, we prove that the usual Eicker-Huber-White heteroskedasticity-robust standard errors suffice for \emph{lag-augmented} local projections. The reason is that, although the regression residuals are serially correlated, the \emph{regression scores} (the product of the residuals and residualized regressor of interest) are serially uncorrelated under weak assumptions. This finding further simplifies local projection inference, as it side-steps the delicate choice of HAR procedure and associated difficult-to-interpret tuning parameters \citep[e.g.,][]{Lazarus2018}.
The robustness properties of lag-augmented local projection inference stand in contrast to the well-known fragility of standard autoregressive procedures. Textbook autoregressive inference methods for impulse responses (such as the delta method) are invalid in some cases with near-unit roots or medium-long to long horizons (e.g., $h_T \propto \sqrt{T})$, as discussed further below. We show that lag-augmented local projection inference is valid when the data has near-unit roots and the horizon sequence satisfies $h_{T}/T \rightarrow 0$. Though the method fails in the case of unit roots and very long horizons $h_T \propto T$, existing VAR-based methods that achieve correct coverage in this case are either highly computationally demanding or result in impractically wide confidence intervals. When the data is stationary and interest centers on short horizons, local projection inference is valid but less efficient than textbook AR inference. Thus, the robustness afforded by our recommended procedure is not a free lunch. We provide a detailed comparison with alternative inference procedures in \cref{sec:comp} below.
Our results rely on assumptions that are similar to those used in the literature on autoregressive inference. In particular, we assume that the true model is a VAR($p$) with possibly conditionally heteroskedastic innovations and known lag length. We discuss the choice of lag length $p$ in \cref{sec:conc}. The key assumption that we require on the innovations is that they are conditionally mean independent of both past and \emph{future} innovations (which is trivially satisfied for i.i.d. innovations). Our strengthening of the usual martingale difference assumption is crucial to avoid HAC inference, but we show that the assumption is satisfied for a large class of conditionally heteroskedastic innovation processes. The robustness property of local projection inference only obtains asymptotically if the researcher controls for all $p$ lags of all of the variables in the VAR system. Thus, our paper highlights the advantages of multivariate modeling even when using single-equation local projections.
To illustrate our theoretical results, we present a small-scale simulation study suggesting that lag-augmented local projection confidence intervals achieve a favorable trade-off between coverage and length. Since local projection estimation is subject to small-sample biases just like VAR estimation \citep{Herbst2020}, we consider a simple and computationally convenient bootstrap implementation of local projection. The simulations suggest that non-augmented autoregressive procedures with delta method standard errors have more severe under-coverage problems than local projection inference, especially at moderate and long horizons. Autoregressive confidence intervals can be meaningfully \emph{shorter} than lag-augmented local projection intervals in \emph{relative} terms, but in \emph{absolute} terms the difference in length is surprisingly modest. Our simulations also indicate that lag-augmented local projections with heteroskedasticity-robust standard errors have better coverage/length properties than more standard \emph{non-augmented} local projections with off-the-shelf HAR standard errors. Finally, although the lag-augmented autoregressive bootstrap procedure of \citet{Inoue2020} achieves good coverage, it yields prohibitively wide confidence intervals at longer horizons when the data is persistent.
\paragraph{Related Literature.}
It is well known that standard autoregressive (AR) inference on impulse responses requires an auxiliary rank condition to rule out super-consistent limit distributions, thus yielding a $\sqrt{T}$-normal limit with strictly positive variance, see Assumption B of \citet{Inoue2020}. When this rank condition holds, the textbook AR impulse response estimator is asymptotically normal even in the presence of (near-)unit roots \citep{Inoue2002}. However, there are two common features of the data that lead to violations of the rank condition. First, the condition can fail when some linear combinations of the variables exhibit no persistence \citep{Benkwitz2000}. Second, in the presence of (near-)unit roots, certain linear combinations of the autoregressive coefficients are necessarily super-consistent \citep{Sims1990}. This compromises textbook AR inference for certain combinations of impulse response horizons and parameter values that typically cannot be ruled out \emph{a priori}, especially in AR(1) or VAR(1) models, but also in higher-order autoregressions (\citealp{Phillips1998}; \citealp[Remark 3, p. 455]{Inoue2020}). In an important paper, \citet{Inoue2020} show that \emph{lag-augmented} autoregressive inference solves the rank problem caused by (near-)unit roots, but data generating processes that lack persistence still need to be ruled out \emph{a priori}. We build on their ideas, which in turn are based on \citet{Toda1995} and \citet{Dolado1996}. As we show, the validity of lag-augmented \emph{local projection} (LP) inference does not hinge on auxiliary rank conditions.
Moreover, the validity of textbook AR inference is also compromised when the length of the impulse response horizon is large \citep{Pesavento2007,Mikusheva2012}. Standard bootstrap methods rectify some of these problems, but not all. Several papers have proposed AR-based methods for impulse response inference at \emph{long} horizons $h=h_T \propto T$ \citep{Wright2000,Gospodinov2004,Pesavento2007,Mikusheva2012,Inoue2020}. With the exception of \citet{Mikusheva2012}, the literature on long-horizon inference has exclusively focused on near-unit root processes as opposed to devising uniformly valid procedures. The \citet{Hansen1999} grid bootstrap analyzed by \citet{Mikusheva2012} is asymptotically valid at short and long horizons. However, it is not valid at intermediate horizons (e.g., $h_T \propto \sqrt{T}$), unlike the LP procedure we analyze. \citeauthor{Mikusheva2012} argues, though, that the grid bootstrap is \emph{close} to being valid at intermediate horizons, although it is much more computationally demanding than our recommended procedure, especially in VAR models with several parameters. \citet{Inoue2020} show that a version of the Efron bootstrap confidence interval, when applied to lag-augmented AR estimators, is valid at long horizons. We show that this procedure delivers impractically wide confidence intervals at moderately long horizons when the data is persistent, unlike lag-augmented LP.
We appear to be the first to prove the \emph{uniform} validity of lag-augmented LP inference. \citet{Mikusheva2007,Mikusheva2012} and \citet{Inoue2020} derive the uniform coverage properties of various AR inference procedures, but they do not consider LP. The \emph{pointwise} properties of LP procedures have been discussed by \citet{Jorda2005}, \citet{Kilian2017}, and \citet{Stock2018}, among others. \citet{Kilian2011} and \citet{Brugnolini2018} present simulation studies comparing AR inference and LP inference. \citet{Brugnolini2018} finds that the lag length in the LP matters, which is consistent with our theoretical results.
Though the theoretical results in this paper appear to be novel, \citet[Section 5]{Dufour2006} and \citet{Breitung2019} have discussed some of the main ideas presented herein. First, both these papers state that lag augmentation in LP avoids unit root asymptotics, but neither paper considers inference at long horizons or derives uniform inference properties. Second, \citet{Breitung2019} further argue that HAC inference in LP can be avoided if the true model is a VAR($p$), although it is not clear from their discussion what are the assumptions needed for this to be true. Neither of these papers provide results concerning the efficiency of lag-augmented LP inference relative to other lag-augmented or non-augmented inference procedures, as we do in \cref{sec:comp}.
Local projections are closely related to multi-step forecasts. \citet{Richardson1989} and \citet{Valkanov2003} develop a non-standard limit distribution theory for long-horizon forecasts. \citet{Chevillon2017} proves a robustness property of direct multi-step inference that involves non-normal asymptotics due to the lack of lag augmentation. \citet{Phillips2013} test the null hypothesis of no long-horizon predictability using a novel approach that requires a choice of tuning parameters, but yields uniformly-over-persistence normal asymptotics. This test is based on an estimator with a faster convergence rate than ours in the non-stationary case. However, to the best of our knowledge, their approach does not carry over immediately to impulse response inference, and it is not obvious whether the procedure is uniformly valid over both short and long horizons.
\paragraph{Outline.}
\cref{sec:ar1} provides a non-technical overview of our results in the context of a simple AR(1) model, including an illustrative simulation study. \cref{sec:comp} provides an in-depth comparison of lag-augmented LP with other inference procedures. \cref{sec:var} presents the formal uniformity result for a general VAR($p$) model. \cref{sec:boot} describes a simple bootstrap implementation of lag-augmented local projection that we recommend for practical use. \cref{sec:conc} concludes. Proofs are relegated to \cref{sec:var_proof_main} and the Online Supplement. \cref{sec:appendix,sec:var_XpX_ar1} contain further simulation and theoretical results. The supplement and a full Matlab code repository are available online.\footnote{\url{https://github.com/jm4474/Lag-augmented_LocalProjections} \label{fn:github}}
\section{Overview of the Results}
\label{sec:ar1}
This section provides an overview of our results in the context of a simple univariate AR(1) model. The discussion here merely intends to illustrate our main points. \cref{sec:var} presents general results for VAR($p$) models.
\subsection{Lag-Augmented Local Projection}
\label{sec:ar1_intuition}
\paragraph{Model.}
Consider the AR(1) model for the data $\lbrace y_t \rbrace$:
\begin{equation} \label{eqn:ar1}
y_t = \rho y_{t-1} + u_t,\quad t=1,2,\dots,T,
\quad y_0 = 0.
\end{equation}
The parameter of interest is a nonlinear transformation of $\rho$, namely the impulse response coefficient at horizon $h \in \mathbb{N}$. We denote this parameter by $\beta(\rho,h) \equiv \rho^h$. In \cref{sec:var} below we argue that the zero initial condition $y_0=0$ is not needed for our results to go through. Our main assumption in the univariate model is:
\begin{asn} \label{asn:u_mds}
$\lbrace u_t \rbrace$ is strictly stationary and satisfies $E(u_t \mid \lbrace u_s \rbrace_{s \neq t})=0$ almost surely.
\end{asn}
\noindent The assumption requires the innovations to be mean independent relative to past and future innovations. This is a slight strengthening of the usual martingale difference assumption on $u_t$. \cref{asn:u_mds} is trivially satisfied if $\lbrace u_t \rbrace$ is i.i.d., but it also allows for stochastic volatility and GARCH-type innovation processes.\footnote{For example, consider processes $u_t = \tau_t \varepsilon_t$, where $\varepsilon_t$ is i.i.d. with $E(\varepsilon_t)=0$, and for which one of the following two sets of conditions hold: (a) $\lbrace \tau_t \rbrace$ and $\lbrace \varepsilon_t \rbrace$ are independent processes; or (b) $\tau_t$ is a function of lagged values of $\varepsilon_t^2$, and the distribution of $\varepsilon_t$ is symmetric. \cref{asn:u_mds} is in principle testable, but that is outside the scope of this paper.}
\paragraph{Local Projections With and Without Lag Augmentation.}
We consider the local projection (LP) approach of \cite{Jorda2005} for conducting inference about the impulse response $\beta(\rho, h)$. A common motivation for this approach is that the AR(1) model \eqref{eqn:ar1} implies
\begin{equation} \label{eqn:ytph_decomp}
y_{t+h} = \beta(\rho,h)y_t + \xi_t(\rho,h),
\end{equation}
where the regression residual (or \emph{multi-step forecast error}),
\[\xi_t(\rho,h) \equiv \sum_{\ell=1}^h \rho^{h-\ell}u_{t+\ell},\]
is generally serially correlated, even if the innovation $u_t$ is i.i.d.
The most straight-forward LP impulse response estimator simply regresses $y_{t+h}$ on $y_t$, as suggested by equation \eqref{eqn:ytph_decomp}, but the validity of this approach is sensitive to the persistence of the data. Specifically, this standard approach leads to a non-normal limiting distribution for the impulse response estimator when $\rho \approx 1$, since the regressor $y_t$ exhibits near-unit-root behavior in this case. Hence, inference based on normal critical values will not be valid uniformly over all values of $\rho \in [-1,1]$ even for fixed forecast horizons $h$. If $\rho$ is safely within the stationary region, then the LP estimator is asymptotically normal, but inference generally requires the use of Heteroskedasticity and Autocorrelation Robust (HAR) standard errors to account for serial correlation in the residual $\xi_t(\rho,h)$.
To robustify and simplify inference, we will instead consider a \emph{lag-augmented} local projection, which uses $y_{t-1}$ as an additional control variable. In the autoregressive literature, ``lag augmentation'' refers to the practice of using more lags for estimation than suggested by the true autoregressive model. Define the covariate vector $x_t \equiv (y_t,y_{t-1})'$. Given any horizon $h \in \mathbb{N}$, the lag-augmented LP estimator $\hat{\beta}(h)$ of $\beta(\rho,h)$ is given by the coefficient on $y_t$ in a regression of $y_{t+h}$ on $y_t$ and $y_{t-1}$:
\begin{equation}
\begin{pmatrix} \label{eqn:LPestimates}
\hat{\beta}(h) \\
\hat{\gamma}(h)
\end{pmatrix} \equiv \left(\sum_{t=1}^{T-h} x_tx_t'\right)^{-1}\sum_{t=1}^{T-h}x_t y_{t+h}.
\end{equation}
Here $\hat{\beta}(h)$ is the impulse response estimator of interest, while $\hat{\gamma}(h)$ is a nuisance coefficient.
The purpose of the lag augmentation is to make the effective regressor of interest stationary even when the data $y_t$ has a unit root. Note that equations \eqref{eqn:ar1}--\eqref{eqn:ytph_decomp} imply
\begin{equation} \label{eqn:ytph_decomp2}
y_{t+h} = \beta(\rho,h)u_t + \beta(\rho,h+1)y_{t-1} + \xi_t(\rho,h).
\end{equation}
If $u_t$ were observed, the above equation suggests regressing $y_{t+h}$ on $u_t$, while controlling for $y_{t-1}$. Intuitively, this will lead to an asymptotically normal estimator of $\beta(\rho,h)$, since the regressor of interest $u_t$ is stationary by \cref{asn:u_mds}, and we control for the term that involves the possibly non-stationary regressor $y_{t-1}$. Fortunately, due to the linear relationship $y_t=\rho y_{t-1}+u_t$, the coefficient $\hat{\beta}(h)$ on $y_t$ in the feasible lag-augmented regression \eqref{eqn:LPestimates} on $(y_t,y_{t-1})$ precisely equals the coefficient on $u_t$ in the desired regression on $(u_t,y_{t-1})$. This argument for why lag-augmented LP can be expected to have a uniformly normal limit distribution even when $\rho \approx 1$ is completely analogous to the reasoning for using lag augmentation in AR inference \citep{Sims1990,Toda1995,Dolado1996,Inoue2002,Inoue2020}. In the LP case, lag augmentation has the additional benefit of simplifying the computation of standard errors, as we now discuss.
\paragraph{Standard Errors.}
We now define the standard errors for the lag-augmented LP estimator. We will show that, contrary to conventional wisdom (e.g., \citealp[p. 166]{Jorda2005}; \citealp[p. 84]{Ramey2016}), HAR standard errors are \emph{not} needed to conduct inference on \emph{lag-augmented} LP, despite the fact that the regression residual $\xi_t(\rho,h)$ is serially correlated. Instead, it suffices to use the usual heteroskedasticity-robust Eicker-Huber-White standard error of $\hat{\beta}(h)$:\footnote{This is computed by the {\tt regress, robust} command in Stata, for example. The usual homoskedastic standard error formula suffices if $u_t$ is assumed to be i.i.d.}
\begin{equation}
\label{eqn:EHW}
\hat{s}(h) \equiv \frac{(\sum_{t=1}^{T-h} \hat{\xi}_t(h)^2 \hat{u}_t(h)^2)^{1/2}}{\sum_{t=1}^{T-h} \hat{u}_t(h)^2},
\end{equation}
where we define the lag-augmented LP residuals
\begin{equation}\label{eqn:estimatedresiduals}
\hat{\xi}_t(h) \equiv y_{t+h} - \hat{\beta}(h)y_t - \hat{\gamma}(h)y_{t-1},\quad t=1,2,\dots,T-h,
\end{equation}
and the residualized regressor of interest
\[\hat{u}_t(h) \equiv y_t - \hat{\rho}(h)y_{t-1},\quad t=1,2,\dots,T-h,\]
\[\hat{\rho}(h) \equiv \frac{\sum_{t=1}^{T-h} y_t y_{t-1}}{\sum_{t=1}^{T-h} y_{t-1}^2}.\]
As mentioned in the introduction, the fact that we may avoid HAR inference simplifies the implementation of LP inference, as there is no need to choose amongst alternative HAR procedures or specify tuning parameters such as bandwidths \citep{Lazarus2018}.
Why is it not necessary to adjust for serial correlation in the residuals? Since lag-augmented LP controls for $y_{t-1}$, equation \eqref{eqn:ytph_decomp2} suggests that the estimator $\hat{\beta}(h)$ is asymptotically equivalent with the coefficient in a linear regression of the (population) residualized outcome $y_{t+h} - \beta(\rho,h+1)y_{t-1}$ on the (population) residualized regressor $u_t = y_t - \rho y_{t-1}$:
\begin{align*}
\hat{\beta}(h) &\approx \frac{\sum_{t=1}^{T-h} \lbrace y_{t+h} - \beta(\rho,h+1)y_{t-1} \rbrace u_t}{\sum_{t=1}^{T-h}u_t^2} \\
&= \beta(\rho,h) + \frac{\sum_{t=1}^{T-h} \xi_t(\rho,h) u_t}{\sum_{t=1}^{T-h}u_t^2}.
\end{align*}
The second term in the decomposition above determines the sampling distribution of the lag-augmented local projection. Although the multi-step regression residual $\xi_t(\rho,h)$ is serially correlated on its own, the \emph{regression score} $\xi_t(\rho,h)u_t$ is serially uncorrelated under \cref{asn:u_mds}.\footnote{\citet{Breitung2019} make this same observation, but they appear to claim that it is sufficient to assume that $\lbrace u_t \rbrace$ is white noise, which is incorrect.} For any $s < t$,
\begin{align}
E[\xi_t(\rho,h)u_t\xi_s(\rho,h)u_s ] &= E[E(\xi_t(\rho,h)u_t\xi_s(\rho,h)u_s \mid u_{s+1},u_{s+2},\dots)] \nonumber \\
&= E[\xi_t(\rho,h) u_t\xi_s(\rho,h)\underbrace{E(u_s \mid u_{s+1},u_{s+2},\dots)}_{=0}]. \label{eqn:scores_uncorr}
\end{align}
Thus, the heteroskedasticity-robust (but not autocorrelation-robust) standard error $\hat{s}(h)$ suffices for doing inference on $\hat{\beta}(h)$.\footnote{\citet[p. 152]{Stock2018} mention a similar conclusion for the distinct case of LP with an instrumental variable, under some conditions on the instrument.} Notice that this result crucially relies on (i) lag-augmenting the local projections and (ii) the strengthening in \cref{asn:u_mds} of the usual martingale difference assumption on $\lbrace u_t \rbrace$ (as remarked above, the strengthening still allows for conditional heteroskedasticity and other plausible features of economic shocks).\footnote{The nuisance coefficient $\hat{\gamma}(h)$ is not interesting \emph{per se}, but note that inference on this coefficient would generally require HAR standard errors, and its limit distribution is in fact non-standard when $\rho \approx 1$.}
Though lag augmentation robustifies and simplifies local projection inference, it is not necessarily a free lunch. We show in \cref{sec:comp} that the relative efficiency of non-augmented and lag-augmented local projection estimators depends on $\rho$ and $h$.
\paragraph{Lag-Augmented Local Projection Inference.}
Define the nominal $100 (1-\alpha)\%$ lag-augmented LP confidence interval for the impulse response at horizon $h$ based on the standard error $\hat{s}(h)$:
\[\hat{C}(h,\alpha) \equiv \left[ \hat{\beta}(h)-z_{1-\alpha/2}\: \hat{s}(h)\:,\: \hat{\beta}(h)+z_{1-\alpha/2}\: \hat{s}(h) \right],\]
where $z_{1-\alpha/2}$ is the $(1-\alpha/2)$ quantile of the standard normal distribution.
Our main result shows that the lag-augmented LP confidence interval above is valid regardless of the persistence of the data, i.e., whether or not the data has a unit root. Crucially, the result does not break down at moderately long horizons $h$. We provide a formal result for VAR($p$) models in \cref{sec:var} and for now just discuss heuristics. Consider any upper bound $\bar{h}_T$ on the horizon which satisfies $\bar{h}_T/T \to 0$. Then \cref{thm:var_lp_inference} below implies that
\begin{equation} \label{eqn:coverage}
\inf_{\rho \in [-1,1]} \inf_{1 \leq h \leq \bar{h}_T} P_{\rho}\left( \beta(\rho,h) \in \hat{C}(h,\alpha) \right ) \to 1-\alpha \quad \text{as } T \to \infty,
\end{equation}
where $P_\rho$ denotes the distribution of the data $\lbrace y_t \rbrace$ under the AR(1) model \eqref{eqn:ar1} with parameter $\rho$. In words, the result states that, for sufficiently large sample sizes, LP inference is valid even under the \emph{worst-case} choices of parameter $\rho \in [-1,1]$ and horizon $h \in [1,\bar{h}_T]$. As is well known, such \emph{uniform} validity is a much stronger result than \emph{pointwise} validity for fixed $\rho$ and $h$. In fact, if we restrict attention to only the stationary region $\rho \in [-1+a,1-a]$, $a \in (0,1)$, then the statement \eqref{eqn:coverage} is true with the upper bound $\bar{h}_T = (1-a)T$ on the horizon. That is, if we know the time series is not close to a unit root, then local projection inference is valid even at long horizons $h$ that are non-negligible fractions of the sample size $T$.
\subsection{Illustrative Simulation Study}
\label{sec:sim_ar1}
We now present a small simulation study to show that lag-augmented LP achieves a favorable trade-off between robustness and efficiency relative to other procedures. For clarity, we continue to assume the simple AR(1) model \eqref{eqn:ar1} with known lag length. Our baseline design considers homoskedastic innovations $u_t \stackrel{i.i.d.}{\sim} N(0,1)$. In \cref{sec:sim_ar1_arch} we present results for ARCH innovations.
We stress that, although we use the AR(1) model for illustration here, the central goal of this paper is to develop a procedure that is feasible even in realistic VAR($p$) models. Thus, we avoid computationally demanding procedures, such as the AR grid bootstrap, which are difficult to implement in applied settings. We provide an extensive theoretical comparison of various inference procedures in \cref{sec:comp}.
\cref{tab:TableMC} displays the coverage and median length of impulse response confidence intervals at various horizons. We consider several versions of AR inference and LP inference, either implemented using the bootstrap or using delta method standard errors. ``LP'' denotes local projection and ``AR'' autoregressive inference. ``LA'' denotes lag augmentation. The subscript ``$b$'' denotes bootstrap confidence intervals constructed from a wild recursive bootstrap design \citep{Goncalves2004}, as described in \cref{sec:boot} (for LP we use the percentile-t confidence interval). Columns without the ``$b$'' subscript use delta method standard errors. For LA-LP, we always use Eicker-Huber-White standard errors as discussed in \cref{sec:ar1_intuition}, whereas non-augmented LP always uses HAR standard errors.\footnote{As an off-the-shelf, state-of-the-art HAR procedure, we choose the Equally Weighted Cosine (EWC) estimator with degrees of freedom as recommended by \citet[equations 4 and 10]{Lazarus2018}. The degrees of freedom depend on the effective sample size $T-h$ and thus differ across horizons $h$.} The column ``AR-LA'' is the Efron bootstrap confidence interval for \emph{lag-augmented} AR estimates developed by \citet{Inoue2020} and discussed further in \cref{sec:comp}.\footnote{We use the \citet{Pope1990} bias-corrected AR estimates to generate the bootstrap samples, as recommended by \citet{Inoue2020}.} All estimation procedures include an intercept. The sample size is $T=240$. We consider data generating processes (DGPs) $\rho \in \lbrace 0,.5,.95,1\rbrace $ and horizons $h$ up to 60 periods (25\% of the sample size, which is not unusual in applied work). The nominal confidence level is 90\%. We use 5,000 Monte Carlo repetitions, with 2,000 bootstrap draws per repetition.
\afterpage{
\begin{landscape}
\begin{table}[p]
\centering
\caption{Monte Carlo results: homoskedastic innovations}
\vspace{0.5\baselineskip}
\begin{tabular}{r|cccccc|cccccc}
& \multicolumn{6}{c|}{Coverage} & \multicolumn{6}{c}{Median length} \\
$h$ & $\text{LP-LA}_b$ & $\text{LP-LA}$ & $\text{LP}_b$ & $\text{LP}$ & $\text{AR-LA}_b$ & $\text{AR}$ & $\text{LP-LA}_b$ & $\text{LP-LA}$ & $\text{LP}_b$ & $\text{LP}$ & $\text{AR-LA}_b$ & $\text{AR}$ \\
\hline
\multicolumn{13}{c}{$\rho = 0.00$} \\
1 & 0.902 & 0.892 & 0.912 & 0.889 & 0.891 & 0.894 & 0.218 & 0.211 & 0.233 & 0.215 & 0.211 & 0.210 \\
6 & 0.908 & 0.899 & 0.916 & 0.898 & 0.000 & 1.000 & 0.219 & 0.214 & 0.233 & 0.220 & 0.000 & 0.000 \\
12 & 0.909 & 0.900 & 0.903 & 0.897 & 0.000 & 1.000 & 0.222 & 0.217 & 0.230 & 0.226 & 0.000 & 0.000 \\
36 & 0.903 & 0.895 & 0.903 & 0.898 & 0.000 & 1.000 & 0.235 & 0.229 & 0.244 & 0.239 & 0.000 & 0.000 \\
60 & 0.898 & 0.886 & 0.894 & 0.889 & 0.000 & 0.979 & 0.252 & 0.244 & 0.261 & 0.255 & 0.000 & 0.000 \\
\multicolumn{13}{c}{$\rho = 0.50$} \\
1 & 0.906 & 0.896 & 0.912 & 0.885 & 0.897 & 0.897 & 0.219 & 0.212 & 0.205 & 0.187 & 0.211 & 0.184 \\
6 & 0.895 & 0.886 & 0.906 & 0.875 & 0.897 & 0.832 & 0.252 & 0.245 & 0.293 & 0.266 & 0.046 & 0.032 \\
12 & 0.906 & 0.894 & 0.903 & 0.887 & 0.897 & 0.766 & 0.255 & 0.248 & 0.293 & 0.280 & 0.002 & 0.001 \\
36 & 0.900 & 0.889 & 0.901 & 0.884 & 0.897 & 0.643 & 0.271 & 0.262 & 0.309 & 0.296 & 0.000 & 0.000 \\
60 & 0.905 & 0.891 & 0.903 & 0.880 & 0.897 & 0.595 & 0.291 & 0.279 & 0.333 & 0.316 & 0.000 & 0.000 \\
\multicolumn{13}{c}{$\rho = 0.95$} \\
1 & 0.892 & 0.878 & 0.842 & 0.827 & 0.882 & 0.850 & 0.220 & 0.212 & 0.076 & 0.072 & 0.212 & 0.075 \\
6 & 0.903 & 0.838 & 0.851 & 0.789 & 0.882 & 0.810 & 0.523 & 0.452 & 0.395 & 0.345 & 1.011 & 0.318 \\
12 & 0.889 & 0.806 & 0.853 & 0.752 & 0.882 & 0.769 & 0.678 & 0.550 & 0.644 & 0.518 & 1.744 & 0.430 \\
36 & 0.885 & 0.814 & 0.865 & 0.674 & 0.882 & 0.656 & 0.728 & 0.625 & 0.859 & 0.612 & 6.567 & 0.272 \\
60 & 0.892 & 0.833 & 0.892 & 0.693 & 0.882 & 0.595 & 0.731 & 0.651 & 0.942 & 0.641 & 23.050 & 0.095 \\
\multicolumn{13}{c}{$\rho = 1.00$} \\
1 & 0.895 & 0.874 & 0.820 & 0.554 & 0.877 & 0.532 & 0.219 & 0.211 & 0.040 & 0.040 & 0.210 & 0.039 \\
6 & 0.875 & 0.777 & 0.836 & 0.503 & 0.877 & 0.494 & 0.564 & 0.498 & 0.243 & 0.222 & 1.206 & 0.214 \\
12 & 0.843 & 0.676 & 0.827 & 0.429 & 0.877 & 0.459 & 0.821 & 0.671 & 0.477 & 0.385 & 2.553 & 0.379 \\
36 & 0.741 & 0.428 & 0.755 & 0.200 & 0.877 & 0.348 & 1.338 & 0.950 & 1.200 & 0.592 & 21.107 & 0.670 \\
60 & 0.642 & 0.276 & 0.712 & 0.156 & 0.877 & 0.295 & 1.434 & 0.978 & 1.667 & 0.637 & 161.250 & 0.731 \\
\end{tabular}
\label{tab:TableMC}
\\
\vspace{0.5\baselineskip}
\begin{minipage}{1.15\textwidth}
{\footnotesize Coverage probability and median length of nominal 90\% confidence intervals at different horizons. AR(1) model with $\rho \in \lbrace 0,.5,.95,1\rbrace $, $T=240$, i.i.d. standard normal innovations. 5,000 Monte Carlo repetitions; 2,000 bootstrap iterations.}
\end{minipage}
\end{table}
\end{landscape}
}
Consistent with our theoretical results, the bootstrap version of lag-augmented local projection (column 1) achieves coverage close to the nominal level in almost all cases, whereas the competing procedures either under-cover or return impractically wide confidence intervals. In contrast, non-augmented LP (columns 3 and 4) exhibits larger coverage distortions in almost all cases. As is well known, textbook AR delta method confidence intervals (column 6) severely under-cover when $\rho>0$ and the horizon is even moderately large.
It is only when both $\rho=1$ and $h \geq 36$ that lag-augmented local projection exhibits serious coverage distortions, again consistent with our theory. However, even in these cases, the coverage distortions are similar to or less pronounced than those for non-augmented LP and for delta method AR inference.
Although the \citet{Inoue2020} lag-augmented AR bootstrap confidence interval (column 5) achieves correct coverage for $\rho>0$ at all horizons, this interval is extremely wide in the problematic cases where $\rho$ is close to 1 and the horizon $h$ is intermediate or long. We explain this fact theoretically in \cref{sec:comp}. Confidence intervals with median width greater than 1 would appear to be of little practical use, since the true impulse response parameter is bounded above by 1 in the AR(1) model.\footnote{In the AR(1) model, we could intersect all confidence intervals with the interval $[-1,1]$. In this case, the median length of the \citet{Inoue2020} confidence interval is close to 1, cf. \cref{sec:comp_details_laarboot}.} Note also that the \citet{Inoue2020} interval severely under-covers when $\rho=0$ at all even (but not odd) horizons $h$, as explained theoretically in \cref{sec:comp}.\footnote{\citet{Inoue2020} assume $\rho \neq 0$ and discuss why this restriction is necessary in their case.}
Although outperformed by bootstrap procedures, the lag-augmented local projection delta method interval (column 2) performs well among the group of delta method procedures. Its coverage distortions are much less severe than textbook AR delta method inference (column 4) and non-augmented LP inference with HAR standard errors (column 6). Recall that the lag-augmented LP confidence interval is at least as easy to compute as these other delta method confidence intervals. The reason why the bootstrap improves on the coverage properties of the delta method procedures is related to the well-known finite-sample bias of AR and LP estimators \citep{Kilian1998, Herbst2020}.\footnote{Our bootstrap implementation of \emph{non-augmented} LP also appears to be quite effective at correcting the most severe coverage distortions of the delta method procedure.}
\cref{tab:TableMC} illustrates the fact that the robustness of lag-augmented local projection inference entails an efficiency loss relative to AR inference when $\rho$ is well below 1, although this loss is not large in absolute terms. In \emph{percentage} terms, local projection confidence intervals are much wider than AR-based confidence intervals when $\rho \ll 1$ and the horizon $h$ is intermediate or long, since AR procedures mechanically impose that the impulse response function tends to 0 geometrically fast with the horizon. Yet, in \emph{absolute} terms, the median length of the LP confidence intervals is not so large as to be a major impediment to applied research. The relative efficiency of lag-augmented LP vs. non-augmented LP cannot be ranked and depends on the DGP and on the horizon. When $\rho$ is close to 1, lag-augmented LP intervals are sometimes (much) narrower than lag-augmented AR intervals. We analytically characterize the various efficiency trade-offs in \cref{sec:comp_details_releff}.
\cref{sec:sim_further} shows that the preceding qualitative conclusions extend to richer models. There we consider a bivariate VAR(4) model with varying degrees of persistence, as well as two empirically calibrated VAR(12) models with four or five observables.
\section{Comparison With Other Inference Procedures}
\label{sec:comp}
The simulations and theoretical results in this paper suggest that lag-augmented local projection is the only known confidence interval procedure that achieves uniformly valid coverage over the DGP and over a wide range of horizons, while preserving reasonable average length and remaining computationally feasible in realistic settings. However, the simulations also suggest that lag-augmented local projection inference is less efficient than standard AR inference when the data is stationary. In this section we discuss in more detail the coverage and length properties of alternative confidence interval procedures for impulse responses. We review the well-known drawbacks of textbook AR inference, provide new results on the relative length of lag-augmented LP vs. non-augmented LP and lag-augmented AR, and discuss the computational challenges of the AR grid bootstrap. We refer the reader back to the small-scale simulation study in \cref{sec:sim_ar1} for illustrations of the following arguments.
\paragraph{Textbook Autoregressive Inference.}
The uniformity result (\ref{eqn:coverage}) for lag-augmented LP stands in stark contrast to textbook AR inference on impulse responses, which suffers from several well-known issues. First, for the standard OLS AR estimator, the usual asymptotic normal limiting theory is invalid when the derivative of the impulse response parameter with respect to the AR coefficients has a singular Jacobian matrix. In the AR(1) model, this occurs at all horizons $h\geq 2$ in the white noise case $\rho=0$ \citep{Benkwitz2000}. Second, as with non-augmented LP, textbook AR inference is not uniformly valid when the data is nearly non-stationary, unless one further restricts the parameter space (\citealp[Remark 2.5]{Phillips1998}; \citealp{Inoue2002}; \citealp[Remark 3, p. 455]{Inoue2020}).\footnote{This is well known in the AR(1) model. In the AR(2) model, a non-normal limit arises at $h=2$ when there is a unit root and the autoregressive coefficients are equal \citep[Remark 3, p. 455]{Inoue2020}.} Third, pre-testing for the presence of a unit root does not yield uniformly valid inference and can lead to poor finite sample performance \citep[e.g.,][p. 1412]{Mikusheva2007}. Fourth, plug-in AR inference with normal critical values must necessarily break down at medium-long horizons $h=h_T \propto T^{1/2}$ and at long horizons $h_T \propto T$, due to the severe nonlinearity of the impulse response transformation at such horizons \citep[Section 4.3]{Mikusheva2012}. \citet{Wright2000} and \citet{Pesavento2006,Pesavento2007} construct confidence intervals for persistent processes at long horizons $h=h_T \propto T$ by inverting the non-standard AR limit distribution, but these tailored procedures do not work uniformly over the parameter space or over the horizon.
The severe under-coverage of delta method AR inference is starkly illustrated in \cref{sec:sim_ar1} (see column 6 of \cref{tab:TableMC}). As discussed in detail by \cite{Inoue2020}, standard bootstrap approaches to AR inference do not solve all the uniformity issues.
We must emphasize, however, that if we restrict attention to stationary processes and short-horizon impulse responses, the standard OLS AR impulse response estimator is more efficient than lag-augmented LP. Hence, there is a trade-off between efficiency in benign settings and robustness to persistence and longer horizons, as is also clear in the simulation results in \cref{sec:sim_ar1}. We expand upon the efficiency properties of the standard AR estimator in \cref{sec:comp_details_releff}.
\paragraph{Lag-Augmented AR Inference.}
The above-mentioned non-uniformity of the textbook AR inference method in the case of near-non-stationary data can be remedied by lag augmentation \citep{Inoue2020}. In the case of an AR(1) model, the lag-augmented AR estimator $\hat{\beta}_\text{ARLA}(h)$ is given by $\hat{\rho}_1^h$, where $(\hat{\rho}_1,\hat{\rho}_2)$ are the OLS coefficients from a regression of $y_t$ on $(y_{t-1},y_{t-2})$ (i.e., we estimate an AR(2) model). The intuition why this guarantees a normal limiting distribution even in the unit root case is the same as in \cref{sec:ar1_intuition}. Lag-augmented AR and lag-augmented LP coincide at horizon $h=1$, but not at longer horizons. Lag augmentation involves a loss of efficiency: The lag-augmented AR estimator is strictly less efficient than the non-augmented AR estimator except when the true process is white noise (see \cref{sec:comp_details_releff}). Note that lag augmentation by itself does not solve the above-mentioned issues that occur when the Jacobian of the impulse response transformation is singular, or when doing inference at medium-long or long horizons.\footnote{The AR(1) simulations in \cref{sec:sim_ar1} show that the coverage of the \citet{Inoue2020} confidence interval is 0 at all \emph{even} horizons when $\rho=0$. This is because the true impulse response is 0, but the bootstrap samples of $\hat{\rho}_1^h$ are all strictly positive. Their procedure achieves uniformly correct coverage at \emph{odd} horizons.}
The bootstrap confidence interval for lag-augmented AR proposed by \citet{Inoue2020} has valid coverage even at long horizons. Specifically, \citet{Inoue2020} show that the \emph{Efron} bootstrap confidence interval---applied to recursive AR bootstrap samples of $\hat{\beta}_\text{ARLA}(h)$---has valid coverage even at long horizons $h=h_T \propto T$, as long as the largest autoregressive root is bounded away from 0.\footnote{For intuition, consider the AR(1) case. The Efron bootstrap preserves monotonic transformations, and the bootstrap transformation $\beta(\rho,h)=\rho^h$ is monotonic (if we restrict attention to $\rho \in (0,1]$ or $\rho \in [-1,0)$). Hence, the Efron confidence interval is valid for $\rho^h$ if it is valid for $\rho$ itself. In more general VAR($p$) models, the same argument can be applied at long horizons, since here only the largest autoregressive root matters for impulse responses (if the roots are well-separated).}
Unfortunately, we show in \cref{sec:comp_details_laarboot} that the expected length of the lag-augmented AR interval is prohibitively large when the data is persistent and the horizon is long. Precisely, in the case of an AR(1) model, $\hat{\beta}_\text{ARLA}(h)=\hat{\rho}_1^h$ is \emph{inconsistent} for sequences of DGPs $\rho=\rho_T$ and horizons $h=h_T$ such that $h_{T} \propto T^{\eta}, \eta \in [1/2,1]$, and $h_T(1-\rho_T) \to a \in [0,\infty)$. The reason is that the lag-augmented coefficient estimator $\hat{\rho}_1$ converges at rate $T^{-1/2}$ even in the unit root case, implying that the estimation error in $\hat{\rho}_1$ is not negligible when raising the estimator to a power of $h=h_T$. This implies that the Efron bootstrap confidence interval is inconsistent (i.e., its length does not shrink to 0 in probability) for such sequences $\rho_T$ and $h_T$. In fact, when $\eta>1/2$, the width of the confidence interval for the $h_{T}$ impulse response is almost equal to the entire positive part of the parameter space $[0,1]$ with probability equal to the nominal confidence level. This contrasts with the lag-augmented LP confidence interval, which is consistent for any sequence $\rho_T \in [-1,1]$ and any sequence $h_T$ such that $h_T/T\to 0$. The large width of the \cite{Inoue2020} interval is illustrated in the simulations in \cref{sec:sim_ar1} (see the second-to-last column in \cref{tab:TableMC}).
Interestingly, \emph{if we restrict attention to stationary processes and short horizons, the relative efficiency of lag-augmented AR and lag-augmented LP inference is ambiguous}. In the context of a stationary, homoskedastic AR(1) model with a fixed horizon $h$ of interest, \cref{fig:se_indiff} shows that lag-augmented AR is more efficient than lag-augmented LP when $\rho$ is small or when the horizon $h$ is large, and vice versa. For any horizon $h$, there exists some cut-off value for $\rho \in (0,1)$, above which lag-augmented LP is more efficient. Intuitively, the nonlinear impulse response transformation $\rho \mapsto \rho^h$ is highly sensitive to values of $\rho$ near 1 whenever $h$ is large, which compounds the effects of estimation error in $\hat{\rho}$, whereas LP is a purely linear procedure.
\begin{figure}[t]
\centering
\includegraphics[width=0.8\linewidth]{se_indiff.eps}
\caption{Efficiency ranking of three different estimators of the fixed impulse response $\beta(\rho,h) = \rho^h$ in the homoskedastic AR(1) model: lag-augmented LP ($\text{LP}_\text{LA}$), non-augmented LP ($\text{LP}_\text{NA}$), and lag-augmented AR ($\text{AR}_\text{LA}$). Gray area: combinations of $(|\rho|,h)$ for which $\text{LP}_\text{LA}$ is more efficient than $\text{AR}_\text{LA}$. Thatched area: $\text{LP}_\text{LA}$ is more efficient than $\text{LP}_\text{NA}$. See \cref{sec:comp_details_releff} for analytical derivations of the indifference curves (thick lines).} \label{fig:se_indiff}
\end{figure}
\paragraph{AR Grid Bootstrap and Projection.}
The grid bootstrap of \citet{Hansen1999} represents a computationally intensive approach to doing valid inference at fixed and long horizons, regardless of persistence, but it is invalid at intermediate horizons, as shown by \citet{Mikusheva2012}. The grid bootstrap is based on test inversion, so it requires running an autoregressive bootstrap on each point in a fine grid of potential values for the impulse response parameter of interest. It also requires estimating a constrained OLS estimator that imposes the hypothesized null on the impulse response at each point in the grid. Recall that lag-augmented LP inference is computationally simple and valid at any horizon $h=h_T$ satisfying $h_T/T \to 0$. However, in the case of unit roots and long horizons $h_T \propto T$, lag-augmented LP inference with normal critical values is not valid, while the grid bootstrap is valid \citep{Mikusheva2012}.
Another computationally intensive approach is to form a uniformly valid confidence set for the AR parameters and then map it into a confidence interval for impulse responses by projection. Although doable in the AR(1) model, this approach would appear to be computationally infeasible and possibly highly conservative in realistic VAR($p$) settings, unlike lag-augmented LP (see \cref{sec:var}).
\paragraph{Other Local Projection Approaches.}
Non-augmented LP is not robust to non-stationarity, as already discussed in \cref{sec:ar1_intuition}. \emph{If the data is stationary and the horizon $h$ is fixed, the relative efficiency of non-augmented LP and lag-augmented LP is generally ambiguous}, as shown in \cref{fig:se_indiff} in the case of a homoskedastic AR(1) model. There are two competing forces. On the one hand, as shown in \cref{sec:ar1_intuition}, non-augmented LP uses the regressor $y_t$, which has higher variance than the effective regressor $u_t$ in the lag-augmented case. By itself, this suggests that non-augmented LP should be more efficient. On the other hand, absent lag augmentation, the LP regression scores are serially correlated and thus have a larger long-run variance. On balance, \cref{sec:comp_details_releff} shows that lag-augmented LP is relatively more efficient the smaller is $\rho$ and the larger is $h$.
In some empirical settings, the researcher may directly observe the autoregressive innovation, or some component of the innovation, for example by constructing narrative measures of economic shocks \citep{Ramey2016}. For concreteness, consider the AR(1) model \eqref{eqn:ar1} and assume we observe the innovation $u_t$. In this case, it is common in empirical practice to simply regress $y_{t+h}$ on $u_t$, without controls. Although this strategy provides consistent impulse response estimates when the data is stationary, it is inefficient relative to lag-augmented LP, since the latter approach additionally controls for the variable $y_{t-1}$, which would otherwise show up in the error term in the representation \eqref{eqn:ytph_decomp2}. Thus, lag augmentation is desirable on robustness and efficiency grounds even if some shocks are directly observed.
\paragraph{Summary.}
Existing and new theoretical results confirm the main message of our simulations in \cref{sec:sim_ar1}: Lag-augmented LP is the only known procedure that is computationally feasible in realistic problems and can be shown to have valid coverage under a wide range of DGPs and horizon lengths, without achieving such valid coverage by returning a confidence interval that is impractically wide. This robustness does come at the cost of a loss of efficiency relative to non-robust AR methods. However, the efficiency loss is large in \emph{relative} terms only in stationary, short-horizon cases, where lag-augmented LP confidence intervals do well in \emph{absolute} terms, as illustrated in \cref{sec:sim_ar1}. Based on these results, we believe that it is only in the case of highly persistent data and very long horizons $h=h_T \propto T$ that the use of alternative robust procedures should be considered, such as the computationally demanding AR grid bootstrap.
\section{General Theory for the VAR(\texorpdfstring{$p$}{p}) Model}
\label{sec:var}
This section presents the inference procedure and theoretical uniformity result for a general VAR($p$) model. In this case, the lag-augmented LP procedure controls for $p$ lags of all the time series that enter into the VAR model. We follow \citet{Mikusheva2012} and \citet{Inoue2020} in assuming that the lag length $p$ is finite and known. We also assume that the VAR process has no deterministic dynamics for simplicity. See \cref{sec:conc} for further discussion of these assumptions.
\subsection{Model and Inference Procedure} \label{sec:var_model}
Consider an $n$-dimensional VAR($p$) model for the data $y_t=(y_{1,t},\dots,y_{n,t})'$:
\begin{equation} \label{eqn:var}
y_{t} = \sum_{\ell=1}^p A_\ell y_{t-\ell} + u_{t}, \quad t=1,2,\dots,T,\quad y_0=\dots=y_{1-p}=0,
\end{equation}
Let $A \equiv (A_1, \ldots, A_{p})$ denote the $n \times np$ matrix collecting all the autoregressive coefficients. The assumption of zero pre-sample initial conditions $y_0=\dots=y_{1-p}=0$ is made for notational simplicity and can be relaxed, as discussed below in the remarks after \cref{thm:var_lp_inference}. As in the AR(1) case, we assume that the $n$-dimensional innovation process $\lbrace u_{t} \rbrace$ satisfies the strengthening of the martingale difference condition in \cref{asn:u_mds} (which from now on will refer to the vector process $\lbrace u_t \rbrace$).
We seek to do inference on a scalar function of the reduced-form impulse responses of the VAR model. Generalizations to \emph{structural} impulse responses and \emph{joint} inference require more notation but are otherwise straight-forward, see \cref{sec:conc}. Let $\beta_{i}(A,h)$ denote the $n \times 1$ vector containing each of variable $i$'s reduced-form impulse responses at horizon $h \geq 0$. Without loss of generality, we focus on the impulse responses of the first variable $y_{1,t}$. Thus, we seek a confidence interval for the scalar parameter $\nu'\beta_{1}(A,h)$, where $\nu \in \mathbb{R}^n \backslash \lbrace 0 \rbrace$ is a user-specified vector. For example, the choice $\nu = e_j$ (the $j$-th unit vector) selects the horizon-$h$ response of $y_{1,t}$ with respect to the $j$-th reduced-form innovation $u_{j,t}$.
Local projection estimators of impulse responses are motivated by the representation
\begin{equation} \label{eqn:var_long_regression}
y_{1,t+h} = \beta_1(A,h)' y_{t} + \sum_{\ell=1}^{p-1} \delta_{1,\ell}(A,h)^{\prime}y_{t-\ell} + \xi_{1,t}(A,h),
\end{equation}
see \citet{Jorda2005} and \citet[Chapter 12.8]{Kilian2017}. Here $\delta_{1,\ell}(A,h)$ is an $n \times 1$ vector of regression coefficients that can be obtained by iterating on the VAR model \eqref{eqn:var}. The model-implied multi-step forecast error in this regression is
\begin{equation}
\xi_{1,t}(A,h) \equiv \sum_{\ell = 1}^{h} \beta_{1}(A,h-\ell)' u_{t+\ell}.
\end{equation}
\paragraph{Multivariate Lag-Augmented Local Projection.}
The lag-augmented LP estimator corresponding to the VAR model \eqref{eqn:var} is motivated by \eqref{eqn:var_long_regression}. We regress $y_{1,t+h}$ on the $n$ variables $y_{t}$, using the $np$ variables $(y'_{t-1}, \ldots, y'_{t-p})$ as additional controls. According to equation \eqref{eqn:var_long_regression}, the population regression coefficients on the last $n$ control variables $y_{t-p}$ equal zero. Thus, we are including one additional lag in the estimation of the impulse response coefficients. Given any horizon $h \in \mathbb{N}$, the lag-augmented LP estimator $\hat{\beta}_1(h)$ of $\beta_1(A,h)$ is given by the vector of coefficients on $y_{t}$ in the regression of $y_{1,t+h}$ on $x_{t} \equiv (y_{t}', y_{t-1}', \ldots, y_{t-p}')'$:
\begin{equation}
\begin{pmatrix} \label{eqn:LPestimates_VAR}
\hat{\beta}_1(h) \\
\hat{\gamma}_1(h)
\end{pmatrix} \equiv \left(\sum_{t=1}^{T-h} x_tx_t'\right)^{-1}\sum_{t=1}^{T-h}x_t y_{1,t+h},
\end{equation}
where $\hat{\beta}_{1}(h)$ is a vector of dimension $n \times 1$.
The usual (Eicker-Huber-White) heteroskedasticity-robust standard error for $\nu'\hat{\beta}_1(h)$ is defined as
\[\hat{s}_1(h,\nu) \equiv \frac{1}{T-h}\left\lbrace \nu'\hat{\Sigma}(h)^{-1} \left(\sum_{t=1}^{T-h} \hat{\xi}_{1,t}(h)^2\hat{u}_t(h)\hat{u}_t(h)' \right) \hat{\Sigma}(h)^{-1}\nu\right\rbrace^{1/2},\]
where
\[\hat{\xi}_{1,t}(h) \equiv y_{1,t+h} - \hat{\beta}_1(h)' y_{t} - \hat{\gamma}_{1}(h)' X_{t}, \quad X_{t} \equiv (y_{t-1}' , \ldots, y_{t-p}')', \]
\[ \hat{u}_{t}(h) \equiv y_{t} -\hat{A}(h)X_{t}, \quad \hat{A}(h) \equiv \left( \sum_{t=1}^{T-h} y_{t}X_{t}' \right) \left ( \sum_{t=1}^{T-h} X_{t} X_{t}' \right)^{-1}, \]
and
\[\hat{\Sigma}(h) \equiv \frac{1}{T-h}\sum_{t=1}^{T-h} \hat{u}_t(h)\hat{u}_t(h)'.\]
The $1-\alpha$ confidence interval for $\nu'\beta_1(A,h)$ is defined as
\[\hat{C}_1(h,\nu,\alpha) \equiv \left[ \nu'\hat{\beta}_1(h)-z_{1-\alpha/2}\: \hat{s}_1(h,\nu)\:,\: \nu'\hat{\beta}_1(h)+z_{1-\alpha/2}\: \hat{s}_1(h,\nu) \right].\]
\paragraph{Parameter Space.}
We consider a class of VAR processes with possibly multiple unit roots combined with arbitrary stationary dynamics. Specifically, we will prove that the confidence interval $\hat{C}_1(h,\nu,\alpha)$ has uniformly valid coverage over the following parameter space. Let $\|M\| \equiv \sqrt{\operatorname*{trace}(M'M)}$ denote the Frobenius matrix norm, and let $I_n$ denote the $n \times n$ identity matrix.
\begin{defn}[VAR parameter space] \label{dfn:param_space}
Given constants $a \in [0,1)$, $C>0$, and $\epsilon \in (0,1)$, let $\mathcal{A}(a,C,\epsilon)$ denote the space of autoregressive coefficients $A=(A_1, \ldots, A_{p})$ such that the associated $p$-dimensional lag-polynomial $A(L)=I_n-\sum_{\ell=1}^p A_\ell L^\ell$ admits the factorization
\begin{equation}
A(L) = B(L) (I_n - \operatorname*{diag}(\rho_1, \ldots, \rho_n) L),
\end{equation}
where $\rho_i \in [a-1,1-a]$ for all $i=1,\dots,n$, and $B(L)$ is a lag polynomial of order $p-1$ with companion matrix $\mathbf{B}$ satisfying $\| \mathbf{B}^\ell \| \leq C (1-\epsilon)^\ell$ for all $\ell = 1,2,\dots$.\footnote{See \cref{sec:var_proof_main} for the standard definition of a companion matrix.}
\end{defn}
This parameter space contains any stationary VAR process (for sufficiently small $a,\epsilon$ and sufficiently large $C$) as well as many---but not all---non-stationary processes. Lag polynomials $A(L)$ in this parameter space imply that the process $\lbrace y_t \rbrace$ can be written in the form $y_t = \operatorname*{diag}(\rho_1,\dots,\rho_n)y_{t-1} + \tilde{y}_t$, where $\tilde{y}_t \equiv B(L)^{-1}u_t$ is a stationary process whose impulse responses at horizon $\ell$ decay at the geometric rate $(1-\epsilon)^\ell$. We allow all the roots $\rho_1,\dots,\rho_n$ to be potentially close to or equal to 1. \citet[Section 4.2]{Mikusheva2012} considers the same class of processes but with $\rho_2=\dots=\rho_n=0$. We are not aware of other uniform inference results that allow multiple near-unit roots. Although the parameter space in \cref{dfn:param_space} appears more restrictive than the local-to-unity framework of \citet[Eqn. 2]{Phillips1988}, we argue below that our uniform coverage result applied to the parameter space $\mathcal{A}(a,C,\epsilon)$ immediately implies an extended result that also covers processes with cointegration among the control variables $y_{2,t},\dots,y_{n,t}$. However, we do impose the restriction that the response variable of interest $y_{1,t}$ has at most one root near unity, as in \citet{Wright2000}, \citet{Pesavento2006}, \citet{Mikusheva2012}, and \citet{Inoue2020}.
\subsection{Additional Assumptions}
\label{sec:var_asn}
Our main result requires two further technical assumptions in addition to \cref{asn:u_mds}. Let $\lambda_{\min}(M)$ denote the smallest eigenvalue of a symmetric positive semidefinite matrix $M$.
\begin{asn} \label{asn:var_u_reg}
\leavevmode
\begin{enumerate}[i)]
\item \label{itm:var_asn_u_bounds} $E(\|u_t\|^8)<\infty$, and there exists $\delta > 0$ such that $\lambda_{\min}(E[u_tu_t' \mid \lbrace u_s \rbrace_{s<t}]) \geq \delta$ almost surely.
\item \label{itm:var_asn_u2_cum} The process $\lbrace u_t \otimes u_t \rbrace$ has absolutely summable cumulants up to order 4.
\end{enumerate}
\end{asn}
\noindent Part (\ref{itm:var_asn_u_bounds}) of \cref{asn:var_u_reg} is a common requirement for consistent estimation of regression standard errors with possibly heteroskedastic residuals. Part (\ref{itm:var_asn_u2_cum}) is a standard weak dependence restriction on the second moments of $u_t$ \citep[Chapter 2.6]{Brillinger2001}.
We will write $\rho(A)=(\rho_1(A),\dots,\rho_n(A))'$ to represent any of the possible vectors of roots $\rho_1,\dots,\rho_n$ corresponding to a collection of autoregressive coefficients $A=(A_1, \ldots, A_{p}) \in \mathcal{A}(0,C,\epsilon)$. This is a slight abuse of notation, since the mapping from $A(L)$ to $\rho_i$'s is one-to-many. Define $g(\rho,h)^2 \equiv \min\lbrace \frac{1}{1-|\rho|},h\rbrace$ and $\rho_i^*(A,\epsilon) \equiv \max\lbrace |\rho_{i}(A)|, 1-\epsilon/2 \rbrace$. Define also the $np \times np$ diagonal matrix $G(A,h,\epsilon) \equiv I_p \otimes \operatorname*{diag}(g(\rho_1^*(A,\epsilon) ,h),\dots,g(\rho_n^*(A,\epsilon),h))$.
\begin{asn} \label{asn:var_XpX}
For any $C>0$ and $\epsilon \in (0,1)$,
\[\lim_{K \to \infty} \lim_{T\to\infty} \inf_{A \in \mathcal{A}(0,C,\epsilon)} P_A\left( \lambda_{\min}\left(G(A,T,\epsilon)^{-1}\left[\frac{1}{T}\sum_{t=1}^T X_tX_t'\right]G(A,T,\epsilon)^{-1} \right) \geq 1/K \right) = 1.\]
\end{asn}
\noindent This high-level assumption ensures that the properly scaled (matrix) ``denominator'' in the VAR OLS estimator $\hat{A}(h)$ is uniformly non-singular asymptotically, so the estimator is uniformly well-defined with high probability in the limit. Hence, the assumption is essentially necessary for our result.
How can \cref{asn:var_XpX} be verified? $G(A_T,T,\epsilon)^{-1}\left[\frac{1}{T}\sum_{t=1}^T X_tX_t'\right]G(A_T,T,\epsilon)^{-1}$ is known to converge in distribution in a \emph{pointwise} sense to an almost surely positive definite (perhaps stochastically degenerate) random matrix under stationary, local-to-unity, or unit root sequences $\lbrace A_T \rbrace$ \citep[e.g.,][]{Phillips1988,Hamilton1994}.\footnote{Note that the diagonal entries of $G(A,T,\epsilon)^{-1}$ are constants for stationary VAR coefficient matrices $A$, whereas these diagonal entries are proportional to $T^{-1/2}$ under local-to-unity or unit root sequences.} \cref{asn:var_XpX} requires that such convergence obtains for \emph{all} possible sequences $\lbrace A_T \rbrace$. In \cref{sec:var_XpX_ar1} we illustrate how the assumption can be verified in the AR(1) model under an additional weak condition on the innovation process.
\subsection{Main Result}
We now state the result that the LP estimator $\nu'\hat{\beta}_1(h)$ is asymptotically normally distributed uniformly over the parameter space in \cref{dfn:param_space}, even at long horizons $h$. Let $P_A$ denote the probability measure of the data $\lbrace y_t \rbrace$ when it is generated by the VAR($p$) model \eqref{eqn:var} with coefficients $A \in \mathcal{A}(a,C,\epsilon)$. The distribution of the innovations $\lbrace u_t \rbrace$ is fixed.
\begin{prop} \label{thm:var_lp_inference}
Let \cref{asn:u_mds,asn:var_u_reg,asn:var_XpX} hold. Let $C>0$ and $\epsilon \in (0,1)$.
\begin{enumerate}[i)]
\item \label{itm:var_lp_inference_stat} Let $a \in (0,1)$. For all $x \in \mathbb{R}$,
\[\sup_{A \in \mathcal{A}(a,C,\epsilon)} \sup_{1 \leq h \leq (1-a)T} \left| P_A\left(\frac{\nu'[\hat{\beta}_1(h)-\beta_1(A,h)]}{\hat{s}_1(h,\nu)} \leq x\right) - \Phi(x) \right| \to 0.\]
\item \label{itm:var_lp_inference_all} Consider any sequence $\lbrace \bar{h}_T \rbrace$ of nonnegative integers such that $\bar{h}_T < T$ for all $T$ and $\bar{h}_T/T \to 0$. Then for all $x \in \mathbb{R}$,
\[\sup_{A \in \mathcal{A}(0,C,\epsilon)} \sup_{1 \leq h \leq \bar{h}_T} \left| P_A\left(\frac{\nu'[\hat{\beta}_1(h)-\beta_1(A,h)]}{\hat{s}_1(h,\nu)} \leq x\right) - \Phi(x) \right| \to 0.\]
\end{enumerate}
\end{prop}
\begin{proof}
See \cref{sec:var_proof_main}.
\end{proof}
The uniform asymptotic normality established above immediately implies that the confidence interval $\hat{C}_1(h,\nu,\alpha)$ has uniformly valid coverage asymptotically. Part (\ref{itm:var_lp_inference_stat}) considers stationary VAR processes whose largest roots are bounded away from 1; then inference is valid even at long horizons $h=h_T \propto T$. Part (\ref{itm:var_lp_inference_all}) allows all or some of the $n$ roots $\rho_1,\dots,\rho_n$ to be near or equal to 1, but then we require $h_T/T \to 0$.
\paragraph{Remarks.}
\begin{enumerate}
\item The proof of \cref{thm:var_lp_inference} shows that the uniform convergence rate of $\hat{\beta}_1(h_T)$ is $O_p((h_T/T)^{1/2})$ if $h_T/T \to 0$. This rate may be slower than that of the possibly super-consistent non-augmented LP estimator, which is the price to pay for uniformity. If we restrict attention to the stationary parameter space $\mathcal{A}(a,C,\epsilon)$, $a>0$, the convergence rate of $\hat{\beta}_1(h_T)$ is $O_p(T^{-1/2})$ provided that $h_T \leq (1-a)T$.
\item There are three main challenges in establishing the uniform validity of local projection inference.
\begin{enumerate}[a)]
\item The variance of the regression residual $\xi_{1,t}(A,h)$ is increasing in the horizon $h$ and also depends on $A$. Thus, the simplest laws of large numbers and central limit theorems for stationary processes do not apply. We instead apply a central limit theorem for martingale difference sequences and derive uniform bounds on moments of relevant variables. The central limit theorem is delicate, since the regression scores $\xi_{1,t}(A,h)u_t$ are not a martingale difference sequence with respect to the natural filtration generated by past $u_t$'s. However, it is possible to ``reverse time'' in a way that makes the scores a martingale difference sequence with respect to an alternative filtration, see the proof of the auxiliary \cref{thm:var_clt}.
\item To handle both unit roots, stationary processes, and everything in between, we must consider various kinds of sequences of drifting parameters $A=A_T$, following the general logic of \citet{Andrews2019}. This is primarily an issue when showing consistency of the standard error $\hat{s}_1(h,\nu)$, which requires deriving the convergence rates of the various estimators along drifting parameter sequences. We do this by explicit calculation of moment bounds that are uniform in the both the DGP and the horizon.
\item Our proof requires bounds on the rate of decay of impulse response functions that are uniform in both the DGP and the horizon. Though the AR(1) case is trivial due to the monotonically decreasing exponential functional form $\beta(\rho,h)=\rho^h$, the bounds for the general VAR($p$) case require more work, see especially \cref{thm:var_bound_for_IRFs_A} in \cref{sec:var_se_proof}. These results may be of independent interest.
\end{enumerate}
\item \cref{thm:var_lp_inference} does not cover the case where $h \propto T$ and some of the roots $\rho_i$ are local-to-unity or equal to unity. Simulation evidence and analytical calculations along the lines of \citet{Hjalmarsson2020} strongly suggest that even in the AR(1) model the asymptotic normality of lag-augmented local projections does \emph{not} go through when $\rho=1$ and $h = \kappa T$ for $\kappa \in (0,1)$. Indeed, in this case the sample variance of the regression scores $\xi_t(\rho,h)u_t$ appears not to converge in probability to a constant, thus violating the conclusion of the key auxiliary \cref{thm:var_se_infeas} below. As discussed in \cref{sec:comp}, the behavior of plug-in autoregressive impulse response estimators is also non-standard when $\rho \approx 1$ and $h \propto T$.
\item A corollary of our main result is that we can allow for cointegrating relationships to exist among the control variables $y_{2,t},\dots,y_{n,t}$. This is because both the LP estimator and the reduced-form impulse responses are equivariant with respect to non-singular linear transformations of these $n-1$ variables. For example, consider a 3-dimensional process $(y_{1,t},y_{2,t},y_{3,t})$ that follows a VAR model in the parameter space in \cref{dfn:param_space} with $\rho_2=1,\rho_3=0$. Now consider the transformed process $(y_{1,t},\tilde{y}_{2,t},\tilde{y}_{3,t}) = (y_{1,t}, y_{2,t} + y_{3,t}, -y_{2,t} + y_{3,t})$. The variables $\tilde{y}_{2,t}$ and $\tilde{y}_{3,t}$ are cointegrated with cointegrating vector $(1,1)'$. Since $(\tilde{y}_{2,t},\tilde{y}_{3,t})$ is a non-singular linear transformation of $(y_{2,t},y_{3,t})$, the conclusions of \cref{thm:var_lp_inference} apply also to the transformed data vector.
\item If the vector of innovations $u_t$ were observed, an alternative estimator would regress $y_{1,t+h}$ onto $u_t$ and $y_{t-1},\dots,y_{t-p}$. As discussed in \cref{sec:ar1_intuition}, this estimator is numerically equivalent with $\hat{\beta}_1(h)$, so the uniformity result carries over.
\item It is easily verified in our proofs that, rather than initializing the process at zero, we can allow the initial conditions $y_0,\dots,y_{1-p}$ to be random variables that are independent of the innovations $\lbrace u_t\rbrace_{t \geq 1}$, as long as $E[\|y_\ell\|^4] < \infty$ for $\ell \leq 0$.
\end{enumerate}
\section{Bootstrap Implementation}
\label{sec:boot}
In this section we describe the bootstrap implementation of lag-augmented local projection that we recommend for practical use. We find in simulations that the bootstrap procedure is effective at correcting small-sample coverage distortions. These distortions arise primarily due to the small-sample bias of local projection, which \citet{Herbst2020} show is analogous to the well-known bias of the VAR OLS estimator \citep{Kilian1998}.
Our baseline algorithm is based on a wild autoregressive bootstrap design, which allows for heteroskedastic VAR innovations \citep{Goncalves2004} as in our theoretical results. Guided by simulation evidence, we construct the bootstrap confidence interval using the equal-tailed percentile-t method, which has a built-in bias correction (\citealp{Kilian1998}; \citealp[Chapter 12.2.6]{Kilian2017}).
The bootstrap procedure for computing a $1-\alpha$ confidence interval proceeds as follows, assuming a VAR($p$) model:
\begin{enumerate}
\item Compute the impulse response estimate of interest $\nu'\hat{\beta}_1(h)$ and its standard error $\hat{s}_1(h,\nu)$ by lag-augmented local projection as in \cref{sec:var_model}.
\item \label{itm:boot_var} Estimate the VAR($p$) model by OLS without lag augmentation. Compute the corresponding VAR residuals $\hat{u}_t$. Bias-adjust the VAR coefficients using the formula in \citet{Pope1990} (this adjustment is optional, but improves finite-sample performance).
\item Compute the impulse response of interest implied by the VAR model estimated in step \ref{itm:boot_var}. Denote this impulse response by $\nu'\hat{\beta}_\text{1,VAR}(h)$.
\item For each bootstrap iteration $b=1,\dots,B$:
\begin{enumerate}[i)]
\item Generate bootstrap residuals $\hat{u}_t^* \equiv U_t \hat{u}_t$, $t=1,\dots,T$, where $U_t \stackrel{i.i.d.}{\sim} N(0,1)$ are computer-generated random variables that are independent of the data.
\item Draw a block of $p$ initial observations $(y_1^*,\dots,y_p^*)$ uniformly at random from the $T-p+1$ blocks of $p$ observations in the original data.
\item Generate bootstrap data $y_t^*$, $t=p+1,\dots,T$, by iterating on the bias-corrected VAR($p$) model estimated in step \ref{itm:boot_var}, using the innovations $\hat{u}_t^*$.
\item Apply the lag-augmented LP estimator to the bootstrap data $\lbrace y_t^* \rbrace$. Denote the impulse response estimate and its standard error by $\nu'\hat{\beta}(h)^*$ and $\hat{s}_1(h,\nu)^*$, respectively.
\item Store $\hat{T}_b^* \equiv (\nu'\hat{\beta}_1(h)^*-\nu'\hat{\beta}_\text{1,VAR}(h))/\hat{s}_1(h,\nu)^*$.\footnote{It is critical that the bootstrap t-statistic $\hat{T}_b^*$ is centered at the VAR-implied impulse response $\nu'\hat{\beta}_\text{1,VAR}(h)$ rather than the LP-estimated impulse response $\nu'\hat{\beta}_1(h)$. This is because the former estimate is the pseudo-true parameter in the recursive bootstrap DGP, and the latter estimate differs from the former by an amount that is not asymptotically negligible.}
\end{enumerate}
\item Compute the $\alpha/2$ and $1-\alpha/2$ quantiles of the $B$ draws of $\hat{T}_b^*$, $b=1,\dots,B$. Denote these by $\hat{Q}_{\alpha/2}$ and $\hat{Q}_{1-\alpha/2}$, respectively.
\item Return the percentile-t confidence interval\footnote{It is not valid to use the Efron bootstrap confidence interval based on the bootstrap quantiles of $\hat{\beta}(h)^*$. This is because the bootstrap samples are asymptotically centered around $\hat{\beta}_\text{VAR}(h)$, not $\hat{\beta}(h)$.}
\[[\nu'\hat{\beta}_1(h)-\hat{s}_1(h,\nu)\hat{Q}_{1-\alpha/2},\nu'\hat{\beta}_1(h)-\hat{s}_1(h,\nu)\hat{Q}_{\alpha/2}].\]
\end{enumerate}
Instead of the above recursive VAR design, it is also possible to use the standard fixed-design pairs bootstrap, as in any linear regression with serially uncorrelated scores.\footnote{This is the bootstrap carried out by Stata's {\tt bootstrap} command with standard settings.} In this case, the usual Efron bootstrap confidence interval is valid, like the percentile-t interval. However, simulations suggest that the pairs bootstrap procedure is less accurate in small samples than the above recursive bootstrap design, mirroring the results in \cite{Goncalves2004} for autoregressive inference.
Our online code repository implements the above recommended bootstrap procedure, as well as several alternative LP- and VAR-based procedures, see \cref{fn:github}.
\section{Conclusion and Directions for Future Research}
\label{sec:conc}
Local projection inference is already popular in the applied macroeconomics literature. The simple nature of local projections has allowed the methods of causal analysis in macroeconomics to connect with the rich toolkit for program evaluation in applied microeconomics; see for example \citet{Angrist2018}, \citet{Nakamura2018}, \citet{Stock2018}, and \citet{rambachan2019econometric}. We hope the novel results in this paper on the statistical properties of local projections may further this convergence.
\paragraph{Recommendations for Applied Practice.}
The simplicity and statistical robustness of \emph{lag-augmented} local projection inference makes it an attractive option relative to existing inference procedures. We recommend that applied researchers conduct inference based on lag-augmented local projections with heteroskedasticity-robust (Eicker-Huber-White) standard errors. This procedure can be implemented using any regression software and has desirable theoretical properties relative to textbook delta method autoregressive inference and to non-augmented local projection methods. In particular, we showed that confidence intervals based on lag-augmented local projections that use robust standard errors with standard normal critical values are uniformly valid over the persistence in the data and for a wide range of horizons. We also suggested a simple bootstrap implementation in \cref{sec:boot}, which seems to achieve even better finite-sample performance.
Conventional VAR-based procedures deliver smaller standard errors than local projections in many cases, but this comes at the cost of fragile coverage, especially at longer horizons. In our opinion, there are only two cases in which the lag-augmented local projection inference method is inferior to competitors: (i) If the data is known to be at most moderately persistent and interest centers on very short impulse response horizons, in which case textbook VAR inference is valid and efficient. (ii) When the data has (near-)unit roots and interest centers on horizons that are a substantial fraction of the sample size, in which case the computationally demanding AR grid bootstrap may be deployed if feasible \citep{Hansen1999,Mikusheva2012}. In all other cases, lag-augmented local projection inference appears to achieve a competitive trade-off between robustness and efficiency.
How should the VAR lag length $p$ be chosen in practice? Naive pre-testing for $p$ causes uniformity issues for subsequent inference \citep{Leeb2005}. Though we leave the development of a formal procedure for future research (see below), our theoretical analysis yields three insights. First, users of local projection should worry about the choice of $p$ in order to obtain robust inference, just as users of VAR methods do. Second, $p$ should be chosen conservatively, as is conventional in VAR analysis \citep[Chapter 2.6.5]{Kilian2017}. In our framework there is no asymptotic efficiency cost of controlling for more than $p_0$ lags if the true model is a VAR($p_0$), and the simulation results in \cref{sec:sim_further} confirm that the cost is also small in finite samples. Third, the logic of \cref{sec:ar1_intuition} suggests that in realistic models where the higher-lag VAR coefficients are relatively small, it is not crucial to get $p$ exactly right: What matters is that we include enough control variables so that the effective regressor of interest approximately satisfies the conditional mean independence condition (\cref{asn:u_mds}).
\paragraph{Directions for Future Research.}
It would be interesting to relax the assumption of a finite lag length $p$ by adopting a VAR($\infty$) framework. We are not aware of existing work on uniform inference in such settings. One possibility would be to base inference on a sieve VAR framework that lets the lag length used for estimation tend to infinity at an appropriate rate as in \citet{Goncalves2007}. A second possibility is to impose \emph{a priori} bounds on the rate of decay of the VAR coefficients, and then take the resulting worst-case bias of finite-$p$ local projection estimators into account when constructing confidence intervals \citep[as in the ``honest inference'' approach of][]{Armstrong_Kolesar:2018}.
Due to space constraints, we leave a proof of the validity of the suggested bootstrap strategy to future work. It appears straight-forward, albeit tedious, to prove its pointwise validity. Proving uniform validity requires extending the already lengthy proof of \cref{thm:var_lp_inference}.
Several extensions of the results in this paper could be pursued by adopting techniques from the VAR literature. First, the results of \citet{PlagborgMoller2019} suggest straight-forward ways to generalize our results on reduced-form impulse response inference to \emph{structural} inference. Second, our assumption of no deterministic dynamics in the VAR model could presumably be relaxed using standard arguments. Third, by considering linear system estimators rather than single-equation OLS, our results on scalar inference could be extended to simultaneous inference on several impulses \citep{IK2016,MontielOlea2019}. Finally, whereas we adopt a frequentist perspective in this paper, it remains an open question whether local projection inference is relevant from a Bayesian perspective.