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.
121,337 characters · 22 sections · 0 citation commands
Simple robust two-stage estimation and inference for generalized impulse responses and multi-horizon causality
\listoftables
\listoffigures
\@starttoc{lth}
\addcontentsline{toc}{section}{List of Definitions, Assumptions, Propositions and Theorems}
\pagenumbering{arabic} \setcounter{section}{0} \setcounter{page}{1}
\pagestyle{headings}
{\thesection} {\thesection} \setcounter{theorem}{0} \setcounter{definition}{0} \setcounter{equation}{0} The concept of causality introduced by \nocite*{wiener1956theory} and \nocite*{granger1969investigating}, often referred to as Granger causality, constitutes a foundational approach for studying dynamic relationships between time series. While \nocite*{granger1969investigating} primarily focused on bivariate causality at a single horizon, the concept was generalized by \nocite*{dufour1998short} to account for multi-horizon causality.\footnote{For related work, see also \nocite*{sims1980macroeconomics}, \nocite*{hsiao1982autoregressive}, and \nocite*{lutkepohl1993testing}.} This extension serves two purposes. First, Granger causality at one horizon may not fully capture the causal dynamics in a multivariate system, as causality at one horizon does not necessarily imply causality at multiple horizons.\footnote{An example from \nocite*{dufour1998short}: In a trivariate system ($x,y,z$), $x_t$ may not Granger-cause $y_t$ at time $t$, but $x_t$ might still help predict $y_t$ several periods ahead through indirect effects mediated by $z_t$. For instance, $x_t$ could predict $z_t$ one period ahead, which in turn affects $y_t$ at a later time.} Second, the concept of impulse response proposed by \nocite*{sims1980macroeconomics} has limitations, as a zero impulse response is neither a necessary nor sufficient condition for non-causality in the sense of Wiener and Granger (see \nocite*{dufour1993relationship}). To address the issue of non-causality between two vectors of variables across multiple horizons in the presence of auxiliary variables, \nocite*{dufour1998short} introduced a necessary and sufficient condition on the coefficients in a multi-horizon linear projection model. These coefficients, referred to as ‘‘generalized impulse responses’’ (GIRs), include Sims’ impulse response as a special case. In this paper, we propose an innovative two-stage estimation method for the GIRs in a multi-horizon linear projection model\footnote{\nocite*{dufour2006short} refer to this model as ‘‘autoregression at horizon $h$’’ or ‘‘($p,h$)-autoregression’’ within a finite lag VAR($p$) framework. The autoregression at horizon $h$ model is also referred to as "Local Projection," which is widely used to estimate impulse response functions due to its straightforward inference [see \nocite*{jorda2005estimation}, \nocite*{montiel2021local}, and \nocite*{wolf2021same}].}, which complements two conventional methods: the recursive method and the Least Squares (LS) projection method.
The recursive method begins by estimating the VAR coefficients, followed by the recursive computation of GIRs, as these represent nonlinear transformations of the underlying VAR coefficients. Statistical tests for non-causality restrictions across multiple horizons (two or more) are expressed as zero constraints on multilinear forms involving the VAR coefficients. However, the Wald-type test criteria, can lead to asymptotically singular covariance matrices, making standard asymptotic theory inapplicable. This issue has been explored by \nocite*{lutkepohl1997modified}, \nocite*{benkwitz2000problems}, \nocite*{inoue2020uniform}, and \nocite*{dufour2023wald}. To address this challenge, \nocite*{pelletier2004problems} and \nocite*{dufour2006short} implemented a multi-horizon linear projection model with LS estimation.\footnote{The multi-horizon linear projection model has gained considerable attention in economics and finance, particularly for forecasting, impulse response estimation, causality testing, and causality measurement. For related work, see \nocite*{fama1988dividend}, \nocite*{campbell1988dividend}, \nocite*{dufour2010short}, \nocite*{dufour2012measuring}, \nocite*{zhang2016exchange}, \nocite*{salamaliki2019transmission}, \nocite*{shi2020causal}.} Appropriate asymptotic inference for the coefficient estimators can be obtained through adjustments for serially correlated residuals, such as using heteroskedasticity and autocorrelation consistent (HAC) estimators. However, LS-based multi-horizon linear projection estimates may suffer from reduced estimation efficiency as the horizon lengthens, and the reliability of confidence intervals constructed with HAC standard error estimates can be fragile. For further discussion, see \nocite*{dufour2006short}, \nocite*{kilian2011reliable}, and \nocite*{montiel2021local}. Although \nocite*{dufour2006short} proposed a closed-form formula for estimating the covariance matrix, they also cautioned that the finite-sample estimate may not always be positive semi-definite. Therefore, in this study, we conduct a comprehensive reevaluation of the estimation method and inference procedures for GIR coefficients, aiming to improve the reliability of covariance matrix estimates and potentially achieve more efficient coefficient estimates.
Our two-stage estimation method for GIRs makes five key contributions to the existing literature on multi-horizon linear projection and causality studies.
First, we introduce an innovative identification method for GIRs that leverages second-moment conditions of observables and VAR innovations. This novel approach offers an alternative to the traditional LS method for estimating multi-horizon linear projection models, drawing on the well-established Two-Stage Least Squares (2SLS) framework but diverging in its application. The method consists of two distinct stages: in the first, VAR residuals are estimated using LS, and in the second, these residuals are employed as instrumental variables in a 2SLS procedure. To address potential issues related to unit roots, we also incorporate a lag-augmented version of the two-stage estimation. Our approach complements the LS method in multi-horizon projection models, enhancing both estimation precision and the robustness of inference.
Second, we compare the estimation efficiency of our two-stage approach with the LS method through both theoretical and empirical analyses. The theoretical discussion is based on an illustrative AR(2) process, where we compare the asymptotic variances of the two methods using explicit formulas and the true values of the underlying data-generating process. This analysis shows that LS estimates outperform the two-stage estimates only when the horizon is very short or the data is highly cyclical. Empirically, we conduct Monte Carlo simulations on a bivariate VAR(2) process, which demonstrate that the two-stage estimates are generally no less efficient than LS estimates. Our findings serve as a reminder to the econometrics community that instrumental variable estimates in time series may not always suffer from efficiency loss compared to LS. Furthermore, our theoretical results show that the two-stage estimates are asymptotically as efficient as the lag-augmented two-stage estimates, contradicting the common belief that lag augmentation in time series typically reduces efficiency.
Third, we develop a simple and robust method for estimating the covariance matrix of GIR estimates, which relies on the long-run variance of the regression score function. This approach eliminates the need to correct for serial correlation in the projection error. Our method removes the necessity for HAC estimates for all coefficients in a multi-horizon linear projection, including those for impulse response functions and lagged coefficients. By bypassing the complications associated with selecting bandwidths and kernel functions in HAC estimation, our approach simplifies the process. Monte Carlo simulations show that our robust covariance estimates lead to more reliable empirical test results compared to HAC estimates.
Fourth, we contribute to econometric theory by deriving the asymptotic normality of two-stage GIR estimates uniformly over the parameter space, allowing the projection horizon to grow with the sample size. Specifically, we consider the persistence of underlying data processes to range from stationary to integrated of order one or order two. The projection horizon is permitted to expand infinitely at varying rates relative to the sample size, depending on the parameter space. This uniformity across the parameter space and allowance for long-range horizons provide robust theoretical justification for empirical macroeconomic applications.
Fifth, our empirical analysis of the GIRs for economic activity in response to economic uncertainty reveals both short-run (1–3 months) and long-run (around 30 months) causal effects. The findings underscore the longer-than-expected impact of uncertainty on economic activity. Additionally, our empirical GIR studies highlight the limitations of impulse response analysis. In the long run, the Wald test on impulse responses is statistically insignificant, while the causality test remains significant, suggesting that impulse responses may fail to capture persistent causal effects over extended horizons.
Relevant Literature--From an economic conceptual perspective, the multi-horizon linear projection model, as examined in this paper, forms a foundational framework in the study of causality and predictability \nocite*{lutkepohl1993testing}, \nocite*{dufour1998short}, \nocite*{dufour2006short}, \nocite*{dufour2010short}, \nocite*{diebold2014network}. The novel identification and estimation method we propose draws from innovation algorithms for (vector) ARMA processes \nocite*{hannan1982recursive}, \nocite*{brockwell1991time}, \nocite*{dufour2022practical}, and is closely related to recent literature on Local Projection-Instrumental Variable (LP-IV) methods, as explored by \nocite*{stock2018identification}. To the best of our knowledge, no existing studies have focused on using estimated VAR residuals as internal instruments to estimate GIRs within a multi-horizon linear projection model.
Regarding the statistical inference of two-stage estimates, the methodology for obviating HAC standard errors, even in the presence of serially correlated residuals, has been explored by \nocite*{montiel2021local}, \nocite*{breitung2023projection}, and \nocite*{xu2023local}. However, these studies primarily focus on Sims; impulse responses rather than all coefficients in the projection model. The concept of uniform inference, which controls the test level in the "worst" case scenario, is thoroughly discussed in \nocite*{dufour2003identification} and \nocite*{andrews2020generic}. Specific applications to time series can be found in \nocite*{mikusheva2007uniform}, \nocite*{mikusheva2012one}, \nocite*{inoue2020uniform}, \nocite*{montiel2021local}, and \nocite*{xu2023local}. Research on long (potentially infinite) projection horizons has also been conducted by \nocite*{montiel2021local} and \nocite*{xu2023local}. Our paper is the first to fully examine uniform inference for all coefficient estimates in a multi-horizon linear projection model, while also addressing restrictions on the projection horizon for stationary, I(1), and I(2) processes.
Outline--Section (ref) introduces the notation and outlines the data generation process. In Section (ref), we present the concept of two-stage identification and estimation. Section (ref) explores robust inference methods that eliminate the need to correct for serial correlation in the projection residuals. The main statistical results on uniform inference for two-stage estimates are detailed in Section (ref). Section (ref) compares the asymptotic efficiency of our two-stage estimates with standard LS estimates. Section (ref) presents the results of Monte Carlo simulations. In Section (ref), we apply the proposed statistical methodology for GIR estimation to analyze the dynamic causality of economic uncertainty on macroeconomic aggregates. Finally, we conclude the paper in Section (ref).
In this paper, we consider a multivariate time series following a Vector Autoregressive (VAR) process,
where ${y}_t$ represents a $K$-dimensional vector of variables, and ${u}_t$ is a vector of innovations. The sequence of ${u}_t$ is defined on a probability space $(\Omega,\mathcal{F},\mathbb{P})$. We define the information set at time $t$, denoted as $\mathcal{F}_t$, which is generated by $ \{y_s\}_{s\leq t}$. We consider ${u}_t$ as a white noise vector process with non-singular covariance matrix $\Sigma_u$. The order $p$ is a finite integer known a priori.\footnote{Two reasons support the assumption that the order $p$ is known a priori. First, the VAR coefficient estimates still have pointwise convergence to a standard Gaussian distribution as long as order $p$ can be consistently estimated by certain types of information criteria (see p62, \nocite*{kilian2017structural}). Second, in practice, order $p$ can be considered a known upper bound for the number of VAR lags. Thus, the assumption that $p$ is known a priori is not as restrictive as it might initially appear.} For the sake of notation simplicity, we omit the intercept and deterministic function in our notation, which generally do not impact the limiting results for the estimated coefficients.\footnote{In practice, researchers may incorporate dummy variables or polynomial functions of time $t$ to “detrend” the data, as discussed in \nocite*{toda1995statistical}.} The initial conditions are set as $y_t=0$ for $-p+1\leq t\leq 0$.
We investigate a broad range of time series processes, including stationary processes, I(1), and I(2) processes. In addition, we account for the presence of a 'local-to-unit' root, aiming to ensure uniform inference across the specified parameter space while allowing the projection horizon to grow infinitely. Our framework also accommodates the potential presence of cointegration within the process, further enhancing its practical applicability. Notably, we highlight that processes conforming to a VARX model can also benefit from our findings. This is because VARX is a specific case of VAR, where the equation for the exogenous variable imposes constraints on the coefficients of the endogenous variables, setting them to zero (see p. 74 in \nocite*{kilian2017structural}).
We build up the following setup for the VAR parameter space definition:
The parameter space encompasses various stationary, I(1), and I(2) VAR processes, with the integer ${\Greekmath 010E}$ serving as a control parameter that represents the degree of data persistence. When all roots are believed to be distant from the unit circle, ${\Greekmath 010E}=0$. If at most one root is close to or on the unit circle, ${\Greekmath 010E}=1$, while in the case where two roots may be unity, ${\Greekmath 010E}=2$. The parameter ${\Greekmath 011A}_{m,i}$, which lies between -1 and 1, can be treated as local-to-unity (e.g., ${\Greekmath 011A}_{m,i}=1-c/T$). In contrast, the roots in $B_{{\Greekmath 010E}}(L)$ are constrained to remain strictly bounded away from the unit circle. Finally, the non-singular rotation matrix $\Pi$ allows the process to incorporate cointegration.
Notably, although I(2) processes are rarely observed in practice, we include them in our framework to enhance the generalizability of our theory. Specifically, we illustrate the parameter space condition where $|{\Greekmath 011A}_{m,2}|=1$ if $|{\Greekmath 011A}_{m,1}|\geq 1-{\Greekmath 010F}/2$. This constraint implies that when $|{\Greekmath 011A}_{m,1}|$ is sufficiently close to the unit circle, $|{\Greekmath 011A}_{m,2}|$ must lie on the unit circle. However, this condition is not as restrictive as it may initially seem, as ${\Greekmath 010F}$ can be chosen arbitrarily small. This setup accommodates a wide range of persistent datasets, including cases where both roots equal unity, one root equals unity while the other is local-to-unity, or one root is local-to-unity while the other is strictly less than unity. The only exception is when both roots are local-to-unity, such as when $|{\Greekmath 011A}_{m,1}|=|{\Greekmath 011A}_{m,2}|=1-c/T$. This consideration arises because two local-to-unit roots offer limited practical meaning and would require additional notations in the derivation of the limiting distribution.
Our purpose is to establish a Gaussian limiting distribution for GIR estimates, covering cases where the data is persistent and integrated up to order two, while allowing the projection horizon to scale with the sample size. We show that as data persistence increases, the constraints on the projection horizon become more restrictive. This is an extension of traditional setups with a finite (fixed) horizon, the variation in the projection residual remains finite, regardless of the integration order.
Under the data generating process of (ref), t he multi-horizon linear projection model is represented as:
where $(\Phi_1^{(h)},\Phi_2^{(h)},\cdots,\Phi_p^{(h)})$ are GIRs defined by \nocite*{dufour1998short}. The recursive formula is detailed as follows,
where the coefficient $\Psi_h$ represents the Sims' impulse response function at horizon $h$. The residual term ${u}_t^{(h)}$ follows an MA($h-1$) process:
Generally, coefficients in (ref) can be estimated through two methodologies. The review of estimation methods is presented in the next subsection.
This subsection reviews two prominent approaches for estimating GIR coefficients: the Least Squares Projection (LS-Proj) method and the Recursive (RC) method. Given that equation (ref) represents a simultaneous equation system with identical regressors across equations, we can, without loss of generality, focus on the first equation. The projection model is expressed as:
where ${\Greekmath 010C}_h$ denotes the $p$-lag coefficients in the first row of the system, where ${\Greekmath 010C}_h = (\Phi_{1\bullet,1}^{(h)},\Phi_{1\bullet,2}^{(h)},\cdots,\Phi_{1\bullet,p}^{(h)})'$. The vector $x_t$ consists of lagged values of $y$, specifically $x_t := (y_t',\cdots,y_{t-p+1}')'$, and $e_{t,h}$ represents the first element of $u_t^{(h)}$.
LS Projection (LS-Proj) method: As demonstrated by \nocite*{dufour2006short}, the orthogonality between the residuals and the $p$-lag regressors allows the LS method to yield consistent coefficient estimates:
However, a key challenge arises in conducting statistical inference due to the presence of serially correlated residuals. It requires estimating the long-run variance (LRV) of the regression score function, $x_t e_{t,h}$,
The summation is truncated at order $h-1$ if $u_t$ satisfies a martingale difference sequence (m.d.s.) condition. When the forecast horizon exceeds one, the presence of serially correlated residuals often leads practitioners to rely on various HAC estimators for constructing confidence intervals or test statistics. \nocite*{dufour2006short} provide an explicit formula for $\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{LRV}(x_t e_{t,h})$. However, they also recognize that in small samples, the resulting estimates may not be positive semi-definite. The explicit covariance matrix estimate proposed by \nocite*{dufour2006short} corresponds to a special case of the HAC estimator with truncated lags and equal weights.
Recursive VAR (RC-VAR) method: The recursive method (RC) involves estimating coefficients using the recursive formula provided in equations (3.7) and (3.8) of \nocite*{dufour1998short}. In this approach, the first stage estimates the VAR slope coefficients, followed by the iterative computation of coefficients across different horizons using the recursive formula. Within the finite-order VAR framework, the RC method is typically implemented through a companion matrix, which transforms a VAR($p$) process into a VAR(1) model. This transformation simplifies the projection equations, as described in equation (2.2.9) of \nocite*{lutkepohl2005new}. Thus, the recursive method can be represented as follows:
where the vector ${\Greekmath 0117}_{1}$ is a unit-size vector with a conformable dimension, where the first element is one and all others are zero. The matrix $\hat{\Phi}$ represents the estimated companion matrix derived from the VAR coefficient estimates. It is important to note that the RC method is closely related to the approach used for computing impulse response functions, as demonstrated by \nocite*{lutkepohl1990asymptotic}.
Compared to the LS projection method, applying the RC method requires caution due to its reliance on the delta method for statistical inference, given the non-linear transformation involved in the estimation process. As noted by \nocite*{dufour2006short}, standard Wald-type test statistics may result in asymptotically singular covariance matrices due to the reduced rank of the Jacobian matrix. This issue can arise, for instance, if the underlying process is white noise, as recognized by \nocite*{benkwitz2000problems}, \nocite*{dufour2015wald}, \nocite*{inoue2020uniform}, and \nocite*{dufour2023wald}, among others. Consequently, conventional asymptotic theory may not apply to such statistics without the appropriate rank condition. This limitation partly explains the preference for the LS projection method, despite its potential inefficiency.
In this section, we present an innovative two-stage method for identifying and estimating the coefficients of multi-horizon linear projections.
We revisit the first row of the multi-horizon linear projection model, as defined in equation (ref):
The conventional identification approach relies on the assumption of weak exogeneity, expressed as $\mathbb{E}[x_t({y}_{1,t+h}-{\Greekmath 010C}_h' x_t)] = 0$, which is typically satisfied when $u_t$ is orthogonal to past observables. As an alternative, we propose using VAR innovations as instruments for identifying the coefficient parameters:
As the underlying process follows a VAR model, the system of equations can be written as:
This formulation provides an alternative identification strategy, utilizing the structure of the VAR model for efficient estimation. where $\overline{\Psi}_p$ is a $pk\times pk$ upper triangular matrix, $\overline{\Psi}_p = [\overline{\Psi}_{ij,p}]_{1\leq i,j\leq p}$, $\overline{\Psi}_{ij,p}=\Psi_{j-i}$ if $j\geq i$ and zero otherwise, $v_t=\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{P}(x_{t}|\mathcal{F}_{t-p})$ is a $pK\times 1$ vector of residuals, $v_t=(v_{1,t}',v_{2,t}',\cdots,v_{p,t}')'$, and $v_{i,t} =y_{t-i+1} - \sum_{j=i}^{p}\Psi_{j-i} u_{t+1-j} $.
The two-equation system described above can essentially be interpreted as a two-stage least squares (2SLS) procedure, where $z_t$ serves as the instrument. Thus, we can replace $x_t$ with its instrumental variable (IV) representation:
Since the instrument $z_t$ comprises innovations from time $t$ to $t-p+1$, it satisfies the weak exogeneity condition: $\mathbb{E}[e_{t,h} z_t] = 0$. Moreover, as each component of $z_t$ corresponds to the innovation term of the related observable, $z_t$ is naturally correlated with the regressor $x_t$. This correlation ensures that the triangular matrix $\overline{\Psi}_p$ is non-singular. Specifically, the covariance matrix $\mathbb{E}[z_t x_t']$ is given by:
The non-singularity of both the covariance matrix $\Sigma_u$ and the triangular matrix $\overline{\Psi}_p$ ensures the validity of the instrument. The coefficient ${\Greekmath 010C}_h$ can then be identified using the following moment conditions:
For instance, in a simple VAR(1) model where $x_t = y_t$ and $z_t = u_t$, we have $\mathbb{E}[z_t x_t'] = \Sigma_u$ and $\mathbb{E}[z_t y_{1,t+h}] = \Sigma_u \Psi_{1\bullet,h}'$. Consequently, ${\Greekmath 010C}_h$ is identified as:
This identification strategy provides an alternative estimation method using the instrument $z_t$, rather than the conventional use of $x_t$.
In this subsection, we propose two types of two-stage estimators: infeasible and feasible. As previously discussed, the identification of ${\Greekmath 010C}_h$ relies on the VAR innovation $u_t$. Infeasible estimates are derived using the actual, though typically unobserved, innovations $u_t$. In contrast, feasible estimates are computed using the LS-estimated residuals, $\hat{u}_t$.
Following the moment-based identification method discussed earlier, and assuming $z_t$ is available, the two-stage estimates can be computed as:
These estimates, $\tilde {\Greekmath 010C}_h^{2S}$, are termed infeasible because the innovation process $u_t$ is generally unobserved. In practice, it is common to replace $u_t$ with the estimated VAR residuals, $\hat{u}_t$, and derive an estimate for $z_t$ as follows:
where $\hat{u}_t = y_t - \sum_{i=1}^{p}\hat{\Phi}_i y_{t-i}$, and $\hat{\Phi}_i$ represents the coefficient matrices estimated through LS on a VAR($p$) model. This leads to the feasible estimator:
A critical step in deriving the asymptotic distribution for the feasible estimator is establishing its asymptotic equivalence with the infeasible estimator. Specifically, it is necessary to demonstrate that incorporating the estimated $z_t$ does not introduce asymptotic bias at the order of the square root of the sample size. This formal result is presented in Proposition (ref). It is important to note that using estimated variables as instruments in two-stage least squares (2SLS) differs fundamentally from employing estimated variables directly as regressors, as in the LS-based VARMA estimation method discussed by \nocite*{hannan1982recursive}. In the latter case, efficiency losses occur when estimated residuals are used as regressors, prompting the development of remedial procedures to improve efficiency. In contrast, our paper demonstrates that using estimated residuals as instruments does not impair efficiency.
With both feasible and infeasible two-stage estimators established, an important question arises: why transition from the well-established, straightforward LS-based estimation to the two-stage method? This question is addressed by highlighting three key advantages.
First, two-stage estimators typically yield more efficient results across a broad spectrum of data-generating processes and projection horizons. A detailed analysis of asymptotic efficiency is provided in Section (ref), where an illustrative AR(2) model offers intuitive insights into the enhanced efficiency of two-stage estimators. Additionally, the section on Monte Carlo simulations examines these efficiency comparisons from an empirical perspective.
Second, the two-stage approach addresses serial correlation in projection residuals by providing robust covariance estimation, thereby eliminating the need for HAC corrections. This leads to more reliable inference, as HAC estimates are known to be unreliable in finite samples. In Section (ref), we introduce a novel method for estimating the long-run variance of the regression score function, with detailed theoretical and methodological discussions.
Third, the two-stage method facilitates the analysis of infinite projection horizons, a significant advantage for macroeconometric research, where the projection horizon often represents a substantial portion of the sample. In contrast, the literature on standard LS estimation remains largely silent on infinite projection horizons. We address this gap by exploring long projection horizons and specifying the conditions under which the Gaussian limiting distribution holds for two-stage estimates.
Our two-stage estimates encounter non-standard convergence issues when the data exhibit a unit root.\footnote{For example, in the case of an AR(1) random walk, the covariance matrix $T^{-1}\sum_{t=1}^T u_t y_t$ does not converge to the variance of $u_t$ but instead to a stochastic process. Although this process has finite mean and variance, it is almost surely bounded by a constant, as guaranteed by Chebyshev's inequality. This ensures consistency, as the sample mean of the score function, $u_t e_{t,h}$, converges in probability to zero. However, this invalidates the use of Wald-type tests with standard critical values.} To address this, we introduce a lag-augmented version of the two-stage estimates, as proposed by \nocite*{toda1995statistical} and \nocite*{dolado1996making}.
The regression model under consideration is augmented with one or two additional lags, represented as:
for ${\Greekmath 010E} = 1, 2$, where the coefficient ${\Greekmath 010D}_{\Greekmath 010E}$ is a $({\Greekmath 010E} K)$-dimensional nuisance parameter associated with the extra lag(s), and ${\Greekmath 010D}_{\Greekmath 010E}$ equals to zero in the population. Here, $x_{t,1} := (x_t', y_{t-p}')'$ and $x_{t,2} := (x_t', y_{t-p}', y_{t-p-1}')'$. Notably, the lag-augmented regression nests the standard regression when ${\Greekmath 010E} = 0$, in which case $x_{t,0} = x_t$ and ${\Greekmath 010D}_0$ is an empty vector.
The inclusion of the extra lag(s) serves as a control variable, allowing the previous regressors to be transformed (e.g., first-differenced) into stationary variables. This transformation facilitates the derivation of the Gaussian limiting distribution. Instead of directly applying the Least Squares (LS) method, as done in the linear projection framework, we employ VAR residuals as instruments with lag augmentation: $z_{t,1}:= (z_t', y_{t-p}')'$ or $z_{t,2}:= (z_t', y_{t-p}', y_{t-p-1}')'$, using a two-stage least squares (2SLS) approach. Here, $z_t$ is defined as in equation (ref). This yields the infeasible lag-augmented estimates as follows:\footnote{When ${\Greekmath 010E}=0$, the estimate $\tilde {\Greekmath 010C}_h^{\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{LA(0)-2S}}$ coincides with the non-lag-augmented two-stage estimate $\tilde {\Greekmath 010C}_h^{\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{2S}}$. Thus, $\tilde {\Greekmath 010C}_h^{2S}$ ($\hat {\Greekmath 010C}_h^{2S}$) and $\tilde {\Greekmath 010C}_h^{LA(0)-2S}$ ($\hat {\Greekmath 010C}_h^{LA(0)-2S}$) are interchangeable in this paper.}:
where $H_{\Greekmath 010E}$ represents a $pK \times (p + {\Greekmath 010E})K$ selection matrix, defined as $H_{\Greekmath 010E} := (I_{pK}, 0_{pK \times {\Greekmath 010E} K})$. Following the rationale for lag augmentation as outlined in \nocite*{dolado1996making} and related literature, the first $p$ regressors in $x_{t,{\Greekmath 010E}}$ can be linearly transformed into a stationary process due to the lag augmentation. This leads to the convergence in probability of the sample covariance. Further details, along with the theoretical proof, are provided in Section (ref).
Since the innovation process $u_t$ is not directly observed, we propose feasible lag-augmented estimates:
where $\hat{z}_{t,{\Greekmath 010E}}$ is obtained by replacing $z_t$ in $z_{t,{\Greekmath 010E}}$ with $\hat{z}_t$, and $\hat{z}_t$ is calculated in the same manner as in the stationary case. The theoretical proof of the asymptotic equivalence between the feasible and infeasible lag-augmented two-stage estimates will be provided in Section (ref).
In this section, we propose an alternative estimation method for the long run variance of the regression score function, which obviates the need to correct the serial correlation of the score function.
In general, the necessity of HAC inference in conventional methods stems from the serial correlation in the regression score function, $x_t e_{t,h}$, particularly when the horizon $h$ exceeds one. This subsection introduces a simple method to obviate the need for HAC inference.
Let the regression score function for the two-stage estimates be denoted as
Under certain regularity conditions on the innovation process, such as strict stationarity and ergodicity, or strong mixing, the scaled summation of the regression score function is expected to converge in law,
where $\Bar{T}=T-h-p+1$, and $\Omega_{s,h}$ denotes the LRV of $s_{t,h}$,
Note the summation would be truncated at order $h-1$ if $u_t$ satisfies a m.d.s. condition. The sum of lead-lag covariance matrices stems from the serial correlation presenting in $e_{t,h}$, given its characterization as an MA($h-1$) process. Notably, the summation is truncated at order $h-1$ due to the property that $\mathbb{E}[s_{t,h} s_{t-k}']=\mathbf{0}$ for all $|k| \geq h$ under the m.d.s. assumption.
Consequently, to conduct statistical tests or establish confidence intervals, it becomes imperative to prove (ref) and propose a consistent estimate of $\Omega_{s,h}$. Typically, due to the lead-lag summation described in (ref) and the need of positive semi-definite LRV estimates in finite samples, researchers commonly resort to utilizing HAC estimators, such as the Newey-West estimator, along with specific bandwidth and weight selections. However, we have recognized that by imposing slightly more stringent assumptions on the innovation process $u_t$, it may be feasible to obviate the need for HAC estimation to achieve a positive semi-definite estimate of $\Omega_{s,h}$.
This subsection derive the limiting distribution for (ref) along with a simple and robust approach to estimate the long-run variance $\Omega_{s,h}$.
We perform an algebraic manipulation of the regression score function and introduce a new series, $s_{t,h}^{*}$, defined as:
Specifically, the construction of $s_{t,h}^{*}$ involves replacing the $j$-th component in $s_{t,h}$ with its $(j-1)$-th leading value. For instance, the second component in $s_{t,h}$, $e_{t,h}u_{t-1}$, is replaced by $e_{t+1,h}u_{t}$. The introduction of $s_{t,h}^{*}$ serves two purposes: (1) it reorders the elements of $s_{t,h}$, ensuring that both series have the same long-run variance; and (2) under certain regularity conditions on the $u_t$ process, $s_{t,h}^{*}$ becomes serially uncorrelated, and its long-run variance is identical to its variance. To clarify this property, we explicitly express $s_{t,h}^{*}$ as:
To avoid the need for HAC estimators and to propose positive semi-definite consistent estimates of $\Omega_{s,h}$, we introduce the following assumption:
There are three main reasons for adopting the mean-independence assumption in our innovation process. First, it provides a sufficient condition for proposing a convenient method to estimate the long-run variance $\Omega_{s,h}$. Typically, when estimating $\Omega_{s,h}$, which involves summing lead-lag covariances, HAC estimators are required to ensure positive semi-definiteness in finite samples. However, under Assumption (ref), $\Omega_{s,h}$ becomes algebraically equivalent to the variance of $s_{t,h}^*$. For illustration, consider the case where $t < {\Greekmath 011C}$, which leads to:
where the second equality holds because the explicit form of $s_{t,h}^*$ reveals that $s_{{\Greekmath 011C},h}^*$ and $(e_{t,h}, e_{t+1,h}, \cdots, e_{t+p-1,h})$ are measurable with respect to the information set ${\Greekmath 011B}(u_{t+1}, u_{t+2}, \cdots)$, since $e_{t+i,h}$ is a linear combination of $u_{t+1+i}, \cdots, u_{t+h+i}$. This framework, previously applied by \nocite*{montiel2021local} in impulse response estimation, is extended in our study to encompass all coefficients within a multi-horizon linear projection model.\footnote{\nocite*{montiel2021local}, footnote 7, notes that inference for lagged coefficients in Local Projection models typically relies on HAC estimators. We address this by demonstrating that using VAR-estimated residuals as instruments yields a standard limiting distribution for these coefficients without HAC estimators.}
It is important to note, as \nocite*{xu2023local} highlights, that the mean-independence assumption is sufficient but not necessary to achieve a zero-correlation result for the regression score function. Moreover, this assumption may not hold in all contexts, especially in the presence of skewed distributions or conditional heteroskedasticity, which could violate it. While we acknowledge these limitations, the second and third motivations for using the mean-independence assumption may outweigh these concerns, making it a valuable assumption in this context.
Second, the mean-independence assumption facilitates the application of the Central Limit Theorem (CLT) for martingale difference sequences to establish convergence in distribution. In conventional frameworks with finite horizons, various versions of the CLT can be employed. The key point here is that a finite horizon guarantees a finite variance for the regression score function, ensuring that ${\Greekmath 0115}_{\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{max}}(\Omega_{s,h}) < c < \infty$. However, when the horizon grows relative to the sample size, as often occurs in macroeconomic applications where the projection horizon constitutes a significant portion of the sample size, the variance may become unbounded, particularly in processes with unit roots. This issue is observed in Monte Carlo simulations, where increasing the horizon degrades the empirical size of statistical tests based on asymptotic variance and critical values from a standard Gaussian distribution, especially in highly persistent processes. Therefore, we rely on the CLT for martingale difference sequences, which accommodates unbounded variance.
Third, while convergence can still be established under weaker conditions, such as assuming a martingale difference sequence and mixing, these conditions impose restrictions on the horizon $h$. \nocite*{xu2023local} proposed an alternative method to derive asymptotic results by restructuring the regression score function in a more complex manner. This approach introduces two asymptotically negligible terms, with an order dependent on the horizon $h$, implying that the theorem would only hold if $h/T \to 0$ in the limit. In our view, this condition applies even for finite-order stationary VAR processes, which are commonly encountered in practice. While this restructuring allows for relaxing the mean-independence assumption to a martingale difference sequence assumption, it limits the horizon $h$ from growing toward infinity at the rate permitted under stationarity. Our results, by contrast, allow the horizon $h$ to represent a non-trivial proportion of the sample size $T$. Therefore, we adopt the mean-independence assumption in this paper.
In light of these considerations, the algebraic manipulation and the result presented in Equation (ref) demonstrate that the variance of $s_{t,h}^*$ is equal to the long-run variance (LRV) of $s_{t,h}$. This insight allows us to bypass the use of HAC standard errors, as the LRV of $s_{t,h}$—which involves summing lead-lag covariances—typically requires HAC estimation to ensure positive semi-definiteness in finite samples. The natural estimate of the variance of $s_{t,h}^*$, however, satisfies this criterion directly.
We formalize these findings in the following lemma:
The proof of Lemma (ref) is provided in Appendix (ref). Lemma (ref) states that $\Omega_{s,h}$ is identical to the variance matrix of $s_{t,h}^*$, offering a straightforward method for estimating the LRV $\Omega_{s,h}$. Empirical researchers need only estimate the sample variance of $s_{t,h}^*$, rather than the LRV of $s_{t,h}$. As the sample variance is naturally positive semi-definite, this approach eliminates the need for HAC correction of the serial correlation.
The asymptotic convergence to normality will be established in two steps. First, we decompose the summation of $s_{t,h}$ into three parts: the primary component, which is the summation of the reordered regression score function $s_{t,h}^*$, and two additional asymptotically negligible terms, $\overline{s}_{1,h}$ and $\overline{s}_{2,h}$. Second, we demonstrate the convergence in distribution of the summation of $s_{t,h}^*$ using the martingale Central Limit Theorem, thereby confirming the convergence of the regression score function. The primary challenge arises from the potentially unbounded variance of $e_{t,h}$ as the horizon $h$ increases, particularly in cases of persistent data.
The summation of $s_{t,h}$ is expanded as follows:
where $\overline{s}_{1,h}$ and $\overline{s}_{2,h}$ are two terms of order ($p-1$), defined as $\overline{s}_{1,h} = \sum_{i=1}^{p-1} s_{i,h}^{**}$ and $\overline{s}_{2,h} = \sum_{i=1}^{p-1} s_{i,h}^{***}$, with $s_{i,h}^{**} = J_{(-i)} s_{i,h}^{*}$ and $s_{i,h}^{***} = J_{(p-i)} s_{\bar{T}+i,h}^{*}$. Here, $\bar{T} = T - h - p + 1$, and $J_{(i)}$ and $J_{(-i)}$ denote the first $iK$ rows and the last $(p-i)K$ rows of an identity matrix of dimension $pK$, respectively, such that $[J_{(i)}', J_{(-i)}'] = I_{pK}$. The terms $\overline{s}_{1,h}$ and $\overline{s}_{2,h}$ arise from the reordering process. When the VAR order $p$ is finite, these terms are bounded by a constant scaled by the standard error of $s_{t,h}^*$, almost surely, as guaranteed by Chebyshev's inequality.
To establish the limiting result, we show that the scaled summations of $\sum_{t=p}^{T-h} s_{t,h}$ and $\sum_{t=p}^{T-h-p+1} s_{t,h}^*$ are asymptotically equivalent.
The proof of Lemma (ref) is provided in Appendix (ref). In essence, Lemma (ref) asserts that the Euclidean norm of $\overline{s}_{1,h} + \overline{s}_{2,h}$ is asymptotically negligible at the scale of $\bar T^{-1/2}$. Although the summation involves a finite number of terms (with finite order $p$), the variance is not guaranteed to be constant when the data exhibits persistence and the horizon $h$ approaches infinity. The proof hinges on showing that the variance of each term becomes negligible relative to $\bar T$, a condition that depends on the boundedness of the horizon $h$ and the persistence of the data.
Building on Lemma (ref), proving the convergence of $\bar T^{-1/2} w' \sum_{t=p}^{T-h} s_{t,h}$ is equivalent to establishing the convergence of $\bar T^{-1/2} w' \sum_{t=p}^{\bar T} s_{t,h}^*$. Next, we impose the following regularity conditions on the innovation process, followed by the result on distributional convergence.
The proof of Proposition (ref) is provided in Appendix (ref). This proposition represents a key step in establishing asymptotic results for our two-stage estimates. It states that, for data processes of orders I(0), I(1), or I(2), and under certain conditions regarding the projection horizon $h$, the scaled summation of the regression score function converges to a Gaussian distribution when normalized by $\bar T^{1/2}$. By combining this result with the earlier finding that the LRV matrix $\Omega_{s,h}$ is equivalent to the covariance matrix of $s_{t,h}^*$, we conclude that $\Omega_{s,h}$ can be consistently estimated using a sample covariance matrix.
Proposition (ref) is one of econometrics contributions to the statistical inference of multi-horizon projection estimates. It demonstrates that, in a general projection (or prediction) model, the estimates can bypass the need for HAC inference, even in the presence of serial correlation in the residuals. Instead, the two-stage estimates can rely on White’s heteroskedasticity-robust inference, reconstructed from the regression score function.
In this section, we present asymptotic uniform inference for our two-stage estimates. Recognizing the potential presence of a unit root in the process and the need to establish the corresponding asymptotic distributional theory, we introduce a square matrix of dimension $pK$, denoted as $G_{{\Greekmath 010E}, T}$. This matrix is a function of both the transformation matrix and the probability scaling matrix, and is defined as follows:
where $\bar P_{\Greekmath 010E}, \bar \Upsilon_{\Greekmath 010E}$ are square matrices of dimension $pK$,
The matrices $\Upsilon_1$ and $\Upsilon_2$ represent $K$-dimensional diagonal probability scaling matrices, defined as $\Upsilon_{{\Greekmath 010E}} = \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{diag}[g_{m,{\Greekmath 010E}}^{-1/2}]_{1 \leq m \leq K}$, where $g_{m,1} = \min(T, \frac{1}{1 - |{\Greekmath 011A}_{m,1}|})$ and $g_{m,2} = g_{m,1}^2 \min(T, \frac{1}{1 - |{\Greekmath 011A}_{m,2}|})$. In the case of ${\Greekmath 010E} = 1$, when ${\Greekmath 011A}_{m,1}$ is bounded away from the unit disk, $g_{m,1}$ is constant. However, when ${\Greekmath 011A}_{m,1}$ is local to unity (e.g., $1 - c/T$) or lies on the unit circle, $g_{m,1}$ becomes proportional to $T$. For ${\Greekmath 010E} = 2$, the situation is more complex. If only one root is local to unity or a unit root, then $\bar{g}_{m,2}$ is proportional to $T$. If two roots are on the unit circle, then $\bar{g}_{m,2}$ scales as $T^3$. This result aligns with the convergence rate of the sample variance required by integrated time series variables, which is $T^2$ and $T^4$ for I(1) and I(2) processes, respectively (see \nocite*{park1988statistical}, \nocite*{phillips1988testing}, \nocite*{toda1995statistical}, among others).
The matrix $\bar{P}_{\Greekmath 010E}$ is a differenced matrix designed to prevent perfect multicollinearity in the limit for cointegration, following standard procedures for autoregressive time series models with unit roots (\nocite*{sims1990inference}). If all diagonal elements of $P_1$ are equal to one, $\bar{P}_1$ simply applies the first difference to $\Pi y_t$. The transformation matrix $G_{{\Greekmath 010E},T}$, combined with the probability scaling matrix, provides the necessary high-level condition for our asymptotic results.
To derive the asymptotic distribution, we propose the following limiting conditions:
Assumption (ref) serves as a commonly required intermediate step in establishing the limiting distribution of OLS estimates in autoregressive models. It states that the scaled covariance matrix is non-singular with probability one. It is a high level assumption seen in Local Projection literature, e.g., see \nocite*{montiel2021local} and \nocite*{xu2023local}.
In this subsection, we first present the asymptotic normality results for the infeasible estimates, which rely on the unobserved innovation process. We then establish the asymptotic equivalence between feasible and infeasible estimates, demonstrating that the two estimators are identical at the scale of $\bar T^{1/2}$.
To illustrate the asymptotic variance of both feasible and infeasible estimates, we define $\Omega_{{\Greekmath 010C},h}$ as follows:
where $\Sigma_{zx} = \mathbb{E}[z_t x_t']$, and $\mathbb{E}[z_t x_t']$ and $\Omega_{s,h}$ are defined in (ref) and (ref), respectively.
The limiting distribution of the infeasible estimates is provided below.
Proposition (ref) provides asymptotic statistical inference for both lag-augmented and none lag-augmented two-stage infeasible estimates. It is important to note that increasingly restrictive conditions on the projection horizon are necessary to account for the persistence in the underlying data process. This restriction primarily stems from the requirements of the martingale Central Limit Theorem and the asymptotic negligibility of the bias introduced by lag-augmentation.
A key insight from Proposition (ref) is that lag-augmented and none lag-augmented two-stage infeasible estimates share the same asymptotic variance, $\Omega_{{\Greekmath 010C},h}$. This reveals a surprising constancy in asymptotic efficiency, regardless of whether lag-augmentation is employed. This finding is a contradiction with the conventional view in standard LS estimation, where asymptotic efficiency is expected to vary with the lag-augmentation.
Moreover, we establish a crucial result confirming the asymptotic equivalence of infeasible and feasible estimates.
The proof of Proposition (ref) is in Appendix (ref). This demonstrates that incorporating estimated residuals as instruments does not introduce asymptotic bias.
This subsection presents the uniform inference for both lag-augmented and non-lag-augmented two-stage estimates, scaled by the estimated covariance matrix. First, we outline the method for estimating the covariance matrix. The estimator for $\Omega_{{\Greekmath 010C},h}$ is given by:
where:
where $\hat{\Sigma}_u = T^{-1} \sum_{t=1}^T \hat{u}_t \hat{u}_t'$, and $\hat{\overline \Psi}_p'$ is constructed from the recursive VAR-based impulse response functions. Here, $\hat{e}_{t,h}$ represents the LS residual of the multi-horizon linear projection with order $p$, and $\hat{u}_t$ is the standard LS residual of the VAR($p$) model.
We deliberately choose the explicit formula for $\Sigma_{zx}$ in (ref) for estimation rather than relying on the sample average, such as $\bar{T}^{-1} \sum_{t=p}^{T-h} \hat{z}_t x_t'$, commonly used in the stationary case. This choice is based on three key considerations: (1) The sample average method loses $h$ observations due to projection. While asymptotically consistent as long as $T-h \rightarrow \infty$, this approach may suffer from efficiency loss and finite sample bias. Recursive VAR-based estimates, in contrast, are well-known for their efficiency. Practitioners can also apply finite sample bias correction for LS VAR slope coefficient estimates (\nocite*{pope1990biases}). (2) The closed-form expression of $\Sigma_{zx}$ reveals that $\overline{\Psi}_p$ is a block lower triangular matrix, allowing for greater precision by setting zeros in the upper triangular blocks and ones on the main diagonal. By contrast, the sample average method estimates the entire matrix directly, with blocks in the upper triangular area expected to be zero in the population but not necessarily so in finite samples. Although practitioners could manually set these blocks to zero to enhance precision, this adds an extra step to the algorithm. (3) The closed-form formula is flexible, applying to both lag-augmented and non-lag-augmented estimates, as they share the same structure. In contrast, the sample average method requires the use of $\hat{z}_{t,{\Greekmath 010E}}$ and $x_{t,{\Greekmath 010E}}$ depending on ${\Greekmath 010E}$. Our approach simplifies implementation in computer algorithms and facilitates theoretical proof of consistency.
We summarize the uniform inference for feasible estimates with the robust covariance matrix estimates in the following proposition.
See the proof of Proposition (ref) in Appendix (ref). It is noteworthy that the condition on the horizon $h$ becomes more restrictive as the underlying process becomes more persistent. Readers may wonder about the intuition behind the term $\max(1, \frac{\bar h}{w'\Omega_{{\Greekmath 010C},h} w})$. This term suggests that the horizon $h$ depends on the scale of $w'\Omega_{{\Greekmath 010C},h} w$: if $w'\Omega_{{\Greekmath 010C},h} w$ is proportional to $h$, the condition on $h$ can be relaxed to $h/T, h^3/T \xrightarrow{p} 0$ for the cases of ${\Greekmath 010E} = 1, 2$, respectively. Conversely, in the worst case where $w'\Omega_{{\Greekmath 010C},h} w$ is constant, the condition on $h$ tightens to $h^2/T, h^4/T \xrightarrow{p} 0$ for ${\Greekmath 010E} = 1, 2$, respectively.
This difference arises from the long-run variance (LRV) of the regression score function, $\Omega_{s,h}$, which is embedded within the asymptotic variance matrix $\Omega_{{\Greekmath 010C},h} = \Sigma_{zx}^{-1} \Omega_{s,h} \Sigma_{zx}^{'-1}$. The maximum eigenvalue of $\Omega_{s,h}$ can be proportional to the horizon $h$ when unit roots are present. However, the minimum eigenvalue is only guaranteed to be bounded below by a constant.\footnote{Consider the following example: in a scalar autoregressive process with $u_t \overset{i.i.d.}{\sim} (0, 1)$, the upper-left $2 \times 2$ block of $\Omega_{s,h}$ has elements $\sum_{i=0}^{h-1} {\Greekmath 0120}_i^2$ on the main diagonal and $\sum_{i=0}^{h-2} {\Greekmath 0120}_i {\Greekmath 0120}_{i+1}$ off the diagonal. If the AR process is a random walk, the impulse response functions are ${\Greekmath 0120}_h = 1$ for all $h$. This results in $\sum_{i=0}^{h-1} {\Greekmath 0120}_i^2 = h$ and $\sum_{i=0}^{h-2} {\Greekmath 0120}_i {\Greekmath 0120}_{i+1} = h-1$. Thus, the maximum eigenvalue of this block equals $2h - 1$, while the minimum eigenvalue remains 1.} Consequently, the scale of $w'\Omega_{{\Greekmath 010C},h} w$ is crucial in determining the condition on $h$ when the data may be non-stationary.
Remarks:
In this section, we aim to compare the asymptotic efficiency of our two-stage GIR estimates with LS GIR estimates using an illustrative AR(2) process:
where ${\Greekmath 011E}_1 = {\Greekmath 011A}_1 + {\Greekmath 011A}_2$ and ${\Greekmath 011E}_2 = -{\Greekmath 011A}_1 {\Greekmath 011A}_2$, with ${\Greekmath 011A}_1$ and ${\Greekmath 011A}_2$ representing the autoregressive roots. The parameter ${\Greekmath 011A}_1$ is varied through a grid search over the interval $[0.01, 0.99]$ in increments of 0.01, while ${\Greekmath 011A}_2$ takes values from the set $\{0.8, 0.5, 0.2, -0.5\}$. We focus on the AR(2) process instead of AR(1) because it captures a wider range of empirical data characteristics, including stationarity, persistence, and cyclicality.
The asymptotic efficiency is evaluated by comparing the asymptotic variance, which is computed using a closed-form formula. For LS projection method estimates in (ref), the asymptotic variance, denoted as $\Omega_{{\Greekmath 010C},h}^{LS} := \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{AsymVar}(\hat{{\Greekmath 010C}}_h^{LS})$, is given by:
The matrix $\mathbb{E}[x_t x_t']$ contains the autocovariance of $y_t$. The truncation of the long-run variance at order $h-1$ reflects the i.i.d. assumption on $u_t$. This assumption also simplifies computation, as $\mathbb{E}[x_t x_{t+k}' e_{t,h} e_{t+k,h}] = \mathbb{E}[x_t x_{t+k}'] \mathbb{E}[e_{t,h} e_{t+k,h}]$. For the two-stage projection estimates, we use the formula for $\Omega_{{\Greekmath 010C},h}$ provided in (ref).
In Figure (ref), we compare the efficiency of two estimators, ${\Greekmath 011E}_1^{(h)}$ and ${\Greekmath 011E}_2^{(h)}$, across horizons $h = 1, 2, \dots, 36$. Green/red coloring is used to indicate whether the two-stage estimate is more/less efficient than the LS projection estimate. Varying shades of green represent the degree to which the two-stage estimate outperforms LS projection estimates in terms of standard error. Standard error, rather than variance, is chosen as it directly reflects the narrowing of the confidence interval. The lightest green indicates that the confidence interval of the two-stage estimate is up to 10% narrower than that of the LS projection estimates. The medium green shows a 10% to 30% narrowing, while the darkest green reflects a confidence interval more than 30% narrower.
Figure (ref) demonstrates that two-stage estimates only exhibit less efficiency than LS estimates for short horizons, mostly at $h = 1$. In this case, the LS projection estimate coincides with the pseudo maximum likelihood estimate, yielding minimal variance. However, as $h$ increases, it becomes increasingly rare for LS projection estimates to outperform two-stage projection estimates. This occurs in only two cases, both at short horizons: (1) when the data is highly persistent, with ${\Greekmath 011A}_1 = 0.8$ and ${\Greekmath 011A}_2$ approaching 1, and (2) when the data is cyclical, with ${\Greekmath 011A}_1 = -0.5$.
In the first scenario, characterized by highly persistent data where both roots approach one, the outperformance of LS projection estimates at short horizons can be attributed to the well-known phenomenon of super-consistency. However, this advantage diminishes quickly as the horizon $h$ increases. This is because the variance matrix of LS projection estimates requires summing $h-1$ lead-lag autocovariances of $e_{t,h}$, resulting in large values due to persistence. In contrast, the variance matrix of two-stage estimates only depends on the variance of $e_{t,h}$, contributing to its greater efficiency.
In the second scenario, where the data exhibits cyclical patterns, the LS projection estimate's efficiency advantage arises from the behavior of the LRV matrix. The summation of $h-1$ lead-lag autocovariances in the LS projection variance matrix is sensitive to the cyclical nature of the data. When the autocovariances of $e_{t,h}$ differ in sign, they can offset one another, leading to a smaller LRV for $e_{t,h}$ compared to cases where both roots have the same sign and autocovariances decay gradually to zero.
In summary, two-stage estimates typically outperform LS projection estimates across a broad range of parameter spaces and projection horizons, particularly when the horizon is moderately large or the data exhibits moderate stationarity. The only cases where LS projection estimates have an advantage are (1) for a horizon of $h = 1$, (2) for extremely persistent data with a short horizon ($h \leq 3$), or (3) for persistent and cyclical data with a short horizon. Overall, two-stage estimates offer significant efficiency gains in a wide range of practical applications.
In this section, we conduct Monte Carlo simulations to assess the performance of our estimation and inference methods for GIR coefficient estimates. The data generating process (DGP) is based on a bivariate VAR(2) model, with initial observations set to $y_1 = y_2 = \mathbf{0}$. The DGP is specified in terms of autoregressive roots, allowing us to control the stationarity or non-stationarity of the process by adjusting the eigenvalues of the root matrices. We consider a range of processes, including white noise, stationary, I(1), and I(2) models, to provide a comprehensive evaluation of the finite sample performance of our two-stage estimation models with and without lag-augmentation.
We examine three versions of our two-stage estimation model: the standard two-stage model without lag-augmentation, the two-stage model with one additional lag, and the two-stage model with two additional lags. As benchmarks, we include the recursive VAR-based estimation method with delta-method inference and the standard least squares (LS) linear projection model with HAC inference. Additionally, we provide $t$-interval bootstrap results for the three versions of our two-stage estimates.\footnote{The algorithm for simulation-based inference is presented in Appendix (ref).} The horizons considered are $h = 1, 3, 6, 12, 24, 36$. The goal of these simulations is to evaluate the finite sample performance of the estimates in terms of bias, root mean squared error (RMSE), empirical test size, confidence interval coverage ratio, and average confidence interval width. For the benchmark LS projection model with HAC inference, we use MATLAB’s "hac" function with a Bartlett kernel and bandwidth equal to $h$ to estimate HAC standard errors. The recursive VAR-based method uses the closed-form formula for covariance matrix estimation as described in \nocite*{lutkepohl1990asymptotic}.
Simulation results for stationary, I(1), and I(2) processes are presented in Tables (ref) and (ref). The abbreviations for each model are as follows: (1) RC-VAR: recursive VAR-based method with delta-method inference, utilizing the explicit formula for the Jacobian matrix reported in Proposition 3.6 (p. 110) of \nocite*{lutkepohl2005new}; (2) LS-Proj: LS projection method with HAC inference, using critical values from the $z$-table (MATLAB: ‘hac’ command, kernel = Bartlett, bandwidth = h); (3) 2S-Proj(${\Greekmath 010E}$): two-step method with ${\Greekmath 010E}$ lag-augmentation, using critical values from the $z$-table; and (4) $\textup{2S-Proj}({\Greekmath 010E})_b$: two-step method with equal-tailed percentile t-interval bootstrap (Wild bootstrap, as per \nocite*{gonccalves2004bootstrapping}). The results for white noise and I(2) processes are shown in Tables (ref) and (ref) in Appendix (ref). Specific parameter values and eigenvalues of the two polynomial root matrices are detailed in the table footnotes. These tables highlight the performance of two coefficients, ${\Greekmath 011E}_{12,1}^{(h)}$ and ${\Greekmath 011E}_{12,2}^{(h)}$, which represent the causality of the second variable on the first variable at horizon $h$.
Tables (ref) and (ref) offer a comprehensive evaluation of the performance of GIR coefficient estimates across different estimation methods. A comparison between the two-stage projection estimates and the LS projection estimates reveals that the two-stage estimates generally exhibit satisfactory finite sample performance. These estimates are characterized by more accurate coverage ratios and, in certain cases, greater efficiency, as reflected in the narrower average width of the confidence intervals. In contrast, the LS projection estimates with Newey-West standard errors show distorted coverage ratios, with the distortion becoming more pronounced for persistent data and at longer horizons.
In terms of efficiency, measured by the average width of confidence intervals, the two-stage estimates generally outperform multi-horizon linear projection estimates, irrespective of whether lag-augmentation is applied. As expected from econometric theory, iterated VAR estimates are more efficient than either LS or two-stage projection methods. However, iterated VAR estimates exhibit the poorest coverage ratios due to reliance on delta-method inference.
In summary, the simulations demonstrate that GIR estimates obtained from the two-stage multi-horizon linear projection, utilizing our robust covariance matrix estimates, show relatively reliable finite sample performance compared to recursive estimates based on the explicit formula of delta-method inference or LS projection estimates with HAC standard error adjustments. While the two-stage estimates experience some efficiency loss due to the multi-horizon projection residual, they generally improve efficiency relative to LS projection estimates.
In this section, we examine the bidirectional Granger causality between economic uncertainty and various macroeconomic indices across multiple horizons. The concept of economic uncertainty has garnered significant attention in recent literature, as highlighted by studies such as \nocite*{baker2016measuring}, \nocite*{jurado2015measuring}, and \nocite*{rossi2015macroeconomic}. Empirical investigations into the causality of economic uncertainty can be found in works like \nocite*{salamaliki2019transmission}, \nocite*{gao2022oil}, and \nocite*{fernandez2023cross}, among others. We implement the Granger causality test over multiple horizons, report the Wald test statistics and corresponding $p$-values, and present both tables and figures for the GIR coefficient estimates, clarifying the causal relationship between economic uncertainty and macroeconomic indices.
The dataset used in this analysis consists of five monthly variables spanning from January 1974 to June 2023, totaling 593 observations. These variables include the Chicago Fed National Activity Index (CFNAI), the economic uncertainty index (sourced from \nocite*{jurado2015measuring}, JLN 3-month Ahead Macroeconomic Uncertainty, hereafter referred to as JLN), the unemployment rate (Unemp), inflation, and the Federal Funds Rate (FFR)\footnote{See data source in Appendix}. The focus of this investigation is the multi-horizon Granger causality from the uncertainty index to the macroeconomic indicators.
The VAR model employed in this analysis includes 12 lags and an intercept term. This lag order selection is based on the information criteria computed using the ‘VARselect’ command in the R package ‘vars’ (version 1.6-1), with lag.max = 20 and type = ‘const’. The results of the information criteria selection are as follows: AIC = 13, HQ = 3, SC = 2, and FPE = 13. Given that our dataset consists of typical monthly macroeconomic aggregates, we select 12 lags. This choice aligns with the argument made by \nocite*{hamilton2004comment} for impulse response estimation in monthly VAR models. However, we acknowledge the potential sensitivity of empirical results to the VAR lag order selection. To address this, we perform a robustness check by extending our estimations to lag orders of 15 and 18, as recommended by \nocite*{kilian2011reliable} in their empirical exercises.
In this study, we apply the Wald test on a set of GIRs to assess causality from one variable to another across specific horizons, as expressed by the null hypothesis: \[ \mathcal{H}_0: \Phi_{ij,k}^{(h)} = 0, \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{ for } k = 1, 2, \dots, 12. \] The coefficient estimates and covariance matrix are derived using the two-stage estimation method outlined in this paper. Table (ref) presents the Granger causality from the row variables to the column variables across multiple horizons. Notably, the cells along the main diagonal demonstrate the predictability of each variable concerning its own future output. It is evident that all five macroeconomic variables exhibit significant persistence.
Of particular interest is the finding that economic uncertainty (JNL) exhibits a causal effect on the CFNAI both in the short run (1–7 months) and the long run (approximately 2.5 years). While the short-term impact of uncertainty on economic activity is expected, its prolonged influence on the economy surpasses conventional expectations. Conversely, economic uncertainty has a sustained effect on the labor market, with JNL exerting causal influence on unemployment over a span of forty months. Surprisingly, JNL does not significantly affect the Federal Funds Rate (FFR), suggesting that the central bank may not adjust interest rates in direct response to fluctuations in economic uncertainty. Instead, the CFNAI shows causality on the FFR, which aligns with the widely held view that central banks use interest rate policy to stimulate or slow down the economy based on economic activity.
Table (ref) also provides empirical evidence that statistical insignificance of impulse responses from zero does not necessarily imply zero causality over multiple horizons, and, conversely, statistically significant impulse responses from zero do not necessarily imply causality. For instance, Sims' impulse response (the first column next to the causality test) is significant for horizons 9, 12, 15, and 18, yet none of these horizons exhibit significant multi-horizon causality tests. This discrepancy may be attributed to multicollinearity. On the other hand, at horizons 27 and 33, the causality test is statistically significant, but Sims' impulse response is not significantly different from zero. Thus, relying solely on Sims' impulse response for assessing multi-horizon causality could yield misleading results.
In our analysis, we specifically investigate the causality from the JNL index to CFNAI by visualizing the Wald test statistics for the null hypothesis of non-causality. Recognizing the potential sensitivity of empirical results to order selection, as highlighted by \nocite*{hamilton2004comment} in the context of \nocite*{bernanke1997systematic}, we conduct tests using three estimation models with lag orders of 12, 15, and 18. We report the minimum Wald test statistic across these models to ensure robustness. This approach ensures that if any Wald test statistic exceeds the critical value threshold in the plot, the results remain robust irrespective of the chosen lag order. The plot clearly illustrates that the uncertainty index does indeed exhibit causality toward economic activity, particularly around two and a half years, across all three sample periods and lag order selections.
This paper introduces a novel estimation and inference method for GIRs within a multi-horizon linear projection model. Our proposed two-stage estimation method with heteroskedasticity-robust inference complements the standard LS approach with HAC inference. The proposed covariance matrix estimation offers a robust alternative, eliminating the need for correcting serial correlation in the projection residuals. This method provides researchers with a valuable tool for applying multi-horizon linear projections in causality testing, impulse response estimation, or predictability analysis, particularly for those concerned about the efficiency of the LS projection method and the finite sample performance of HAC covariance estimates.
We derive uniform inference for two-stage estimates across the parameter space and extensively highlight the advantages of two-stage estimates over LS estimates from two key perspectives: (1) two-stage estimates generally exhibit greater efficiency than LS estimates across a wide range of parameter spaces and projection horizons, and (2) two-stage estimates can be implemented without relying on HAC standard error estimates. Instead, our explicit formula for the covariance matrix provides a positive semi-definite, heteroskedasticity-robust estimate. By allowing the projection horizon to grow with the sample size, we emphasize the crucial dependence of asymptotic normality of projection coefficient estimates on both the projection horizon and data persistence.
Our research on GIRs within a linear projection model underscores the importance of multi-horizon Granger causality in economics and finance for capturing the full dynamics of causal relationships. It also highlights the limitations of relying solely on Granger causality at horizon one or impulse response analysis when examining multi-horizon causality.
For future research, extending the finite-order VAR to an infinite-order VAR remains a promising area of interest from both theoretical and empirical perspectives. Inference for infinite-lag VAR models may require imposing restrictions on the norm of the lagged VAR coefficient matrices and carefully selecting a lag length that grows with the sample size, as discussed by \nocite*{lewis1985prediction}.
Another potential extension involves exploring high-dimensional VAR models instead of finite-dimensional ones, as investigated by \nocite*{basu2015regularized}, \nocite*{adamek2023lasso}, \nocite*{hecq2023granger}, and others. With datasets expanding exponentially, the exploration of high-dimensionality has gained considerable attention, particularly regarding the need for robust inference on regularized coefficients. Developing a high-dimensional version of the two-stage estimation method would thus be of significant interest for future work.