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.
62,117 characters · 17 sections · 32 citation commands
\def\spacingset#1 \spacingset{1}
\pagenumbering{arabic}\setcounter{page}{1}
\if11 {
\footnotetext[1]{D\'epartement de sciences \'economiques, Universit\'e de Montr\'eal, 3150 rue Jean-Brillant, Montr\'eal, H3T 1N8, Canada. Email: [email removed]. Web: \url{https://eugenedettaa.github.io/}.}
\footnotetext[2]{Department of Economics, University of Mannheim, L7, 3--5, Room 124, 68161 Mannheim, Germany. Email: [email removed]. Web: \url{http://www.endongwang.com}.}
} \fi
\if01 {
} \fi
{\it Keywords:} High-dimensional inference; Local projections; Neyman orthogonalization; Network connectedness; Robust testing.
{\it JEL Classification:} C13, C32, C36, C38, G12.
\spacingset{1.8}
This paper studies high-dimensional local projections (LPs) and develops Gaussian and robust inference for horizon-specific Granger-causal parameters in large dynamic systems. Local projections provide a transparent representation of multi-step dynamics and underpin a broad class of multi-horizon objects, including Granger-causality measures and network connectedness statistics\footnote{See, e.g., lutkepohl1993testing,dufour1998short,jorda2005estimation,dufour2006short,dufour2010short,dufour2012measuring,zhang2016exchange,diebold2014network,salamaliki2019transmission,diebold2023past.}. Multi-horizon Granger causality is particularly relevant for connectedness analysis because spillovers may arise with delay and persist across horizons; focusing solely on one-step-ahead causality can therefore understate systemic interdependence, even though formal Granger-causality-based connectedness measures are largely confined to a single horizon billio2012econometric,brunetti2019interconnectedness,franch2024temporal.
We depart from the existing literature along two dimensions. First, we impose sparsity only on the underlying data-generating VAR rather than separately at each projection horizon. Although LP estimands are numerically equivalent to those implied by the recursive VAR representation dufour1998short,wolf2021same, the mapping from VAR slope parameters to horizon-specific LP coefficients is highly nonlinear and involves iterated dynamic propagation; consequently, even sparse VARs may induce dense LP coefficients at longer horizons, rendering horizon-by-horizon sparsity assumptions internally inconsistent. Second, we dispense with HAC-type corrections for serial correlation in projection errors. While standard LP inference relies on such corrections, HAC variance estimators are fragile and particularly unreliable at long horizons montiel2021local,xu2023local. We instead develop inference procedures that remain valid without correcting for serial correlation.
Motivated by these considerations, we propose a two-step estimation framework for high-dimensional LPs. In the first step, we obtain regularized estimates of the underlying VAR that accommodate approximate sparsity. In the second step, we exploit closed-form mappings from VAR parameters to horizon-specific targets and apply a de-biasing correction to remove regularization-induced bias. This construction enables valid inference on multi-horizon parameters without imposing sparsity restrictions on the LP equations themselves. We illustrate the empirical relevance of the proposed methodology through multi-horizon Granger-causality tests in large financial systems, where horizon-specific restrictions yield a natural measure of system-wide connectedness.
This paper contributes to the literature on high-dimensional time series and multi-horizon Granger causality. To the best of our knowledge, we provide the first unified econometric framework for estimation and inference on horizon-specific Granger-causal effects in sparse high-dimensional VAR models.
First, building on the two-stage method of wang2024 developed for fixed-dimensional systems, we extend this framework to sparse high-dimensional VARs. Our approach relies on sparsity only at the level of the underlying VAR and does not impose sparsity on LP equations at each horizon. While existing results establish asymptotic properties in low-dimensional settings, their validity under high dimensionality has not been formally characterized. We close this gap by deriving large-sample theory for two-stage estimators when the system dimension may exceed the sample size.
Second, we develop a robust inference theory for the resulting high-dimensional procedures. Under suitable regularity conditions, we establish asymptotic Gaussian inference and construct heteroskedasticity-consistent (HC) variance estimators that remain valid without HAC corrections. This approach avoids the finite-sample distortions associated with kernel and bandwidth selection in Newey--West-type inference and is computationally tractable in high-dimensional environments.
Third, we propose a de-biasing strategy for the two-stage estimators that builds on the de-sparsification framework of van2014asymptotically and its extension to dynamic systems in adamek2023lasso. The procedure combines regularized VAR estimation with Neyman-orthogonal score constructions chernozhukov2018double, facilitating valid Gaussian inference for low-dimensional targets in high-dimensional settings.
We assess the finite-sample performance of Wald tests based on the proposed estimators through Monte Carlo experiments calibrated to empirically relevant designs. The results show that the two-stage procedure combined with HC variance estimation achieves accurate size control across all horizons, including long horizons where conventional HAC-based methods and post-double-selection approaches exhibit substantial distortions.
Finally, we illustrate the empirical relevance of our approach in an application to financial network connectedness. Focusing on major crisis episodes, the Lehman Brothers collapse, the Flash Crash, and the U.S. sovereign credit downgrade, we show that networks based on higher-horizon Granger-causality measures exhibit stronger and more persistent connectedness than networks constructed from one-step-ahead causality, highlighting the importance of longer-horizon transmission channels in assessing systemic risk.
Relevant literature. This paper contributes to the literature on regularized estimation in high-dimensional time series, including basu2015regularized,medeiros2016l1,wong2020lasso,masini2022regularized,adamek2023lasso. We depart from this line of work by developing de-biased estimation and inference procedures for parameters in LP equations constructed from regularized VAR estimators. Our analysis builds on the desparsified estimation literature belloni2012sparse,van2014asymptotically,chernozhukov2018double,krampe2023structural and complements recent work on single-horizon Granger causality in large systems hecq2023granger,babii2024high. In contrast to adamek2024local, who study de-biased inference for Sims impulse responses in high-dimensional LP models, we focus on multi-horizon Granger-causal parameters, following the conceptual distinction of dufour1998short. We further extend the literature on robust inference in LP models montiel2021local,breitung2023projection,xu2024new,wang2024 to high-dimensional settings.
Outlines. Section (ref) outlines the econometric framework. Section (ref) presents the de-biased two-stage estimation method. Asymptotic Gaussian and robust inference results are derived in Section (ref). Monte Carlo evidence is reported in Section (ref), and an empirical application is provided in Section (ref). Section (ref) concludes. Proofs are collected in the Online Appendix.
Notation. Throughout, $C$ denotes a generic positive constant whose value may change from line to line, and $v_j$ denotes a conformable selection vector with one in the $j$-th position and zeros elsewhere. For a vector $x$, $\|x\|_1$ and $\|x\|_2$ denote the $\ell_1$ and $\ell_2$ norms, respectively. For a matrix $B$, $\|B\|_1$, $\|B\|_2$, and $\|B\|_\infty$ denote the induced matrix norms, and $\|B\|_{\max}$ the elementwise maximum norm. For a symmetric matrix, $\lambda_{\min}(\cdot)$ and $\lambda_{\max}(\cdot)$ denote its smallest and largest eigenvalues.
Granger causality assesses whether the past of one variable improves the linear prediction of another beyond its own history granger1969investigating,geweke1984inference. Extending the concept to multi-step forecasts yields multi-horizon Granger causality lutkepohl1993testing,dufour1998short. Let $w_t=(y_t,x_t,q_t')'$, where $x_t$ and $y_t$ are scalar variables of interest and $q_t$ collects a possibly high-dimensional set of controls. Variable $x$ does not Granger-cause $y$ at horizon $h$ if
where $\textup{P}_L(\cdot\mid\cdot)$ denotes linear projection, $\mathcal{F}_{t}=\sigma(w_t,w_{t-1},\ldots)$, and $\mathcal{F}_{-x,t}$ excludes the history of $x$.
Suppose that $w_t$ follows a high-dimensional $\operatorname{VAR}(p)$, \[ w_t = A_1 w_{t-1} + \cdots + A_p w_{t-p} + u_t, \] where $\{u_t\}$ is a martingale difference sequence with zero mean and nonsingular covariance matrix $\Sigma_u$, satisfying $\lambda_{\min}(\Sigma_u)>0$. The lag order $p$ is fixed, while the dimension $d$ may grow with the sample size $n$.
Under the VAR, the linear projection of $y_{t+h}$ on $\mathcal F_t$ admits the local projection representation\footnote{See dufour1998short,wolf2021same.} $\textup{P}L(y_{t+h}\mid\mathcal F_t)=\beta_h' W_t$, where $W_t=(w_t',w_{t-1}',\ldots,w_{t-p+1}')'$. The coefficient vector $\beta_h$ corresponds to the row of the $h$-th power of the VAR companion matrix $\mathbf A$ associated with $y$. Here, $\mathbf A$ denotes the $dp\times dp$ companion matrix constructed from the VAR slope matrices ${A_1,\ldots,A_p}$, so multi-horizon predictability is governed by dynamic propagation through $\mathbf A^h$ (see, e.g., kilian2017structural).
Partition $W_t=(W_{1,t}',W_{2,t}')'$, where $W_{1,t}=(x_t,x_{t-1},\ldots,x_{t-p+1})'$ and $W_{2,t}$ contains all remaining predictors. Then
where $e_{t,h}$ is the projection error, which may be serially correlated. The null hypothesis of Granger non-causality at horizon $h$ is $\mathcal H_0:\beta_{1,h}=0$. Valid inference therefore requires estimating $\beta_{1,h}$ in the presence of the high-dimensional nuisance parameter $\beta_{2,h}$, even though sparsity of the VAR need not imply sparsity of (ref).
At horizon $h=1$, post--double-selection methods can be used to construct de-biased estimators under the sparsity of the VAR and auxiliary regressions; see, for example, hecq2023granger. At longer horizons, however, feasibility typically requires approximate sparsity of the local projection coefficients themselves babii2024high, an assumption that is generally untenable in dynamic systems.
The key difficulty is that $\beta_h$ is determined by powers of the VAR companion matrix $\mathbf A$, which are nonlinear functions of the primitive VAR coefficients. As a result, sparse short-run VAR dynamics can generate dense horizon-specific local projection coefficients as indirect transmission paths accumulate over time. To illustrate this mechanism, consider a $d$-dimensional VAR(1) with coefficient matrix $A$ satisfying
Then $A$ is a sparse band matrix and is stable since $\rho(A)=|\alpha|<1$. However, the horizon-$h$ coefficients correspond to $A^h$, whose support expands monotonically with $h$ as longer transmission paths become active. In particular,
where $\|\cdot\|_0$ denotes the cardinality of the support of a vector. Thus, even when the VAR is sparse, the implied local projection coefficients become dense as $h$ grows. This example demonstrates that sparsity of the VAR does not propagate to local projections, rendering horizon-by-horizon sparsity assumptions inappropriate.
Existing Granger-causality-based measures of network connectedness with formal testing are largely confined to a single horizon billio2012econometric,brunetti2019interconnectedness. By contrast, diebold2014network study multi-horizon connectedness using forecast error variance decompositions without explicit inference. A multi-horizon Granger-causality framework with formal testing therefore bridges these two strands of the literature, while allowing for economically meaningful horizon choices that range from short-term risk monitoring to longer-run spillovers.
Let $p_{ij}(h)$ denote the $p$-value for testing the null that variable $j$ does not Granger-cause variable $i$ at horizon $h$. For a given confidence level $\alpha\in(0,1)$, define
and summarize system-wide connectedness by the Degree of Granger Causality,
which measures the density of statistically significant predictive links and takes values in $[0,1]$. Aggregating pairwise tests in this way, DGC is not interpreted as a familywise error--controlled count of causal links, but as a relative measure of network connectedness across horizons and subsamples. The use of high confidence thresholds mitigates concerns about false discovery accumulation, and the qualitative patterns are stable across alternative thresholds, indicating systematic volatility spillovers rather than noise-driven link formation.
This section proposes a two-stage procedure for estimating the horizon-$h$ Granger causal parameters.
Suppose the reduced-form innovations $u_t$ satisfy a weak exogeneity condition and the innovation covariance matrix $\Sigma_u$ is nonsingular. Let $U_t=(u_t',u_{t-1}',\ldots,u_{t-p+1}')'$ denote the stacked innovation vector. Under these conditions, the orthogonality restriction $\mathrm P_L(y_{t+h}-\beta_h' W_t \mid U_t)=0$ holds, implying that $U_t$ is a valid instrument set for $W_t$. Identification therefore follows from the moment condition
where the population moments admit closed-form expressions, $\mathbb E[U_t W_t']=(I_p\otimes\Sigma_u)\Psi(p)$ and $\mathbb E[U_t y_{t+h}] =\mathbb E[U_t w_{t+h}']v_1 =(I_p\otimes\Sigma_u) (\Psi_h',\Psi_{h+1}',\ldots,\Psi_{h+p-1}')'v_1$.\footnote{Since $y_t$ is the first element of $w_t$, $y_t=v_1'w_t$.} Here $\Psi(p)$ is a $p\times p$ block matrix whose $(i,j)$ block equals $\Psi_{i-j}'$ for $i\ge j$ and zero otherwise, and $\Psi_h=J\mathbf A^h J'$. Since $\Sigma_u$ is full rank, $\mathbb E[U_t W_t']$ is nonsingular, ruling out under-identification, as in the static IV case.
To isolate the coefficient of interest, partition the instruments conformably with the regressors. Define $U_{1,t}=R_1 U_t$ and $U_{2,t}=R_2 U_t$, where the selection matrices $R_1$ and $R_2$ coincide with those used to define $W_{1,t}=R_1 W_t=(x_t,x_{t-1},\ldots,x_{t-p+1})'$, and $W_{2,t}=R_2 W_t$ includes the rest of variables. Consequently, $U_{1,t}$ collects the VAR innovations corresponding to $W_{1,t}$, and analogously for $U_{2,t}$ and $W_{2,t}$. The Frisch--Waugh--Lovell theorem yields
where the residualized instrument $U_{1,t}^{\perp}=U_{1,t}-\Gamma U_{2,t}$, with $\Gamma=\mathbb E[U_{1,t} W_{2,t}'] \mathbb E[U_{2,t} W_{2,t}']^{-1}$, satisfies the orthogonality condition $\mathrm P_L(U_{1,t}^{\perp}\mid W_{2,t})=0$. By construction, $U_{1,t}^{\perp}$ remains relevant for $W_{1,t}$ while being orthogonal to the control regressors $W_{2,t}$.
This section develops a de-biased two-stage estimator for the multi-horizon Granger-causal coefficient $\beta_{1,h}$ in the local-projection framework. The construction is grounded in the population representation in (ref). Identification is achieved through the residualized instrument $U_{1,t}^{\perp}$, which partials out high-dimensional controls from the raw instrument stack and isolates variation orthogonal to the nuisance component.
Let $\Sigma_{UW}:=\mathbb E[U_t W_t']$ denote the population cross-moment matrix between the instrument vector $U_t$ and the regressor vector $W_t$. This matrix captures the linear dependence between innovations and lagged regressors and enters the definition of the partialling-out coefficient $\Gamma$. Using the population expression for $\Gamma$, the residualized instrument admits the closed-form representation
which shows that $U_{1,t}^{\perp}$ depends only on second moments of $(U_t,W_t)$.
To implement (ref), we exploit that $\Sigma_{UW}$ admits a closed-form representation in terms of the reduced-form impulse responses $\{\Psi_h\}_{h\ge 0}$ and the innovation covariance matrix $\Sigma_u$. We estimate the VAR using a regularized estimator $\widehat{\mathbf A}_{1:p}^{(\mathrm{re})}$, obtained by applying thresholding to Lasso-type estimates of the VAR slope coefficients, $\widehat{\mathbf A}_i^{(\mathrm{re})} = T_{\tau_A}(\widehat{\mathbf A}_i^{(\mathrm{lasso})})$, with $\tau_A\asymp(\log(dp)/T)^{1/2}$, which ensures boundedness of the column-wise $\ell_1$ norm as required by Assumption (ref)(iii). Based on these estimates, we form the fitted innovations $\hat u_t := w_t-\widehat{\mathbf A}_{1:p}^{(\mathrm{re})}W_{t-1}$ and the sample covariance estimator $\widehat{\Sigma}_u := n^{-1}\sum_{t=p+1}^{n}\hat u_t\hat u_t'$. To stabilize covariance estimation in high dimensions under sparsity of $\Sigma_u$ (Assumption (ref)), we apply entrywise thresholding and define $\widehat{\Sigma}_u^{(\mathrm{re})}:=T_{\tau_u}(\widehat{\Sigma}_u)$, where $T_{\tau_u}(\cdot)$ sets entries with absolute value below $\tau_u$ to zero. In theory, $\tau_u$ is chosen on the order of $\|\widehat{\Sigma}_u-\Sigma_u\|_{\max}$ (see Section (ref)); in applications, we set $\tau_u\asymp(\log d/T)^{1/2}$.
Impulse responses are estimated from the regularized VAR. Let $\widehat{\mathbf A}^{(\mathrm{re})}$ denote the associated companion matrix and define $\hat\Psi_h := J(\widehat{\mathbf A}^{(\mathrm{re})})^h J'$. Substituting $\{\hat\Psi_h\}$ and $\widehat{\Sigma}_u^{(\mathrm{re})}$ into the population representation yields the plug-in estimator
Given $\widehat{\Sigma}_{UW}^{(\mathrm{re})}$, the two-stage estimator of $\beta_{1,h}$ is obtained by instrumenting $W_{1,t}$ with residualized instruments $\hat U_{1,t}^\perp$,
where $\hat U_{1,t}^{\perp} = \bigl(R_1(\widehat{\Sigma}_{UW}^{(\mathrm{re})})^{-1}R_1'\bigr)^{-1} R_1(\widehat{\Sigma}_{UW}^{(\mathrm{re})})^{-1}\hat U_t, $ $ \hat U_t=(\hat u_t',\ldots,\hat u_{t-p+1}')'.$ In high-dimensional settings, $\hat \beta_{1,h}^{(2S)}$ is generally not $\sqrt n$-consistent, since the moment conditions are not Neyman-orthogonal with respect to the high-dimensional nuisance parameter $\beta_{2,h}$.
We therefore apply a de-biasing correction and define
where $\hat \beta_{2,h}$ is a consistent estimator of $\beta_{2,h}$ with controlled bias. A convenient choice is the plug-in estimator from the regularized VAR, $\hat \beta_{2,h}=J(\widehat{\mathbf A}^{(\mathrm{re})})^h$. This correction removes the leading bias induced by projection onto $W_{2,t}$ and restores first-order orthogonality, enabling valid inference.
\noindentRemarks.
We next derive feasible inference for the de-biased two-stage estimator by characterizing its asymptotic distribution. The argument proceeds by isolating the leading stochastic component in the $\sqrt{n}$-normalized estimation error and verifying that all remaining terms are asymptotically negligible. We will show that the asymptotic behavior of the debiased two-stage estimator is driven by a martingale-type sum involving the instrument $U_{t}$ and the projection error $e_{t,h}$. Importantly, all effects stemming from first-stage regularization enter only through higher-order remainder terms and therefore do not affect the limiting distribution.
The asymptotic variance of $\hat{\beta}_{1,h}^{(\mathrm{de\text{-}2S})}$ is given by
where $\Omega_{U,h} := \lim_{n\to\infty} \operatorname{Var}\!\left( n^{-1/2}\sum_t U_{t} e_{t,h} \right)$ denotes the long-run variance of the score process.
For practical implementation, we construct a feasible estimator of the asymptotic variance using a HAC procedure,
where $\hat\Omega_{U,h}^{(hac)}$ is a consistent HAC estimator of $\Omega_{U,h}$, such as the Newey--West estimator. In practice, the latent instrument $U_{t}$ is replaced by its sample analogue obtained from the sample analogue defined in (ref). Moreover, since the projection error $e_{t,h}$ is unobservable, we replace it by the residual $\hat e_{t,h}=y_{t+h}-\hat\beta_h W_t$, where $\hat\beta_h=v_1'(\widehat{\mathbf A}^{(\mathrm{re})})^h$. The resulting plug-in score $\hat U_{t}\hat e_{t,h}$ is then used to construct the HAC covariance estimator. Under standard regularity conditions, estimation error from the first-stage regularization and the residualization step enters only through higher-order terms and therefore does not affect first-order asymptotics. This ensures the validity of the proposed HAC-based inference.
HAC-type covariance estimators are well known to perform poorly in finite samples, particularly in long-horizon local projection settings, where they often induce substantial size distortions. Their implementation further depends on nontrivial kernel and bandwidth choices. These considerations motivate alternative covariance estimators based solely on heteroskedasticity-robust methods.
Replacing HAC estimation requires conditions under which serial dependence in the score process can be eliminated. HAC estimators are primarily used to obtain a positive semidefinite estimate of the long-run variance, a property that is generally not preserved by naive aggregation of autocovariances. An alternative approach is to construct a transformation of the score process that is serially uncorrelated, so that the long-run variance coincides with the contemporaneous covariance matrix. The long-run variance of the original score process is the summation of all lead-lag autocovariances,
The score process $U_t e_{t,h}$, where $U_t=(u_t',u_{t-1}',\ldots,u_{t-p+1}')'$, is serially correlated, which motivates the use of HAC corrections in standard inference. Following wang2024, we instead construct an alternative score sequence by re-indexing and stacking the score components as
This transformation aligns all score components with the contemporaneous innovation $u_t$ and shifts serial dependence forward in time.
Under Assumption (ref), the transformed process $\{s_t\}$ is serially uncorrelated. Consequently, the long-run variance of the original score process coincides with the contemporaneous covariance matrix of $s_t$, namely
This identification ensures that the asymptotic variance appearing in the limit distribution is consistently estimated by a heteroskedasticity-robust covariance estimator, thereby obviating the need for HAC corrections.
Assumption (ref)(i) imposes a martingale-difference structure on the innovation process. Assumption (ref)(ii) is a fourth-order orthogonality condition: it requires that, for any $t<\tau$ and any forward offset $k>0$, the cross-product $u_tu_\tau'$ is orthogonal (in the unconditional fourth moment) to the future quadratic form $u_{\tau+k}u_{\tau+k}'$. This restriction is imposed for variance estimation: it rules out precisely the type of intertemporal fourth-moment dependence that would generate nonzero autocovariances in the transformed score sequence $\{s_t\}$, so that a heteroskedasticity-consistent variance estimator based only on contemporaneous score variation is valid.
Importantly, Assumption (ref)(ii) does not require serial independence of $\{u_t\}$ and does not rule out conditional heteroskedasticity or volatility clustering. It allows the conditional covariance matrix $\mathbb E[u_tu_t' \mid \{u_s\}_{s<t}]$ to vary over time, as long as this variation does not induce the above cross--fourth-moment dependence across nonoverlapping time blocks that would transmit serial dependence to $\{s_t\}$. The condition is satisfied by a broad class of disturbance processes, including i.i.d.\ shocks, mean-independent disturbances, conditionally homoskedastic martingale differences, and conditionally Gaussian ARCH-type models (including Gaussian GARCH specifications) in which the innovation is conditionally symmetric and the relevant cross-products are orthogonal to future quadratic terms.
By contrast, Assumption (ref)(ii) can be violated in environments where second-moment dynamics feed back into future quadratic variation, e.g., through leverage-type effects, which may induce serial dependence in the score sequence and thus require HAC-type corrections. In the present framework, Assumption (ref) is imposed to rule out this source of score dependence. In particular, it implies that the transformed score $s_t$ is serially uncorrelated. The following proposition formalizes this implication.
See Online Appendix for the proof. Accordingly, a heteroskedasticity-robust estimator of the asymptotic covariance matrix for the de-biased two-stage estimator is given by
where $\widehat{\mathrm{Var}}(\hat s_t) = \frac{1}{n-h}\sum_t \hat s_t \hat s_t'$, $\hat s_t = (\hat e_{t,h},\hat e_{t+1,h},\ldots,\hat e_{t+p-1,h})\otimes\hat u_t$, $\hat u_t = w_t-\hat\Phi_{1:p}^{(re)}W_{t-1}$. Here, $\hat e_{t,h}$ denotes the local-projection residual defined in (ref). Under standard regularity conditions, the estimator in (ref) is consistent for the asymptotic variance of $\sqrt n\,\hat\beta_{1,h}^{(\mathrm{de\text{-}2S})}$.
This section studies the large-sample behavior of the proposed de-biased two-stage (de-2S) estimator. We proceed in two steps. First, we impose high-level assumptions for the regularized VAR estimator. Second, we establish asymptotic normality of the de-2S estimator and consistency of feasible variance estimators. Both HAC and heteroskedasticity-robust standard errors are covered.
We begin with conditions used in the preliminary consistency results. To justify Lasso-type regularization for the VAR slope matrices $\widehat{\mathbf A}^{(\mathrm{re})}_j$, $j=1,\ldots,p$, and their stacked form $\widehat{\mathbf A}^{(\mathrm{re})}$, we impose approximate sparsity on the VAR companion matrix. Following krampe2023structural and bickel2008covariance, define the row-wise approximately sparse class
This class includes exact sparsity as $\mu=0$ (interpreting $\sum_{j=1}^s |b_{ij}|^\mu$ as the number of nonzero entries in row $i$) and allows many small coefficients when $\mu\in(0,1)$, while controlling effective row complexity through $k$.\footnote{In our application, $k$ may depend on $(d,p)$ and is allowed to grow with $n$.}
Rates for $\ell_1$-regularized estimators in high-dimensional time series and VARs are available under conditions of this type; see, e.g., basu2015regularized,adamek2023lasso. We summarize the required inputs in a high-level assumption.
Assumption (ref)(i)--(ii) imposes approximate sparsity and stability of the $dp\times dp$ companion matrix. The sparsity index $k_A$ may grow with $(d,p)$, while stability yields geometric decay of $\mathbf A^m$ and standard weak-dependence properties for the stacked state.
Assumption (ref)(iii) is a high-level rate for the regularized estimator in $\|\cdot\|_1$ and $\|\cdot\|_\infty$, which are convenient for bounding remainder terms. The factor $\nu_n$ collects dimension and tail effects driving regularization error; for instance, one may take $\nu_n=\log(dp)$ under sub-Gaussian tails, while heavier tails typically lead to larger $\nu_n$.
Assumption (ref)(iv) imposes a uniform concentration condition on quadratic forms of the sample covariance estimator of $u_t$. Equivalently, this requirement corresponds to an operator-norm concentration bound for $\widehat{\Sigma}_u-\Sigma_u$. The condition does not impose sparsity on $\Sigma_u$ itself and its effective dimensional complexity is summarized by the factor $\tilde{\nu}_n$. Assumption (ref)(v) provides a mixing and moment condition sufficient for central limit and law-of-large-numbers arguments for the quantities driving the de-biasing step. Assumption (ref)(vi) ensures that $\Sigma_u$ is well conditioned and approximately sparse, which supports thresholding-type estimation of $\Sigma_u$.
Assumption (ref)(vii) controls the growth of $|\Sigma_{UW}^{-1}|\infty$, which governs the sensitivity of the de-biasing correction to estimation error in $\Sigma{UW}$. This condition is imposed purely for inferential stability and is standard in high-dimensional de-biasing arguments; see, for example, Assumption 2(v) in krampe2023structural. Finally, Assumption (ref)(viii) ensures finite fourth moments and uniform conditioning of the long-run variance $\Omega_{U,h}$; together with (v), it underpins consistency of both HAC and heteroskedasticity-robust variance estimators used for feasible inference.
This subsection establishes asymptotic normality and feasible inference for the de-biased two-stage estimator. The argument has two ingredients. First, we derive high-dimensional consistency rates for the covariance and cross-covariance estimators that enter the de-biasing correction. Second, we impose growth conditions under which the resulting higher-order remainder terms are asymptotically negligible, so that the studentized statistic admits a standard Gaussian limit.
Under Assumption (ref), we obtain the following bounds for $\widehat{\Sigma}_u$, $\widehat{\Sigma}_u^{(\mathrm{re})}$ and $\widehat{\Sigma}_{UW}^{(\mathrm{re})}$.
See Online Appendix for the proof. The rate $\delta_n$ collects two components. The first is the baseline high-dimensional sampling error $\big\|n^{-1}\sum_{t=1}^n u_t u_t' - \Sigma_u\big\|_{\max} = O_p(\sqrt{\tilde{\nu}_n/n})$ from Assumption (ref)(iv). The second is the plug-in error induced by estimating the VAR dynamics via $\widehat{\mathbf A}^{(\mathrm{re})}$, which affects fitted innovations (and hence $\widehat{\Sigma}_u^{(\mathrm{re})}$) and also propagates into the estimated impulse responses used to form $\widehat{\Sigma}_{UW}^{(\mathrm{re})}$. The second is sampling variability in covariance estimation, which would remain even if the innovations $\{u_t\}$ were observed; this component is captured by Assumption (ref)(iv). In particular, $\delta_n$ collects both (a) the bias induced by using $\hat u_t$ in place of $u_t$ and (b) the high-dimensional sampling error of the sample covariance of $\{u_t\}$. The bound in part (iii) then combines the regularization error in $\widehat{\Sigma}_u^{(\mathrm{re})}$ with the approximation error in the impulse response coefficients $\hat{\Psi}_h$ constructed from powers of the regularized transition matrix $(\widehat{\mathbf A}^{(\mathrm{re})})^h$.
As a benchmark, suppose $k_A$ and $k_U$ are fixed, the VAR slope matrices and $\Sigma_u$ are exactly sparse (i.e., $\mu=\mu_u=0$), and the innovation process has sub-Gaussian tails so that $\nu_n$ and $\tilde{\nu}_n$ are of order $\log(d)$. Then both $\|\widehat{\Sigma}_u^{(\mathrm{re})}-\Sigma_u\|_{\infty}$ and $\|\widehat{\Sigma}_{UW}^{(\mathrm{re})}-\Sigma_{UW}\|_{\infty}$ achieve the rate $O_p\big(\sqrt{\log(d)/n}\big)$.
Condition (ref) summarizes the growth restrictions needed for (a) the $\sqrt{n}$-normalized de-biasing representation to be asymptotically linear and (b) the plug-in standard error to be consistent. Part (i) enforces that the dominant first-stage nuisance error entering the de-biasing correction vanishes after scaling by $k_W$ and the stochastic fluctuation level $\tilde\nu_n^{1/2}$; equivalently, it implies $k_W \tilde\nu_n^{1/2}\|\widehat{\Sigma}_{UW}^{(\mathrm{re})}-\Sigma_{UW}\|_{\infty} =o_p(1)$ up to the rate in Lemma (ref)(iii). Parts (ii) and (iii) ensure that the estimation error in the asymptotic variance is negligible. These restrictions control the effect of using $\widehat{\Sigma}_{UW}^{(\mathrm{re})}$ and the feasible long-run variance estimators $\widehat{\Omega}_{U,h}^{(hac)}$ or $\widehat{\Omega}_{U,h}^{(hc)}$ in place of their population counterparts. In the benchmark case discussed above with fixed $k_W$, the conditions reduce to standard requirements such as $\log(d)/\sqrt{n}\to 0$ (and analogous rate restrictions implied by the remaining terms).
We now state the main inference result.
See Online Appendix for the proof. The HAC standard error is robust to general serial dependence in the innovation process under the mixing and moment conditions in Assumption (ref). As an alternative, the heteroskedasticity-robust standard error in (ref) avoids long-run variance estimation but requires the stronger conditions on $\{u_t\}$ summarized in Assumptions (ref).
This section reports a Monte Carlo study assessing the finite-sample performance of the proposed de-biased two-stage estimation and inference procedures for high-dimensional local projections. The experiments focus on the size properties of Wald tests for multi-horizon Granger-causal parameters under empirically relevant departures from sparsity.
We consider vector autoregressive models of order $p=2$. Stationarity is imposed via a factorization of the characteristic polynomial: two root matrices $\{\Lambda_k\}_{k=1}^2$ are constructed and the VAR slope matrices are recovered as $A_1=\Lambda_1+\Lambda_2$ and $A_2=-\Lambda_1\Lambda_2$, ensuring that all eigenvalues of the companion matrix lie strictly inside the unit circle. The innovations $\{u_t\}$ are i.i.d.\ Gaussian, $u_t\sim\mathcal N(0,\Sigma_u)$, with $\Sigma_{u,ij}=0.5^{|i-j|}$. We consider $(d,T)\in\{(60,120),(60,240)\}$. Under the VAR$(2)$ specification, the number of regressors in the associated local projection equations is of the same order as the sample size, rendering conventional OLS-based inference invalid.
Two designs for the root matrices are considered. In the first design, the roots are tridiagonal, with $\Lambda_{ij,k}=\rho^{|i-j|+1}$ for $|i-j|\le q$, where $q=3$ and $\rho=0.3$, and zero otherwise; the largest eigenvalue of the VAR companion matrix equals $0.549$. In the second design, the roots are upper triangular, with $\Lambda_{ij,k}=0.1$ if $0\le j-i\le q$ for $q=3$, and zero otherwise; the corresponding largest eigenvalue of the VAR companion matrix equals $0.147$.
Both designs are sparse at horizon one, but sparsity deteriorates rapidly as the projection horizon increases due to dynamic propagation. In the tridiagonal design, repeated matrix multiplication activates an expanding set of small but non-negligible coefficients, rendering local projection equations increasingly dense. In the upper-triangular design, the one-sided propagation structure induces even faster support expansion, with the number of nonzero coefficients in the $h$-step companion matrix growing at rate proportional to $qh$. As a result, approximate sparsity breaks down at moderate horizons in both designs.
We conduct $1{,}000$ Monte Carlo replications. In each replication and for each horizon $h=1,\ldots,24$, we compare four estimators: post--double-selection LASSO with HAC inference; de-biased LASSO with HAC inference; the proposed de-biased two-stage estimator with HAC inference; and the proposed de-biased two-stage estimator with heteroskedasticity-robust inference. For all procedures relying on HAC inference, long-run variance matrices are estimated using a Bartlett kernel with bandwidth equal to the projection horizon $h$. Regularization parameters are selected using BIC-type criteria.
Figures (ref)--(ref) report empirical rejection frequencies of 5% Wald tests for the two designs, computed as one minus the empirical coverage probabilities of the corresponding 95% confidence sets. Left (right) panels correspond to $(d,T)=(60,120)$ ($(60,240)$). Ejection frequencies are based on the Wald statistic
where $\hat{\boldsymbol\beta}_{1,h}$ denotes the estimated vector of Granger-causal coefficients at horizon $h$. The covariance matrix $\hat{\Omega}_h$ is the corresponding asymptotic variance estimates. Under the null hypothesis, $W_h \xrightarrow{d} \chi^2_q$, where $q=2$ is the number of restrictions in our VAR(2) simulations.
Both de-biased LASSO and post--double-selection LASSO exhibit substantial size distortions across horizons, reflecting violations of the sparsity conditions required for valid de-biasing and variable selection. Although the VAR dynamics are sparse at short horizons, local projection equations become progressively dense as the horizon increases, placing these methods outside their intended high-dimensional regimes.
Across both designs and sample sizes, the proposed de-biased two-stage estimator combined with heteroskedasticity-robust inference delivers accurate and stable size control, with rejection frequencies remaining close to the nominal level uniformly across horizons. In contrast, procedures relying on HAC variance estimation display increasing distortions at longer horizons, consistent with documented finite-sample failures of HAC estimators in long-horizon local projections montiel2021local,xu2023local,wang2024. Increasing the sample size mitigates distortions for all methods. Overall, the results indicate that the proposed de-biased two-stage procedure provides reliable finite-sample inference for multi-horizon Granger non-causality testing in high-dimensional VARs.
We apply the proposed methodology to study volatility transmission and network connectedness in U.S.\ equity markets. A large literature documents that shocks to firm-level volatility propagate across assets and sectors, generating system-wide risk interdependence; see, among others, diebold2014network,mcaleer2008realized,hecq2023granger,miao2023high. Our contribution is to characterize how such spillovers unfold across multiple forecast horizons in a high-dimensional setting.
We use daily realized variance data for 30 large-cap U.S.\ equities, corresponding to constituents of the Dow Jones Industrial Average and spanning the main sectors of the U.S.\ economy.\if11 {\footnote{We thank a colleague for providing the 10-minute realized variance data used in hecq2023granger.} } \fi \if01 {\footnote{We thank a colleague for providing the data used in this study.} } \fi Realized variances are constructed from intraday 10-minute returns and log-transformed to stabilize dispersion. The sample spans March 2008 to February 2017, comprising 2{,}236 trading days.
We construct directed networks of predictive relations using multi-horizon Granger-causality tests. For each horizon $h\in\{1,5,10,20\}$ and each ordered pair $(i,j)$ with $i\neq j$, we test the null that asset $j$ does not Granger-cause asset $i$ at horizon $h$. A directed link is recorded whenever the null is rejected at a fixed significance level. The resulting adjacency matrix captures the extensive margin of connectedness. System-wide connectedness is summarized by the degree of Granger causality $\mathrm{DGC}(h;\alpha)$ defined in (ref), evaluated at $\alpha=0.99$ and $\alpha=0.999$.
Figures (ref) and (ref) report rolling-window estimates of $\mathrm{DGC}(h;\alpha)$ based on 100-day windows. We focus on the period from August 2008 to February 2012, covering three major market events: the Lehman Brothers bankruptcy (2008--09--15), the Flash Crash (2010--05--06), and the S&P sovereign downgrade (2011--08--05), indicated by vertical lines. Within each window, we estimate a sparse VAR(4) via adaptive LASSO and conduct horizon-specific Granger-causality tests.
Two main patterns emerge. First, volatility connectedness is highly state-dependent and rises sharply during periods of market stress. All three events coincide with pronounced spikes in $\mathrm{DGC}(h;\alpha)$, indicating a rapid densification of the volatility transmission network. Tightening the significance level from $\alpha=0.99$ to $\alpha=0.999$ sharpens these spikes, isolating episodes of acute systemic connectedness, while looser thresholds capture more persistent but weaker dependence.
Second, the horizon dimension is central to identifying economically meaningful spillovers. One-day-ahead connectedness responds only modestly around crisis events, whereas multi-day horizons exhibit large and persistent increases. This suggests that short-horizon measures understate the importance of slower-moving propagation mechanisms, while longer horizons capture cumulative volatility transmission. Notably, a pronounced one-day spike in late 2010 is not associated with a major systemic event and disappears at longer horizons, underscoring the value of multi-horizon analysis for filtering spurious signals.
To visualize these dynamics, Figures (ref)--(ref) plot directed networks before and after each event at the one-day and two-week horizons. Firms are grouped by sector and color-coded accordingly. While one-day networks change little around events, two-week networks become markedly denser following each shock, indicating broad-based amplification of volatility spillovers across sectors. Our results highlight that equity-market volatility transmission operates significantly through multi-day propagation channels. Crisis episodes are characterized by sharp increases in predictive interdependence, with individual stocks becoming more informative about future volatility dynamics of others. Monitoring horizon-specific Granger-causal linkages therefore provides a financially meaningful and complementary measure of market connectedness.
We develop estimation and inference methods for multi-horizon Granger causality in high-dimensional dynamic systems. A key insight is that sparsity of the underlying VAR representation does not imply sparsity of local-projection coefficients at horizons $h>1$, because powers of the VAR transition matrix typically produce dense horizon-specific LP coefficients. This disconnect undermines direct application of high-dimensional regularization and post-selection methods to LP regressions.
We exploit sparsity at the VAR level and derive explicit analytic mappings to recover horizon-specific LP coefficients. Based on Neyman-orthogonal scores, we construct a de-biased two-stage estimator and establish asymptotic normality for low-dimensional multi-horizon Granger-causal parameters with growing system dimension. Our Wald tests are valid under approximate sparsity and weak dependence even as the number of variables increases with the sample size.
We also propose a heteroskedasticity-robust inference approach that obviates long-run variance estimation. By re-indexing score equations and imposing mild structural conditions on the innovation process, this route eliminates HAC corrections while preserving first-order validity. Monte Carlo evidence shows stable size control and improved finite-sample performance at long horizons relative to conventional HAC inference.
An empirical application to realized-volatility networks demonstrates that horizon-specific connectedness measures reveal propagation patterns that one-step analyses miss, particularly around market stress episodes. The findings underscore the value of multi-horizon causality measures for characterizing dynamic interdependence in large financial systems.
\if11
We are extremely grateful to Marine Carrasco, Benoit Perron, Mathieu Marcoux, Jean-Marie Dufour, Carsten Trenkler, and Victoria Zinde-Walsh for helpful discussions and guidance. We also thank Ren\'e Garcia, Prosper Dovonon, and participants in the CIREQ Econometrics Conference in Honor of Eric Ghysels, the 2024 NBER--NSF Time Series Conference, the International Association for Applied Econometrics 2025 Conference, and the Dagenais Econometrics Seminars for constructive comments. This work was supported by the Fonds de recherche du Qu\'ebec -- Soci\'et\'e et culture (FRQSC) and the research funds provided by the University of Mannheim. \fi
\pagenumbering{arabic}\setcounter{page}{1}