EconBase
← Back to paper

Inference for High-Dimensional Local Projection

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.

88,069 characters · 11 sections · 51 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

{8pt} {8pt}

titlepage\begin{center} \begin{spacing}{1.5} { Inference for High-Dimensional Local Projection} \end{spacing} {\sc Jiti Gao$^{\dag}$, Fei Liu$^\sharp$ and Bin Peng$^{\dag}$} $^\dag$Monash University and $^\sharp$University of Bath \today \end{center} \begin{abstract} This paper rigorously analyzes the properties of the local projection (LP) methodology within a high-dimensional (HD) framework, with a central focus on achieving robust long-horizon inference. We integrate a general dependence structure into $h$-step ahead forecasting models via a flexible specification of the residual terms. Additionally, we study the corresponding HD covariance matrix estimation, explicitly addressing the complexity arising from the long-horizon setting. Extensive Monte Carlo simulations are conducted to substantiate the derived theoretical findings. In the empirical study, we utilize the proposed HD LP framework to study the impact of business news attention on U.S. industry-level stock volatility. Keywords: high-dimensional local projection; long-horizon analysis; $h$-step ahead forecasting models; covariance matrix estimation JEL Classification: C32, C53, C55 \end{abstract}

Introduction

The local projection (LP) method, introduced in the seminal work of Oscar_2005, has garnered significant attention in macroeconomics and econometrics. We refer interested readers to JT2025 for the latest comprehensive literature review.

A key advantage of the LP method, which has led to its widespread adoption and numerous extensions, is its simplicity: impulse responses can be recovered via a set of simple least squares regressions, bypassing traditional short-run, long-run, or sign restrictions. Consequently, researchers have intuitively applied different generalizations of the least squares regression to the LP method. The literature has thus expanded to include extensions concerning nonlinear structures, robustness estimation, and such models associated with high-dimensional (HD) regressors (covariates), just to name a few. Given the vastness of this literature, we acknowledge the limitation of the following literature review, and summarize only the key articles that directly motivate our research, after which we situate our specific contribution within this field.

The first stream we review concerns the robustness and inference of the LP method. For example, montiel2021local establish long-horizon inference and provide a practical bootstrap method for conducting it. xu2023local extends this work, demonstrating the efficiency of the LP method within a vector autoregression of infinite order (VAR($\infty$)) framework. Furthermore, both herbst2024bias and mei2023nickell investigate the biases associated with the LP method under different model specifications, drawing attention to its application in panel data studies. Given the inherent connection between panel and HD data, a group of studies focuses on the application of the LP method to HD settings. Recently, ASW2024 apply the desparsified least absolute shrinkage and selection operator (LASSO) to carry on the HD LP approach, while leaving the impulse response parameter of interest unpenalized; cha2024 considers a similar setting without specifically relying on sparsity assumptions, and establishes some useful concentration inequalities. Methodologically, basu2015regularized and MIAO2023155 provide useful insights by analyzing specific HD--VAR models; these serve as essential benchmarks for investigating the LP approach in certain HD settings.

Up to this point, we emphasize the necessity of integrating the LP approach within a HD framework, recognizing the intrinsic relationship between LP and VAR models. To illustrate this, consider $N$ time series observed over $T$ time periods. A critical challenge arises in the HD--LP setup: the effective sample size is only of the order $T$, yet the number of unknown parameters can scale as $N$ or even $N^2$ (e.g., as shown in (ref) of Example (ref) and (ref) of Section (ref) respectively). In the latter case, this means even a modest cross-sectional dimension, such as $N=10$, yields at least $100\cdot p$ unknown parameters in which $p$ is the number of lags to be specified. It then calls for a toolkit to understand and implement the LP approach robustly in HD settings, thereby enabling accurate estimation of impulse response functions when the cross-sectional dimension is large -- a common scenario in empirical macroeconomics and finance. To the best of our knowledge, very limited research has investigated long-horizon inference for the LP approach within the HD setting. In this regard, we believe that we make a significant contribution by deriving a set of basic results (e.g., concentration inequality and HD covariance estimation) for HD $h$-step ahead forecasting models in the context of long-horizon analysis.

The second main stream that we review applies the LP method to nonlinear models. For instance, BB2019 consider the LP approach based on a varying coefficient framework. GONCALVES2021107 introduce a semi-parametric model and examine the validity of the LP method, while INOUE2024105726 further introduce a nonparametric time-varying framework in order to model local instabilities. In sharp contrast to this existing literature, we motivate and model the nonlinearity via the residual terms as proposed in the relevant literature, such as in Assumption 2 of WP2009, and Assumption A.3 of DG2018. This approach has direct and critical impacts on the long-horizon inference established in our paper, thereby offering a distinct path to further generalize the aforementioned nonlinear models.

The third stream that is relevant to our work considers simultaneous equations and IV approaches. For a comprehensive survey of the simultaneous equations, we refer interested readers to GONCALVES2024105702 for details, wherein they also refer to this line of research as “state-dependent local projection”. PW2021 investigate the IV based approach, and further bridge the LP method and the VAR literature. Although we do not have any specific IV included, our framework can be considered as a set of simultaneous equations. While one could impose additional identification restrictions, similar to those often applied in the VAR literature, these are generally application-driven; hence, we do not pursue this avenue further in the current paper.

Having reviewed the relevant literature, we summarize our contributions below:

enumerate[leftmargin=24pt, parsep=2pt, topsep=2pt] • We establish long-horizon inference specifically for the LP approach in a class of HD moving average infinity--HDMA($\infty$) models, which is considerably different from such models associated with high--dimensional regressors (covariants) discussed in the relevant literature. Our approach relies only on the underlying DGP, circumventing the need for restrictive assumptions like mixing conditions or near-epoch dependence imposed in the relevant literature. • While BL2008 lay the foundation for HD covariance matrix estimation, our study offers the first application of this approach to long--horizon inference within the literature. Based on this result, we establish some asymptotic normality extending the long-horizon inference analysis (montiel2021local) to our class of HD--MA($\infty$) models. • Motivated by Assumption 2 of WP2009 and Assumption A.3 of DG2018, we model the nonlinearity via the residual terms, and offer a distinct path to further generalize the existing literature of nonlinear models. The main challenge is how to integrate this general structure into the $h$-step ahead forecasting models and work out its theoretical properties for such a class of HDMA($\infty$) models. • In addition, we conduct extensive simulations to examine the theoretical results. In the empirical study, we utilize the proposed HD--LP framework to revisit the impact of business news attention on U.S. industry-level stock volatility. Though the structural VAR models, as a powerful toolkit, have been widely used in the relevant literature diebold2009measuring, diebold2014network to trace how volatility shocks propagate across markets or industries, such methods impose rigid assumptions on system dynamics. Our method, on the other hand, allows for more flexibility in capturing the effect of news shocks and the volatility spillover effects among industries. Finally, we compute impulse-response functions at different horizons, and summarize how news-driven volatility shocks evolve over time.

The remainder of the paper is organized as follows. Section (ref) introduces the model setup and methodology, and establishes the asymptotic results. In Section (ref), we conduct extensive simulations to corroborate the theoretical findings. Section (ref) presents the empirical study on the impact of business news attention on U.S. industry-level stock volatility. Section (ref) concludes. Appendix (ref) provides supplementary results, while Appendix (ref) contains the technical proofs and preliminary lemmas.

Before proceeding further, we introduce some notations which will be repeatedly used in the paper. For a positive integer $L$, we let $[L]\coloneqq \{1,\ldots,N\}$. Let $\mathbf{e}_i $ be a $N\times 1$ selection vector with the $i^{th}$ element being 1 and the others being 0, and $i\in [N]$. For a vector $\mathbf{v} =(v_1,\ldots, v_N)^\top$, we let $\|\mathbf{v} \|_1\coloneqq \sum_{i=1}^N |v_i|$, and $\|\mathbf{v} \|_2\coloneqq \{\sum_{i=1}^N |v_i|^2\}^{1/2}$. For a matrix $\mathbf{B}\coloneqq \{b_{ij}\}$, we let $\|\mathbf{B}\|$, $\|\mathbf{B}\|_2$, $\|\mathbf{B}\|_{\max}$, $\|\mathbf{B}\|_1$, and $\|\mathbf{B}\|_\infty$ define its Frobenius norm, spectral norm, max norm, column norm, and row norm respectively. Additionally, we define the following operator for a symmetric square matrix $\mathbf{B}$:

eqnarray[eqnarray omitted — 134 chars of source]

where $\eta$ is a user-specified threshold parameter. For a random variable $z$, let $\|z\|_p\coloneqq (E|z|_{p}^p )^{1/p}$. For $a\in \mathbb{R}$, let $\lfloor a \rfloor$ be the largest positive integer less than or equal to $a$. For two numbers, we let $a\lesssim b$ stand for $a=O(b)$; let $a\asymp b$ stands for $a\lesssim b$ and $b\lesssim a$; $a\wedge b$ and $a\vee b$ stand for $\min(a,b)$ and $\max(a,b)$ respectively. For two random variables (say, $a$ and $b$ again for simplicity), we let $a\lesssim b$ stand for $a=O_P(b)$, and accordingly define $a\asymp b$.

Setup and Asymptotic Properties

In this section, we first present the DGP and justify its generality (Section (ref)). Next, Section (ref) introduces the setup for $h$-step ahead forecasting models and their useful properties. We then detail the estimation procedure and establish its corresponding asymptotics under minimal conditions in Section (ref). Finally, Section (ref) imposes more structural assumptions, specifically considering a HD covariance matrix estimation, and establishes the asymptotic distribution for practical analysis. Secondary results are relegated to Appendices (ref)-(ref) for the sake of space.

Data Generating Process

Suppose that the following panel dataset is observable.

eqnarray[eqnarray omitted — 85 chars of source]

where $i$ indexes the number of individuals, $t$ indexes the time periods, and $p \, (\ge 1)$ is introduced to accommodate the lag terms in the dynamic structure to be specified and estimated later. In an AR(1) model, $p=1$ ensures the observable dataset is $\{x_{it}\mid i\in [N], \, t \in [T]\}$. The parameter $p$ is preserved in (ref) solely to simplify the notation during theoretical derivations.

Additionally, suppose the data generating process (DGP) of $\{x_{it}\}$ admits the following HDMA($\infty$) process:

eqnarray[eqnarray omitted — 115 chars of source]

where, for $\forall \ell\ge 0$,

eqnarray*[eqnarray* omitted — 174 chars of source]

is an $N\times N$ matrix, $\mathbf{x}_t=(x_{1t},\ldots, x_{Nt})^\top$, and $\pmb{\varepsilon}_{t} =(\varepsilon_{1t},\ldots, \varepsilon_{Nt})^\top$ with $\varepsilon_{it}$ being independent and identically distributed (i.i.d.) over both dimensions.

The popularity and importance of MA($\infty$) processes is governed by the Wold decomposition (fan2003nonlinear). Its recent extensions to fixed multi--dimensional settings include gpy2024a and ygp2025. Equation (ref) is a substantial extension of such MA($\infty$) processes to a HDMA($\infty$) setting.

Throughout, we suppose that $E[\pmb{\varepsilon}_{t}]=\mathbf{0}$ and $E[\pmb{\varepsilon}_{t}\pmb{\varepsilon}_{t}^\top]=\mathbf{I}_N$ without loss of generality. Otherwise, one always encounters the identification issue due to the fact that

eqnarray*[eqnarray* omitted — 137 chars of source]

for any conformable matrix $\mathbf{W}$ with invertibility.

Our goal is to infer the impulse responses defined below via the LP approach.

eqnarray[eqnarray omitted — 235 chars of source]

where $h\ge 1$, $j\in [N]$, and $\pmb{\varepsilon}_{t\mid j}(a) \coloneqq (\varepsilon_{1t},\ldots,\varepsilon_{j-1,t}, a, \varepsilon_{j+1,t},\ldots, \varepsilon_{Nt})^\top$. Using (ref) and (ref), simple algebra shows that

eqnarray[eqnarray omitted — 70 chars of source]

so the key question is how to recover $\mathbf{B}_h$ for $h\ge 1$.

Before proceeding further, we show the flexibility of the above setup.

example\normalfont If the DGP of (ref) follows a VAR(1) process $\mathbf{x}_t = \mathbf{a}_1\mathbf{x}_{t-1}+\pmb{\varepsilon}_t$, it then yields the following $h$-step ahead forecasting model for $h\ge 1$ \begin{eqnarray} \mathbf{x}_{t+h} =\mathbf{a}_1^{h} \mathbf{x}_t+\mathbf{u}_{t,h}, \end{eqnarray} where $\mathbf{u}_{t,h} \coloneqq \mathbf{a}_1^{h-1}\pmb{\varepsilon}_{t+1}+ \cdots+\mathbf{a}_1\pmb{\varepsilon}_{t+h-1}+\pmb{\varepsilon}_{t+h}.$ Notably, the following two facts hold: (i). $\mathbf{u}_{t,h}$ is $(h-1)$-dependent along the time dimension (i.e., $\mathbf{u}_{t,h}$ and $\mathbf{u}_{s,h}$ are independent for $|t-s|\ge h$); (ii). for $\forall t$, $\{\mathbf{u}_{t,h}\mid h\ge 1 \}$ and $\{\mathbf{x}_s \mid s\le t\}$ are independent of each other. One can further extend the above DGP to include (conditional) heteroskedasticity. For example, consider an autoregressive conditional heteroskedasticity (ARCH) structure: \begin{eqnarray*} \mathbf{x}_t = \mathbf{a}_1\mathbf{x}_{t-1}+\tilde{\pmb{\varepsilon}}_t,\quad \tilde{\pmb{\varepsilon}}_t =\pmb{\sigma}_t^{1/2}\pmb{\varepsilon}_t, \quadand\quad \pmb{\sigma}_t = \mathbf{b}_0+ \mathbf{b}_1^\top \pmb{\sigma}_{t-1} \mathbf{b}_1, \end{eqnarray*} where $\mathbf{b}_0$ and $\mathbf{b}_1$ need to fulfil certain conditions. In this case, the $h$-step ahead forecasting model is almost the same as (ref) with minor modification: \begin{eqnarray} \mathbf{x}_{t+h} =\mathbf{a}_1^{h} \mathbf{x}_t+\tilde{\mathbf{u}}_{t,h}, \end{eqnarray} where $\tilde{\mathbf{u}}_{t,h} \coloneqq (\mathbf{a}_1^{h-1}\pmb{\sigma}_{t+1}^{1/2},\ldots, \mathbf{a}_1 \pmb{\sigma}_{t+h-1}^{1/2}, \pmb{\sigma}_{t+h}^{1/2}) \mathbf{v}_{t,h}$ with $\mathbf{v}_{t,h}\coloneqq (\pmb{\varepsilon}_{t+1}^\top,\ldots,\pmb{\varepsilon}_{t+h}^\top)^\top$. In (ref), $\mathbf{v}_{t,h}$ is $(h-1)$-dependent, and for $\forall t$, $\{\mathbf{v}_{t,h}\mid h\ge 1 \}$ and $\{\mathbf{x}_s \mid s\le t\}$ are still independent of each other.
example\normalfont If the DGP of (ref) follows a VAR($p$) process, we have \begin{eqnarray} \begin{pmatrix} \mathbf{x}_{t} \\ \vdots \\ \mathbf{x}_{t-p+1} \end{pmatrix} =\begin{pmatrix} \mathbf{a}_{1:(p-1)}&\mathbf{a}_p\\ \begin{matrix} \mathbf{I} \end{matrix}& \mathbf{0} \end{pmatrix} \begin{pmatrix} \mathbf{x}_{t-1}\\ \vdots \\ \mathbf{x}_{t-p} \end{pmatrix} + \begin{pmatrix} \pmb{\varepsilon}_t \\ \mathbf{0} \end{pmatrix}\eqqcolon \tilde{\mathbf{A}}_p \mathbf{X}_{t-1} + \pmb{\mathcal{E}}_t, \end{eqnarray} where $\mathbf{a}_{1:(p-1)}=(\mathbf{a}_1 , \ldots , \mathbf{a}_{p-1})$, $\mathbf{I}$ and $\mathbf{0}$ are conformable matrix and vector respectively, and the definitions of $\tilde{\mathbf{A}}_p$, $\mathbf{X}_t$ and $\pmb{\mathcal{E}}_t$ are evident. According to Appendix (ref), we can further obtain that \begin{eqnarray*} \mathbf{x}_{t+h}=\mathbf{S}_{p}\tilde{\mathbf{A}}_p^{h}\mathbf{X}_t + \mathbf{u}_{t,h}, \end{eqnarray*} where $\mathbf{S}_{p} =(\mathbf{I}_N, \mathbf{0}_{N\times N(p-1)} )$, $\mathbf{u}_{t,h} \coloneqq \tilde{\mathbf{a}}_p^{h-1}\pmb{\varepsilon}_{t+1}+ \cdots+ \tilde{\mathbf{a}}_p\pmb{\varepsilon}_{t+h-1}+\pmb{\varepsilon}_{t+h}$, and $\tilde{\mathbf{a}}_p^\ell\coloneqq\mathbf{S}_{p} \tilde{\mathbf{A}}_p^{\ell}\mathbf{S}_{p}^\top$. Still, $\mathbf{u}_{t,h}$ is $(h-1)$-dependent, and for $\forall t$, $\{\mathbf{u}_{t,h}\mid h\ge 1 \}$ and $\{\mathbf{x}_s \mid s\le t\}$ are independent of each other. Again, one can extend (ref) to capture (conditional) hetroskedasticity, and we will no longer discuss it due to similarity.
example\normalfont We now provide an example to connect our study with some popular examples of the literature. For notational simplicity we focus on the case with one lag only, i.e., Example (ref) again. The extension with multiple lags such as Example (ref) should be straightforward in view of the following justification. Firstly, we partition $\mathbf{x}_t$ of Example (ref) into two parts ($\mathbf{x}_t^*$ and $\mathbf{x}_t^\dag$) as follows: \begin{eqnarray} \begin{pmatrix} \mathbf{x}_t^* \\ \mathbf{x}_t^\dag \end{pmatrix} = \begin{pmatrix} \mathbf{a}_1^* & \mathbf{a}_1^\dag\\ \mathbf{a}_2^* & \mathbf{a}_2^\dag \end{pmatrix}\begin{pmatrix} \mathbf{x}_{t-1}^* \\ \mathbf{x}_{t-1}^\dag \end{pmatrix}+\begin{pmatrix} \pmb{\varepsilon}_t^* \\ \pmb{\varepsilon}_t^\dag \end{pmatrix}, \end{eqnarray} wherein $\{\pmb{\varepsilon}_t^*\}$ and $\{\pmb{\varepsilon}_t^\dag\}$ are independent of each other. It yields that \begin{eqnarray*} \left\{ \begin{array}{l} \mathbf{x}_t^* = \mathbf{a}_1^* \mathbf{x}_{t-1}^* + \mathbf{a}_1^\dag \mathbf{x}_{t-1}^\dag +\pmb{\varepsilon}_t^*\\ \mathbf{x}_t^\dag =\mathbf{a}_2^* \mathbf{x}_{t-1}^* + \mathbf{a}_2^\dag \mathbf{x}_{t-1}^\dag +\pmb{\varepsilon}_t^\dag \end{array} \right. , \end{eqnarray*} which is then similar to Eq. (2) of GONCALVES2024105702 with constant parameters. One may further regulate $\{\mathbf{x}_t^\dag \}$ to be exogenous regressors (i.e., $\mathbf{a}_2^*\equiv \mathbf{0}$), so (ref) reduces to \begin{eqnarray*} \begin{pmatrix} \mathbf{x}_t^* \\ \mathbf{x}_t^\dag \end{pmatrix} = \begin{pmatrix} \mathbf{a}_1^* & \mathbf{a}_1^\dag\\ \mathbf{0} & \mathbf{a}_2^\dag \end{pmatrix}\begin{pmatrix} \mathbf{x}_{t-1}^* \\ \mathbf{x}_{t-1}^\dag \end{pmatrix}+\begin{pmatrix} \pmb{\varepsilon}_t^* \\ \pmb{\varepsilon}_t^\dag \end{pmatrix}, \end{eqnarray*} In this case, if we let $\mathbf{a}_2^\dag \equiv \mathbf{0}$, $\{\mathbf{x}_t^\dag \ (\equiv \pmb{\varepsilon}_t^\dag) \}$ become purely white noises. The setting therefore is in the same sprit as in INOUE2024105726 and GONCALVES2024105702. Additionally, we may focus exclusively on modeling $\mathbf{x}_t^*$. The first equation of (ref) then represents a standard dynamic model incorporating a HD set of exogenous regressors: \begin{eqnarray} \mathbf{x}_t^* = \mathbf{a}_1^* \mathbf{x}_{t-1}^* +\mathbf{a}_1^\dag\mathbf{x}_{t-1}^{\dag} + \pmb{\varepsilon}_t^*. \end{eqnarray} In the limiting case where $\mathbf{x}_t^*$ is a scalar, $\mathbf{a}_1^\dag$ becomes a HD row vector. In this configuration, (ref) aligns with the frameworks established by ASW2024 and cha2024.

The Setup

Up to this point, we have discussed our goal and the flexibility of the DGP process. However, (ref) has infinite parameters, so we have to further impose certain structure in order to get some meaningful results practically while maintaining the assumptions as flexible as possible. In Appendix (ref), we provide a Proposition (ref) to briefly answer the question that how far we can go in approximating (ref) without knowing the underlying DGP, which however is not the main focus of this paper.

Having said that, we suppose that $\{\mathbf{x}_t\}$ also admit the following set of regression models:

eqnarray[eqnarray omitted — 260 chars of source]

where $p$ is a fixed positive constant, $h \ (\ge 1)$ may diverge, $\tilde{\mathbf{A}}=\{\tilde{A}_{ij} \}_{N\times Np}\coloneqq (\mathbf{A}_1,\ldots, \mathbf{A}_p)$, and $\mathbf{u}_{t,h}=(u_{1t,h},\ldots, u_{Nt,h})^\top$ has an $(h-1)$-dependent structure to be specified below. It should be understood that $(\mathbf{A}_1,\ldots, \mathbf{A}_p)$ vary with respect to $h$, which is suppressed in the sub-indices for notational simplicity when no misunderstanding arises. The set of regression models given in (ref) is similar to Eq. (2) of Oscar_2005 with a focus on long-horizon inference under the HD setting.

To facilitate development, we impose the following assumptions.

assumption$\{\varepsilon_{it}\}$ are i.i.d. over both $i$ and $t$, and satisfy that $E[\varepsilon_{11}]=0$, $E[\varepsilon_{11}^2]=1$, and $E|\varepsilon_{11}|^J<\infty$ with a fixed $J \, (\ge 4)$.

Assumption (ref) only regulates the random components of (ref). We will impose more restrictions on the deterministic matrices involved whenever necessary.

assumption\begin{enumerate}[leftmargin=24pt, parsep=2pt, topsep=2pt] • Assume that $g(\cdot)$ is a smooth function such that $\mathbf{u}_{t,h}=\mathbf{g}(\pmb{\varepsilon}_{t+h},\ldots, \pmb{\varepsilon}_{t+1})$ satisfies that $E[\mathbf{u}_{t,h}]=\mathbf{0}$, $E[\mathbf{u}_{t,h}\mathbf{u}_{t,h}^\top]=\pmb{\Sigma}_h$, and $\max_{i,h}\|\mathbf{e}_i^\top\mathbf{u}_{t,h} \|_J<\infty$, where $J$ is given in Assumption (ref). • Suppose that (i). $\max_{j}\sum_{\ell = 0}^{\infty} \ell^{\frac{J}{2(J+1)}} \|\pmb{\beta}_{\ell, j} \|_2^{\frac{J}{J+1}}<\infty$ with $\{\pmb{\beta}_{\ell, j}\}$ being defined below ((ref)); (ii). $\frac{d_{\pmb{\beta}}^4\log N}{T}\to 0$ with $d_{\pmb{\beta}}\coloneqq \sum_{\ell=0}^\infty \|\mathbf{B}_\ell\|_{\infty} $; (iii). $h/T\to 0$ as $(h, T)\to (\infty,\infty)$. \end{enumerate}

Assumption (ref) is readily justified in light of Examples (ref) and (ref). Specifically, Assumption (ref).2.(i) requires a reasonably slow decay rate for the norm $\|\pmb{\beta}_{\ell, j} \|_2$. At this stage, the sample size $N$ can be large in view of the condition in Assumption (ref).2.(ii). Assumption (ref).2.(ii) allows that $N$ and $T$ can be proportional to each other. Assumption (ref).2.(iii) permits that $h$ can go to $\infty$ as long as $\frac{h}{T}\rightarrow 0$.

With the above setup, we present the first result of this paper.

propositionSuppose that $\mathbf{x}_t$ admits both representations (ref) and (ref), Assumptions (ref) and (ref).1 hold, and $\mathbf{B}_0 =\mathbf{I}_N$. For $\forall h\ge 1$, we obtain that $\mathbf{A}_{1} = \mathbf{B}_{h} = ({\normalfont \textbf{IR}_{h,1}},\ldots, {\normalfont \textbf{IR}_{h,N}})$.

The condition $\mathbf{B}_0 =\mathbf{I}$ serves a standard identification purpose, aligning with Oscar_2005. Without this normalization, we can only determine the relationship $\mathbf{A}_1 \mathbf{B}_0 = \mathbf{B}_h$. Given the definition of the impulse responses in (ref), this lack of identification means the LP method would only recover the impulse responses up to a rotation matrix. The matrices $\mathbf{A}_j$ for $j\ge 2$, while being essential components of the underlying model (ref), are not directly connected to the impulse responses.

To conclude our model setup, we establish the following concentration inequality for the $h$-step ahead forecasting models based on (ref) and (ref), which will facilitate the estimation of the parameters of interest in the next subsection.

theoremUnder Assumptions (ref)-(ref), there exit positive constants $c_1,c_2,c_3, c_4$ such that for $\forall i,j\in [N]$, $\forall\ell \in 0\cup [p-1]$, and $\delta >c_4\sqrt{T}\mu_{h}^{1+1/J}$ \begin{eqnarray*} \Pr\left(\max_{t\in [T-h]}\left|\sum_{s=1}^{t}u_{is,h} x_{j,s-\ell}\right|\ge \delta\right) \le c_1 \frac{T}{\delta^J} \mu_h^{J+1}+c_2\exp\left(-\frac{c_3\delta^2}{T \mu_h^{2+2/J}}\right), \end{eqnarray*} where $\displaystyle\mu_h= \max_{j\in [N]}\sum_{\ell=0}^\infty [(\ell+h)^{\frac{J}{2}-1} \|\pmb{\beta}_{\ell, j} \|_2^J]^{\frac{1}{J+1}}\lesssim h^{\frac{J-2}{2(J+1)}}$.

Theorem (ref) establishes a concentration inequality for the $h$-step ahead forecasting models of (ref), contributing to the literature on long-horizon inference in a HD setting. Notably, the bound contains the term $\mu_h$, which diverges as the horizon $h$ increases. The term $\mu_h$ is simply due to lack of structure in $\mathbf{u}_{t,h}$ and (ref). Since $J\ge 4$ as stated in Assumption (ref), we obtain that $h^{\frac{J-2}{2(J+1)}}\in [h^{1/5}, h^{1/2}]$. Therefore, in the worst case scenario, $\mu_h$ diverges at the rate $h^{1/2}$. From a modeling perspective, this result explains the inherent loss of accuracy in forecasting models as $h$ becomes large.

Furthermore, the derivation of this theorem relies only on the underlying DGP, circumventing the need for restrictive assumptions like mixing conditions or near-epoch dependence. In this regard, it substantially improves upon the findings in Section 4 of cha2024.

Estimation

We are now ready to consider the estimation in this subsection. Firstly, we assume $p$ is known, and discuss how to select the optimal lag later.

According to (ref), we define the following objective function:

eqnarray[eqnarray omitted — 219 chars of source]

where $\tilde{\mathbf{a}}$ is a generic $N\times Np$ matrix. However, due to the large amount of elements included in $\tilde{\mathbf{A}}=\{\tilde{A}_{ij} \}_{N\times Np}$, one cannot fully recover everything from (ref) unless the effective sample size $T-h$ is extremely large practically. Consequently, it might not always be feasible for real data analysis. In order to find a balance between feasibility and generality, we impose sparsity on $\tilde{\mathbf{A}}$ of (ref), and define some new notation to facilitate the investigation. Let

eqnarray*[eqnarray* omitted — 202 chars of source]

Apparently, $d_{\mathscr{A}} \coloneqq \|\mathbf{S}_{\mathscr{A}}\|^2$ gives the total number of nonzero elements in $\tilde{\mathbf{A}}$. It is noteworthy that we allow for the case that $\tilde{\mathbf{A}}=\mathbf{0}$, i.e., $d_{\mathscr{A}}=0$. For this extreme case, $\mathbf{x}_t$ reduces to a HD white noise.

Using (ref), we then introduce the following estimation procedure to recover $\tilde{\mathbf{A}}$ and its sparsity.

enumerate[leftmargin=48pt, parsep=2pt, topsep=2pt] • Initial estimation via LASSO: \begin{eqnarray} \hat{\tilde{\mathbf{a}}} =\operatorname*{\arg\!\min}_{\tilde{\mathbf{a}}} (Q_0(\tilde{\mathbf{a}})+\gamma\|\operatorname*{\normalfontvec}(\tilde{\mathbf{a}})\|_1), \end{eqnarray} where $\gamma$ is a tuning parameter. • Refined estimation via adaptive LASSO: \begin{eqnarray} \hat{\tilde{\mathbf{a}}}_{\phi} =\operatorname*{\arg\!\min}_{\tilde{\mathbf{a}}}( Q_0(\tilde{\mathbf{a}})+ \gamma\|\operatorname*{\normalfontvec}(\tilde{\mathbf{a}})\circ \pmb{\phi}\|_1), \end{eqnarray} where $\pmb{\phi}=(\phi_1,\cdots, \phi_{N^2p})^\top$ is a vector of predetermined weights.

The above estimation procedure connects the adaptive LASSO developed in Zou2006 with the LP approach in a HD framework. For the adaptive LASSO, we refer interested readers to Zou2006 for extensive theoretical investigation and to mcilhagga2016penalized for detailed numerical implementation. In the second step, we may let $\phi_\ell=|\hat{\tilde{a}}_{\ell}|^{-\zeta}$ with $\hat{\tilde{a}}_{\ell}$ being the $\ell^{th}$ element of $\hat{\tilde{\mathbf{a}}}$, in which $\zeta>0$ is an arbitrary positive constant, as long as certain conditions to be specified below are fulfilled.

To proceed, we introduce some additional notation and assumptions. Let $\phi_{\mathscr{A},\ell}$ and $\phi_{\bar{\mathscr{A}},\ell}$ denote the $\ell^{th}$ elements of $\pmb{\phi}_{\mathscr{A}}$ and $\pmb{\phi}_{\bar{\mathscr{A}}}$ respectively, where $\pmb{\phi}_{\mathscr{A}}$ and $\pmb{\phi}_{\bar{\mathscr{A}}}$ contain the elements of $\pmb{\phi}$ that correspond to the non--zero elements of $\operatorname*{\normalfont\textrm{vec}}(\mathbf{S}_{\mathscr{A}})$ and $\operatorname*{\normalfont\textrm{vec}}(\mathbf{S}_{\bar{\mathscr{A}}})$. Let $\pmb{\Sigma}_{\mathbf{B}} \coloneqq \{\pmb{\Sigma}_{\mathbf{B}, ji}\}$ with $0\le i,j\le p-1$ and

eqnarray*[eqnarray* omitted — 342 chars of source]
assumption\begin{enumerate}[leftmargin=24pt, parsep=2pt, topsep=2pt] • Suppose that $\operatorname*{\normalfont\textrm{vec}}(\mathbf{a})^\top (\pmb{\Sigma}_{\mathbf{B}} \otimes \mathbf{I}_N)\operatorname*{\normalfont\textrm{vec}}(\mathbf{a})\ge \alpha \|\operatorname*{\normalfont\textrm{vec}}(\mathbf{a})\|_2^2$ for $\forall\mathbf{a}\in \mathbb{A}(\mathbf{S}_{\mathscr{A}})$ and $\mathbf{a}\ne \mathbf{0}$ uniformly in $N$, where $\alpha>0$ is a fixed positive constant, and \begin{eqnarray*} \mathbb{A}(\mathbf{S}_{\mathscr{A}} ) &\coloneqq &\{\mathbf{a}=\{a_{ij}\}_{N\times Np}\mid \|\operatorname*{\normalfontvec}(\mathbf{S}_{\bar{\mathscr{A}}}\circ \mathbf{a} ) \|_1\le 3 \|\operatorname*{\normalfontvec}(\mathbf{S}_{\mathscr{A}}\circ \mathbf{a})\|_1\}. \end{eqnarray*} • Suppose that $d_{\mathscr{A}} d_{\pmb{\beta}}^2\frac{\sqrt{\log N}}{\sqrt{T}}\to 0$ and $\frac{N^2 }{T^{J-1}\log N}\to 0$, where $J$ is given in Assumption (ref). \end{enumerate}

The first condition of Assumption (ref).1 is the so-called restricted eigenvalue condition, which has been fully discussed in the literature. See BRT2009 and raskutti10a for example. In the definition of $\mathbb{A}(\mathbf{S}_{\mathscr{A}} )$, the targeted set $\|\operatorname*{\normalfont\textrm{vec}}(\mathbf{S}_{\bar{\mathscr{A}}}\circ \mathbf{a} ) \|_1\le 3 \|\operatorname*{\normalfont\textrm{vec}}(\mathbf{S}_{\mathscr{A}}\circ \mathbf{a})\|_1$ is obvious in view of the development given in (ref) of the appendix. The second condition further regulates the sparsity and the sample size involved.

With these in hand, we are able to achieve the following results for the 2-step procedure.

theoremLet Assumptions (ref)-(ref) hold and $\gamma \asymp \frac{\mu_{h}^{1+1/J}\sqrt{\log (N)}}{\sqrt{T}}$, where $\mu_h$ is the same as in Theorem 1. \begin{enumerate}[leftmargin=24pt, parsep=2pt, topsep=2pt] • For (ref), the following results hold: \begin{enumerate}[leftmargin=24pt, parsep=2pt, topsep=2pt] • $\|\operatorname*{\normalfont\textrm{vec}}(\tilde{\mathbf{A}}-\hat{\tilde{\mathbf{a}}}) \|_2=O_P \big( \frac{\mu_{h}^{1+1/J}\sqrt{d_{\mathscr{A}}\log N}}{\sqrt{T}} \big);$$\|\operatorname*{\normalfont\textrm{vec}}(\tilde{\mathbf{A}}-\hat{\tilde{\mathbf{a}}}) \|_1 =O_P \big( \frac{\mu_{h}^{1+1/J}d_{\mathscr{A}}\sqrt{\log N}}{\sqrt{T}} \big)$. \end{enumerate} • For (ref), let the elements of $\pmb{\phi}$ satisfy \begin{enumerate}[leftmargin=24pt, parsep=2pt, topsep=2pt] • $\min_{\ell \in [d_\mathscr{A}]}|\beta_{0,\mathscr{A},\ell}| \gtrsim \frac{\mu_{h}^{1+1/J}\sqrt{d_{\mathscr{A}}\log N}}{\sqrt{T}} \cdot \max_{\ell\in[d_{\mathscr{A}}]}\phi_{\mathscr{A},\ell}$; • $ \min_{\ell\in[d_{\bar{\mathscr{A}}}]}\phi_{\bar{\mathscr{A}},\ell} \gtrsim d_{\mathscr{A}} \max_{\ell\in[d_{\mathscr{A}}]}\phi_{\mathscr{A},\ell}$. \end{enumerate} Then $\operatorname*{\normalfont\textrm{sgn}}(\hat{\tilde{\mathbf{a}}}_{\phi}) =\operatorname*{\normalfont\textrm{sgn}} (\tilde{\mathbf{A}})$ with probability approaching one. \end{enumerate}

The first result of Theorem (ref) explains how the $h$-step ahead forecasting framework affects the estimation results when $h$ diverges. The second result proves the sign consistency, i.e., the identification of 0's of $\tilde{\mathbf{A}}$, wherein the additional restrictions can be easily fulfilled in view of the discussion under (ref).

Selection of the lags \ With Theorem (ref) in hand, we can then choose the lag terms involved. The following selection criterion is based on the first step of the estimation procedure (i.e., (ref)).

Define the following information criterion:

eqnarray*[eqnarray* omitted — 246 chars of source]

where $\hat{\tilde{\mathbf{a}}}_{\mathfrak{p}}$ is obtained via (ref) using $\mathfrak{p}$ lags, and $\xi$ is a tuning parameter satisfying certain condition to be specified. The optimal lag is then estimated by

eqnarray[eqnarray omitted — 129 chars of source]

where $\mathfrak{p}^*$ is a user-specified large and fixed integer.

theoremLet Assumptions (ref)-(ref) hold, $\frac{\mu_h^{2+2/J} d_{\mathscr{A}} \log N}{T\xi}\to 0$, and $\xi\to 0$. Then $\Pr(\widehat{p} = p)\to 1$.

In order to select the optimal number of lags practically, we make two comments:

enumerate[leftmargin=24pt, parsep=2pt, topsep=2pt] • An intuitive choice of $\xi$ is to let $\xi\asymp \gamma$, so the condition $\frac{\mu_h^{2+2/J} d_{\mathscr{A}} \log N}{T\xi}\to 0$ reduces to $\frac{\mu_h^{1+1/J} d_{\mathscr{A}} \sqrt{\log N}}{\sqrt{T} }\to 0$, which is required by Theorem (ref).1 already. Therefore, such a choice does not create any additional restriction. • In fact, one can further simplify the selection process. It is worth pointing out that Oscar_2005 defines the LP approach as a series of regressions such as those in (ref), where the number of lags is independent of the horizon $h$. Theorem (ref) establishes a result that holds for $\forall h$. Therefore, from a practical standpoint, we suggest that using a small horizon (e.g., $h=1$ or $2$) is preferable for selecting the optimal lag, as this significantly simplifies the required restrictions. For instance, the term $\mu_h$ no longer plays a role in such cases. In the simulation, we shall further examine this point.

Asymptotic Distribution

It is worth mentioning that we have not imposed too many restrictions on $\mathbf{u}_{t,h}$ so far. However, to achieve the HD covariance estimation and establish the asymptotic normality in what follows, we need to add more structures to the residuals (i.e., $\mathbf{u}_{t,h}$'s).

Firstly, we consider the estimation of the covariance matrix. To be precise, our goal is to estimate $\pmb{\Omega}_h \coloneqq E[\tilde{\pmb{\Omega}}_h]$, where

eqnarray[eqnarray omitted — 218 chars of source]

Note that if $\tilde{\mathbf{A}}$ is sparse, it will naturally pass the sparsity to $\pmb{\Sigma}_h$, and consequently pass the sparsity to $\pmb{\Omega}_h$. To see this point, consider an extreme case using Example (ref).

Example 1 (Cont.) If $\mathbf{a}_1=\mathbf{0}$, then $\pmb{\Sigma}_h=\mathbf{I}_N$ and $\pmb{\Omega}_h =\mathbf{I}_{N^2}$ for $\forall h\ge 1$.

The justification of the above statement should be obvious, so we omit the details. Therefore, without loss of generality, we suppose that $\pmb{\Omega}_h\in \mathbb{U}(c_a, c_N, \bar{c})$, where

eqnarray[eqnarray omitted — 203 chars of source]

$0<c_a<1$, and $\bar{c}$ is a fixed number. The restriction (ref) is the so-called densest sparse condition. We refer interested readers to BL2008 and extensive extensions since then for discussion on (ref).

We now move to estimate $\pmb{\Omega}_h$. Recall the the operator $\mathscr{G}_\eta(\cdot)$ defined in Section (ref), and let

eqnarray[eqnarray omitted — 233 chars of source]

where $\hat{\mathbf{u}}_{t,h}\coloneqq \mathbf{x}_{t+h}-(\mathbf{X}_t^\top \otimes \mathbf{I}_N) \operatorname*{\normalfont\textrm{vec}}(\hat{\tilde{\mathbf{a}}}_{\phi})$. Putting these together, the estimator of $\pmb{\Omega}_h$ is given by

eqnarray[eqnarray omitted — 83 chars of source]

where $\eta \asymp \sqrt{h\log N /T}$.

In order to investigate (ref), we impose some extra conditions in Assumption (ref).

assumption\begin{enumerate}[leftmargin=24pt, parsep=2pt, topsep=2pt] • There exists a constant $c_\delta>2$ such that for $l=1,\ldots, h-1$ with $h>1$ \begin{eqnarray*} \max_i\|\mathbf{e}_i^\top (\mathbf{u}_{t,h} - \mathbf{u}_{t,h,l}^*)\|_{J} \eqqcolon \max_{i}\delta_{J,i}(h,l) =O(l^{-c_\delta}), \end{eqnarray*} where $\mathbf{u}_{t,h,l}^*$ is the coupled version of $\mathbf{u}_{t,h}$ by replacing $\pmb{\varepsilon}_{t+h-l}$ with an independent copy $\pmb{\varepsilon}_{t+h-l}^*$. • There exists a constant $c_b>2$, such that $\max_j\sum_{\ell=m+1}^{\infty} \|\pmb{\beta}_{\ell,j}\|_2^{2}=O (m^{1-2c_b})$ as $m\to \infty$. • Suppose that for $J>4$, \begin{enumerate} • $\mu_{h}^{2+2/J}d_{\mathscr{A}}^2\sqrt{h\log N}/\sqrt{T}\rightarrow 0$, • $N^4T^{1-J/4}(\log N)^{-J/4}(\log T) ( T^{(-J/4)/c_0} + h^{- J/4-1} )\rightarrow 0,$ \end{enumerate} where $c_0 = \big(\frac{2c_b-1}{2}\big)\big(\frac{c_\delta\wedge c_b-1}{c_\delta\wedge c_b}\big)$. \end{enumerate}

Note that $h$ is involved in a few places of Assumption (ref). The first condition essentially regulates the decay rate of the time series dependence of $\{\mathbf{u}_{t,h}\}$, while the second condition imposes a further restriction on the matrices involved in (ref). To ensure the third condition holds, $J$ has to be larger than 4.

theoremUnder Assumptions (ref)-(ref), \begin{enumerate}[leftmargin=24pt, parsep=2pt, topsep=2pt] • $\max_{i,j}|\hat{\Omega}_{ij}-\Omega_{ij}|=O_P (\sqrt{h\log N /T} )$, • $ \|\mathscr{G}_\eta(\hat{\pmb{\Omega}}_h)- \pmb{\Omega}_h \|_2=O_P((h\log N /T)^{(1-c_a)/2} c_{N})$, \end{enumerate} where $\hat{\Omega}_{ij}$ and $\Omega_{ij}$ are the $(i,j)^{th}$ elements of $\hat{\pmb{\Omega}}_h$ and $\pmb{\Omega}_h$, respectively.

Based on Theorem (ref), we can then derive a central limit theory for the purpose of inference, which further extends the long--horizon inference of montiel2021local to a HD setting.

Specifically, we adopt the node-wise LASSO method to construct the debiased estimator. For notational simplicity, define $\mathbf{Z}_t\coloneqq \mathbf{X}_t\otimes \mathbf{I}_N$. Accordingly, for $i\in[N^2p]$, let $\mathbf{Z}_{i,t}$ and $\mathbf{Z}_{-i,t}$ be the $i^{th}$ row of $\mathbf{Z}_{t}$ and the submatrix of $\mathbf{Z}_{t}$ constructed by removing its $i^{th}$ row, respectively. Additionally, define $\pmb{\rho}_{\mathbf{a}}=(\rho_1,\cdots, \rho_{N^2p})^\top$ as the selection vector such that $\|\pmb{\rho}_{\mathbf{a}}\|_1<\infty$. For each $i$, let $s_{\mathscr{A},i} \coloneqq \|\{I(|\Sigma_{s,ij}|>0)\}_{Np\times 1}\|^2$ denote the sparsity of the $i^{th}$ row of $\pmb{\Sigma}_{\mathbf{B}}^{-1} $, where $\Sigma_{s,ij}$ denotes its $(i,j)^{th}$ element of $\pmb{\Sigma}_{\mathbf{B}}^{-1} $.

Then, we conduct the node-wise LASSO estimation:

eqnarray[eqnarray omitted — 239 chars of source]

where $\tilde{\gamma}_i\asymp \frac{\mu_{h}^{1+1/J}\sqrt{\log (N)}}{\sqrt{T}}$. Let $\hat{\tau}_i^2=\frac{1}{T-h}\sum_{t=1}^{T-h}\|\mathbf{Z}_{i,t}^\top-\mathbf{Z}_{-i,t}^\top\hat{\mathbf{b}}_i \|_2^2+\tilde{\gamma}_i\|\hat{\mathbf{b}}_i\|_1$, and we then define the debiased estimator \(\hat{\tilde{\mathbf{a}}}_{\mathrm{bc}}\) as

eqnarray[eqnarray omitted — 372 chars of source]

where \(\hat{\pmb{\Omega}}_z=\hat{\mathbf{T}}_z^{-1}\hat{\mathbf{C}}_z\), \(\hat{\mathbf{T}}_z=\mathrm{diag}(\hat{\tau}_1^2,\ldots,\hat{\tau}_{N^2p}^2)\), and \(\hat{\mathbf{C}}_z=(\hat{\mathbf{C}}_{z,1},\ldots,\hat{\mathbf{C}}_{z,N^2p})^\top\) is an \(N^2p\times N^2p\) matrix. Each row \(\hat{\mathbf{C}}_{z,i}^\top\) is constructed by placing a 1 in the \(i^{th}\) position and filling the remaining \(N^2p-1\) entries with the elements of \(-\hat{\mathbf{b}}_i\) in their natural order.

Finally, it is noteworthy that under the HD--LP framework, the node-wise LASSO estimation can be much simplified numerically in practical analysis. For the sake of space, we provide the details in Appendix (ref).

We are now ready to establish the following CLT for the debiased estimator.

theoremLet Assumptions (ref)-(ref) hold and $ (d_{\mathscr{A}}+\max_i s_{\mathscr{A},i})(\log N) T^{-1/2}\mu_{h}^{2+2/J}\rightarrow0 $. Then \begin{eqnarray*} \sqrt{T}\,(\pmb{\rho}_{\mathbf{a}}^\top\hat{\pmb{\Omega}}_z\mathscr{G}_\eta(\hat{\pmb{\Omega}}_h)\hat{\pmb{\Omega}}_z\pmb{\rho}_{\mathbf{a}})^{-1/2}\pmb{\rho}_{\mathbf{a}}^\top\operatorname*{\normalfontvec}(\hat{\tilde{\mathbf{a}}}_{\rm bc}-\tilde{\mathbf{A}})\xrightarrow{ D }\mathcal{N}(0,1). \end{eqnarray*}

We have emphasized that the preceding framework is not limited to the large-sample $N\to \infty$ regime; it also encompasses such settings associated with fixed $N$ as special cases. A corollary detailing the corresponding asymptotic distribution for the fixed---$N$ case is provided in Appendix (ref) below.

Up to this point, we have completed our theoretical investigation. It is worth mentioning again that some additional secondary results are also presented in Appendix (ref) for the sake of space. In what follows, we examine these theoretical findings via numerical analyses.

Simulation

Without loss of generality, let the true DGP be as follows:

eqnarray[eqnarray omitted — 128 chars of source]

where $\pmb{\varepsilon}_t\sim N(\mathbf{0},\mathbf{I}_N)$. By (ref), it is obvious that $p=2$. The coefficient matrices $\mathbf{a}_1=\{a_{1,ij}\}_{N\times N}$ and $\mathbf{a}_2=\{a_{2,ij}\}_{N\times N}$ are designed as follows:

eqnarray*[eqnarray* omitted — 390 chars of source]

Also, we let $t\in \{-T_B, -T_B+1,\ldots, T \}$, in which $T_B$ is a sufficiently large positive number corresponding to a burn-in period. According to Example (ref), the model (ref) admits an HDMA($\infty$) representation $\mathbf{x}_t = \sum_{\ell =0}^{\infty}\mathbf{B}_\ell \pmb{\varepsilon}_{t-\ell},$ in which

eqnarray*[eqnarray* omitted — 253 chars of source]

For each generated dataset, we fit the observations $\{\mathbf{x}_t\mid, t=2-p, 3-p,\ldots, T\}$ to the following set of $h$-step ahead forecasting models:

eqnarray[eqnarray omitted — 128 chars of source]

where $h\ge 1$. We suppress the subindex $h$ in $\mathbf{A}_1$ and $\mathbf{A}_2$ when no misunderstanding arises. According to Proposition (ref), we have $\mathbf{A}_1=\mathbf{B}_h$. Additionally, for $h=1$, (ref) and (ref) infer that $\mathbf{a}_1=\mathbf{A}_1=\mathbf{B}_1$.

In what follows, we consider the cases with $N\in \{20, 30, 40\}$ and $T\in \{300, 400, 500\}$, so the number of parameters involved ranges from $800$ to $3200$. Based on the above settings, we firstly present Figure (ref) to demonstrate the sparsity involved in $\mathbf{B}_h$. It is clear that the sparsity exists in all $\mathbf{B}_h$'s, and also the pattern varies with respect to $h$. For example, when $h=10$, $\mathbf{B}_h$ has more non-zero elements, but many of them are close to zero and may become negligible.

figure[figure omitted — 150 chars of source]

Based on the specified DGP, we perform $R$ replications to evaluate our proposed estimation procedure. Without loss of generality, we let $R=500$ throughout. The complete numerical algorithm is detailed in Appendix (ref). We assess the estimation accuracy of $\mathbf{B}_h =\{ b_{h,ij}\}$ after identifying (1). $\hat{p}$, and (2). the sparsity.

We start from examining the selection of the optimal lag, and emphasize that two facts:

enumerate[leftmargin=24pt, parsep=2pt, topsep=2pt] • Over-selection is generally permissible, as it still yields reasonably reliable (consistent) estimates, though it may result in a decrease in estimation efficiency; • As explained under Theorem (ref), we only need to select the optimal $p$ once for the group of regressions such as those in (ref) practically. To minimize the restrictions required, we therefore consider $h=1,2$ only.

To quantify performance, we introduce the following measures:

equation*[equation* omitted — 209 chars of source]

where $\ell\in \{1,2\}$ corresponds to the lag orders $h\in \{1,2\}$ when implementing (ref). These indices represent the empirical frequencies of under-selection ($S_{\ell}^-$), correct selection ($S_{\ell}$), and over-selection ($S_{\ell}^+$), where $\hat{p}_{\ell, r}$ denotes the estimated lag in the $r^{th}$ replication for $\forall \ell\in \{1,2\}$.

The selection results are reported in Panel A of Table (ref). For $\forall \ell\in \{1,2 \}$, we observe that the values of $S_\ell$ are close to 1, which aligns with theoretical expectations. When $\ell =2$, we have a small proportion over-section (i.e., $S_{\ell}^+ >0$), and $S_{\ell}$ moves towards 1 as $T$ increases. As explained above, over-section is not ideal but permissible. It is not surprising that $S_2$ decreases slightly as $N$ goes up. As we have seen in Figure (ref), the coefficients have the strongest signal for $h=1$ when implementing LASSO procedure. Hence, it is reasonable to expect $h=1$ offers better finite sample performance from different perspectives such as those presented by Panel A of Table (ref).

For each $\hat{p}_\ell$ with $\ell\in \{ 1,2\}$, we carry on the second step of the numerical algorithm of Appendix (ref). To measure the performance, we quantify the accuracy of sparsity detection as:

eqnarray*[eqnarray* omitted — 142 chars of source]

where SL stands for the selection, the subscript $\mid \hat{p}_\ell$ indicates the selection results are calculated based on the estimated lag order $\hat{p}_\ell$, $\hat{\mathbf{A}}_{a} =\{ \hat{a}_{a,ij}\}$ defines the adaptive LASSO estimate, and $r$ indexes the $r^{th}$ replication. The SL measure focuses on the overall accuracy of zero/non-zero element classification rather than distinguishing between false negative (incorrectly identified zero elements) and false positive (incorrectly identified non-zero elements).

We observe in Panel A of Table (ref) that the SL value tends to become relatively large as the horizon $h$ increases. This phenomenon is primarily attributable to a source of false negative, as illustrated in Figure (ref). Specifically, as $h$ increases, the number of non-zero elements in the true matrix $\mathbf{B}_h$ may rise; however, many of these elements become negligible (i.e., close to zero) due to weak signal propagation. The adaptive LASSO procedure may incorrectly detect these negligible elements as zero, thereby leading to a higher overall SL score.

Finally, based on $\hat{p}_\ell$ and the estimated sparsity, we introduce the following measure:

eqnarray*[eqnarray* omitted — 154 chars of source]

where AD stands for the averaged distance, $\star \in \{a,d\}$, the subscript $\mid \hat{p}_\ell$ again indicates the results are calculated based on the estimated lag order $\hat{p}_\ell$, and $\hat{\mathbf{S}}_r\coloneqq\{ I( \hat{a}_{a,ij, r}\ne 0)\}$. The definition of $\text{AD}_{\mathbf{B}_h}$ infers that for the debiased LASSO estimator, we only consider those elements which are not identified as 0 by the adaptive LASSO.

Panel B of Table (ref) summarizes the relevant results. The results conditional on $\hat{p}_1$ and $\hat{p}_2$ are very similar, which should be expected. Despite the relatively large number of parameters, which ranges from 800 to 3200, both the adaptive LASSO estimation error ($\text{AD}_{a}$) and the debiased LASSO estimation error ($\text{AD}_{d}$) are observed to be small. Also, both $\text{AD}_{a}$ and $\text{AD}_{d}$ move towards 0 as $T$ increases for all $(N,h)$. Overall, the debiased LASSO exhibits superior finite sample performance, which is consistent with the theoretical expectation that debiasing mitigates the shrinkage bias inherent in the standard LASSO procedure.

table[table omitted — 3,761 chars of source]

Empirical Study

In the era of big data, textual sources such as news and policy announcements are increasingly used in economic research. For example, baker2016measuring develop an economic policy uncertainty index from newspaper articles and find that higher news-based uncertainty is associated with greater stock price volatility. Similarly, bybee2024business construct 180 news attention indices from thousands of Wall Street Journal (WSJ) articles and show that attention topics like recession news have economically large predictive power for future output variables. Their news-attention measures closely track a wide range of macroeconomic and financial series. In particular, it is shown that shifts in WSJ coverage of specific themes coincide closely with observed industry-level market volatility.

In this paper, we utilize the proposed high-dimensional local projection framework to revisit the impact of business news attention on U.S. industry-level stock volatility. Additionally, we explicitly include lagged volatilities from all industries as regressors to model cross-sector spillovers. The spillover of volatility across assets and sectors has long been a focus in finance diebold2009measuring, diebold2014network. Though the structural VAR models, as a powerful toolkit, have been widely used in the relevant literature to trace how volatility shocks propagate across markets or industries, such methods impose rigid assumptions on system dynamics. Our method, on the other hand, allows for more flexibility in capturing the effect of news shocks and the volatility spillover effects among industries. Finally, we compute impulse-response functions at different horizons to trace the short-, medium-, and long-run effects of a news-attention shock on industry volatility. By deriving these responses, we summarize how news-driven volatility shocks evolve over time.

Data and Variables

To construct U.S. industry-level volatility measures, we collect daily industry portfolio returns for 37 industry categories from Kenneth French's data library.\footnote{The data are available at \url{https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html}. We exclude the sector labeled “Other” in the original dataset, as it lacks a clear economic interpretation.} Monthly industry volatility is computed as the standard deviation of daily returns within each calendar month, which is a standard approach in the empirical asset pricing and volatility literature. Details on the industry classifications are provided in Table (ref) and the descriptive statistics for the resulting industry-level volatility series are reported in Table (ref). To assess the time-series properties of the data, we conduct augmented Dickey--Fuller (ADF) unit root tests. As shown in Table (ref), all industry volatility series are found to be stationary at conventional significance levels.

For the news attention data, we rely on the monthly topic-level attention indices constructed by bybee2024business, which cover 180 distinct news topics over the period 1984-2017.\footnote{The data are available at \url{www.structureofnews.com}.} From this universe, we select 23 topics that are most closely related to economic growth and financial markets. Prior evidence in bybee2024business highlights the {\it Recession} topic (REC) as particularly informative for macro-financial dynamics. Accordingly, in our empirical analysis, we progressively incorporate the recession-related news attention series, followed by the remaining seven topics in the economic growth category and fifteen topics in the financial market category, into the local projection framework to examine their effects on industry-level volatilities. A detailed description of the selected topics is provided in Table (ref).

By merging the industry volatility dataset with the news attention indices, we obtain a multivariate time-series dataset spanning January 1984 to June 2017 ($T=402$), comprising a total of 60 variables: 37 industry-level volatility series and 23 news topic attention measures. All volatility series and news attention indices are subsequently standardized to have zero mean and a variance of one. We consider three local projection specifications that differ in the set of regressors included. These models are summarized as follows:

itemize• Model 1: 37 industry-level volatility variables and the recession-related news attention index; • Model 2: 37 industry-level volatility variables and 8 news attention indices related to economic growth; • Model 3: 37 industry-level volatility variables, 8 news attention indices related to economic growth, and 15 news attention indices related to financial markets.
table[table omitted — 5,357 chars of source]
table[table omitted — 2,138 chars of source]

Estimation Results

The numerical implementation is the same as that documented in Appendix (ref). In constructing the asymptotic covariance estimator, we set the threshold parameter to \(\eta = 2\sqrt{h \log N / T}\). To determine the optimal lag order in the local projection regressions, we employ the information criterion proposed in Section (ref). As discussed earlier, it is sufficient to consider relatively small forecasting horizons when selecting the lag length \(p\). Accordingly, we set \(h=1,2,3\) and, for each horizon, compute the values of \(\text{IC}(p)\) for lag orders up to a maximum of ten. The resulting information criterion values are reported in Table (ref). The results indicate that the criterion favors \(\widehat{p}=5\) when \(h=1\), and \(\widehat{p}=3\) for both \(h=2\) and \(h=3\). Based on these findings, we adopt \(\widehat{p}=3\) as the common lag order in the subsequent analysis.

The estimated local projection coefficients for different horizons ($h=1,6,12,24$) using Model 1 are reported in Figure (ref). In the short run, the recession-attention news index exerts a broadly positive effect on industry-level return volatility. This result coincides with the expectation that recession-related news can raise perceived macroeconomic uncertainty in the market, which can in turn increase the cross-sectional dispersion of firms' economic activities (e.g., investment) and eventually increases asset-return volatility. This pattern is consistent with a large literature confirming that positive contribution of the market uncertainty to the asset return volatility Bloom2009. The effect turns weak at medium horizons (e.g., $h=6$) and becomes negligible at longer horizons (e.g., $h=24$), indicating that the volatility impulse triggered by recession news is predominantly transitory. Similar patterns are observed in Figures (ref) and (ref), which present the estimated coefficients for Models 2 and 3. This indicates that the positive impact of recession-related news attention on volatility remains robust to the inclusion of additional news indices. To document cross-industry heterogeneity more precisely, Table (ref) reports the estimated short-run ($h=1$) impulse responses and their standard errors for each industry; several industries exhibit estimates exceeding 0.2, whereas many others remain indistinguishable from zero. These heterogeneous patterns are economically plausible: the short-run responses are strongest in cyclical and finance-sensitive sectors (for example, construction, transportation equipment, and finance sectors), whereas defensive sectors (such as mining, food or tobacco products) display statistically insignificant responses. Table (ref) also presents the short-run impulse-response estimates from Models 2 and 3. Even after adding attention series for other news topics that are related to economic growth or financial markets, the recession topic maintains positive impacts on volatility across most industries, though several estimated effects decline in magnitude.

To assess the temporal evolution of recession-driven volatility shocks, we compute impulse responses at horizons spanning one month to two years for each industry. Owing to page constraints, Figure (ref) displays the estimated impulse responses and their 90% confidence intervals for six representative industries (AG, OG, ST, CN, MF, and FI). The results confirm positive, statistically significant short-run effects of recession-related news attention. For example, the value for ST industry remains significantly positive through horizon $h=4$, while the responses for CN and OG industries are significantly positive at horizons up to 3. At medium and long horizons, on the other hand, the effects become statistically insignificant as confidence bands widen. These findings corroborate those reported in Figure (ref).

In addition to analyzing the effects of recession-related news on volatility, Models 2 and 3 allow us to examine the influence of other economics- and finance-related news topics. As shown in Figures (ref) and (ref), these topics exert divisive effects on industry-level volatilities. In particular, the topic {\it European Sovereign Debt} (ESD) exhibits a strong short-run effect in increasing volatility, and this effect is broadly uniform across industries. This finding is not surprising, as news related to sovereign debt crises, which is similar to the recession-related news, tends to heighten market perceptions of macroeconomic and financial uncertainty, thereby amplifying volatility. By contrast, the {\it Federal Reserve} (FED) topic displays a pervasive volatility-reducing effect across industries. One possible explanation is that Federal Reserve–related news often conveys policy guidance and information about the future path of monetary policy, which can reduce uncertainty by anchoring market expectations.

We then turn to examine volatility spillovers across industries. In our estimates, the Finance, Insurance, and Real Estate (FI) sector emerges as a prominent source of volatility spillovers to other industries at both short and medium horizons, as shown in Figure (ref). This result is consistent with the broader literature, which finds that financial sectors often act as net transmitters of volatility due to their central role in credit intermediation and risk sharing diebold2015financial. It also extends the existing literature that documents risk spillovers among financial institutions themselves tobias2016covar, and further underscores the pivotal role of the financial sector in propagating volatility shocks across the entire economic network. By contrast, the transportation industry appears to exert a restraining effect on volatility spillovers, which indicates sectors that are more closely tied to final demand and production may absorb volatility rather than amplify it.

As a robustness check, we replace historical volatility with realized volatility. Specifically, monthly realized volatility is constructed as the square root of the sum of squared daily returns for each industry. We then re-estimate the model using the same local projection framework. Due to page constraints, the corresponding results are reported in Figures (ref), (ref), and (ref) in Appendix (ref). The main empirical findings remain largely unchanged. Recession-related news attention continues to exert a positive effect on volatility across most industries, and the presence of volatility spillovers from the financial sector to other sectors is again confirmed.

table[table omitted — 637 chars of source]
table[table omitted — 4,976 chars of source]
sidewaysfigure[htp!] \subfloat[$h=1$ ] \subfloat[$h=6$ ] \\ \subfloat[$h=12$ ] \subfloat[$h=24$ ] \\ \caption{Estimated High-Dimensional Local Projection Coefficients for Model 1. Notes: Vertical axis shows the local projection response of the dependent variable at horizon \(t+h\); horizontal axis lists the predictor variable. Responses are shown separately for horizons \(h\in\{1,6,12,24\}\). }
sidewaysfigure[htp!] \subfloat[$h=1$ ] \subfloat[$h=6$ ] \\ \subfloat[$h=12$ ] \subfloat[$h=24$ ] \\ \caption{Estimated High-Dimensional Local Projection Coefficients for Model 2. Notes: Vertical axis shows the local projection response of the dependent variable at horizon \(t+h\); horizontal axis lists the predictor variable. Responses are shown separately for horizons \(h\in\{1,6,12,24\}\). }
sidewaysfigure[htp!] \subfloat[$h=1$ ] \subfloat[$h=6$ ] \\ \subfloat[$h=12$ ] \subfloat[$h=24$ ] \\ \caption{Estimated High-Dimensional Local Projection Coefficients for Model 3. Notes: Vertical axis shows the local projection response of the dependent variable at horizon \(t+h\); horizontal axis lists the predictor variable. Responses are shown separately for horizons \(h\in\{1,6,12,24\}\). For brevity, estimates for spillovers within the volatility series are omitted; only the responses of volatility to news-attention indices are reported. }
figure[figure omitted — 875 chars of source]

Conclusion

In this paper, we rigorously analyze the properties of the LP methodology within a HD framework, with a central focus on achieving robust long-horizon inference. We integrate a general dependence structure into $h$-step ahead forecasting models via a flexible specification of the residual terms. Additionally, we study the corresponding HD covariance matrix estimation, explicitly addressing the complexity arising from the long-horizon setting. Extensive Monte Carlo simulations are conducted to substantiate the derived theoretical findings.

In the empirical study, we utilize the proposed HD LP framework to revisit the impact of business news attention on U.S. industry-level stock volatility. Though the structural VAR models, as a powerful toolkit, have been widely used in the relevant literature diebold2009measuring, diebold2014network to trace how volatility shocks propagate across markets or industries, such methods impose rigid assumptions on system dynamics. Our method, on the other hand, allows for more flexibility in capturing the effect of news shocks and the volatility spillover effects among industries. Finally, we compute impulse-response functions at different horizons, and summarize how news-driven volatility shocks evolve over time.

{