#  Discovering the Network Granger Causality in Large Vector Autoregressive
Models

Yoshimasa Uematsu Correspondence: Yoshimasa Uematsu, Department of Social Data
Science, Hitotsubashi University, 2-1 Naka, Kunitachi, Tokyo 186-8601, Japan
(E-mail: yoshimasa.uematsu@r.hit-u.ac.jp).  Takashi Yamagata†

###### Abstract

This paper proposes novel inferential procedures for the network Granger
causality in high-dimensional vector autoregressive models. In particular, we
offer two multiple testing procedures designed to control discovered networks’
false discovery rate (FDR). The first procedure is based on the limiting
normal distribution of the $t$-statistics constructed by the debiased lasso
estimator. The second procedure is based on the bootstrap distributions of the
$t$-statistics made by imposing the null hypotheses. Their theoretical
properties, including FDR control and power guarantee, are investigated. The
finite sample evidence suggests that both procedures can successfully control
the FDR while maintaining high power. Finally, the proposed methods are
applied to discovering the network Granger causality in a large number of
macroeconomic variables and regional house prices in the UK.

Keywords. Multiple testing, FDR and power, Debiased lasso, Bootstrap.

##  1 Introduction

Revealing a dynamic interrelationship among variables is a critical challenge
in economics, finance, neuroscience, genomics, and so forth. In particular,
identifying the ability of a time series variable to predict the future values
of other time series variables has been of great interest in various fields.
Following the vast literature published after Granger (1969), when the past
values of a time series $x:=\\{x_{t}\\}$ can predict the future values of
another time series $y$, this will be expressed as “$x$ is Granger-causal for
$y$” in this article. This Granger causality has conventionally been discussed
in the bivariate relationship, typically with (bivariate) vector
autoregressive (VAR) models; see Sims (1972) and Hosoya (1977), among many
others.

###  1.1 Network Granger causality

Suppose that an $N$-dimensional time series $\mathbf{y}$ follows the
stationary VAR($K$) model:

$\displaystyle\mathbf{y}_{t}=\bm{\Phi}_{1}\mathbf{y}_{t-1}+\cdots+\bm{\Phi}_{K}\mathbf{y}_{t-K}+\mathbf{u}_{t},$

where $\bm{\Phi}_{k}=(\phi_{ij,k})$ is the $k$th coefficient matrix and
$\mathbf{u}_{t}$ is the error vector. The predictability relationship among
the $N$ time series manifests as the sparsity pattern of the coefficient
matrices in the VAR($K$) model. For this VAR model, the Granger-causal network
is defined as the graph $G=(V,E)$ with vertex set $V=\\{1,\dots,N\\}=:[N]$ and
edge set $E$ such that for distinct $i,j\in V$, $i\to j\in E$ if and only if
$\phi_{ji,k}\not=0$ for some $k\in[K]$. The same definition has been adopted
by Basu et al. (2015) and Eichler (2007, 2012b),111He includes contemporaneous
connections via non-zero error covariance in his causal graph called the path
diagram, but we do not consider them in this article. for instance. See also
Shojaie and Fox (2021) for recent advances in the analysis of network Granger
causality. As an illustration, we consider the 4-dimensional VAR(2) model with
the sparsity pattern of the coefficient matrices shown in Figure 1(a), where
the non-zero and zero elements in $\bm{\Phi}_{1}$ and $\bm{\Phi}_{2}$ are
indicated as the gray and white cells, respectively. From this sparse
structure, we readily obtain the associated Granger-causal network in Figure
1(b).

$\bm{\Phi}_{1}$$\bm{\Phi}_{2}$ (a) Sparsity pattern of the coefficients

$y_{1}$$y_{2}$$y_{3}$$y_{4}$ (b) Granger-causal network

Figure 1: Sparse VAR$(2)$ coefficients and associated Granger-causal network

Although the network Granger causality has been simply defined as the direct
causal network as above, we have to remark the existence of further indirect
causalities over the lagged variables. For the aforementioned VAR(2) model,
Figure 2 depicts the causal chain that takes all the lags into account
(Eichler, 2012a, Sec. 3), where the chain continues to the infinite past and
future due to the stationarity. In this figure, we can identify many indirect
causalities; for instance, $y_{3,t-2}$ is indirectly causal to $y_{2t}$ via
$y_{1,t-1}$. To refine the concept of such direct and indirect causalities,
Dufour and Renault (1998) introduced short- and long-run causality called the
$h$-step ($h=1,2,\dots$) causality. In this context, the $h$-step non-
causality is characterized by zero restrictions to a part of the first $N$
rows of the first to $h$th powers of the companion coefficient matrices
(Lütkepohl, 2005, p.50). Indeed, in Figure 2, it can be shown that the one-
step (i.e., direct) causality is expressed by the thick arrows in Figure 2
while the aforementioned indirect effects will be included in $h$-step
causalities for $h=2,3,\dots$. Apparently, this one-step causality corresponds
to the network Granger causality we have defined.

$y_{1}^{t+2}$$y_{2}^{t+2}$$y_{3}^{t+2}$$y_{4}^{t+2}$$y_{1}^{t+1}$$y_{2}^{t+1}$$y_{3}^{t+1}$$y_{4}^{t+1}$$y_{1}^{t}$$y_{2}^{t}$$y_{3}^{t}$$y_{4}^{t}$$y_{1}^{t-1}$$y_{2}^{t-1}$$y_{3}^{t-1}$$y_{4}^{t-1}$$y_{1}^{t-2}$$y_{2}^{t-2}$$y_{3}^{t-2}$$y_{4}^{t-2}$$y_{1}^{t-3}$$y_{2}^{t-3}$$y_{3}^{t-3}$$y_{4}^{t-3}$$y_{1}^{t-4}$$y_{2}^{t-4}$$y_{3}^{t-4}$$y_{4}^{t-4}$
Figure 2: Causal chain of the VAR(2)

The Granger-causal network can be discovered just by estimating the sparse
$\bm{\Phi}_{k}$’s in an appropriate manner. In early empirical studies,
including Fujita et al. (2007) and Lozano et al. (2009) among others, such
statistical methods were already employed and a network Granger causality was
estimated. Recently the theoretical properties of these methods have begun to
be investigated. Basu and Michailidis (2015), Kock and Callot (2015), and Basu
et al. (2019) study properties of the lasso (Tibshirani, 1996) and its
variants in high-dimensional VAR models. Han et al. (2015) extends the Dantzig
selector of Candès and Tao (2007) to be applicable to (weakly) sparse VAR
models. Kock and Callot (2015) and Barigozzi and Brownlees (2019) propose the
adaptive lasso estimator. Davis et al. (2016) consider an alternative two-step
estimator, which uses estimates of the partial spectral coherence.

###  1.2 Classical Granger causality test

Estimation-based network detection such as the lasso selection is appealing
because it is simple to use. However, it is quite unstable, and this, along
with the difficulty of evaluating the type I error of the selection result,
may cause a lack of reproducibility. Thus, this article discusses the
statistical inference of the Granger causality network. Indeed, in Granger
(1969), statistical inference of Granger causality was the central issue, as
Sir Clive Granger himself put it: “the problem is how to devise a definition
of causality and feedback and to test for their existence” (Granger, 1969,
p.428).

Regardless of the model dimensionality, it is possible to test the null of
Granger non-causality if an asymptotically normal estimator is available. For
the purpose of illustration, consider an $N$-dimensional VAR(1) model. Given
$\mathcal{H}\subset[N]\times[N]$, we may test $H_{0}:\phi_{ij}=0$ for all
$(i,j)\in\mathcal{H}$ versus $H_{1}:\phi_{ij}\not=0$ for some
$(i,j)\in\mathcal{H}$ using the asymptotic Wald statistic constructed by the
debiased lasso estimator (Zheng and Raskutti, 2019; Zhu and Liu, 2020; Babii
et al., 2021). This is effective when a very small $\mathcal{H}$ is of our
interest. Unfortunately, however, this approach is not suitable for
discovering the virtually high-dimensional Granger causality networks in
$\mathcal{H}$. This is because, when $\mathcal{H}$ is large, the rejection of
this null hypothesis tells us only that there may exist Granger causality
between some of the $N$ variables. This is not informative for our purpose.

###  1.3 Discovering the network with FDR control

The network discovering is characterized as the multiple testing for the
sequence of hypotheses, $H_{0}^{(i,j)}:\phi_{ij}=0$ versus
$H_{1}^{(i,j)}:\phi_{ij}\not=0$ for each $(i,j)\in\mathcal{H}$. The false
discovery rate (FDR), which is considered as a measure of type I error in
multiple tests, is defined as
$\operatorname{FDR}=\operatorname{\mathbb{E}}[{|\hat{\mathcal{S}}\cap\mathcal{S}^{c}|}/{(|\hat{\mathcal{S}}|\vee
1)}]$ with $\mathcal{S}=\\{(i,j)\in\mathcal{H}:\phi_{ij}\not=0\\}$ and
$\hat{\mathcal{S}}=\\{(i,j)\in\mathcal{H}:H_{0}^{(i,j)}\text{ is rejected}\\}$
(Benjamini and Hochberg, 1995). The FDR controlled multiple testing is
expected to exhibit higher power
(Power=$\operatorname{\mathbb{E}}[{|\hat{\mathcal{S}}\cap\mathcal{S}|}/{|\mathcal{S}|}]$)
than that controlling the family-wise error rate
($\operatorname{FWER}=\operatorname{\mathbb{P}}(|\hat{\mathcal{S}}\cap\mathcal{S}^{c}|\geq
1)$), typically by the so-called Bonferroni correction (Bonferroni, 1935;
Holm, 1979). This is because the FDR is smaller than or equal to the FWER; the
FDR control allows more “aggressive” selection.

In this paper, we propose two novel inferential methods for discovering the
network Granger causality in high-dimensional VAR models, based on asymptotic
and bootstrap $t$-statistics, respectively, with the debiased lasso estimator.
The key aspect is that they are designed to control the FDR, which will lead
to a more stable discovery. Their theoretical properties, including the FDR
control and power guarantee, are investigated. To the best of our knowledge,
this is the first paper to propose such inferential procedures for high-
dimensional VAR models and to investigate their theoretical properties. The
finite sample evidence suggests that both of the procedures can successfully
control the FDR while maintaining the high power. The proposed methods are
applied to a large dataset of macroeconomic time series and regional house
prices in the UK.

Finally, there are a few studies that consider Granger causality between
groups of variables, typically using group lasso; see Basu et al. (2015), Lin
and Michailidis (2017) and Basu et al. (2019). This approach requires the
researcher to know the group members. Guðmundsson and Brownlees (2021) develop
methods to statistically identify groups among the variables. We also notice
that the FDR control for time series is essentially quite difficult, and there
are very few studies. One exception is Chi et al. (2021), but this does not
appear readily applicable to our purpose.

###  1.4 Organization and notation

The paper is organized as follows. Section 2 formally defines the VAR model.
Section 3 proposes two methods for discovering the networks based on a
multiple testing. Section 4 explores the statistical theory for the FDR
control and power guarantee of our methods. Section 5 confirms the finite
sample validity via Monte Carlo experiments. Section 6 applies our methods to
large datasets. Section 7 concludes. All the proofs of our theoretical results
and supplemental analyses are collected in Supplementary Materials.

Notation. For any matrix $\mathbf{M}=(m_{ti})\in\mathbb{R}^{T\times N}$,
denote by $\|\mathbf{M}\|_{\mathrm{F}}$, $\|\mathbf{M}\|_{2}$,
$\|\mathbf{M}\|_{1}$, $\|\mathbf{M}\|_{\max}$, and $\|\mathbf{M}\|_{\infty}$
the Frobenius norm, induced $\ell_{2}$ (spectral) norm, entrywise
$\ell_{1}$-norm, entrywise $\ell_{\infty}$-norm, and induced
$\ell_{\infty}$-norm, respectively. Specifically, they are defined by
$\|\mathbf{M}\|_{\mathrm{F}}=(\sum_{t,i}m_{ti}^{2})^{1/2}$,
$\|\mathbf{M}\|_{2}=\lambda_{1}^{1/2}(\mathbf{M}^{\prime}\mathbf{M})$,
$\|\mathbf{M}\|_{1}=\sum_{t,i}|m_{ti}|$,
$\|\mathbf{M}\|_{\max}=\max_{t,i}|m_{ti}|$, and
$\|\mathbf{M}\|_{\infty}=\max_{t}\sum_{i}|m_{ti}|$, where
$\lambda_{i}(\mathbf{S})$ refers to the $i$th largest eigenvalue of a
symmetric matrix $\mathbf{S}$. Denote by $\mathbf{I}_{N}$ and
$\mathbf{0}_{T\times N}$ the $N\times N$ identity matrix and $T\times N$
matrix with all the entries being zero, respectively. We use $\lesssim$
($\gtrsim$) to represent $\leq$ ($\geq$) up to a positive constant factor. For
any positive sequence $a_{n}$ and $b_{n}$, we write $a_{n}\asymp b_{n}$ if
$a_{n}\lesssim b_{n}$ and $a_{n}\gtrsim b_{n}$. Moreover, denote by $a_{n}\sim
b_{n}$ if $a_{n}/b_{n}\to 1$. For any positive values $a$ and $b$, $a\vee b$
and $a\wedge b$ stand for $\max(a,b)$ and $\min(a,b)$, respectively. The
indicator function is denoted by $1\\{\cdot\\}$. $\operatorname{sgn}$ is the
sign function defined as $\operatorname{sgn}(x)=x/|x|$ for $x\not=0$ and
$\operatorname{sgn}(0)=0$.

##  2 Model

Suppose that the $N$-dimensional vector of stationary time series,
$\mathbf{y}_{t}=(y_{1t},\dots,y_{Nt})^{\prime}$, is generated from the
VAR($K$) model:

$\displaystyle\mathbf{y}_{t}$
$\displaystyle=\bm{\Phi}_{1}\mathbf{y}_{t-1}+\dots+\bm{\Phi}_{K}\mathbf{y}_{t-K}+\mathbf{u}_{t}=\bm{\Phi}\mathbf{x}_{t}+\mathbf{u}_{t},$
(1)

where $\bm{\Phi}_{k}$’s are the $N\times N$ sparse coefficient matrices with
$\bm{\Phi}_{k}=(\phi_{ij,k})$,
$\mathbf{x}_{t}=(\mathbf{y}_{t-1}^{\prime},\dots,\mathbf{y}_{t-K}^{\prime})^{\prime}$
is the $KN$-dimensional vector of lagged variables, and
$\mathbf{u}_{t}=(u_{1t},\dots,u_{Nt})^{\prime}$ is an error vector with mean
zero and finite positive definite covariance matrix
$\bm{\Sigma}_{u}=(\sigma_{ij})=\operatorname{\mathbb{E}}\mathbf{u}_{t}\mathbf{u}_{t}^{\prime}$.
Let
$\bm{\Sigma}_{x}=\operatorname{\mathbb{E}}\mathbf{x}_{t}\mathbf{x}_{t}^{\prime}$
denote the finite positive definite covariance matrix. Moreover, define
$\bm{\Gamma}_{y}(k-1)=\operatorname{\mathbb{E}}[\mathbf{y}_{t-1}\mathbf{y}_{t-k}^{\prime}]$
with $\bm{\Gamma}_{y}(k-1)=\bm{\Gamma}_{y}(-k+1)^{\prime}$ for $k\in[K]$. Then
$\bm{\Sigma}_{x}$ is composed of submatrices $\bm{\Gamma}_{y}(k-1)$. Moreover,
define the precision matrix of $\mathbf{x}_{t}$ as
$\bm{\Omega}=(\omega_{ij})=\bm{\Sigma}_{x}^{-1}$. We also denote
$\sigma_{i}^{2}=\sigma_{ii}$ and $\omega_{i}^{2}=\omega_{ii}$. Denote by
$\bm{\omega}_{j}$ the $j$th column vector of $\bm{\Omega}$ for $j\in[KN]$.
Throughout the paper, we allow $N$ to be as large as or possibly larger than
$T$, and $N\wedge T$ tends to infinity while $K$ is fixed.

Stacking the $T$ observations in columns such that
$\mathbf{Y}=(\mathbf{y}_{1},\dots,\mathbf{y}_{T})\in\mathbb{R}^{N\times T}$,
$\mathbf{X}=(\mathbf{x}_{1},\dots,\mathbf{x}_{T})\in\mathbb{R}^{KN\times T}$,
and $\mathbf{U}=(\mathbf{u}_{1},\dots,\mathbf{u}_{T})\in\mathbb{R}^{N\times
T}$, we have the matrix form

$\displaystyle\mathbf{Y}=\bm{\Phi}\mathbf{X}+\mathbf{U},$ (2)

where $\bm{\Phi}=\left(\bm{\Phi}_{1},\dots,\bm{\Phi}_{K}\right)$. Denote by
$\mathbf{a}_{\cdot i}$ the $i$th row vector of matrix $\mathbf{A}$. Then, by
extracting the $i$th row of (2), the model is also written as

$\displaystyle\mathbf{y}_{i\cdot}=\bm{\phi}_{\cdot
i}\mathbf{X}+\mathbf{u}_{i\cdot},$ (3)

which is $1\times T$.

To describe the sparsity pattern of $\bm{\Phi}$, define
$\mathcal{S}=\\{(i,j)\in[N]\times[KN]:\phi_{ij}\not=0\\}$ with cardinality
$s=|\mathcal{S}|$. Similarly, the sparsity pattern of $\bm{\phi}_{i\cdot}$ is
described as $\mathcal{S}_{i}=\\{j\in[KN]:\phi_{ij}\not=0\\}$ with
$s_{i}=|\mathcal{S}_{i}|$ and $\bar{s}=\max_{i\in[N]}s_{i}$. There is a one-
to-one correspondence between $\mathcal{S}$ and the Granger-causal network as
mentioned in Section 1.1. Our goal is to discover the network controlling the
FDR.

##  3 Inferential Methodology

We propose two multiple testing procedures that can control the FDR of
discovering the Granger-causal networks, which are based on asymptotic and
bootstrap $t$-statistics, respectively. The associated statistical theory is
developed in Section 4.

###  3.1 Debiased lasso estimator

We start with constructing the row-wise lasso estimator,
$\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}$, defined as

$\displaystyle\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}=\operatorname*{\arg\min}_{\bm{\phi}_{i\cdot}\in\mathbb{R}^{1\times
KN}}~{}(2T)^{-1}\|\mathbf{y}_{i\cdot}-\bm{\phi}_{i\cdot}\mathbf{X}\|_{2}^{2}+\lambda\|\bm{\phi}_{i\cdot}\|_{1},$
(4)

where $\lambda>0$ is a regularization parameter. It is well-known that the
lasso estimator has a bias caused by this regularization. Following Javanmard
and Montanari (2014), we remove the bias, which leads to the asymptotic
normality of each element; see also van de Geer et al. (2014) and Zhang and
Zhang (2014) for a related discussion. Define the debiased lasso estimator as

$\displaystyle\hat{\bm{\phi}}_{i\cdot}=\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}+(\mathbf{y}_{i\cdot}-\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}\mathbf{X})\mathbf{X}^{\prime}\hat{\bm{\Omega}}/T,$
(5)

where $\hat{\bm{\Omega}}$ is a consistent estimator of the precision matrix
$\bm{\Omega}$. Let $\hat{\bm{\Sigma}}_{x}=\mathbf{X}\mathbf{X}^{\prime}/T$.
Then by a simple calculation, (5) is equivalently written as

$\displaystyle\sqrt{T}(\hat{\bm{\phi}}_{i\cdot}-\bm{\phi}_{i\cdot})=\mathbf{z}_{i\cdot}+\mathbf{r}_{i\cdot},$
(6)

where

$\displaystyle\mathbf{z}_{i\cdot}=\mathbf{u}_{i\cdot}\mathbf{X}^{\prime}\bm{\Omega}/\sqrt{T}~{}~{}~{}\mbox{and}~{}~{}~{}\mathbf{r}_{i\cdot}=\mathbf{u}_{i\cdot}\mathbf{X}^{\prime}(\hat{\bm{\Omega}}-\bm{\Omega})/\sqrt{T}-\sqrt{T}(\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}-\bm{\phi}_{i\cdot})(\hat{\bm{\Sigma}}_{x}\hat{\bm{\Omega}}-\mathbf{I}_{KN}).$

For each $i\in[N]$, it is expected that each element of $\mathbf{z}_{i\cdot}$
can be asymptotically normal while $\|\mathbf{r}_{i\cdot}\|_{\max}$ becomes
negligible for a relevant choice of the estimator, $\hat{\bm{\Omega}}$. The
asymptotic behavior of this estimator is formally investigated in Section 4.1.

###  3.2 Multiple test

For $\mathcal{H}\subset[N]\times[KN]$, consider discovering the underlying
Granger causal network in $\mathcal{H}$. This problem is understood as the
multiple test for the sequence of hypotheses:

$\displaystyle
H_{0}^{(i,j)}:\phi_{ij}=0~{}~{}~{}\mbox{versus}~{}~{}~{}H_{1}^{(i,j)}:\phi_{ij}\not=0~{}~{}\text{
for each }~{}~{}(i,j)\in\mathcal{H}.$ (7)

As observed in Section 3.1, under regularity conditions, each entry of the
debiased lasso estimator will be expressed as
$\hat{\phi}_{ij}=z_{ij}+o_{p}(1)$, where the leading term $z_{ij}$ is expected
to be asymptotically normal with the variance,

$\displaystyle\operatorname{Var}(z_{ij})=\sigma_{i}^{2}\bm{\omega}_{j}^{\prime}\bm{\Sigma}_{x}{\bm{\omega}}_{j}=\sigma_{i}^{2}\omega_{j}^{2}.$
(8)

Due to the expression of (8), the $t$-test for each pair of hypotheses in (7)
is performed with the $t$-statistic, either

$\displaystyle\texttt{T}_{ij}=\frac{\sqrt{T}\hat{\phi}_{ij}}{\hat{\sigma}_{i}\sqrt{\hat{\bm{\omega}}^{\prime}_{j}\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}}}\text{~{}~{}or~{}~{}}\frac{\sqrt{T}\hat{\phi}_{ij}}{\hat{\sigma}_{i}\hat{\omega}_{j}},$
(9)

where the denominator (standard error, s.e.) is composed of some consistent
estimators. Hereafter, denote by $\hat{m}_{ij}$ either
$\hat{\sigma}_{i}\sqrt{\hat{\bm{\omega}}^{\prime}_{j}\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}}$
or $\hat{\sigma}_{i}\hat{\omega}_{j}$. Repeating the $t$-test over
$(i,j)\in\mathcal{H}$ with a critical value t leads to a set of discoveries,
$\hat{\mathcal{S}}(\texttt{t}):=\\{(i,j)\in\mathcal{H}:|\texttt{T}_{ij}|\geq\texttt{t}\\}$.
Here, the choice of t is critical; we propose two procedures to determine t
that achieves controlling the directional FDR of
$\hat{\mathcal{S}}(\texttt{t})$,

$\displaystyle\operatorname{dFDR}=\operatorname{\mathbb{E}}\left[\operatorname{dFDP}\right],~{}~{}~{}\operatorname{dFDP}=\frac{|\\{(i,j)\in\hat{\mathcal{S}}(\texttt{t}):\operatorname{sgn}(\hat{\phi}_{ij})\not=\operatorname{sgn}(\phi_{ij})\\}|}{|\hat{\mathcal{S}}(\texttt{t})|\vee
1},$

below a preassigned level. Notice that the $\operatorname{dFWER}$ and
$\operatorname{dPower}$ are also defined in a similar manner.

There are several ways to construct consistent estimators of the nuisance
parameters, $\sigma_{i}^{2}$ and $\bm{\omega}_{j}$. In this paper, we use
$\hat{\sigma}_{i}^{2}=\sum_{t=1}\hat{u}_{it}^{2}/(T-d_{i})$ with
$\hat{u}_{it}$ the $(i,t)$th element of the lasso residual matrix
$\hat{\mathbf{U}}$ and the CLIME estimator (Cai et al., 2011)
$\hat{\bm{\Omega}}=(\hat{\bm{\omega}}_{1},\dots,\hat{\bm{\omega}}_{KN})$,
respectively. Here, $d_{i}=o(T)$ is a positive number for degrees of freedom
adjustment. A typical choice is $d_{i}=\hat{s}_{i}$, where
$\hat{s}_{i}=|\hat{\mathcal{S}}_{i}^{{\textsf{L}}}|$ with
$\hat{\mathcal{S}}_{i}^{{\textsf{L}}}=\operatorname{supp}(\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}})$.

####  3.2.1 First procedure: Limiting normal distribution

We construct the set of discoveries as follows.

######  Procedure 1.

  1. 1.

Set $\bar{\texttt{t}}=\sqrt{2\log(|\mathcal{H}|)-a\log\log(|\mathcal{H}|)}$
for given $\mathcal{H}\subset[N]\times[KN]$, where $a>3$ is an arbitrary fixed
constant.

  2. 2.

For any level $q\in[0,1]$, compute

$\displaystyle\texttt{t}_{0}=\inf\left\\{\texttt{t}\in[0,\bar{\texttt{t}}]:\frac{2|\mathcal{H}|Q(\texttt{t})}{|\hat{\mathcal{S}}(\texttt{t})|\vee
1}\leq q\right\\},$ (10)

where
$|\hat{\mathcal{S}}(\texttt{t})|=\sum_{(i,j)\in\mathcal{H}}1\\{|\texttt{T}_{ij}|\geq\texttt{t}\\}$
is the total number of discoveries in $\mathcal{H}$ with critical value t, and
$Q(\texttt{t})=\operatorname{\mathbb{P}}(\mathcal{Z}>\texttt{t})$ with
$\mathcal{Z}$ a standard normal random variable. If (10) does not exist, set

$\displaystyle\texttt{t}_{0}=\sqrt{2\log(|\mathcal{H}|)}.$ (11)

  3. 3.

For each $(i,j)\in\mathcal{H}$, reject $H_{0}^{(i,j)}$ in (7) if
$|\texttt{T}_{ij}|\geq\texttt{t}_{0}$, and obtain the set of discoveries,
$\hat{\mathcal{S}}(\texttt{t}_{0})=\\{(i,j)\in\mathcal{H}:|\texttt{T}_{ij}|\geq\texttt{t}_{0}\\}$.

This procedure is designed to asymptotically control the dFDR of
$\hat{\mathcal{S}}(\texttt{t}_{0})$ to be less than or equal to $q$. A similar
procedure can be found in Liu (2013), Javanmard and Javadi (2019), and Uematsu
and Yamagata (2021), for instance. The $\operatorname{dFDR}$ control always
implies the $\operatorname{FDR}$ control since
$\operatorname{FDR}\leq\operatorname{dFDR}$. Threshold (11) will even control
the dFWER (and hence the dFDR), which occurs if it is large enough not to
exist in $[0,\bar{\texttt{t}}]$. A theoretical justification for the
$\operatorname{dFDR}$ control and $\operatorname{dPower}$ guarantee is given
in Section 4.2.

######  Remark 1.

  1. (a)

The constant $a>3$ is required for technical reasons. In practice, we can
choose an arbitrary value that is slightly larger than three.

  2. (b)

A choice of $d_{i}$ effects the FDR control to some extent, but
$d_{i}=\hat{s}_{i}$ seems to bring good results by simulation studies in
Section 5. For other choices of $d_{i}$ and constructions of consistent
estimator of $\sigma_{i}^{2}$, see Reid et al. (2016).

####  3.2.2 Second procedure: Bootstrapped distribution

Even in a low-dimensional setting, Brüggemann et al. (2016) point out that the
finite sample properties of asymptotic VAR inference can be rather poor, and
the use of bootstrap methods is often advocated. See also Kilian (1999),
Gonçalves and Kilian (2004), Hafner and Herwartz (2009), indicating that
bootstrap methods can be very beneficial in (V)ARs.

To improve the performance, we propose a second method based on the fixed-
design wild bootstrap (FWB). The bootstrapped multiple $t$-test we propose is
different from a conventional bootstrapped $t$-test for a single hypothesis.
Roughly speaking, it attempts to replicate a bootstrapped distribution under
the null by the FWB, which is substituted for the limit normal distribution
$Q$ in Procedure 1.

Let $\tilde{\mathcal{S}}\subset[N]\times[KN]$ denote an index set of
discoveries that is expected to satisfy
$\operatorname{\mathbb{P}}(\tilde{\mathcal{S}}\supset\mathcal{S})\to 1$. The
next procedure is available for $\mathcal{H}$ such that
$\tilde{\mathcal{S}}^{c}\cap\mathcal{H}\not=\emptyset$.

######  Procedure 2.

  1. 1.

Obtain bootstrap version of the $t$-statistics,
$\\{\texttt{T}_{ij}^{*(b)}:(i,j)\in\tilde{\mathcal{S}}^{c}\cap\mathcal{H}\\}_{b=1}^{B}$,
by repeating (a)–(f) $B$ times:

    1. (a)

Generate $\\{\zeta_{t}\\}$, a sequence of i.i.d. random variables with mean
zero and variance one.

    2. (b)

Obtain $\\{\mathbf{u}_{t}^{*}\\}_{t=1}^{T}$ by the FWB;
$\mathbf{u}_{t}^{*}=\hat{\mathbf{u}}_{t}\zeta_{t}$, where
$\hat{\mathbf{u}}_{t}=\mathbf{y}_{t}-\hat{\bm{\Phi}}^{{\textsf{L}}}\mathbf{x}_{t}$.

    3. (c)

Generate $\\{\mathbf{y}_{t}^{*}\\}_{t=1}^{T}$ by
$\mathbf{y}_{t}^{*}=\hat{\bm{\Phi}}^{{\textsf{L}}}\mathbf{x}_{t}+\mathbf{u}_{t}^{*}$,
and set $\mathbf{Y}^{*}=(\mathbf{y}_{1}^{*},\dots,\mathbf{y}_{T}^{*})$.

    4. (d)

Compute a bootstrap version of the lasso estimate
$\hat{\bm{\Phi}}^{{\textsf{L}}*}$ by $\mathbf{Y}^{*}$ and $\mathbf{X}$.

    5. (e)

Construct a bootstrap version of the debiased lasso estimate

$\displaystyle\hat{\bm{\Phi}}^{*}=\hat{\bm{\Phi}}^{{\textsf{L}}*}+\left(\mathbf{Y}^{*}-\hat{\bm{\Phi}}^{{\textsf{L}}*}\mathbf{X}\right)\mathbf{X}^{\prime}\hat{\bm{\Omega}}/T.$

    6. (f)

Construct a bootstrap version of the $t$-statistics, either

$\displaystyle\texttt{T}_{ij}^{*}=\frac{\sqrt{T}\hat{\phi}_{ij}^{*}}{\hat{\sigma}^{*}_{i}\sqrt{\hat{\bm{\omega}}^{\prime}_{j}\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}}}\text{~{}~{}or~{}~{}}\frac{\sqrt{T}\hat{\phi}_{ij}^{*}}{\hat{\sigma}^{*}_{i}\hat{\omega}_{j}}~{}~{}~{}\text{for}~{}~{}~{}(i,j)\in\tilde{\mathcal{S}}^{c}\cap\mathcal{H},$
(12)

where $\hat{\sigma}^{*2}_{i}$ is given by
$\hat{\sigma}_{i}^{*2}=\sum_{t=1}^{T}\hat{u}_{it}^{*2}/(T-\hat{s}_{i})$ with
$\hat{\mathbf{U}}^{*}=\mathbf{U}^{*}-(\hat{\bm{\Phi}}^{{\textsf{L}}*}-\hat{\bm{\Phi}}^{{\textsf{L}}})\mathbf{X}$
and $\mathbf{U}^{*}=(\mathbf{u}_{1}^{*},\dots,\mathbf{u}_{T}^{*})$.

  2. 2.

Compute the empirical distribution

$\displaystyle\mathbb{Q}_{B}^{*}(\texttt{t})$
$\displaystyle=\frac{1}{|\tilde{\mathcal{S}}^{c}\cap\mathcal{H}|}\sum_{(i,j)\in\tilde{\mathcal{S}}^{c}\cap\mathcal{H}}\left[\frac{1}{B}\sum_{b=1}^{B}1\left\\{\texttt{T}_{ij}^{*(b)}>\texttt{t}\right\\}\right].$

  3. 3.

Run Procedure 1 with (13) replacing (10):

$\displaystyle\texttt{t}_{0}=\inf\left\\{\texttt{t}\in[0,\bar{\texttt{t}}]:\frac{|\mathcal{H}|\left\\{\mathbb{Q}_{B}^{*}(\texttt{t})+1-\mathbb{Q}_{B}^{*}(-\texttt{t})\right\\}}{|\hat{\mathcal{S}}(\texttt{t})|\vee
1}\leq q\right\\}.$ (13)

######  Remark 2.

  1. (a)

We may adopt several distributions for $\zeta_{t}$. For example, Mammen (1993)
suggests using $\zeta_{t}$ that takes values $\mp(\sqrt{5}\mp 1)/2$ with
probability $(\sqrt{5}\pm 1)/(2\sqrt{5})$, respectively, while Davidson and
Flachaire (2008) propose using the Rademacher random variable, $\zeta_{t}=\pm
1$ with probability $1/2$.

  2. (b)

A typical choice of $\tilde{\mathcal{S}}$ is the lasso-selected variables,
$\hat{\mathcal{S}}_{{\textsf{L}}}$. We may also use
$\hat{\mathcal{S}}(\texttt{t}_{0})$ obtained by Procedure 1.

  3. (c)

We do not need to assign large $B$ because many
($|\tilde{\mathcal{S}}^{c}\cap\mathcal{H}|$) $t$-values are used to form the
null distribution $\mathbb{Q}_{B}^{*}(\texttt{t})$.

  4. (d)

Unlike the construction of $\hat{\sigma}_{i}^{2}$ in Procedure 1, the degree
of freedom adjustment in $\hat{\sigma}_{i}^{*2}$ is not sensitive to the
results of resulting FDR control.

##  4 Statistical Theory

We develop a formal statistical theory for the inferential methodology
proposed in Section 3. Throughout this section, we set
$\mathcal{H}=[N]\times[KN]$, $d_{i}=0$, and
$\tilde{\mathcal{S}}=\hat{\mathcal{S}}_{{\textsf{L}}}$ to alleviate
unnecessary technical complications. We suppose $N$ is possibly larger than
$T$ at most polynomials of $T$.

######  Condition 1.

The error term, $\\{\mathbf{u}_{t}\\}$, is a sequence of i.i.d. sub-Gaussian
random vectors with mean zero and covariance matrix $\bm{\Sigma}_{u}$; for
every $N>0$, there exists some constant $c_{u}>0$ such that for all $x>0$,
$\max_{i\in[N]}\operatorname{\mathbb{P}}\left(|u_{it}|>x\right)\leq
2\exp(-x^{2}/c_{u})$.

######  Condition 2.

All the eigenvalues of the companion matrix of
$(\bm{\Phi}_{1},\dots,\bm{\Phi}_{K})$ are strictly less than one in modulus
uniformly in $N\in\mathbb{N}$. (Hereafter, denote the companion matrix by
$\mathbf{A}$.)

######  Condition 3.

For all $N\in\mathbb{N}$, there exists some constant $\gamma>0$ such that
$\gamma\leq\lambda_{\min}(\bm{\Sigma}_{x})\leq\lambda_{\max}(\bm{\Sigma}_{x})\leq
1/\gamma$.

Conditions 1 and 2 are commonly used in the literature. Condition 2 guarantees
that the VAR($K$) model is stable and stationary, so that it is inverted to
the VMA($\infty$) model:

$\displaystyle\mathbf{y}_{t}=\sum_{\ell=0}^{\infty}\mathbf{B}_{\ell}\mathbf{u}_{t-\ell},~{}~{}~{}~{}~{}b:=\sum_{\ell=0}^{\infty}\|\mathbf{B}_{\ell}\|_{\infty}<\infty$
(14)

uniformly in $N$, where $\mathbf{B}_{0}=\mathbf{I}_{N}$ and
$\mathbf{B}_{\ell}=\mathbf{J}^{\prime}\mathbf{A}^{\ell}\mathbf{J}$ with
$\mathbf{J}^{\prime}=(\mathbf{I}_{N},\mathbf{0}_{N\times(KN-N)})$ for
$\ell=1,2,\dots$. The summability condition in (14) is verified by Lemma 2 of
Supplementary Materials. Throughout this section, we set the lasso
regularization parameter in (4) to be

$\displaystyle\lambda=8bc_{uu}\sqrt{2(\nu+5)^{3}T^{-1}\log^{3}(N\vee T)},$

where $\nu>0$ is a fixed constant.

###  4.1 Theory for the debiased lasso estimator

The first proposition states that the lasso estimator,
$\hat{\bm{\Phi}}^{{\textsf{L}}}$, has nonasymptotic error bounds.

######  Proposition 1 (Nonasymptotic error bounds for the lasso).

If Conditions 1–3 are assumed, then the lasso estimator defined in (4)
satisfies the following inequalities with probability at least $1-O((N\vee
T)^{-\nu})$:

$\displaystyle(a)~{}$
$\displaystyle~{}\left\|\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}-\bm{\phi}_{i\cdot}\right\|_{2}\leq\frac{12\sqrt{s_{i}}\lambda}{\gamma-8bs_{i}\lambda},$
$\displaystyle(b)~{}$
$\displaystyle~{}\left\|\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}-\bm{\phi}_{i\cdot}\right\|_{1}\leq\frac{48s_{i}\lambda}{\gamma-8bs_{i}\lambda},$
$\displaystyle(c)~{}$
$\displaystyle~{}\frac{1}{T}\left\|(\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}-\bm{\phi}_{i\cdot})\mathbf{X}\right\|_{2}^{2}\leq\frac{144s_{i}\lambda^{2}}{\gamma-8bs_{i}\lambda}$

for all $i\in[N]$ such that $8bs_{i}\lambda<\gamma$.

Using Proposition 1, we next show the asymptotic linearity of the debiased
lasso estimator, $\hat{\bm{\Phi}}$. This requires some conditions on the
precision matrix, $\bm{\Omega}$.

######  Condition 4.

There exist some positive numbers $s_{\omega}$ and $M_{\omega}$ and some
constant $r\in[0,1)$ such that
$\max_{i\in[KN]}\sum_{j=1}^{KN}|\omega_{ij}|^{r}\leq s_{\omega}$ and
$\max_{j\in[KN]}\|{\bm{\omega}}_{j}\|_{1}\leq M_{\omega}$, where $s_{\omega}$
and $M_{\omega}$ can diverge as $N,T\to\infty$.

Condition 4 has frequently been used to derive the rate of convergence of the
CLIME estimator $\hat{\bm{\Omega}}$; see Cai et al. (2011) and Shu and Nan
(2019), for instance. This condition requires that $\bm{\Omega}$ is
(approximately) sparse. Especially when $r=0$, the condition reduces to the
exact sparsity assumption. Denote by $\\{\mathbf{e}_{j}\\}$ the standard basis
of $\mathbb{R}^{KN}$.

######  Proposition 2 (Asymptotic linearity of the debiased lasso).

The debiased lasso estimator defined in (5) has the representation
$\sqrt{T}(\hat{\bm{\Phi}}-\bm{\Phi})=\mathbf{Z}+\mathbf{R}$, where the
$(i,j)$th elements of $\mathbf{Z}$ and $\mathbf{R}$ are respectively given by

$\displaystyle
z_{ij}=T^{-1/2}\sum_{t=1}^{T}u_{it}\mathbf{x}_{t}^{\prime}{\bm{\omega}}_{j},~{}~{}~{}r_{ij}=\mathbf{u}_{i\cdot}\mathbf{X}^{\prime}(\hat{\bm{\omega}}_{j}-{\bm{\omega}}_{j})/\sqrt{T}-\sqrt{T}(\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}-\bm{\phi}_{i\cdot})(\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}-\mathbf{e}_{j}).$

Furthermore, if Conditions 1–4 are assumed, then the following inequality
holds with probability at least $1-O((N\vee T)^{-\nu})$:

$\displaystyle\max_{j\in[KN]}|r_{ij}|\lesssim\left(s_{\omega}M_{\omega}^{2-2r}\lambda^{1-r}+M_{\omega}s_{i}\lambda\right)\sqrt{\log^{3}(N\vee
T)}=:\bar{r}_{i}$

for all $i\in[N]$ such that $s_{i}\lambda=o(1)$.

To achieve the asymptotic normality of $\sqrt{T}(\hat{\phi}_{ij}-\phi_{ij})$,
we need $\bar{r}_{i}=o(1)$. Then $z_{ij}$ becomes the dominating term, which
can converge in distribution to a Gaussian random variable.

###  4.2 Theory for the multiple testing

We are now ready to develop a statistical theory of the multiple testing of
(7) by Procedures 1 and 2 in Section 3.2. We first establish the asymptotic
normality of the $t$-statistics in (9).

######  Condition 5.

For all $N\in\mathbb{N}$, there exists some constant $\gamma>0$ such that
$\gamma\leq\min_{i\in[N]}\sigma_{i}^{2}\leq\max_{i\in[N]}\sigma_{i}^{2}\leq
1/\gamma$ and
$\gamma\leq\min_{j\in[KN]}\omega_{j}^{2}\leq\max_{j\in[KN]}\omega_{j}^{2}\leq
1/\gamma$.

Condition 5 complements Condition 3, and is required to deal with the standard
errors (see Lemma 7 in Supplementary Materials).

######  Theorem 1 (Asymptotic normality of $t$-statistic).

If Conditions 1–5 are assumed, then both $t$-statistics $\texttt{T}_{ij}$
defined in (9) satisfy that $\texttt{T}_{ij}-\sqrt{T}\phi_{ij}/\hat{m}_{ij}$
converge in distribution to $\mathcal{N}(0,1)$ for all $(i,j)\in\mathcal{H}$
such that $\bar{v}_{i}:=\bar{r}_{i}+M_{\omega}^{2}\lambda=o(1)$.

In the theorem, $\bar{v}_{i}=o(1)$ entails that both $\bar{r}_{i}$ in
Proposition 2 and the estimation error of the standard error are
asymptotically negligible.

####  4.2.1 Theory for the FDR control

For
$\left((i,j),(k,\ell)\right)\in\mathcal{H}\times\mathcal{H}=:\mathcal{H}^{2}$,
define the correlation between $z_{ij}$ and $z_{k\ell}$ as

$\displaystyle\rho_{(i,j),(k,\ell)}=\frac{\sigma_{ik}\omega_{j\ell}}{\sigma_{i}\sigma_{k}\omega_{j}\omega_{\ell}},$
(15)

which becomes $\rho_{(i,j),(k,\ell)}\asymp\sigma_{ik}\omega_{j\ell}$ under
Condition 5. Obviously $\rho_{(i,j),(i,j)}=1$ for all $(i,j)\in\mathcal{H}$,
but the “off-diagonal” correlations in
$\mathcal{H}_{\text{off}}^{2}:=\left\\{((i,j),(k,\ell))\in\mathcal{H}^{2}:(i,j)\not=(k,\ell)\right\\}$,
where $|\mathcal{H}_{\text{off}}^{2}|=|\mathcal{H}^{2}|-|\mathcal{H}|$, take
non-trivial values. To manipulate them, impose the following condition.

######  Condition 6.

There exists a partition of $\mathcal{H}_{\text{off}}^{2}$ denoted as
$\mathcal{H}_{\text{off}}^{2}=\mathcal{H}_{w}^{2}\cup\mathcal{H}_{s}^{2}$ such
that for some constants $c>0$ and $\bar{\rho}\in(0,1)$,

$\displaystyle|\rho_{(i,j),(k,\ell)}|\in\begin{cases}\left[0,c/\log^{2}(KN^{2})\right]&\text{for
}~{}((i,j),(k,\ell))\in\mathcal{H}_{w}^{2}~{}~{}~{}\text{(weak
correlations)},\\\ \left(c/\log^{2}(KN^{2}),\bar{\rho}\right]&\text{for
}~{}((i,j),(k,\ell))\in\mathcal{H}_{s}^{2}~{}~{}~{}\text{(strong
correlations)},\end{cases}$

where
$|\mathcal{H}_{w}^{2}|=|\mathcal{H}_{\text{off}}^{2}|-|\mathcal{H}_{s}^{2}|$
and
$|\mathcal{H}_{s}^{2}|=O\left(|\mathcal{H}_{\text{off}}^{2}|/\log^{2}(N)\right)$.

Condition 6 says that most of the correlations are “weak” so that they are
close to zero, while some are allowed to be “strong” enough to be non-
vanishing. Because $\texttt{T}_{ij}$’s are asymptotically normal as shown in
Theorem 1, Condition 6 ensures that most of $\texttt{T}_{ij}$’s are
asymptotically independent. This property leads to the main theoretical
results. Let $\bar{v}=\max_{i\in[N]}\bar{v}_{i}$.

######  Theorem 2 (FDR control: Limiting normal distribution).

Suppose $\nu>4$ and $\bar{v}=O(T^{-\kappa_{1}})$ for some constant
$\kappa_{1}\in(0,1/2)$. If Conditions 1–6 are assumed, then for any
predetermined level $q\in[0,1]$, Procedure 1 with either $t$-statistics in (9)
achieves the following: If $\texttt{t}_{0}$ is given by (10), the obtained
$\hat{\mathcal{S}}(\texttt{t}_{0})$ satisfies

$\displaystyle\limsup_{N,T\to\infty}\operatorname{dFDR}\leq q\text{ and
}\lim_{N,T\to\infty}\operatorname{\mathbb{P}}\left(\operatorname{dFDP}\leq
q+\varepsilon\right)=1\text{ for any $\varepsilon>0$}.$

If $\texttt{t}_{0}$ is given by (11), the obtained
$\hat{\mathcal{S}}(\texttt{t}_{0})$ satisfies

$\displaystyle\limsup_{N,T\to\infty}\operatorname{dFWER}\leq q.$

We also show the directional FDR control by Procedure 2 with setting
$\tilde{\mathcal{S}}=\hat{\mathcal{S}}_{{\textsf{L}}}$. First, we need a
distributional assumption on the wild bootstrap.

######  Condition 7.

$\\{\zeta_{t}\\}$ is a sequence of i.i.d. sub-gaussian random variables with
$\operatorname{\mathbb{E}}\zeta_{t}=0$ and
$\operatorname{\mathbb{E}}\zeta_{t}^{2}=1$, where for every $T>0$, there
exists some constant $c_{\zeta}>0$ such that for all $x>0$,
$\max_{t\in[T]}\operatorname{\mathbb{P}}\left(|\zeta_{t}|>x\right)\leq
2\exp(-x^{2}/c_{\zeta})$.

Let

$\displaystyle\mathbb{Q}^{*}(\texttt{t})$
$\displaystyle=\frac{1}{KN^{2}-\hat{s}}\sum_{(i,j)\in\hat{\mathcal{S}}_{{\textsf{L}}}^{c}}\operatorname{\mathbb{P}}^{*}\left(\texttt{T}_{ij}^{*}>\texttt{t}\right),$

where
$\operatorname{\mathbb{P}}^{*}(\texttt{T}_{ij}^{*}>\texttt{t}):=p\lim_{B\to\infty}B^{-1}\sum_{b=1}^{B}1\\{\texttt{T}_{ij}^{*(b)}>\texttt{t}\\}$.
It suffices to show that for computing $\texttt{t}_{0}$ using $\mathbb{Q}^{*}$
is asymptotically the same as using $Q$; more precisely, we prove the event,
$|\mathbb{Q}^{*}(\texttt{t})/Q(\texttt{t})-1|$ and
$|\\{1-\mathbb{Q}^{*}(-\texttt{t})\\}/Q(\texttt{t})-1|$ converge to zero
uniformly in $\texttt{t}\in[0,\bar{\texttt{t}}]$, occurs with high
probability. Define
$\bar{\mu}=\bar{v}+s_{\omega}M_{\omega}^{3-2r}\lambda^{1-r}\log^{2}(N\vee T)$.

######  Theorem 3 (FDR control: Bootstrap).

Suppose $\nu>4$, $\bar{\mu}=O(T^{-\kappa_{1}})$ for some constant
$\kappa_{1}\in(0,1/2)$, and
$\mathcal{S}\subset\hat{\mathcal{S}}_{{\textsf{L}}}$ with high probability. If
Conditions 1–7 are assumed, then for any predetermined level $q\in[0,1]$,
Procedure 2 with either $t$-statistics in (12) achieves the same results as
Theorem 2.

The proof focuses on the $t$-statistic that has the “sandwich” s.e., and will
complete in two steps. First, show that “sandwich” $\texttt{T}_{ij}^{*}$ can
be approximated by the self-normalized sum. Then, verify that it can be
uniformly normally approximated in the relative error using the Cramér-type
large deviation theory of Jing et al. (2003). The additional cost
$s_{\omega}M_{\omega}^{3-2r}\lambda^{1-r}\log^{2}(N\vee T)$ in $\bar{\mu}$ is
due to the approximation error of the “sandwich” $t$-statistic to the self-
normalized sum,
$\sum_{t=1}^{T}\hat{u}_{it}^{*}\mathbf{x}_{t}^{\prime}\hat{\bm{\omega}}_{j}/\sum_{t=1}^{T}(\hat{u}_{it}^{*}\mathbf{x}_{t}^{\prime}\hat{\bm{\omega}}_{j})^{2}$.

####  4.2.2 Theory for the power guarantee

We next investigate the asymptotic power of our Procedures 1 and 2. A
condition on the signal strength is required to distinguish the nonzero
elements from zeros.

######  Condition 8.

For $\mathcal{S}\subset\mathcal{H}$, it holds that

$\displaystyle\min_{(i,j)\in\mathcal{S}}\frac{|\phi_{ij}|}{\sigma_{i}\omega_{j}}\geq
4\sqrt{\frac{2\log(KN^{2})}{T}}.$

######  Theorem 4 (Power guarantee).

Suppose $\nu>4$ and $\bar{v}=O(T^{-\kappa_{1}})$ for some constant
$\kappa_{1}\in(0,1/2)$. If Conditions 1–5 and 8 are assumed, then both
Procedures 1 and 2 achieve $\lim_{N,T\to\infty}\operatorname{dPower}=1$.

This theorem guarantees that both Procedures 1 and 2 do not asymptotically
miss any important relation between variables, whichever threshold of either
(10) or (11) is selected. The minimum signal of Condition 8 is strong enough
to achieve this.

##  5 Monte Carlo Experiments

We investigate the finite sample behavior of the proposed procedures by means
of Monte Carlo experiments.

###  5.1 First experiment: Stability

Before the comprehensive experiment, we carried out a small experiment to
illustrate the superiority of the selection result by our FDR-controlled
approach to those by the lasso and adaptive lasso selections. An $N\times 1$
vector ${\mathbf{y}}_{t}=$ for $t=-50,\dots,T$ is generated from the
stationary VAR(1) model for $N=20$ and $T=100$, with
$\mathbf{u}_{t}\sim\text{i.i.d.}N(\mathbf{0},\mathbf{I}_{N})$. The sparse
$\bm{\Phi}$ is generated as described in Section 5.2.1 in such a way that the
elements in the diagonal and $j$-diagonals for $j=-2,-1,1,2$, are non-zero.
Figure 3 shows heatmaps of the frequencies with which the ($i,j$)th element of
$\bm{\Phi}$ is selected as nonzero by different methods in 1000 replications.
Figure 3(a) shows the true sparsity pattern in $\bm{\Phi}$, and Figures 3(b),
(c), and (d) report the results using the lasso, adaptive lasso, and Procedure
1, respectively, with assigning the target FDR level $q=0.2$. The overall
frequency of false positives by the lasso is very high, and for certain
elements unreasonably so (shown by dark blue). The adaptive lasso reduces
overall false positives, but cannot do so sufficiently for the elements
wrongly selected in large numbers by lasso. In addition, some nonzero
$\phi_{ij}$’s are selected too infrequently by lasso, which is largely
inherited by adaptive lasso (shown by very weak red). Meanwhile, when the
multiple testing procedure is used, the frequency of false positives is evenly
spread among the elements and well controlled. Furthermore, nonzero elements
are selected more evenly than the (adaptive) lasso. These suggest that the
proposed inferential methods can provide more stable and reproducible results
for discovering the network Granger causality than existing popular estimation
methods.

(a) true Sparsity

(b) lasso

(c) adaptive lasso

(d) multiple test ($q=0.2$)

Figure 3: Heatmap of selection frequencies in $\bm{\Phi}_{1}$

###  5.2 Second experiment: Performance of our methods

####  5.2.1 Design

We consider the following Data Generating Process (DGP). The $N\times 1$
vector is generated as
$\mathbf{y}_{t}=\sum_{\ell=1}^{K}{\bm{\Phi}}_{\ell}\mathbf{y}_{t-\ell}+\mathbf{u}_{t}$
for $t=-(50+K),\dots,T$ with $\mathbf{y}_{-(50+K)}=\mathbf{0}$, and the set
$\\{\mathbf{y}_{-(50+K)},\dots,\mathbf{y}_{-K}\\}$ is discarded. Define
$\mathbf{u}_{t}={\bm{\Sigma}}_{u}^{1/2}{\bm{\varepsilon}}_{t}$, where
${\bm{\Sigma}}_{u}=\operatorname{diag}(\sigma_{1}^{2},\dots,\sigma_{N}^{2})$
with $\sigma_{i}^{2}\sim\text{i.i.d.}U(0.5,1.5)$. Another construction of
$\bm{\Sigma}_{u}$ with non-zero off-diagonals is considered in Supplementary
Materials. Two different distributions of $\varepsilon_{ti}$ are considered:

  1. (i)

Standard normal, $\varepsilon_{ti}\sim\text{i.i.d.}N\left(0,1\right)$.

  2. (ii)

Standardized mixture normal,
$\varepsilon_{ti}=(\eta_{ti}-\mu_{\eta})/\sigma_{\eta}$ with
$\eta_{ti}=q_{ti}\xi_{ti}+\left(1-q_{ti}\right)\zeta_{ti}$,
$\mu_{\eta}=\mathbb{E}\eta_{ti}$, and
$\sigma_{\eta}^{2}=\mathbb{E}\eta_{ti}^{2}-\mu_{\eta}^{2}$, where
$q_{ti}\sim\text{i.i.d.}Ber(\pi)$,
$\xi_{ti}\sim\text{i.i.d.}N(\mu_{\xi},\sigma_{\xi}^{2})$, and
$\zeta_{ti}\sim\text{i.i.d.}N(\mu_{\zeta},\sigma_{\zeta}^{2})$. We have chosen
the parameter values $\mu_{\xi}=2$, $\sigma_{\xi}=2$, $\mu_{\zeta}=4$,
$\sigma_{\zeta}=10$, and $\pi=0.9$, which give $\mu_{\eta}=0.4$ and
$\sigma_{\eta}=3.88$. The resulting error variable, $\varepsilon_{ti}$, is
unimodal yet it has skewness 1.86 and kurtosis 3.53.

We focus on the model with $K=1$. The coefficient matrix
${\bm{\Phi}}=(\phi_{ij})$ is constructed as follows. First generate
$\mathbf{\Psi}=(\psi_{ij})$ with $\psi_{ij}=\rho^{1+|i-j|/4}1\\{|i-j|\leq
m\\}$. Now form an $N\times N$ sign matrix
$\bm{\Upsilon}^{(r)}=(\upsilon_{ij}^{(r)})$ for $r=1,2,\dots$, where
$\upsilon_{ij}^{(r)}=2\varphi_{ij}^{(r)}-1$ with
$\varphi_{ij}^{(r)}\sim\text{i.i.d.}Ber(0.5)$. Repeatedly compute
${\bm{\Phi}}^{(r)}={\bm{\Psi}}\circ\bm{\Upsilon}^{(r)}$ for $r=1,2,\dots$
until $\lambda_{\max}({\bm{\Phi}}^{(r)})\leq 0.96$. If the process stops at
the $R$th repetition, set ${\bm{\Phi}}={\bm{\Phi}}^{(R)}$. We consider
$\rho=0.4$ and $m\in\\{2,4,7\\}$. Observe that, except for the first and last
rows, we have $s_{i}=5,9,15$ for $m=2,4,7$, respectively. Set $q=0.1$ and
consider all the combinations of $N\in\\{50,100,200,300\\}$ and
$T\in\\{200,300\\}$.

The model is estimated for each $i$th row using the R package, “glmnet.” In
the bootstrap procedure, the same values of the lasso tuning parameters at the
estimation stage are used. The precision matrix estimator is constructed by
the CLIME using the R package, “fastclime.” All the results are based on 1000
replications and 100 bootstrap samples.

####  5.2.2 Results

Table 1 summarizes the experimental results with
${\bm{\Sigma}}_{u}=\operatorname{diag}(\sigma_{11},\dots,\sigma_{NN})$, based
on $\texttt{T}_{ij}$ and its bootstrap version with
$\hat{m}_{ij}=\hat{\sigma}_{i}\sqrt{\hat{\bm{\omega}}^{\prime}_{j}\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}}$.
It contains three panels for $m=2,4,7$. The larger the value of $m$, the
larger the number of nonzero elements in ${\bm{\Phi}}$. The results with error
variance matrix with non-zero off-diagonals are qualitatively very similar, so
that they are reported in Table E2 in Supplementary Materials.

As can be seen, both the asymptotic and the bootstrap thresholds control the
dFDR to be around the predetermined level $q=0.1$, maintaining the high power.
In particular, the good performance for large $N\geq T$ is very encouraging.
The asymptotic threshold tends to produce slightly larger dFDR’s than $q$
while the bootstrap one tends to provide a conservative selection. The mixture
normal error moderately amplifies these tendencies. Despite this
conservativeness, the power based on the bootstrap threshold is almost
identical to that based on the asymptotic one. The exaggeration and the
conservativeness by both procedures are mitigated if $T$ is getting large. As
expected, when the nonzero elements in $\bm{\Phi}$ increase (i.e. $m$ rises),
the power goes down. However, it quickly rises as $T$ increases.

The results based on $\hat{m}_{ij}=\hat{\sigma}_{i}\hat{\omega}_{j}$ are
reported in Table E1. Bootstrap dFDR is virtually identical to bootstrap dFDR
based on
$\hat{m}_{ij}=\hat{\sigma}_{i}\sqrt{\hat{\bm{\omega}}^{\prime}_{j}\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}}$,
whilst asymptotic dFDR often becomes more conservative than bootstrap dFDR. As
a result, the asymptotic power is often lower than the bootstrap power. These
results suggest that the bootstrap procedure is expected to provide stable
performance regardless of the choice of $\hat{m}_{ij}$.

Table 1: Directional FDR and power using asymptotic and bootstrap thresholds
for $q=0.1$, with cross-sectionally uncorrelated errors $m=2$ ($\max s_{i}=5$)
---
|  |  | $T=200$ |  | $T=300$
|  |  | asymptotic |  | bootstrap |  | asymptotic |  | bootstrap
|  |  | dFDR | PWR |  | dFDR | PWR |  | dFDR | PWR |  | dFDR | PWR
Standard Normal Error
| $N=50$ |  | 9.3 | 97.5 |  | 6.7 | 96.8 |  | 8.8 | 99.8 |  | 7.0 | 99.7
| $N=100$ |  | 10.7 | 94.7 |  | 6.9 | 93.2 |  | 9.9 | 99.4 |  | 7.4 | 99.2
| $N=200$ |  | 10.7 | 91.9 |  | 7.7 | 90.5 |  | 10.3 | 99.0 |  | 8.4 | 98.8
| $N=300$ |  | 10.4 | 89.6 |  | 6.8 | 87.3 |  | 9.7 | 98.5 |  | 8.5 | 98.4
Mixture Normal Error
| $N=50$ |  | 10.4 | 94.2 |  | 5.7 | 91.9 |  | 9.5 | 99.0 |  | 6.3 | 98.7
| $N=100$ |  | 12.6 | 89.5 |  | 5.5 | 84.5 |  | 11.5 | 97.8 |  | 6.5 | 96.8
| $N=200$ |  | 13.1 | 87.0 |  | 6.6 | 82.5 |  | 11.7 | 97.2 |  | 7.6 | 96.4
| $N=300$ |  | 14.0 | 84.4 |  | 5.5 | 77.4 |  | 12.6 | 96.2 |  | 7.7 | 95.1
$m=4$ ($\max s_{i}=9$)
|  |  | $T=200$ |  | $T=300$
|  |  | asymptotic |  | bootstrap |  | asymptotic |  | bootstrap
|  |  | dFDR | PWR |  | dFDR | PWR |  | dFDR | PWR |  | dFDR | PWR
Standard Normal Error
| $N=50$ |  | 8.3 | 89.1 |  | 6.7 | 87.9 |  | 7.9 | 96.7 |  | 6.9 | 96.4
| $N=100$ |  | 8.5 | 83.7 |  | 8.5 | 83.6 |  | 8.3 | 94.7 |  | 8.5 | 94.8
| $N=200$ |  | 9.3 | 78.1 |  | 9.6 | 78.4 |  | 8.8 | 91.9 |  | 9.6 | 92.3
| $N=300$ |  | 10.3 | 73.0 |  | 8.8 | 71.8 |  | 9.4 | 89.3 |  | 9.9 | 89.5
Mixture Normal Error
| $N=50$ |  | 8.5 | 85.1 |  | 5.8 | 82.5 |  | 7.8 | 94.8 |  | 6.1 | 93.9
| $N=100$ |  | 8.9 | 79.9 |  | 7.2 | 78.2 |  | 8.7 | 92.4 |  | 7.6 | 91.9
| $N=200$ |  | 11.5 | 74.7 |  | 8.2 | 71.7 |  | 10.6 | 88.9 |  | 8.8 | 88.0
| $N=300$ |  | 12.0 | 69.9 |  | 6.9 | 65.3 |  | 11.0 | 86.1 |  | 8.9 | 85.1
$m=7$ ($\max s_{i}=15$)
|  |  | $T=200$ |  | $T=300$
|  |  | asymptotic |  | bootstrap |  | asymptotic |  | bootstrap
|  |  | dFDR | PWR |  | dFDR | PWR |  | dFDR | PWR |  | dFDR | PWR
Standard Normal Error
| $N=50$ |  | 6.8 | 68.6 |  | 6.0 | 67.5 |  | 6.7 | 80.4 |  | 6.1 | 79.8
| $N=100$ |  | 7.6 | 60.6 |  | 8.0 | 61.0 |  | 7.6 | 74.2 |  | 8.0 | 74.6
| $N=200$ |  | 9.5 | 53.0 |  | 9.0 | 52.5 |  | 8.6 | 67.8 |  | 9.1 | 68.2
| $N=300$ |  | 9.4 | 50.1 |  | 9.3 | 50.1 |  | 8.5 | 65.3 |  | 10.1 | 66.4
Mixture Normal Error
| $N=50$ |  | 6.4 | 65.5 |  | 4.9 | 63.2 |  | 6.1 | 77.9 |  | 5.1 | 76.6
| $N=100$ |  | 7.9 | 58.1 |  | 6.8 | 56.7 |  | 7.9 | 71.8 |  | 7.2 | 71.1
| $N=200$ |  | 9.4 | 51.1 |  | 7.5 | 49.3 |  | 9.0 | 65.6 |  | 8.3 | 65.0
| $N=300$ |  | 10.5 | 48.9 |  | 7.5 | 46.5 |  | 9.6 | 63.4 |  | 9.1 | 62.9

##  6 Two Empirical Applications

We apply our proposed methods to two large datasets to discover the underlying
Granger-causal networks. Section 6.1 investigates the large macroeconomic and
financial variables. Section 6.2 analyzes the regional house price growths in
the UK.222The results reported in this section are based on
$\hat{m}_{ij}=\hat{\sigma}_{i}\sqrt{\hat{\bm{\omega}}^{\prime}_{j}\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}}$.

###  6.1 Large macroeconomic variables

The FRED-MD macroeconomic and financial data file of May 2019 is obtained from
McCracken’s website, and the variables are transformed as instructed by
McCracken and Ng (2016). The data consists of a balanced panel of 128 monthly
series spanning the period from June 1999 to May 2019. All series are
standardized before the analysis. Following McCracken and Ng (2016), the
series are categorized into eight groups: G1, Output and Income; G2, Labour
Market; G3, Consumption, Orders and Inventories; G4, Housing; G5, Interest and
Exchange Rate; G6, Prices; G7, Money and Credit; G8, Stock Market. (The group
order is different from McCracken and Ng, 2016.) The variables are numbered
from 1 to 128, and the descriptions for all the variables are reported in
Table F1 in Supplementary Materials.

We estimate a VAR(1) model. The asymptotic and bootstrap thresholds for the
$t$-ratios with $q=0.05$ were 2.88 and 4.41, respectively. Here we summarize
the result with the bootstrap threshold in Figure 4 as a network Granger
causality diagram. The result with the asymptotic threshold is reported in
Supplementary Materials. The nodes represent the variables, and their colors
show the eight categories. The size of a node indicates the number of
variables it significantly predicts; the larger the node size, the more
variables it can predict. The arrows show the direction of the Granger
causality; if $y_{1}$ can significantly predicts $y_{2}$, $y_{1}\rightarrow
y_{2}$. The self-lag effects are excluded in the figure. The main dynamic
inter-relations are clustered within the groups, yet interesting interlinkages
between the variable groups are observed. In particular, Price variables are
clustered together, and eleven price variables are caused by variable 112,
Real M2 Money Supply. Price variable 94, Crude Oil, also Granger-causes seven
other price variables. Finally, Real Manufactures and Trade Industries Sales
(variable 49) Granger-cause three producer price indices (variables 90, 91,
92). These findings make a lot of sense from an economic point of view. In
addition, it is easy to identify the variables which cause many other
variables (many edges come out from the node) and those which are caused by
many other variables (the node surrounded by many pointing arrows); see the
Housing variable cluster, for example.

Figure 4: Network Granger-causality: 128 macroeconomic variables, bootstrap
$t_{0}$, $q=0.05$

###  6.2 UK regional house price growths

We obtained the monthly average house prices at the local authority district
level, published in November 2021 by HM Land Registry in the UK. Before
analysis, we seasonally adjusted the prices, then deflated them by the UK
consumer price index (CPI).333The CPI index (D7BT, not seasonally adjusted) is
obtained from Office for National Statistics, UK. The house prices and the CPI
are seasonally adjusted using the R package “seasonal.”  Denoting by $HP_{it}$
the seasonally adjusted real house price of the district $i$ at month $t$, the
monthly house price growth is computed as $\Delta
hp_{it}=\log(HP_{it}/HP_{it-1})$. In this analysis, we choose variables of 86
districts of Scotland, Wales, and the London area, spanning 209 months, from
February 2004 to June 2021. The variables are numbered from 1 to 86, and the
full list of the district names is reported in Table F2 in Supplementary
Materials. All the series are demeaned before the analysis.

We estimate a VAR(1) model. The asymptotic and bootstrap thresholds for the
$t$-ratios with $q=0.05$ were 2.73 and 4.22, respectively. Following the
previous subsection, we summarize the result with the bootstrap threshold in
Figure 5 as a network Granger causality diagram. The result with the
asymptotic threshold is reported in Supplementary Materials.

Figure 5: Network Granger causality: 86 UK regional house prices, bootstrap
$t_{0}$, $q=0.05$

The results show that house price growth causality networks are more or less
clustered in each of the three regions. However, there are interesting inter-
regional Granger-causal relationships. The London network Granger-causes
variables 21 and 24, which are the house price growth in the City of Edinburgh
and West Lothian. These two regions comprise a large part of Edinburgh, the
capital of Scotland. Furthermore, the London regional network is shown to
Granger-cause and also to be Granger-caused by Cardiff, the capital of Wales,
variable 78. These results suggest that house price fluctuations in London
interact dynamically with those in the capitals of other countries in the UK.

There is one unique inter-regional Granger Causality. Falkirk in Scotland
(Variable 6) has a significant concentration of arrows from London, Scotland,
and Wales. A little searching leads to the BBC News of 21 July 2020 reporting
that the Scottish and UK governments have pledged 90 million pounds in “Growth
Deal” funds to stimulate the economy around
Falkirk.444https://www.bbc.co.uk/news/uk-scotland-tayside-central-53471904 The
deal was signed off by the UK, Scottish governments, and Falkirk council on 21
December 2021.555https://www.bbc.co.uk/news/uk-scotland-scotland-
business-59734937 This large deal may have attracted investment to the Falkirk
property market and was statistically identified as a Granger Causality for
house price fluctuations in Falkirk.

##  7 Conclusion

This paper has proposed two multiple testing procedures that control the FDR
for discovering the network Granger causality in high-dimensional VAR models.
The first procedure is based on the limiting normal distribution of the
$t$-statistics constructed by the debiased lasso estimator, and the second
procedure is based on the bootstrap distributions of the $t$-statistics. Their
theoretical properties, including FDR control and power guarantee, have been
investigated. The finite sample evidence suggests that both procedures can
successfully control the FDR while maintaining high power. The proposed method
has been applied to a number of macroeconomic variables as well as UK regional
house prices, and interesting economically rational causal networks have been
uncovered.

The main focus of this article is statistical inference on network Granger
causality, which embodies the predictability relationship among a large number
of variables. Our FDR controlling procedures are fairly robust against
contemporaneous correlations among the variables of interest; see Condition 6.
As considered in Barigozzi and Brownlees (2019), extending the proposed method
to discover contemporaneous correlation networks in addition to network
Granger causality could be an interesting future research direction.

It has been pointed out in the literature that Granger causality found by
statistical methods can be spurious. For example, Stokes and Purdon (2017)
argue that if the lag order ($K$) is too small, the estimator will be biased,
and if $K$ is too large, the bias will be reduced, but the variance will
increase, both leading to spurious Granger causality findings. Our proposed
method is not expected to suffer much from this problem, so long as a
sufficiently (but not too) large value of $K$ is chosen. However,
investigating the ‘optimal’ choice of $K$ can be another interesting topic for
future research.

## Supplementary Material

Section A: Proof of Theorems, Section B: Proof of Propositions, Section C:
Lemmas and their proofs, Section D: Precision Matrix Estimation, Section E:
Additional Empirical Results, Section F: List of Variable Names in Empirical
Applications, Section G: Additional Results of Empirical Applications

## Acknowledgments

The authors thank …

## Funding

This work was supported by …

## References

  * Babii et al. (2021) Babii, A., E. Ghysels, and J. Striaukas (2021).  High-dimensional granger causality tests with an application to vix and news.  SSRN: https://ssrn.com/abstract=3615718.
  * Barigozzi and Brownlees (2019) Barigozzi, M. and C. Brownlees (2019).  NETS: Network estimation for time series.  Journal of Applied Econometrics 34, 347–364.
  * Basu et al. (2019) Basu, S., X. Li, and G. Michailidis (2019).  Low rank and structured modeling of high-dimensional vector autoregressions.  IEEE Transactions on Signal Processing 67, 1207–1222.
  * Basu and Michailidis (2015) Basu, S. and G. Michailidis (2015).  Regularized estimation in sparse high-dimensional time series models.  Annals of Statistics 43, 1535–1567.
  * Basu et al. (2015) Basu, S., A. Shojaie, and G. Michailidis (2015).  Network granger causality with inherent grouping structure.  Journal of Machine Learning Research 16, 417–453.
  * Benjamini and Hochberg (1995) Benjamini, Y. and Y. Hochberg (1995).  Controlling the false discovery rate: a practical and powerful approach to multiple testing.  Journal of the Royal Statistical Society Series B 57, 289–300.
  * Bercu et al. (2015) Bercu, B., B. Delyon, and E. Rio (2015).  Concentration Inequalities for Sums and Martingales (1st ed.).  Springer.
  * Bonferroni (1935) Bonferroni, C. E. (1935).  Il calcolo delle assicurazioni su gruppi di teste.  Studi in Onore del Professore Salvatore Ortu Carboni, 13–60.
  * Brüggemann et al. (2016) Brüggemann, R., C. Jentsch, and C. Trenkler (2016).  Inference in vars with conditional heteroskedasticity of unknown form.  Journal of Econometrics 191, 69–85.
  * Cai et al. (2011) Cai, T., W. Liu, and X. Luo (2011).  A constrained $\ell_{1}$ minimization approach to sparse precision matrix estimation.  Journal of the Amerian Statistical Association 106, 594–607.
  * Candès and Tao (2007) Candès, E. and T. Tao (2007).  The dantzig selector: statistical estimation when $p$ is much larger than $n$.  Annals of Statistics 35, 2313–2351.
  * Chi et al. (2021) Chi, C.-M., Y. Fan, C.-K. Ing, and J. Lv (2021).  High-dimensional knockoffs inference for time series data.  arXiv:2112.09851.
  * Davidson and Flachaire (2008) Davidson, R. and E. Flachaire (2008).  The wild bootstrap, tamed at last.  Journal of Econometrics 146, 162–169.
  * Davis et al. (2016) Davis, R. A., P. Zang, and T. Zheng (2016).  Sparse vector autoregressive modeling.  Journal of Computational and Graphical Statistics 25, 1077–1096.
  * Dufour and Renault (1998) Dufour, J.-M. and E. Renault (1998).  Short run and long run causality in time series: Theory.  Econometrica 66, 1099––1125.
  * Eichler (2007) Eichler, M. (2007).  Granger causality and path diagrams for multivariate time series.  Journal of Econometrics 137, 334–353.
  * Eichler (2012a) Eichler, M. (2012a).  Causal inference in time series analysis.  In C. Berzuini, P. Dawid, and L. Bernardinelli (Eds.), Causality: Statistical Perspectives and Applications, pp. 327–354. Wiley.
  * Eichler (2012b) Eichler, M. (2012b).  Graphical modelling of multivariate time series.  Probability Theory and Related Fields 153, 233–268.
  * Fujita et al. (2007) Fujita, A., J. R. Sato, H. M. Garay-Malpartida, R. Yamaguchi, S. Miyano, M. C. Sogayar, and C. E. Ferreira (2007).  Modeling gene expression regulatory networks with the sparse vector autoregressive model.  BMC Systems Biology 1, 1–11.
  * Gonçalves and Kilian (2004) Gonçalves, S. and L. Kilian (2004).  Bootstrapping autoregressions with conditional heteroskedasticity of unknown form.  Journal of Econometrics 123, 89–120.
  * Granger (1969) Granger, C. W. J. (1969).  Investigating causal relations by econometric models and cross-spectral methods.  Econometrica 37, 424–438.
  * Guðmundsson and Brownlees (2021) Guðmundsson, G. S. and C. Brownlees (2021).  Detecting groups in large vector autoregressions.  Journal of Econometrics forthcoming.
  * Hafner and Herwartz (2009) Hafner, C. M. and H. Herwartz (2009).  Testing for linear vector autoregressive dynamics under multivariate generalized autoregressive heteroskedasticity.  Statistica Neerlandica 63, 294–323.
  * Han et al. (2015) Han, F., H. Lu, and H. Liu (2015).  A direct estimation of high dimensional stationary vector autoregressions.  Journal of Machine Learning Research 16, 3115–3150.
  * Holm (1979) Holm, S. (1979).  A simple sequentially rejective multiple test procedure.  Scandinavian Journal of Statistics 6, 65–70.
  * Horn and Johnson (2012) Horn, R. A. and C. R. Johnson (2012).  Matrix Analysis (2nd ed.).  Cambridge University Press.
  * Hosoya (1977) Hosoya, Y. (1977).  On the Granger condition for non-causality.  Econometrica 45, 1735–1736.
  * Javanmard and Javadi (2019) Javanmard, A. and H. Javadi (2019).  False discovery rate control via debiased lasso.  Electronic Journal of Statistics 13, 1212–1253.
  * Javanmard and Montanari (2014) Javanmard, A. and A. Montanari (2014).  Confidence intervals and hypothesis testing for high-dimensional regression.  Journal of Machine Learning Research 15, 2869–2909.
  * Jing et al. (2003) Jing, B.-Y., Q.-M. Shao, and Q. Wang (2003).  Self-normalized cramér-type large deviations for independent random variables.  Annals of Probability 31, 2167–2215.
  * Kifer (2013) Kifer, Y. (2013).  Strong approximations for nonconventional sums and almost sure limit theorems.  Stochastic Processes and their Applications 123, 2286–2302.
  * Kilian (1999) Kilian, L. (1999).  Finite-sample properties of percentile and percentile-t bootstrap confidence intervals for impulse responses.  Review of Economics and Statistics 81, 652–660.
  * Kock and Callot (2015) Kock, A. B. and L. Callot (2015).  Oracle inequalities for high dimensional vector autoregressions.  Journal of Econometrics 186, 325–344.
  * Lin and Michailidis (2017) Lin, J. and G. Michailidis (2017).  Regularized estimation and testing for high-dimensional multi-block vector-autoregressive models.  Journal of Machine Learning Research 18, 1–49.
  * Liu (2013) Liu, W. (2013).  Gaussian graphical model estimation with false discovery rate control.  Annals of Statistics 41, 2948–2978.
  * Lozano et al. (2009) Lozano, A. C., N. Abe, Y. Liu, and S. Rosset (2009).  Grouped graphical granger modeling for gene expression regulatory networks discovery.  Bioinformatics 25, i110–i118.
  * Lütkepohl (2005) Lütkepohl, H. (2005).  New Introduction to Multiple Time Series Analysis.  Springer Berlin Heidelberg.
  * Mammen (1993) Mammen, E. (1993).  Bootstrap and wild bootstrap for high dimensional linear models.  Annals of Statistics 21, 255–285.
  * McCracken and Ng (2016) McCracken, M. W. and S. Ng (2016).  Fred-md: A monthly database for macroeconomic research.  Journal of Business & Economic Statistics 34(4), 574–589.
  * Reid et al. (2016) Reid, S., R. Tibshirani, and J. Friedman (2016).  A study of error variance estimation in lasso regression.  Statistica Sinica 26, 35–67.
  * Shojaie and Fox (2021) Shojaie, A. and E. B. Fox (2021).  Granger causality: A review and recent advances.  arXiv:2105.02675.
  * Shu and Nan (2019) Shu, H. and B. Nan (2019).  Estimation of large covariance and precision matrices from temporally dependent observations.  Annals of Statistics 47, 1321–1350.
  * Sims (1972) Sims, C. A. (1972).  Money, income, and causality.  American Economic Review 62, 540–552.
  * Stokes and Purdon (2017) Stokes, P. A. and P. L. Purdon (2017).  A study of problems encountered in granger causality analysis from a neuroscience perspective.  Proceedings of the National Academy of Sciences 114(34), E7063–E7072.
  * Szarek and Werner (1999) Szarek, S. J. and E. Werner (1999).  A nonsymmetric correlation inequality for gaussian measure.  Journal of Multivariate Analysis 68, 193–211.
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society Series B, 267–288.
  * Uematsu and Yamagata (2021) Uematsu, Y. and T. Yamagata (2021).  Inference in sparsity-induced weak factor models.  Journal of Business & Economic Statistics.
  * van de Geer et al. (2014) van de Geer, S., P. Bühlmann, Y. Ritov, and R. Dezeure (2014).  On asymptotically optimal confidence regions and tests for high-dimensional models.  Annals of Statistics 42, 1166–1202.
  * Vershynin (2018) Vershynin, R. (2018).  High-Dimensional Probability: An Introduction with Applications in Data Science.  Cambridge University Press.
  * Zhang and Zhang (2014) Zhang, C.-H. and S. S. Zhang (2014).  Confidence intervals for low dimensional parameters in high dimensional linear models.  Journal of Royal Statistical Society Series B 76, 217–242.
  * Zheng and Raskutti (2019) Zheng, L. and G. Raskutti (2019).  Testing for high-dimensional network parameters in auto-regressive models.  Electronic Journal of Statistics 13, 4977–5043.
  * Zhu and Liu (2020) Zhu, K. and H. Liu (2020).  Confidence intervals for parameters in high-dimensional sparse vector autoregression.  arXiv:2009.09462.

Supplementary Material for

Discovering the Network Granger Causality
in Large Vector Autoregressive Models

Yoshimasa Uematsu∗ and Takashi Yamagata†

*Center for the Promotion of Social Data Science Education and Research, Hitotsubashi University

*Department of Economics and Management, Tohoku University

$\dagger$Department of Economics and Related Studies, University of York

$\dagger$Institute of Social Economic Research, Osaka University

##  Appendix A Proofs of Theorems

###  A.1 Proof of Theorem 1

###### Proof.

Fix any $(i,j)\in\mathcal{H}$ such that $\bar{v}_{i}=o(1)$. Denote
$m_{ij}=\sqrt{\sigma_{i}^{2}\omega_{jj}}=\sqrt{\sigma_{i}^{2}{\bm{\omega}}_{j}^{\prime}\bm{\Sigma}_{x}{\bm{\omega}}_{j}}$
and
$\hat{m}_{ij}=\hat{\sigma}_{i}\sqrt{\hat{\bm{\omega}}_{j}^{\prime}\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}}$
or $\hat{\sigma}_{i}\hat{\omega}_{j}$. By the construction of the debiased
lasso estimator with Proposition 2 and Lemma 7, the $t$-statistic is written
as

$\displaystyle\texttt{T}_{ij}-\frac{\sqrt{T}\phi_{ij}}{\hat{m}_{ij}}$
$\displaystyle=\frac{\sqrt{T}(\hat{\phi}_{ij}-\phi_{ij})}{\hat{m}_{ij}}=\frac{z_{ij}+r_{ij}}{m_{ij}}\left\\{1+\left(\frac{m_{ij}}{\hat{m}_{ij}}-1\right)\right\\}$
$\displaystyle=\left\\{\frac{z_{ij}}{m_{ij}}+O\left(\bar{r}_{i}\right)\right\\}\left\\{1+O\left(M_{\omega}^{2}\lambda\right)\right\\},$

which holds with probability at least $1-O((N\vee T)^{-\nu})$. Here,
$\bar{r}_{i}$ and $M_{\omega}^{2}\lambda$ are asymptotically negligible (with
high probability) by the assumed condition.

Next, prove the asymptotic normality of $z_{ij}/m_{ij}$ for each
$(i,j)\in\mathcal{H}$. Recall that

$\displaystyle\frac{z_{ij}}{m_{ij}}=\sum_{t=1}^{T}\xi_{Tt}^{(i,j)},~{}~{}~{}~{}~{}\xi_{Tt}^{(i,j)}=\frac{u_{it}\mathbf{x}_{t}^{\prime}{\bm{\omega}}_{j}}{\sqrt{T}m_{ij}}.$

By a simple calculation, we have $\operatorname{Var}(\xi_{Tt}^{(i,j)})=1/T$
and $\sum_{t=1}^{T}(\xi_{Tt}^{(i,j)})^{2}\to_{p}1$. Furthermore, by an
application of Lemma 3, we obtain

$\displaystyle\max_{t}|\xi_{Tt}^{(i,j)}|\leq
T^{-1/2}\max_{t}\|u_{it}\mathbf{x}_{t}\|_{\infty}\|{\bm{\omega}}_{j}\|_{1}/m_{ij}\lesssim
M_{\omega}T^{-1/2}\sqrt{\log^{3}(N\vee T)}$

with high probability, which implies $\max_{t}|\xi_{Tt}^{(i,j)}|\to_{p}0$ when
$M_{\omega}^{2}\lambda=o(1)$. Thus McLeish’s central limit theorem is
applicable to achieve
$\texttt{T}_{ij}-\sqrt{T}\phi_{ij}/\hat{m}_{ij}\to_{d}N(0,1)$ for each
$(i,j)\in\mathcal{H}$ such that $\bar{r}_{i}+M_{\omega}^{2}\lambda=o(1)$. This
completes the proof. ∎

###  A.2 Proof of Theorem 2

###### Proof.

Throughout this proof, let $p=|\mathcal{H}|=KN^{2}$ denote the total number of
parameters in the (augmented) coefficient matrix, $\bm{\Phi}$. Define

$\displaystyle\mathcal{H}_{\leq 0}=\mathcal{H}\cap\\{(i,j):\phi_{ij}\leq
0\\},~{}~{}~{}\mathcal{H}_{\geq 0}=\mathcal{H}\cap\\{(i,j):\phi_{ij}\geq
0\\},$
$\displaystyle\mathcal{S}_{<0}=\mathcal{H}\cap\\{(i,j):\phi_{ij}<0\\},~{}~{}~{}\mathcal{S}_{>0}=\mathcal{H}\cap\\{(i,j):\phi_{ij}>0\\}.$

By the definition of
$s=|\mathcal{S}|=|\mathcal{H}\cap\\{(i,j):\phi_{ij}\not=0\\}|$, there is some
sequence $\pi=\pi_{np}\in[0,1]$ such that $|\mathcal{S}_{<0}|=\pi s$ and
$|\mathcal{S}_{>0}|=(1-\pi)s$. We also have $|\mathcal{H}_{\leq
0}|=p-(1-\pi)s$ and $|\mathcal{H}_{\geq 0}|=p-\pi s$.

Case 1. Consider when (10) does not exist and $\texttt{t}_{0}=\sqrt{2\log p}$.
First, we observe that

$\displaystyle\operatorname{dFDR}(\texttt{t}_{0})$
$\displaystyle\leq\operatorname{dFWER}=\operatorname{\mathbb{P}}\left(\sum_{(i,j)\in\hat{\mathcal{S}}(\texttt{t}_{0})}1\\{\operatorname{sgn}(\hat{\phi}_{ij})\not=\operatorname{sgn}(\phi_{ij})\\}\geq
1\right)$
$\displaystyle\leq\operatorname{\mathbb{P}}\left(\sum_{(i,j)\in\mathcal{H}_{\leq
0}}1\\{\texttt{T}_{ij}\geq\texttt{t}_{0}\\}\geq
1\right)+\operatorname{\mathbb{P}}\left(\sum_{(i,j)\in\mathcal{H}_{\geq
0}}1\\{\texttt{T}_{ij}\leq-\texttt{t}_{0}\\}\geq 1\right).$ (A.1)

The first probability of (A.1) is further bounded as

$\displaystyle\operatorname{\mathbb{P}}\left(\sum_{(i,j)\in\mathcal{H}_{\leq
0}}1\\{\texttt{T}_{ij}\geq\texttt{t}_{0}\\}\geq
1\right)\leq\sum_{(i,j)\in\mathcal{H}_{\leq
0}}\operatorname{\mathbb{P}}\left(\texttt{T}_{ij}\geq\texttt{t}_{0}\right)$
$\displaystyle\quad\leq\sum_{(i,j)\in\mathcal{S}^{c}}\operatorname{\mathbb{P}}\left(\mathcal{Z}_{ij}\geq\texttt{t}_{0}-\delta_{1}\right)+\sum_{(i,j)\in\mathcal{S}_{<0}}\operatorname{\mathbb{P}}\left(\mathcal{Z}_{ij}+\sqrt{T}\phi_{ij}/{\hat{m}_{ij}}\geq\texttt{t}_{0}-\delta_{1}\right)+pO((N\vee
T)^{-\nu+2})$
$\displaystyle\quad\leq\sum_{(i,j)\in\mathcal{S}^{c}}\operatorname{\mathbb{P}}\left(\mathcal{Z}_{ij}\geq\texttt{t}_{0}-\delta_{1}\right)+\sum_{(i,j)\in\mathcal{S}_{<0}}\operatorname{\mathbb{P}}\left(\mathcal{Z}_{ij}\geq\texttt{t}_{0}-\delta_{1}\right)+pO((N\vee
T)^{-\nu+2})$ $\displaystyle\quad=(p-s+\pi
s)Q(\texttt{t}_{0}-\delta_{1})+pO((N\vee T)^{-\nu+2})$ $\displaystyle\quad\leq
pQ(\texttt{t}_{0}-\delta_{1})+pO((N\vee T)^{-\nu+2}),$

where the second inequality holds for some positive sequence
$\delta_{1}=O(T^{-\kappa})$ for some $\kappa\in(0,\kappa_{1}]$ by Lemma 8(a)
and $Q(\texttt{t})=\operatorname{\mathbb{P}}(\mathcal{Z}>\texttt{t})$ is the
upper tail probability of a standard normal random variable $\mathcal{Z}$.
Because $Q(\texttt{t})\leq(\texttt{t}\sqrt{2\pi})^{-1}\exp(-\texttt{t}^{2}/2)$
and $\texttt{t}_{0}-\delta_{1}=\sqrt{2\log p}+o(1)$, we have

$\displaystyle pQ(\texttt{t}_{0}-\delta_{1})\lesssim p(\sqrt{\log
p})^{-1}\exp(-\log p+o(1))=O(1/\sqrt{\log p}).$

We also obtain $pO((N\vee T)^{-\nu+2})=o(1)$ for $\nu>4$. The second
probability of (A.1) is bounded in the same way. We thus conclude
$\operatorname{dFDR}(\texttt{t}_{0})\leq\operatorname{dFWER}=o(1)$.

Case 2. Consider when $\texttt{t}_{0}$ is given by (10). Define

$\displaystyle
V=\sup_{\texttt{t}\in[0,\bar{\texttt{t}}]}\left|\frac{\sum_{(i,j)\in\mathcal{H}_{\leq
0}}\left[1\\{\texttt{T}_{ij}\geq\texttt{t}\\}-Q(\texttt{t})\right]+\sum_{(i,j)\in\mathcal{H}_{\geq
0}}\left[1\\{\texttt{T}_{ij}\leq-\texttt{t}\\}-Q(\texttt{t})\right]}{2pQ(\texttt{t})}\right|.$

Then we have

$\displaystyle\operatorname{dFDP}(\texttt{t}_{0})$
$\displaystyle=\frac{|\\{(i,j)\in\hat{\mathcal{S}}(\texttt{t}_{0}):\operatorname{sgn}(\hat{\phi}_{ij})\not=\operatorname{sgn}(\phi_{ij})\\}|}{|\hat{\mathcal{S}}(\texttt{t}_{0})|\vee
1}$
$\displaystyle=\frac{2pQ(\texttt{t}_{0})}{|\hat{\mathcal{S}}(\texttt{t}_{0})|\vee
1}\frac{\sum_{(i,j)\in\mathcal{H}_{\leq
0}}1\\{\texttt{T}_{ij}\geq\texttt{t}_{0}\\}+\sum_{(i,j)\in\mathcal{H}_{\geq
0}}1\\{\texttt{T}_{ij}\leq-\texttt{t}_{0}\\}}{2pQ(\texttt{t}_{0})}$
$\displaystyle\leq q\left[\frac{\sum_{(i,j)\in\mathcal{H}_{\leq
0}}1\\{\texttt{T}_{ij}\geq\texttt{t}_{0}\\}+\sum_{(i,j)\in\mathcal{H}_{\geq
0}}1\\{\texttt{T}_{ij}\leq-\texttt{t}_{0}\\}}{2pQ(\texttt{t}_{0})}\right]$
$\displaystyle\leq q\left[V+\frac{(|\mathcal{H}_{\leq 0}|+|\mathcal{H}_{\geq
0}|)Q(\texttt{t}_{0})}{2pQ(\texttt{t}_{0})}\right]\leq q\left(V+1\right).$

Thus the dFDR and dFDP are controlled if we prove $V=o_{p}(1)$ in view of
Fatou’s lemma and Markov’s inequality, respectively. Note that

$\displaystyle
V\leq\sup_{\texttt{t}\in[0,\bar{\texttt{t}}]}\left|\frac{\sum_{(i,j)\in\mathcal{H}_{\leq
0}}\left[1\\{\texttt{T}_{ij}\geq\texttt{t}\\}-Q(\texttt{t})\right]}{2pQ(\texttt{t})}\right|+\sup_{\texttt{t}\in[0,\bar{\texttt{t}}]}\left|\frac{\sum_{(i,j)\in\mathcal{H}_{\geq
0}}\left[1\\{\texttt{T}_{ij}\leq-\texttt{t}\\}-Q(\texttt{t})\right]}{2pQ(\texttt{t})}\right|.$

We only prove that the first term is $o_{p}(1)$ by symmetry.

To this end, consider discretization. That is, we partition
$[0,\bar{\texttt{t}}]$ into small intervals,
$0=\texttt{t}_{0}<\texttt{t}_{1}<\dots<\texttt{t}_{h}=\bar{\texttt{t}}=(2\log
p-a\log\log p)^{1/2}$, such that $\texttt{t}_{m}-\texttt{t}_{m-1}=v_{p}$ for
$m\in\\{1,\dots,h-1\\}$ and $\texttt{t}_{h}-\texttt{t}_{h-1}\leq v_{p}$, where
$v_{p}=(\log p\log\log p)^{-1/2}$. Then a simple calculation gives $1/h\leq
v_{p}/\bar{\texttt{t}}=O(1/(\log p\sqrt{\log\log p}))$. Fix arbitrary
$m\in\\{1,\dots,h\\}$. For any
$\texttt{t}\in[\texttt{t}_{m-1},\texttt{t}_{m}]$, we have

$\displaystyle\frac{\sum_{(i,j)\in\mathcal{H}_{\leq
0}}1\\{\texttt{T}_{ij}\geq\texttt{t}\\}}{2pQ(\texttt{t})}\leq\frac{\sum_{(i,j)\in\mathcal{H}_{\leq
0}}1\\{\texttt{T}_{ij}\geq\texttt{t}_{m-1}\\}}{2pQ(\texttt{t}_{m-1})}\frac{Q(\texttt{t}_{m-1})}{Q(\texttt{t}_{m})}$

and

$\displaystyle\frac{\sum_{(i,j)\in\mathcal{H}_{\leq
0}}1\\{\texttt{T}_{ij}\geq\texttt{t}\\}}{2pQ(\texttt{t})}\geq\frac{\sum_{(i,j)\in\mathcal{H}_{\leq
0}}1\\{\texttt{T}_{ij}\geq\texttt{t}_{m}\\}}{2pQ(\texttt{t}_{m})}\frac{Q(\texttt{t}_{m})}{Q(\texttt{t}_{m-1})}.$

Lemma 7.2 of Javanmard and Javadi (2019) gives

$\displaystyle\frac{Q(\texttt{t}_{m-1})}{Q(\texttt{t}_{m})}\leq\frac{Q(\texttt{t}_{m}-v_{p})}{Q(\texttt{t}_{m})}=1+O(v_{p}+v_{p}\bar{\texttt{t}})=1+o(1)$

uniformly in $m\in\\{1,\dots,h\\}$. Thus the proof completes if the following
is true:

$\displaystyle\tilde{V}:=\max_{m\in\\{1,\dots,h\\}}\left|\frac{\sum_{(i,j)\in\mathcal{H}_{\leq
0}}\left[1\\{\texttt{T}_{ij}\geq\texttt{t}_{m}\\}-Q(\texttt{t}_{m})\right]}{2pQ(\texttt{t}_{m})}\right|=o_{p}(1).$

Prove $\tilde{V}=o_{p}(1)$. Fix arbitrary $\varepsilon>0$. Then we have

$\displaystyle\operatorname{\mathbb{P}}\left(\tilde{V}>\varepsilon\right)$
$\displaystyle\leq
h\max_{m\in\\{1,\dots,h\\}}\operatorname{\mathbb{P}}\left(\left|\frac{\sum_{(i,j)\in\mathcal{H}_{\leq
0}}\left[1\\{\texttt{T}_{ij}\geq\texttt{t}_{m}\\}-Q(\texttt{t}_{m})\right]}{2pQ(\texttt{t}_{m})}\right|>\varepsilon\right)$
$\displaystyle\leq
h\max_{m\in\\{1,\dots,h\\}}\operatorname{\mathbb{E}}\left|\frac{\sum_{(i,j)\in\mathcal{H}_{\leq
0}}\left[1\\{\texttt{T}_{ij}\geq\texttt{t}_{m}\\}-Q(\texttt{t}_{m})\right]}{2pQ(\texttt{t}_{m})}\right|^{2}/\varepsilon^{2},$

where $h\lesssim\log p\sqrt{\log\log p}$. Consider bounding the expectation
uniformly in $m\in\\{1,\dots,b\\}$. Denote by $\sum$ and $\sum\sum$ the
summations over $(i,j)\in\mathcal{H}_{\leq 0}$ and
$(i,j),(k,\ell)\in\mathcal{H}_{\leq 0}$, respectively. By a simple calculus
and Lemma 8(a)(b), we obtain

$\displaystyle\operatorname{\mathbb{E}}\left[\frac{\sum\sum\left[1\\{\texttt{T}_{ij}\geq\texttt{t}_{m}\\}-Q(\texttt{t}_{m})\right]\left[1\\{\texttt{T}_{k\ell}\geq\texttt{t}_{m}\\}-Q(\texttt{t}_{m})\right]}{4p^{2}Q(\texttt{t}_{m})^{2}}\right]$
$\displaystyle=\frac{\sum\sum\operatorname{\mathbb{P}}(\texttt{T}_{ij}\geq\texttt{t}_{m},\texttt{T}_{k\ell}\geq\texttt{t}_{m})}{4p^{2}Q(\texttt{t}_{m})^{2}}-\frac{|\mathcal{H}_{\leq
0}|\sum\operatorname{\mathbb{P}}(\texttt{T}_{ij}\geq\texttt{t}_{m})}{2p^{2}Q(\texttt{t}_{m})}+\frac{|\mathcal{H}_{\leq
0}|^{2}}{4p^{2}}$
$\displaystyle\leq\frac{\sum\sum\operatorname{\mathbb{P}}(\mathcal{Z}_{ij}\geq\texttt{t}_{m}-\delta_{1},\mathcal{Z}_{k\ell}\geq\texttt{t}_{m}-\delta_{1})}{4p^{2}Q(\texttt{t}_{m})^{2}}-\frac{(1+O(s/p))Q(\texttt{t}_{m}+\delta_{1})}{2Q(\texttt{t}_{m})}+\frac{1}{4}$
$\displaystyle=:(i)+(ii)+1/4,$

where $(\mathcal{Z}_{ij},\mathcal{Z}_{k\ell})$ is a bivariate standard normal
random vector and $\delta_{1}=O(T^{-\kappa})$ for some constant $\kappa>0$.
Thus we will conclude $\tilde{V}=o_{p}(1)$ if we show that $(i)\leq
1/4+o(1/h)$ and $(ii)\leq-1/2+o(1/h)$, where $1/h=O(1/(\log p\sqrt{\log\log
p}))$.

First consider $(ii)$. Expand $Q(\texttt{t}_{m}+\delta_{1})$ around
$\delta_{1}=0$ by the mean value theorem. Then there exists $\delta_{1}^{*}$
between $0$ and $\delta_{1}$ such that

$\displaystyle(ii)$
$\displaystyle=-\frac{Q(\texttt{t}_{m}+\delta_{1})}{2Q(\texttt{t}_{m})}(1+o(1))=-\frac{Q(\texttt{t}_{m})+Q^{\prime}(\texttt{t}_{m}+\delta_{1}^{*})\delta_{1}}{2Q(\texttt{t}_{m})}(1+o(1))$
$\displaystyle\leq-1/2-\frac{(2\pi)^{-1/2}\exp\left\\{-(\texttt{t}_{m}+\delta_{1}^{*})^{2}/2\right\\}\delta_{1}}{2(2\pi)^{-1/2}\texttt{t}_{m}^{-1}\exp\left\\{-\texttt{t}_{m}^{2}/2\right\\}}(1+o(1))$
$\displaystyle=-1/2-\texttt{t}_{m}\exp\left\\{-\delta_{1}^{*}(\texttt{t}_{m}+\delta_{1}^{*}/2)\right\\}\delta_{1}(1+o(1))/2$
$\displaystyle=-1/2+o(1/h),$ (A.2)

where the last equality holds uniformly in $m\in\\{1,\dots,h\\}$ because
$\delta_{1}$ is polynomially decaying while $\texttt{t}_{m}$ is the
logarithmic function for all $m\in\\{1,\dots,h\\}$.

Next consider $(i)$ by decomposing the summation into two parts,
$(i,j)=(k,\ell)$ and $(i,j)\not=(k,\ell)$. First we see the summation over
$(i,j)=(k,\ell)$, which has $p$ entries. We have

$\displaystyle\frac{\sum\sum\operatorname{\mathbb{P}}(\mathcal{Z}_{ij}\geq\texttt{t}_{m}-\delta_{1},\mathcal{Z}_{k\ell}\geq\texttt{t}_{m}-\delta_{1})1\\{(i,j)=(k,\ell)\\}}{4p^{2}Q(\texttt{t}_{m})^{2}}$
$\displaystyle=\frac{|\mathcal{H}_{\geq
0}|}{4p}\frac{1}{pQ(\texttt{t}_{m})}\frac{Q(\texttt{t}_{m}-\delta_{1})}{Q(\texttt{t}_{m})}=o(1/h).$
(A.3)

The last estimate is true because we have $|\mathcal{H}_{\geq 0}|/p=O(1)$ and
$Q(\texttt{t}_{m}-\delta_{1})/Q(\texttt{t}_{m})=1+o(1/h)$ by the same reason
as above, and by Szarek and Werner (1999),

$\displaystyle\frac{1}{pQ(\texttt{t}_{m})}$
$\displaystyle\leq\frac{\bar{\texttt{t}}+(\bar{\texttt{t}}^{2}+4)^{1/2}}{p2(2\pi)^{-1/2}\exp\left\\{-\bar{\texttt{t}}^{2}/2\right\\}}\lesssim\frac{\bar{\texttt{t}}}{p\exp\left\\{-\bar{\texttt{t}}^{2}/2\right\\}}$
$\displaystyle\lesssim\log^{1/2}p\cdot\exp\\{-\log\log^{a/2}p\\}=O(\log^{1/2-a/2}p)=o(1/h)$
(A.4)

for any $a>3$.

Finally we bound $(i)$ with summation over $(i,j)\not=(k,\ell)$, which
contains $p^{2}-p$ entries. We have

$\displaystyle\frac{\sum\sum\operatorname{\mathbb{P}}(\mathcal{Z}_{ij}\geq\texttt{t}_{m}-\delta_{1},\mathcal{Z}_{k\ell}\geq\texttt{t}_{m}-\delta_{1})}{4p^{2}Q(\texttt{t}_{m})^{2}}$
$\displaystyle\leq\frac{\sum\sum\operatorname{\mathbb{P}}(\mathcal{Z}_{ij}\geq\texttt{t}_{m}-\delta_{1},\mathcal{Z}_{k\ell}\geq\texttt{t}_{m}-\delta_{1})}{4p^{2}Q(\texttt{t}_{m}-\delta_{1})^{2}}\frac{Q(\texttt{t}_{m}-\delta_{1})^{2}}{Q(\texttt{t}_{m})^{2}},$
(A.5)

where $Q(\texttt{t}_{m}-\delta_{1})^{2}/Q(\texttt{t}_{m})^{2}=1+o(1/h)$ as
above. By Lemma 9 and the inequality $1/\sqrt{1-x^{2}}\leq
1+|x|/\sqrt{1-x^{2}}$ for any $x\in(-1,1)$, we obtain

$\displaystyle\frac{\sum\sum\operatorname{\mathbb{P}}(\mathcal{Z}_{ij}\geq\texttt{t}_{m}-\delta_{1},\mathcal{Z}_{k\ell}\geq\texttt{t}_{m}-\delta_{1})}{4p^{2}Q(\texttt{t}_{m}-\delta_{1})^{2}}\leq\frac{1}{4p^{2}}\sum\sum\frac{1}{\sqrt{1-\rho_{(i,j),(k,\ell)}^{2}}}$
$\displaystyle=1/4+o(1/h)+\frac{1}{4p^{2}}\sum\sum\frac{|\rho_{(i,j),(k,\ell)}|}{\sqrt{1-\rho_{(i,j),(k,\ell)}^{2}}}.$
(A.6)

We evaluate the sum using Condition 6. Denote
$\mathcal{H}_{w}^{2}=\mathcal{H}_{w1}\times\mathcal{H}_{w2}$ and
$\mathcal{H}_{s}^{2}=\mathcal{H}_{s1}\times\mathcal{H}_{s2}$. Then it is
decomposed as

$\displaystyle\sum_{(i,j)\in\mathcal{H}_{\leq
0}}\sum_{(k,\ell)\in\mathcal{H}_{\leq 0}}=\sum_{(i,j)\in\mathcal{H}_{\leq
0}\cap\mathcal{H}_{w1}}\sum_{(k,\ell)\in\mathcal{H}_{\leq
0}\cap\mathcal{H}_{w2}}+\sum_{(i,j)\in\mathcal{H}_{\leq
0}\cap\mathcal{H}_{s1}}\sum_{(k,\ell)\in\mathcal{H}_{\leq
0}\cap\mathcal{H}_{s2}}.$

They are bounded by

$\displaystyle\frac{1}{4p^{2}}\sum_{(i,j)\in\mathcal{H}_{\leq
0}\cap\mathcal{H}_{w1}}\sum_{(k,\ell)\in\mathcal{H}_{\leq
0}\cap\mathcal{H}_{w2}}\frac{|\rho_{(i,j),(k,\ell)}|}{\sqrt{1-\rho_{(i,j),(k,\ell)}^{2}}}$
$\displaystyle\leq\frac{1}{4p^{2}}\sum_{(i,j)\in\mathcal{H}_{w1}}\sum_{(k,\ell)\in\mathcal{H}_{w2}}\frac{|\rho_{(i,j),(k,\ell)}|}{\sqrt{1-\rho_{(i,j),(k,\ell)}^{2}}}$
$\displaystyle\leq\frac{p^{2}-p-|\mathcal{H}_{s1}\times\mathcal{H}_{s2}|}{4p^{2}}\frac{c/\log^{2}p}{\sqrt{1-c^{2}/\log^{4}p}}=O(1)O(1/\log^{2}p)=o(1/h)$
(A.7)

and

$\displaystyle\frac{1}{4p^{2}}\sum_{(i,j)\in\mathcal{H}_{\leq
0}\cap\mathcal{H}_{s1}}\sum_{(k,\ell)\in\mathcal{H}_{\leq
0}\cap\mathcal{H}_{s2}}\frac{|\rho_{(i,j),(k,\ell)}|}{\sqrt{1-\rho_{(i,j),(k,\ell)}^{2}}}$
$\displaystyle\leq\frac{1}{4p^{2}}\sum_{(i,j)\in\mathcal{H}_{s1}}\sum_{(k,\ell)\in\mathcal{H}_{s2}}\frac{|\rho_{(i,j),(k,\ell)}|}{\sqrt{1-\rho_{(i,j),(k,\ell)}^{2}}}$
$\displaystyle\leq\frac{|\mathcal{H}_{s1}\times\mathcal{H}_{s2}|}{4p^{2}}\frac{|\bar{\rho}|}{\sqrt{1-\bar{\rho}^{2}}}=O(1/\log^{2}p)O(1)=o(1/h).$
(A.8)

From (A.3) and (A.5)–(A.8), we have $(i)=1/4+o(1/h)$. This completes the
proof. ∎

###  A.3 Proof of Theorem 3

Before proceeding the proof, we recall the notation. In the proof, we mainly
use the $t$-statistics with the “sandwich” s.e., defined as

$\displaystyle\texttt{T}_{ij}^{*}=\frac{\sqrt{T}\hat{\phi}_{ij}^{*}}{\hat{\sigma}^{*}_{i}\sqrt{\hat{\bm{\omega}}^{\prime}_{j}\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}}}$

for $(i,j)\in\hat{\mathcal{S}}_{{\textsf{L}}}^{c}$, where
$\hat{\sigma}_{i}^{*2}=(T-\hat{s}_{i})^{-1}\sum_{t=1}^{T}\hat{u}_{it}^{*2}$.
The bootstrap debiased lasso estimator is also defined as

$\displaystyle\hat{\bm{\Phi}}^{*}$
$\displaystyle=\hat{\bm{\Phi}}^{{\textsf{L}}*}+\left(\mathbf{Y}^{*}-\hat{\bm{\Phi}}^{{\textsf{L}}*}\mathbf{X}\right)\mathbf{X}^{\prime}\hat{\bm{\Omega}}/T$
$\displaystyle=\hat{\bm{\Phi}}^{{\textsf{L}}}+\mathbf{U}^{*}\mathbf{X}^{\prime}\bm{\Omega}/T+\mathbf{U}^{*}\mathbf{X}^{\prime}(\hat{\bm{\Omega}}-\bm{\Omega})/T+\left(\hat{\bm{\Phi}}^{{\textsf{L}}}-\hat{\bm{\Phi}}^{{\textsf{L}}*}\right)(\hat{\bm{\Sigma}}_{x}\hat{\bm{\Omega}}-\mathbf{I}).$

By the definition, we have $\hat{\phi}_{ij}^{{\textsf{L}}}=0$ for all
$(i,j)\in\hat{\mathcal{S}}_{{\textsf{L}}}^{c}$, and hence

$\displaystyle\hat{\phi}_{ij}^{*}=\mathbf{u}_{i}^{*}\mathbf{X}^{\prime}\bm{\omega}_{j}/T+\mathbf{u}_{i}^{*}\mathbf{X}^{\prime}(\hat{\bm{\omega}}_{j}-\bm{\omega}_{j})/T+\left(\hat{\bm{\phi}}_{i}^{{\textsf{L}}}-\hat{\bm{\phi}}_{i}^{{\textsf{L}}*}\right)(\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}-\mathbf{e}_{j})$

for all $(i,j)\in\hat{\mathcal{S}}_{{\textsf{L}}}^{c}$. We have also defined

$\displaystyle\mathbb{Q}^{*}(\texttt{t})$
$\displaystyle=\frac{1}{p-\hat{s}}\sum_{(i,j)\in\hat{\mathcal{S}}_{{\textsf{L}}}^{c}}\operatorname{\mathbb{P}}^{*}\left(\texttt{T}_{ij}^{*}>\texttt{t}\right)~{}~{}\text{(conditional
on the observations)},$

where $p=KN^{2}$, and $\operatorname{\mathbb{P}}^{*}$ indicates the
probability measure induced by the bootstrap.

###### Proof.

It is sufficient to consider the case when $\texttt{t}_{0}$ is computed by
(13). From the beginning of Case 2 in the proof of Theorem 2 with the same
argument, the FDR is decomposed as

$\displaystyle\operatorname{dFDP}(\texttt{t}_{0})$
$\displaystyle=\frac{p\left\\{\mathbb{Q}^{*}(\texttt{t}_{0})+1-\mathbb{Q}^{*}(-\texttt{t}_{0})\right\\}}{|\hat{\mathcal{S}}(\texttt{t}_{0})|\vee
1}\cdot\frac{2Q(\texttt{t}_{0})}{\mathbb{Q}^{*}(\texttt{t}_{0})+1-\mathbb{Q}^{*}(-\texttt{t}_{0})}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\cdot\frac{\sum_{(i,j)\in\mathcal{I}_{\leq
0}}1\\{\texttt{T}_{ij}\geq\texttt{t}_{0}\\}+\sum_{(i,j)\in\mathcal{I}_{\geq
0}}1\\{\texttt{T}_{ij}\leq-\texttt{t}_{0}\\}}{2pQ(\texttt{t}_{0})}$
$\displaystyle\leq
q\cdot\frac{2}{\mathbb{Q}^{*}(\texttt{t}_{0})/Q(\texttt{t}_{0})+\\{1-\mathbb{Q}^{*}(-\texttt{t}_{0})\\}/Q(\texttt{t}_{0})}\left\\{V+1+o(1)\right\\},$

where $V$ is defined in the proof of Theorem 2 and has been shown
$V=o_{p}(1)$. Hence, it suffices to prove that the event,
$|\mathbb{Q}^{*}(\texttt{t}_{0})/Q(\texttt{t}_{0})-1|$ and
$|\\{1-\mathbb{Q}^{*}(-\texttt{t}_{0})\\}/Q(\texttt{t}_{0})-1|$ converge to
zero uniformly in $\texttt{t}_{0}\in[0,\bar{\texttt{t}}]$, occurs with high
probability.

To this end, define

$\displaystyle S_{ij}^{*}$
$\displaystyle=T^{-1/2}\sum_{t=1}^{T}{u}_{it}^{*}\mathbf{x}_{t}^{\prime}\hat{\bm{\omega}}_{j},~{}~{}R_{ij}^{*}=T^{-1/2}\mathbf{u}_{i}^{*}\mathbf{X}^{\prime}(\hat{\bm{\omega}}_{j}-{\bm{\omega}}_{j})+\sqrt{T}(\hat{\bm{\phi}}_{i}^{L}-\hat{\bm{\phi}}_{i}^{L*})(\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}-\mathbf{e}_{j}),$
$\displaystyle\hat{m}_{ij}^{*2}$
$\displaystyle=\hat{\sigma}_{i}^{*2}\hat{\bm{\omega}}_{j}^{\prime}\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j},~{}~{}\tilde{m}_{ij}^{*2}=T^{-1}\sum_{t=1}^{T}{u}_{it}^{*2}\hat{\bm{\omega}}_{j}^{\prime}\mathbf{x}_{t}\mathbf{x}_{t}^{\prime}\hat{\bm{\omega}}_{j},$

where
$\hat{\sigma}_{i}^{*2}=(T-\hat{s}_{i})^{-1}\sum_{t=1}^{T}\hat{u}_{it}^{*2}$
with
$\hat{u}_{it}^{*2}={u}_{it}^{*2}-2{u}_{it}^{*}\mathbf{x}_{t}^{\prime}\bm{\delta}_{i}^{*}+{\bm{\delta}_{i}^{*}}^{\prime}\mathbf{x}_{t}\mathbf{x}_{t}^{\prime}\bm{\delta}_{i}^{*}$
and ${u}_{it}^{*2}=\hat{u}_{it}^{2}\zeta_{t}^{2}$. Then by the construction,
we have $\texttt{T}_{ij}^{*}=(S_{ij}^{*}+R_{ij}^{*})/\hat{m}_{ij}^{*}$. We
first check if $\texttt{T}_{ij}^{*}$ is asymptotically equivalent to the self-
normalized sum, $S_{ij}^{*}/\tilde{m}_{ij}^{*}$, and then verify that it can
be uniformly normally approximated in the relative error.

Define the events,

$\displaystyle
A_{1}^{*}=\left\\{\max_{i\in[N]}\max_{j\in[KN]}\left|\frac{\tilde{m}_{ij}^{*}}{\hat{m}_{ij}^{*}}-1\right|\lesssim\bar{\mu}_{1}\right\\},~{}~{}~{}A_{2}^{*}=\left\\{\max_{i\in[N]}\max_{j\in[KN]}\left|\frac{R_{ij}^{*}}{\hat{m}_{ij}^{*}}\right|\lesssim\bar{\mu}_{2}\right\\},$

where

$\displaystyle\bar{\mu}_{1}$
$\displaystyle=\max\left\\{\bar{\tau}_{1},\bar{\tau}_{2},\bar{\tau}_{3}\right\\},$
$\displaystyle\bar{\mu}_{2}$
$\displaystyle=\max\left\\{\lambda^{2-r}M_{\omega}^{2-2r}s_{\omega},\
\bar{s}\lambda M_{\omega}\log^{3/2}(N\vee T)\right\\}$

with $\bar{\tau}_{1}$, $\bar{\tau}_{2}$, and $\bar{\tau}_{3}$ defined in
Lemmas 13, 14, and 15, respectively. Then the event,

$\displaystyle
A:=\left\\{\hat{\mathcal{S}}_{{\textsf{L}}}\supset\mathcal{S}\right\\}\cap\left\\{\operatorname{\mathbb{P}}^{*}\left(A_{1}^{*c}\right)=O((NT)^{-\nu})\right\\}\cap\left\\{\operatorname{\mathbb{P}}^{*}\left(A_{2}^{*c}\right)=O((NT)^{-\nu})\right\\},$

occurs with high probability by the assumed condition on
$\hat{\mathcal{S}}_{L}$ and Lemmas 16 and 17. By a simple computation, we have
$\bar{\mu}_{1}\vee\bar{\mu}_{2}\asymp\bar{\mu}$ with

$\displaystyle\bar{\mu}$
$\displaystyle=\max\left\\{M_{\omega}^{3-2r}\lambda^{1-r}s_{\omega}\log^{2}(N\vee
T),~{}M_{\omega}\bar{s}\lambda\log^{3/2}(N\vee
T),~{}M_{\omega}^{2}\lambda\log(N\vee T)\right\\}.$

Because $\bar{\mu}$ tends to zero polynomially, we observe that conditional on
$A$,

$\displaystyle\mathbb{Q}^{*}(\texttt{t})=\frac{1}{|\hat{\mathcal{S}}_{{\textsf{L}}}^{c}|}\sum_{(i,j)\in\hat{\mathcal{S}}_{{\textsf{L}}}^{c}}\operatorname{\mathbb{P}}^{*}\left(\frac{S_{ij}^{*}+R_{ij}^{*}}{\hat{m}_{ij}^{*}}>\texttt{t}\right)$
$\displaystyle\leq\frac{1}{|\hat{\mathcal{S}}_{{\textsf{L}}}^{c}|}\sum_{(i,j)\in\hat{\mathcal{S}}_{{\textsf{L}}}^{c}}\operatorname{\mathbb{P}}^{*}\left(\frac{S_{ij}^{*}}{\tilde{m}_{ij}^{*}}\left\\{1+\left(\frac{\tilde{m}_{ij}^{*}}{\hat{m}_{ij}^{*}}-1\right)\right\\}>\texttt{t}-\frac{R_{ij}^{*}}{\hat{m}_{ij}^{*}},~{}A_{1}^{*}\cap
A_{2}^{*}\right)+\operatorname{\mathbb{P}}^{*}\left(A_{1}^{*c}\cup
A_{2}^{*c}\right)$
$\displaystyle\leq\frac{1}{|\hat{\mathcal{S}}_{{\textsf{L}}}^{c}|}\sum_{(i,j)\in\hat{\mathcal{S}}_{{\textsf{L}}}^{c}}\operatorname{\mathbb{P}}^{*}\left(\frac{S_{ij}^{*}}{\tilde{m}_{ij}^{*}}>\frac{\texttt{t}-O(\bar{\mu}_{2})}{1+O(\bar{\mu}_{1})}\right)+O((N\vee
T)^{-\nu})$ (A.9)

with high probability, where the terms, $O(\bar{\mu}_{1})$,
$O(\bar{\mu}_{2})$, and $O((N\vee T)^{-\nu})$, depend neither on t nor
$(i,j)$. Since $(i,j)\in\hat{\mathcal{S}}_{{\textsf{L}}}^{c}$ implies
$(i,j)\in\mathcal{S}^{c}$ on event $A$, Lemma 18 entails the normal
approximation of the self-normalized sum,

$\displaystyle\operatorname{\mathbb{P}}^{*}\left(\frac{S_{ij}^{*}}{\tilde{m}_{ij}^{*}}>\frac{\texttt{t}-O(\bar{\mu}_{2})}{1+O(\bar{\mu}_{1})}\right)$
$\displaystyle=Q\left(\frac{\texttt{t}-O(\bar{\mu}_{2})}{1+O(\bar{\mu}_{1})}\right)\left\\{1+O\left(\frac{M_{\omega}\log
NT}{T^{1/2}}\right)\left(1+\frac{\texttt{t}-O(\bar{\mu}_{2})}{1+O(\bar{\mu}_{1})}\right)^{3}\right\\},$
(A.10)

with high probability. Combining (A.9) and (A.10) and dividing the both sides
by $Q(\texttt{t})$ yield

$\displaystyle\frac{\mathbb{Q}^{*}(\texttt{t})}{Q(\texttt{t})}$
$\displaystyle=\frac{Q\left(\left\\{\texttt{t}-O(\bar{\mu}_{2})\right\\}/\left\\{1+O(\bar{\mu}_{1})\right\\}\right)}{Q(\texttt{t})}\left\\{1+O\left(\frac{M_{\omega}\log
NT}{T^{1/2}}\right)\left(1+\frac{\texttt{t}-O(\bar{\mu}_{2})}{1+O(\bar{\mu}_{1})}\right)^{3}\right\\}$
$\displaystyle+\frac{O((N\vee T)^{-\nu})}{Q(\texttt{t})}.$

By the assumed condition, $\bar{\mu}_{1}$, $\bar{\mu}_{2}$ and
$M_{\omega}/T^{1/2}$ decay polynomially while
$\texttt{t}\in[0,\bar{\texttt{t}}]$ at most diverges logarithmically. Thus we
obtain uniformly in $\texttt{t}\in[0,\bar{\texttt{t}}]$,

$\displaystyle 1+O\left(\frac{M_{\omega}\log
NT}{T^{1/2}}\right)\left(1+\frac{\texttt{t}-O(\bar{\mu}_{2})}{1+O(\bar{\mu}_{1})}\right)^{3}=1+o(1)$

and, by an application of Lemma 7.2 of Javanmard and Javadi (2019),

$\displaystyle\frac{Q\left(\\{\texttt{t}-O(\bar{\mu}_{2})\\}/\\{1+O(\bar{\mu}_{1})\\}\right)}{Q(\texttt{t})}=1+O(\bar{\mu}_{2})(1+\texttt{t})+O(\bar{\mu}_{1})(1+\texttt{t}^{2})=1+o(1)$

uniformly in $\texttt{t}\in[0,\bar{\texttt{t}}]$. By (A.4), we have

$\displaystyle\sup_{\texttt{t}\in[0,\bar{\texttt{t}}]}\frac{O((N\vee
T)^{-\nu})}{Q(\texttt{t})}=\frac{O((N\vee
T)^{-\nu})}{Q(\bar{\texttt{t}})}=O((N\vee T)^{-\nu+2}).$

Consequently, it holds that

$\displaystyle\sup_{\texttt{t}\in[0,\bar{\texttt{t}}]}\left|\frac{\mathbb{Q}^{*}(\texttt{t})}{Q(\texttt{t})}-1\right|=o(1)$

with high probability for any $\nu>2$. The same argument is applied to showing
$|\\{1-\mathbb{Q}^{*}(-\texttt{t})\\}/Q(\texttt{t})-1|=o(1)$ uniformly in
$\texttt{t}\in[0,\bar{\texttt{t}}]$. This completes the proof. ∎

###  A.4 Proof of Theorem 4

###### Proof.

We use the same notation as in the proof of Theorem 2. Let
$\bar{\texttt{t}}_{0}=\sqrt{2\log p}$ denote the upper bound of the critical
value, $\texttt{t}_{0}$. By the definition of power and monotonicity of
probability, we have

$\displaystyle\operatorname{Power}(\texttt{t}_{0})$
$\displaystyle=\operatorname{\mathbb{E}}\left[\frac{|\\{(i,j)\in\hat{\mathcal{S}}(\texttt{t}_{0}):\operatorname{sgn}(\hat{\phi}_{ij})=\operatorname{sgn}(\phi_{ij})\\}|}{s\vee
1}\right]$
$\displaystyle=\frac{1}{s}\sum_{(i,j)\in\mathcal{S}_{<0}}\operatorname{\mathbb{P}}\left(\texttt{T}_{ij}\leq-\texttt{t}_{0}\right)+\frac{1}{s}\sum_{(i,j)\in\mathcal{S}_{>0}}\operatorname{\mathbb{P}}\left(\texttt{T}_{ij}\geq\texttt{t}_{0}\right)$
$\displaystyle\geq\frac{1}{s}\sum_{(i,j)\in\mathcal{S}_{<0}}\operatorname{\mathbb{P}}\left(\texttt{T}_{ij}\leq-\bar{\texttt{t}}_{0}\right)+\frac{1}{s}\sum_{(i,j)\in\mathcal{S}_{>0}}\operatorname{\mathbb{P}}\left(\texttt{T}_{ij}\geq\bar{\texttt{t}}_{0}\right).$

Consider the second probability in the lower bound. Proposition 2 gives

$\displaystyle\texttt{T}_{ij}=\sqrt{T}\hat{\phi}_{ij}/\hat{m}_{ij}=(\sqrt{T}\phi_{ij}+z_{ij}+r_{ij})/\hat{m}_{ij},$

where $\phi_{ij}=0$ if and only if $(i,j)\in\mathcal{S}^{c}$. It holds that

$\displaystyle\max_{(i,j)\in\mathcal{S}_{>0}}\operatorname{\mathbb{P}}\left(\texttt{T}_{ij}\leq\bar{\texttt{t}}_{0}\right)$
$\displaystyle\leq\max_{(i,j)\in\mathcal{S}_{>0}}\operatorname{\mathbb{P}}\left(\frac{z_{ij}+r_{ij}}{\hat{m}_{ij}}+\frac{m_{ij}}{\hat{m}_{ij}}\frac{\sqrt{T}\phi_{ij}^{0}}{m_{ij}}\leq\bar{\texttt{t}}_{0},~{}\frac{m_{ij}}{\hat{m}_{ij}}>\frac{1}{2}\right)+\max_{(i,j)\in\mathcal{S}_{>0}}\operatorname{\mathbb{P}}\left(\frac{m_{ij}}{\hat{m}_{ij}}\leq\frac{1}{2}\right)$
$\displaystyle\leq\max_{(i,j)\in\mathcal{S}_{>0}}\operatorname{\mathbb{P}}\left(\mathcal{Z}_{ij}\leq\bar{\texttt{t}}_{0}-\min_{(i,j)\in\mathcal{S}_{>0}}\frac{\sqrt{T}\phi_{ij}^{0}}{2m_{ij}}+\delta_{1}\right)+O\left((N\vee
T)^{-\nu}\right)$ $\displaystyle\leq\Phi\left(-\sqrt{2\log
p}+\delta_{1}\right)+O\left((N\vee T)^{-\nu}\right),$

where the second inequality follows from Lemmas 7 and 8, and the third
inequality holds by $\bar{\texttt{t}}_{0}=\sqrt{2\log p}$ and Condition 8.
Since $\delta_{1}\sqrt{\log p}=o(1)$ by the assumption, the Gaussian
probability in the upper bound tends to zero. The same result is obtained for
the first probability. Thus the power goes to unity because
$|\mathcal{S}_{<0}|+|\mathcal{S}_{>0}|=s$. This completes the proof. ∎

##  Appendix B Proofs of Propositions

###  B.1 Proof of Proposition 1

###### Proof.

Derive the non-asymptotic error bound for the Lasso estimator. First define
two events:

$\displaystyle\mathcal{E}_{1}=\left\\{\left\|T^{-1}\mathbf{X}\mathbf{U}^{\prime}\right\|_{\max}\leq\lambda/2\right\\},~{}~{}~{}\mathcal{E}_{2}=\left\\{\left\|\hat{\bm{\Sigma}}_{x}-\bm{\Sigma}_{x}\right\|_{\max}\leq
b\lambda/2\right\\},$

where $\lambda=8bc_{uu}\sqrt{2(\nu+5)^{3}\log^{3}(N\vee T)/T}$ with any
positive constant $\nu$. We work on event $\mathcal{E}_{1}\cap\mathcal{E}_{2}$
since Lemmas 3 and 4 guarantee that $\mathcal{E}_{1}\cap\mathcal{E}_{2}$
occurs with probability at least $1-O((N\vee T)^{-\nu})$.

Define
$\bm{\delta}_{i\cdot}=\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}-\bm{\phi}_{i\cdot}$.
Because $\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}$ minimizes the objective
function, we have

$\displaystyle(2T)^{-1}\|\mathbf{y}_{i\cdot}-\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}\mathbf{X}\|_{2}^{2}+\lambda\|\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}\|_{1}$
$\displaystyle\leq(2T)^{-1}\|\mathbf{y}_{i\cdot}-\bm{\phi}_{i\cdot}\mathbf{X}\|_{2}^{2}+\lambda\|\bm{\phi}_{i\cdot}\|_{1},$

which is equivalently written as

$\displaystyle(2T)^{-1}\|\bm{\delta}_{i\cdot}\mathbf{X}\|_{2}^{2}$
$\displaystyle\leq
T^{-1}\mathbf{u}_{i\cdot}\mathbf{X}^{\prime}\bm{\delta}_{i\cdot}^{\prime}+\lambda\|\bm{\phi}_{i\cdot}\|_{1}-\lambda\|\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}\|_{1}.$

By Hölder’s inequality and the triangle inequality, we have

$\displaystyle(2T)^{-1}\|\bm{\delta}_{i\cdot}\mathbf{X}\|_{2}^{2}$
$\displaystyle\leq\left|T^{-1}\mathbf{u}_{i\cdot}\mathbf{X}^{\prime}\bm{\delta}_{i\cdot}^{\prime}\right|+\lambda\|\bm{\phi}_{i\cdot}\|_{1}-\lambda\|\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}\|_{1}$
$\displaystyle\leq\|T^{-1}\mathbf{u}_{i\cdot}\mathbf{X}^{\prime}\|_{\infty}\|\bm{\delta}_{i\cdot}\|_{1}+\lambda\|\bm{\delta}_{i\cdot}\|_{1}.$

On event $\mathcal{E}_{1}$, we thus obtain the upper bound of
$\|\bm{\delta}_{i\cdot}\mathbf{X}\|_{2}^{2}$ as

$\displaystyle T^{-1}\|\bm{\delta}_{i\cdot}\mathbf{X}\|_{2}^{2}\leq
3\lambda\|\bm{\delta}_{i\cdot}\|_{1}.$ (A.11)

Next we bound $\|\bm{\delta}_{i\cdot}\mathbf{X}\|_{2}^{2}$ from below. Lemma 5
states that $\bm{\delta}_{i\cdot}$ lies in
$\mathcal{D}=\\{\bm{\delta}_{i\cdot}\in\mathbb{R}^{1\times
KN}:\|\bm{\delta}_{i\mathcal{S}_{i}^{c}}\|_{1}\leq
3\|\bm{\delta}_{i\mathcal{S}_{i}}\|_{1}\\}$ on $\mathcal{E}_{1}$. Hence under
$\mathcal{E}_{1}$ and $\mathcal{E}_{2}$, Lemma 6 entails

$\displaystyle
T^{-1}\|\bm{\delta}_{i\cdot}\mathbf{X}\|_{2}^{2}/\|\bm{\delta}_{i\cdot}\|_{2}^{2}\geq\gamma-16s_{i}\|\hat{\bm{\Sigma}}_{x}-\bm{\Sigma}_{x}\|_{\max}\geq\gamma-8bs_{i}\lambda.$
(A.12)

By (A.11), (A.12), Lemma 5, and the Cauchy–Schwarz inequality, we have

$\displaystyle(\gamma-8bs_{i}\lambda)\|\bm{\delta}_{i\cdot}\|_{2}^{2}$
$\displaystyle\leq
3\lambda\|\bm{\delta}_{i\cdot}\|_{1}=3\lambda\|\bm{\delta}_{i\mathcal{S}_{i}}\|_{1}+3\lambda\|\bm{\delta}_{i\mathcal{S}_{i}^{c}}\|_{1}$
$\displaystyle\leq 12\lambda\|\bm{\delta}_{i\mathcal{S}_{i}}\|_{1}\leq
12s_{i}^{1/2}\lambda\|\bm{\delta}_{i\mathcal{S}_{i}}\|_{2}\leq
12s_{i}^{1/2}\lambda\|\bm{\delta}_{i\cdot}\|_{2}.$

This concludes that

$\displaystyle\|\bm{\delta}_{i\cdot}\|_{2}\leq\frac{12s_{i}^{1/2}\lambda}{\gamma-8bs_{i}\lambda}.$

Next we derive the prediction error bound. By Lemma 5 and the Cauchy–Schwarz
inequality again, the error bound in the element-wise $\ell_{1}$-norm is given
by

$\displaystyle\|\bm{\delta}_{i\cdot}\|_{1}=\|\bm{\delta}_{i\mathcal{S}_{i}}\|_{1}+\|\bm{\delta}_{i\mathcal{S}_{i}^{c}}\|_{1}\leq
4\|\bm{\delta}_{i\mathcal{S}_{i}}\|_{1}\leq
4s_{i}^{1/2}\|\bm{\delta}_{i\mathcal{S}_{i}}\|_{2}\leq
4s_{i}^{1/2}\|\bm{\delta}_{i\cdot}\|_{2}\leq\frac{48s_{i}\lambda}{\gamma-8bs_{i}\lambda}.$
(A.13)

From (A.11) and (A.13), the prediction error bound is obtained as

$\displaystyle T^{-1}\|\bm{\delta}_{i\cdot}\mathbf{X}\|_{2}^{2}\leq
3\lambda\|\bm{\delta}_{i\cdot}\|_{1}\leq\frac{144s_{i}\lambda^{2}}{\gamma-8bs_{i}\lambda}.$

This completes the proof. ∎

###  B.2 Proof of Proposition 2

###### Proof.

The first assertion is trivial by the construction of the debiased lasso
estimator. We derive the upper bound of $|r_{ij}|$. Observe that for each
$i\in[N]$,

$\displaystyle|r_{ij}|$
$\displaystyle\leq\left\|T^{-1/2}\mathbf{u}_{i\cdot}\mathbf{X}^{\prime}\right\|_{\infty}\left\|\hat{\bm{\omega}}_{j}-{\bm{\omega}}_{j}\right\|_{1}+\sqrt{T}\left\|\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}-\bm{\phi}_{i\cdot}\right\|_{1}\left\|\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}-\mathbf{e}_{j}\right\|_{\infty}$
$\displaystyle\leq\left\|T^{-1/2}\mathbf{U}\mathbf{X}^{\prime}\right\|_{\max}\max_{j}\left\|\hat{\bm{\omega}}_{j}-{\bm{\omega}}_{j}\right\|_{1}+\sqrt{T}\left\|\hat{\bm{\phi}}_{i\cdot}^{{\textsf{L}}}-\bm{\phi}_{i\cdot}\right\|_{1}\left\|\hat{\bm{\Sigma}}_{x}\hat{\bm{\Omega}}-\mathbf{I}_{KN}\right\|_{\max}$

uniformly in $j\in[KN]$. Let
$\lambda_{1}=b\|\bm{\Omega}\|_{\ell_{1}}\lambda/2$, and consider the event

$\displaystyle\\{\|\mathbf{X}\mathbf{U}^{\prime}/T\|_{\max}\leq\lambda/2\\}\cap\\{\|\hat{\bm{\Sigma}}_{x}-\bm{\Sigma}_{x}\|_{\max}\leq\lambda_{1}/\|\bm{\Omega}\|_{\ell_{1}}\\}.$

Then this occurs with probability at least $1-O((N\vee T)^{-\nu})$ by Lemmas 3
and 4 with the proof of Proposition 1. On this event, the proof of Theorem 6
in Cai et al. (2011) establishes the bounds,
$\max_{j}\|\hat{\bm{\omega}}_{j}-{\bm{\omega}}_{j}\|_{1}\lesssim(\|\bm{\Omega}\|_{\ell_{1}}\lambda_{1})^{1-r}s_{\omega}$
and
$\|\hat{\bm{\Sigma}}_{x}\hat{\bm{\Omega}}-\mathbf{I}_{KN}\|_{\max}\leq\lambda_{1}$,
under Condition 4. Therefore, together with the Lasso bound derived in
Proposition 1 and $\lambda_{1}\lesssim M_{\omega}\lambda$, we obtain

$\displaystyle\max_{j\in[KN]}|r_{ij}|$
$\displaystyle\lesssim\sqrt{T}\lambda(\|\bm{\Omega}\|_{\ell_{1}}\lambda_{1})^{1-r}s_{\omega}/2+\frac{48s_{i}\sqrt{T}\lambda\lambda_{1}}{\gamma-8bs_{i}\lambda}$
$\displaystyle\lesssim\sqrt{T}\lambda
M_{\omega}^{1-r}(M_{\omega}\lambda)^{1-r}s_{\omega}+s_{i}\sqrt{T}\lambda^{2}M_{\omega}$
$\displaystyle\lesssim\left(s_{\omega}M_{\omega}^{2-2r}\lambda^{1-r}+s_{i}M_{\omega}\lambda\right)\sqrt{\log^{3}(N\vee
T)}~{}(=:\bar{r}_{i})$

for each $i\in[N]$ such that $s_{i}\lambda=o(1)$. This completes the proof. ∎

##  Appendix C Lemmas and their Proofs

######  Lemma 1.

If Condition 1 is true, then for any $i,j\in[N]$ and $s,t\in[T]$, there exist
constants $c_{uu},c_{m}>0$ such that

$\displaystyle(a)~{}~{}$
$\displaystyle\operatorname{\mathbb{P}}\left(\left|\frac{1}{T}\sum_{t=1}^{T}(u_{it}u_{jt}-\operatorname{\mathbb{E}}[u_{it}u_{jt}])\right|>x\right)\leq
2\exp\left\\{-\frac{1}{2}\left(\frac{Tx^{2}}{c_{uu}^{2}}\wedge\frac{Tx}{c_{uu}}\right)\right\\},$
$\displaystyle(b)~{}~{}$
$\displaystyle\operatorname{\mathbb{E}}\left[\max_{i\in[N]}\max_{t\in[T]}|u_{it}|^{m}\right]\leq
c_{m}\log^{m/2}(N\vee T),$

where $m\in\mathbb{N}$ is an arbitrary fixed constant.

###### Proof of Lemma 1.

Prove $(a)$. Since $u_{it}$ is sub-Gaussian by Condition 1, Lemma 2.7.7 of
Vershynin (2018) entails that for any $i,j\in[N]$,
$\\{u_{it}u_{jt}-\operatorname{\mathbb{E}}[u_{it}u_{jt}]\\}_{t}$ is a sequence
of i.i.d. (centered) sub-exponential random variables. Thus the result
directly follows by Bernstein’s inequality (Vershynin, 2018, Theorem 2.8.1).

Prove $(b)$. For any fixed $m\in\mathbb{N}$, the sub-Gaussian tail property
for $u_{it}$ in Condition 1 implies for all $x>0$,

$\displaystyle\operatorname{\mathbb{P}}(|u_{it}|^{m}>x)=\operatorname{\mathbb{P}}(|u_{it}|>x^{1/m})\leq
2\exp(-x^{2/m}/c_{u}).$

Using this tail probability with the union bound, we have

$\displaystyle\operatorname{\mathbb{E}}\left[\max_{i,t\in\mathbb{N}}\frac{|u_{it}|^{m}}{(1+\log
it)^{m/2}}\right]\leq\int_{0}^{\infty}\operatorname{\mathbb{P}}\left(\max_{i,t\in\mathbb{N}}\frac{|u_{it}|^{m}}{(1+\log
it)^{m/2}}>x\right)\mathop{}\\!\mathrm{d}x$
$\displaystyle\quad\leq\int_{0}^{(2c_{u})^{m/2}}\mathop{}\\!\mathrm{d}x+2\sum_{i=1}^{\infty}\sum_{t=1}^{\infty}\int_{(2c_{u})^{m/2}}^{\infty}\exp\left(-\frac{x^{2/m}(1+\log
i+\log t)}{c_{u}}\right)\mathop{}\\!\mathrm{d}x$
$\displaystyle\quad\leq(2c_{u})^{m/2}+2\sum_{i=1}^{\infty}i^{-2}\sum_{t=1}^{\infty}t^{-2}\int_{(2c_{u})^{m/2}}^{\infty}\exp\left(-\frac{x^{2/m}}{c_{u}}\right)\mathop{}\\!\mathrm{d}x,$

where the upper bound is further bounded by some positive constants. Thus
there exists some constant $M>0$ such that

$\displaystyle M$
$\displaystyle\geq\operatorname{\mathbb{E}}\left[\max_{i,t\in\mathbb{N}}\frac{|u_{it}|^{m}}{(1+\log
it)^{m/2}}\right]$
$\displaystyle\geq\operatorname{\mathbb{E}}\left[\max_{i\in[N],t\in[T]}\frac{|u_{it}|^{m}}{(1+\log
it)^{m/2}}\right]\geq\operatorname{\mathbb{E}}\left[\max_{i\in[N],t\in[T]}|u_{it}|^{m}\right]\frac{1}{(1+\log
NT)^{m/2}}.$

Replacing the constant factor appropriately gives the result. This completes
the proof. ∎

######  Lemma 2.

If Condition 2 is true, then the summability condition in (14) holds.
Furthermore, there exists a constant $\delta\in(0,1)$ such that for any
monotonically divergent sequence $r_{T}>0$ with sufficiently large $T$,

$\displaystyle\sum_{\ell=r_{T}}^{\infty}\|\mathbf{B}_{\ell}\|_{\infty}\leq\frac{\delta^{r_{T}}}{1-\delta}.$
(A.14)

In particular, the summability condition in (14) holds.

###### Proof of Lemma 2.

The VAR($K$) model is written as the VAR(1) companion form with coefficient

$\displaystyle\mathbf{A}=\begin{pmatrix}\bm{\Phi}_{1}&\bm{\Phi}_{2}&\dots&\bm{\Phi}_{K-1}&\bm{\Phi}_{K}\\\
\mathbf{I}_{N}&\mathbf{0}&\dots&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}_{N}&\dots&\mathbf{0}&\mathbf{0}\\\
\vdots&\vdots&&\vdots&\vdots\\\
\mathbf{0}&\mathbf{0}&\dots&\mathbf{I}_{N}&\mathbf{0}\end{pmatrix}.$

Condition 2 entails that the spectral radius of $\mathbf{A}$, defined as
$\rho(\mathbf{A})=\max_{j}|\lambda_{j}(\mathbf{A})|$, is strictly less than
one uniformly in $N$. This implies that the VAR(1) model is invertible to the
VMA($\infty$). Therefore, we have the representation
$\mathbf{y}_{t}=\sum_{\ell=0}^{\infty}\mathbf{B}_{\ell}\mathbf{u}_{t-\ell}$,
where $\mathbf{B}_{0}=\mathbf{I}_{N}$ and
$\mathbf{B}_{\ell}=\mathbf{J}^{\prime}\mathbf{A}^{\ell}\mathbf{J}$ with
$\mathbf{J}^{\prime}=(\mathbf{I}_{N},\mathbf{0}_{N\times(KN-N)})$ for
$\ell=1,2,\dots$. Note that

$\displaystyle\sum_{\ell=0}^{\infty}\|\mathbf{B}_{\ell}\|_{\infty}\leq\sum_{\ell=0}^{\infty}\|\mathbf{J}\|_{\infty}^{2}\|\mathbf{A}^{\ell}\|_{\infty}=\sum_{\ell=0}^{\infty}\|\mathbf{A}^{\ell}\|_{\infty}.$

Again, since we have $\rho(\mathbf{A})<1$ uniformly in $N$ by Condition 2, we
can always pick up a small constant $\varepsilon>0$ such that
$\delta:=\rho(\mathbf{A})+\varepsilon$ lies in $(0,1)$. We apply the Gelfand
formula (Horn and Johnson, 2012); for this choice of $\varepsilon$, there
exists $T^{*}$ such that $\|\mathbf{A}^{r_{T}}\|_{\infty}\leq\delta^{r_{T}}$
for all $T\geq T^{*}$, where $r_{T}>0$ is any monotonically diverging sequence
as $T\to\infty$. Therefore, for any $T\geq T^{*}$ we have

$\displaystyle\sum_{\ell=r_{T}}^{\infty}\|\mathbf{A}^{\ell}\|_{\infty}\leq\sum_{\ell=r_{T}}^{\infty}\delta^{\ell}=\frac{\delta^{r_{T}}}{1-\delta},$

which proves (A.14). Finally, because
$\limsup_{\ell\to\infty}\|\mathbf{A}^{\ell}\|_{\infty}^{1/\ell}\leq\delta<1$,
the summability condition in (14) holds by the root test. This completes the
proof. ∎

######  Lemma 3.

Let $b=\sum_{\ell=0}^{\infty}\left\|\mathbf{B}_{\ell}\right\|_{\infty}$. If
Conditions 1 and 2 are true, then for any $\bar{c}_{1}>0$ and arbitrary
divergent sequence $r_{T}>0$ with sufficiently large $T$, we have for all
$x>0$,

$\displaystyle\operatorname{\mathbb{P}}\left(\left\|T^{-1}\mathbf{X}\mathbf{U}^{\prime}\right\|_{\max}>x\right)$
$\displaystyle\leq
2KN^{2}\exp\left(-\frac{x^{2}T}{2\bar{c}_{1}^{2}}\right)+2r_{T}KN^{2}T\exp\left(-\frac{\bar{c}_{1}}{4bc_{uu}}\right)+2c_{1}^{2}KT\frac{\delta^{r_{T}}\log
N}{\bar{c}_{1}(1-\delta)},$

where constants $\delta\in(0,1)$ and $c_{uu},c_{1}>0$ are given by Lemmas 1
and 2, respectively. In particular, if we set
$x=\sqrt{2\bar{c}_{1}^{2}(\nu+2)\log(N\vee T)/T}$,
$\bar{c}_{1}=4bc_{uu}(\nu+4)\log(N\vee T)$, and $r_{T}=T$ with any constant
$\nu>0$, then the upper bound becomes $O((N\vee T)^{-\nu})$.

###### Proof of Lemma 3.

Following the definition, we have for any $x,\bar{c}_{1}>0$,

$\displaystyle\operatorname{\mathbb{P}}\left(\left\|T^{-1}\mathbf{X}\mathbf{U}^{\prime}\right\|_{\max}>x\right)=\operatorname{\mathbb{P}}\left(\max_{k\in[K]}\left\|T^{-1}\sum_{t=1}^{T}\mathbf{y}_{t-k}\mathbf{u}_{t}^{\prime}\right\|_{\max}>x\right)$
$\displaystyle\leq\operatorname{\mathbb{P}}\left(\max_{k\in[K]}\left\|T^{-1}\sum_{t=1}^{T}\mathbf{y}_{t-k}\mathbf{u}_{t}^{\prime}\right\|_{\max}>x\mid\max_{k\in[K]}\max_{t\in[T]}\|\mathbf{y}_{t-k}\mathbf{u}_{t}^{\prime}\|_{\max}\leq\bar{c}_{1}\right)$
$\displaystyle\qquad+\operatorname{\mathbb{P}}\left(\max_{k\in[K]}\max_{t\in[T]}\|\mathbf{y}_{t-k}\mathbf{u}_{t}^{\prime}\|_{\max}>\bar{c}_{1}\right).$

Consider the first term. Note for any $k\in[K]$ that
$\\{\mathbf{y}_{t-k}\mathbf{u}_{t}^{\prime}\\}$ forms a martingale difference
sequence with respect to the filtration,
$\mathcal{F}_{t}=\sigma\\{\mathbf{u}_{t-s}:s=0,1,\dots\\}$. Thus the union
bound and the Azuma-Hoeffding inequality (Bercu et al., 2015) give

$\displaystyle\operatorname{\mathbb{P}}\left(\max_{k\in[K]}\left\|T^{-1}\sum_{t=1}^{T}\mathbf{y}_{t-k}\mathbf{u}_{t}^{\prime}\right\|_{\max}>x\mid\max_{k\in[K]}\max_{t\in[T]}\|\mathbf{y}_{t-k}\mathbf{u}_{t}^{\prime}\|_{\max}\leq\bar{c}_{1}\right)$
$\displaystyle\leq
KN^{2}\max_{k\in[K]}\max_{i,j\in[N]}\operatorname{\mathbb{P}}\left(\left|T^{-1}\sum_{t=1}^{T}y_{i,t-k}u_{jt}\right|>x\mid\max_{k\in[K]}\max_{t\in[T]}\|\mathbf{y}_{t-k}\mathbf{u}_{t}^{\prime}\|_{\max}\leq\bar{c}_{1}\right)$
$\displaystyle\leq 2KN^{2}\exp\left(-\frac{x^{2}T}{2\bar{c}_{1}^{2}}\right).$

For the second term, applying the triangle and Hölder’s inequalities yields

$\displaystyle\left\|\mathbf{y}_{t-k}\mathbf{u}_{t}^{\prime}\right\|_{\max}=\left\|\sum_{\ell=0}^{\infty}\mathbf{B}_{\ell}\mathbf{u}_{t-k-\ell}\mathbf{u}_{t}^{\prime}\right\|_{\max}\leq\left(\sum_{\ell=0}^{r-1}+\sum_{\ell=r}^{\infty}\right)\left\|\mathbf{B}_{\ell}\right\|_{\infty}\left\|\mathbf{u}_{t-k-\ell}\mathbf{u}_{t}^{\prime}\right\|_{\max},$

and hence,

$\displaystyle\operatorname{\mathbb{P}}\left(\max_{k\in[K]}\max_{t\in[T]}\|\mathbf{y}_{t-k}\mathbf{u}_{t}^{\prime}\|_{\max}>\bar{c}_{1}\right)$
$\displaystyle\leq\operatorname{\mathbb{P}}\left(\max_{k\in[K]}\max_{t\in[T]}\sum_{\ell=0}^{r-1}\left\|\mathbf{B}_{\ell}\right\|_{\infty}\left\|\mathbf{u}_{t-k-\ell}\mathbf{u}_{t}^{\prime}\right\|_{\max}>\bar{c}_{1}/2\right)$
$\displaystyle\qquad+\operatorname{\mathbb{P}}\left(\max_{k\in[K]}\max_{t\in[T]}\sum_{\ell=r}^{\infty}\left\|\mathbf{B}_{\ell}\right\|_{\infty}\left\|\mathbf{u}_{t-k-\ell}\mathbf{u}_{t}^{\prime}\right\|_{\max}>\bar{c}_{1}/2\right).$

Consider the probability with $\sum_{\ell=0}^{r-1}$. The union bound and the
summability condition in Lemma 2 yield

$\displaystyle\operatorname{\mathbb{P}}\left(\max_{k\in[K]}\max_{t\in[T]}\sum_{\ell=0}^{r-1}\left\|\mathbf{B}_{\ell}\right\|_{\infty}\left\|\mathbf{u}_{t-k-\ell}\mathbf{u}_{t}^{\prime}\right\|_{\max}>\bar{c}_{1}/2\right)$
$\displaystyle\leq
rKT\max_{k\in[K]}\max_{t\in[T]}\max_{\ell=0,\dots,r-1}\operatorname{\mathbb{P}}\left(\left\|\mathbf{u}_{t-k-\ell}\mathbf{u}_{t}^{\prime}\right\|_{\max}>\frac{\bar{c}_{1}}{2\sum_{\ell=0}^{r-1}\left\|\mathbf{B}_{\ell}\right\|_{\infty}}\right)$
$\displaystyle\leq
rKN^{2}T\max_{k\in[K+r-1]}\max_{t\in[T]}\max_{i,j\in[N]}\operatorname{\mathbb{P}}\left(\left|u_{i,t-k}u_{jt}\right|>\frac{\bar{c}_{1}}{2b}\right)$
$\displaystyle\leq 2rKN^{2}T\exp\left(-\frac{\bar{c}_{1}}{4bc_{uu}}\right),$

where the last inequality holds by Lemma 1$(a)$ since $u_{i,t-k}u_{jt}$ is
sub-exponential (Vershynin, 2018, Proposition 2.7.1 and Lemma 2.7.7). Consider
the probability with $\sum_{\ell=r}^{\infty}$. By the union bound and the
Markov inequality, we have

$\displaystyle\operatorname{\mathbb{P}}\left(\max_{k\in[K]}\max_{t\in[T]}\sum_{\ell=r}^{\infty}\left\|\mathbf{B}_{\ell}\right\|_{\infty}\left\|\mathbf{u}_{t-k-\ell}\mathbf{u}_{t}^{\prime}\right\|_{\max}>\bar{c}_{1}/2\right)$
$\displaystyle\leq
2KT\sum_{\ell=r}^{\infty}\left\|\mathbf{B}_{\ell}\right\|_{\infty}\operatorname{\mathbb{E}}\left[\left\|\mathbf{u}_{t-k-\ell}\mathbf{u}_{t}^{\prime}\right\|_{\max}\right]/\bar{c}_{1}$
$\displaystyle\leq
2KT\frac{\delta^{r}}{\bar{c}_{1}(1-\delta)}\operatorname{\mathbb{E}}\left[\left\|\mathbf{u}_{t}\right\|_{\max}\right]^{2}\leq
2c_{1}^{2}KT\frac{\delta^{r}\log N}{\bar{c}_{1}(1-\delta)},$

where the last inequality holds by Lemma 1(b). Combining these upper bounds
completes the proof. ∎

######  Lemma 4.

Let $b=\sum_{\ell=0}^{\infty}\left\|\mathbf{B}_{\ell}\right\|_{\infty}$. If
Conditions 1 and 2 are true, then for any $\bar{c}_{2}\geq c_{uu}$ and
arbitrary divergent sequence $r_{T}>0$ with sufficiently large $T$, we have
for all $x>0$,

$\displaystyle\operatorname{\mathbb{P}}\left(\left\|T^{-1}\mathbf{X}\mathbf{X}^{\prime}-\operatorname{\mathbb{E}}[T^{-1}\mathbf{X}\mathbf{X}^{\prime}]\right\|_{\max}>x\right)$
$\displaystyle\leq
2K^{2}r_{T}^{2}N^{2}\exp\left(-\frac{x^{2}T}{8b^{4}\bar{c}_{2}^{2}}\right)+2K^{2}r_{T}^{2}N^{2}T\exp\left(-\frac{\bar{c}_{2}}{2c_{uu}}\right)+6K^{2}bc_{2}\frac{\delta^{r_{T}}\log
N}{x(1-\delta)},$

where constants $\delta\in(0,1)$ and $c_{uu},c_{2}>0$ are given by Lemmas 1
and 2, respectively. In particular, if we set
$x=\sqrt{8(\nu+4)b^{4}\bar{c}_{2}^{2}\log(N\vee T)/T}$,
$\bar{c}_{2}=2(\nu+5)c_{uu}\log(N\vee T)$, and $r_{T}=T$ with any constant
$\nu>0$, the upper bound becomes $O((N\vee T)^{-\nu})$.

###### Proof of Lemma 4.

For $g,h\in[K]$, define

$\displaystyle\mathbf{W}_{g,h}=T^{-1}\sum_{t=1}^{T}\left(\mathbf{u}_{t-g}\mathbf{u}_{t-h}^{\prime}-\operatorname{\mathbb{E}}[\mathbf{u}_{t-g}\mathbf{u}_{t-h}^{\prime}]\right).$

By the VMA($\infty$) representation in Lemma 2, we have

$\displaystyle\mathbf{y}_{t-g}\mathbf{y}_{t-h}^{\prime}=\sum_{\ell=0}^{\infty}\sum_{m=0}^{\infty}\mathbf{B}_{\ell}\mathbf{u}_{t-g-\ell}\mathbf{u}_{t-h-m}^{\prime}\mathbf{B}_{m}^{\prime},$

so that

$\displaystyle\left\|T^{-1}\mathbf{X}\mathbf{X}^{\prime}-\operatorname{\mathbb{E}}[T^{-1}\mathbf{X}\mathbf{X}^{\prime}]\right\|_{\max}\leq\max_{g,h\in[K]}\left\|T^{-1}\sum_{t=1}^{T}\left(\mathbf{y}_{t-g}\mathbf{y}_{t-h}^{\prime}-\operatorname{\mathbb{E}}[\mathbf{y}_{t-g}\mathbf{y}_{t-h}^{\prime}]\right)\right\|_{\max}$
$\displaystyle\leq\max_{g,h\in[K]}\sum_{\ell=0}^{\infty}\sum_{m=0}^{\infty}\left\|\mathbf{B}_{\ell}\mathbf{W}_{g+\ell,h+m}\mathbf{B}_{m}^{\prime}\right\|_{\max}$
$\displaystyle\leq\max_{g,h\in[K]}\left(\sum_{\ell=0}^{r-1}+\sum_{\ell=r}^{\infty}\right)\left(\sum_{m=0}^{r-1}+\sum_{m=r}^{\infty}\right)\|\mathbf{B}_{\ell}\|_{\infty}\left\|\mathbf{W}_{g+\ell,h+m}\mathbf{B}_{m}^{\prime}\right\|_{\max}$
$\displaystyle\leq\max_{g,h\in[K]}\left(\sum_{\ell=0}^{r-1}\sum_{m=0}^{r-1}+3\sum_{\ell=0}^{\infty}\sum_{m=r}^{\infty}\right)\|\mathbf{B}_{\ell}\|_{\infty}\|\mathbf{B}_{m}\|_{\infty}\left\|\mathbf{W}_{g+\ell,h+m}\right\|_{\max}.$

Therefore, we obtain

$\displaystyle\operatorname{\mathbb{P}}\left(\left\|T^{-1}\mathbf{X}\mathbf{X}^{\prime}-\operatorname{\mathbb{E}}[T^{-1}\mathbf{X}\mathbf{X}^{\prime}]\right\|_{\max}>x\right)$
$\displaystyle\leq\operatorname{\mathbb{P}}\left(\max_{g,h\in[K]}\sum_{\ell=0}^{r-1}\sum_{m=0}^{r-1}\|\mathbf{B}_{\ell}\|_{\infty}\|\mathbf{B}_{m}\|_{\infty}\left\|\mathbf{W}_{g+\ell,h+m}\right\|_{\max}>x/2\right)$
$\displaystyle\quad+\operatorname{\mathbb{P}}\left(\max_{g,h\in[K]}3\sum_{\ell=0}^{\infty}\sum_{m=r}^{\infty}\|\mathbf{B}_{\ell}\|_{\infty}\|\mathbf{B}_{m}\|_{\infty}\left\|\mathbf{W}_{g+\ell,h+m}\right\|_{\max}>x/2\right).$

We consider the first probability. The union bound gives

$\displaystyle\operatorname{\mathbb{P}}\left(\max_{g,h\in[K]}\sum_{\ell=0}^{r-1}\sum_{m=0}^{r-1}\|\mathbf{B}_{\ell}\|_{\infty}\|\mathbf{B}_{m}\|_{\infty}\left\|\mathbf{W}_{g+\ell,h+m}\right\|_{\max}>x/2\right)$
$\displaystyle\leq
K^{2}\max_{g,h\in[K]}\operatorname{\mathbb{P}}\left(\max_{\ell,m=0,\dots,r-1}\left\|\mathbf{W}_{g+\ell,h+m}\right\|_{\max}\sum_{\ell=0}^{r-1}\sum_{m=0}^{r-1}\|\mathbf{B}_{\ell}\|_{\infty}\|\mathbf{B}_{m}\|_{\infty}>x/2\right)$
$\displaystyle\leq
K^{2}r^{2}\max_{g,h\in[K]}\max_{\ell,m=0,\dots,r-1}\operatorname{\mathbb{P}}\left(\left\|\mathbf{W}_{g+\ell,h+m}\right\|_{\max}>x/(2b^{2})\right)$
$\displaystyle=K^{2}r^{2}\max_{k\in[K+r-1]}\operatorname{\mathbb{P}}\left(\left\|\mathbf{W}_{k,1}\right\|_{\max}>x/(2b^{2})\right).$

If $k=1$, each component of $\mathbf{W}_{1,1}$ is a sample average of the
i.i.d. random variables,
$\\{u_{it}u_{jt}-\operatorname{\mathbb{E}}[u_{it}u_{jt}]\\}_{t}$. Clearly this
is a martingale difference sequence with respect to filtration
$\mathcal{F}_{t}=\sigma\\{u_{i,t-s},u_{j,t-s}:s=0,1,\dots\\}$. If $k\geq 2$,
we have $\operatorname{\mathbb{E}}[u_{i,t-k}u_{j,t-1}]=0$ and the sequence,
$\\{u_{i,t-k+1}u_{jt}\\}_{t}$, is also martingale difference with respect to
$\mathcal{F}_{t}$. Therefore, the Azuma-Hoeffding inequality (Bercu et al.,
2015) with the conditioning argument and the union bound as in the proof of
Lemma 3, we have

$\displaystyle\max_{k\in[K+r-1]}\operatorname{\mathbb{P}}\left(\left\|\mathbf{W}_{k,1}\right\|_{\max}>x/(2b^{2})\right)$
$\displaystyle\leq\max_{k\in[K+r-1]}\operatorname{\mathbb{P}}\left(\left\|\mathbf{W}_{k,1}\right\|_{\max}>x/(2b^{2})\mid\max_{t\in[T]}\left\|\mathbf{u}_{t-k}\mathbf{u}_{t-1}^{\prime}-\operatorname{\mathbb{E}}[\mathbf{u}_{t-k}\mathbf{u}_{t-1}^{\prime}]\right\|_{\max}\leq\bar{c}_{2}\right)$
$\displaystyle\quad+\max_{k\in[K+r-1]}\operatorname{\mathbb{P}}\left(\max_{t\in[T]}\left\|\mathbf{u}_{t-k}\mathbf{u}_{t-1}^{\prime}-\operatorname{\mathbb{E}}[\mathbf{u}_{t-k}\mathbf{u}_{t-1}^{\prime}]\right\|_{\max}>\bar{c}_{2}\right)$
$\displaystyle\leq
2N^{2}\exp\left(-\frac{x^{2}T}{8b^{4}\bar{c}_{2}^{2}}\right)+2N^{2}T\exp\left(-\frac{\bar{c}_{2}}{2c_{uu}}\right).$

We next consider the second probability. The union bound and the Markov
inequality with Lemma 2 give

$\displaystyle\operatorname{\mathbb{P}}\left(3\max_{g,h\in[K]}\sum_{\ell=0}^{\infty}\sum_{m=r}^{\infty}\|\mathbf{B}_{\ell}\|_{\infty}\|\mathbf{B}_{m}\|_{\infty}\left\|\mathbf{W}_{g+\ell,h+m}\right\|_{\max}>x/2\right)$
$\displaystyle\leq
6K^{2}\max_{g,h\in[K]}\sum_{\ell=0}^{\infty}\sum_{m=r}^{\infty}\|\mathbf{B}_{\ell}\|_{\infty}\|\mathbf{B}_{m}\|_{\infty}\operatorname{\mathbb{E}}\left\|\mathbf{W}_{g+\ell,h+m}\right\|_{\max}/x$
$\displaystyle\leq
6K^{2}\frac{b\delta^{r}}{x(1-\delta)}\max_{k\in\\{1,2\\}}\operatorname{\mathbb{E}}\left[\left\|\mathbf{u}_{t-k}\mathbf{u}_{t-1}^{\prime}-\operatorname{\mathbb{E}}[\mathbf{u}_{t-k}\mathbf{u}_{t-1}^{\prime}]\right\|_{\max}\right].$

By Lemma 1, the last expectation is evaluated as

$\displaystyle\max_{k\in\\{1,2\\}}\operatorname{\mathbb{E}}\left[\left\|\mathbf{u}_{t-k}\mathbf{u}_{t-1}^{\prime}-\operatorname{\mathbb{E}}\mathbf{u}_{t-k}\mathbf{u}_{t-1}^{\prime}\right\|_{\max}\right]\leq
c_{2}\log N.$

Combining the obtained results so far yields the desired inequality. This
completes the proof. ∎

######  Lemma 5.

Let $\bm{\delta}_{i\cdot}^{\prime}\in\mathbb{R}^{KN}$ denote the $i$-th column
vector of $\bm{\Delta}^{\prime}\in\mathbb{R}^{KN\times N}$ with
$\bm{\Delta}=\hat{\bm{\Phi}}-\bm{\Phi}^{0}$. If inequality (A.11) is true,
then it holds that

$\displaystyle\|\bm{\delta}_{i\mathcal{S}_{i}^{c}}\|_{1}\leq
3\|\bm{\delta}_{i\mathcal{S}_{i}}\|_{1}.$

###### Proof of Lemma 5.

Note that $\bm{\phi}=\bm{\phi}_{\mathcal{S}}$ and
$\mathbf{v}=\mathbf{v}_{\mathcal{S}}+\mathbf{v}_{\mathcal{S}^{c}}$ for any
$KN$-dimensional vector $\mathbf{v}$. Hence the lower bound of
$\|\hat{\bm{\phi}}_{i\cdot}\|_{1}$ is computed as

$\displaystyle\|\hat{\bm{\phi}}_{i\cdot}\|_{1}$
$\displaystyle=\|\bm{\phi}_{i\cdot}+\bm{\delta}_{i\cdot}\|_{1}=\|\bm{\phi}_{i\mathcal{S}_{i}}+\bm{\delta}_{i\mathcal{S}_{i}}+\bm{\delta}_{i\mathcal{S}_{i}^{c}}\|_{1}\geq\|\bm{\phi}_{i\mathcal{S}_{i}}+\bm{\delta}_{i\mathcal{S}_{i}^{c}}\|_{1}-\|\bm{\delta}_{i\mathcal{S}_{i}}\|_{1}$
$\displaystyle=\|\bm{\phi}_{i\mathcal{S}_{i}}\|_{1}+\|\bm{\delta}_{i\mathcal{S}_{i}^{c}}\|_{1}-\|\bm{\delta}_{i\mathcal{S}_{i}}\|_{1}=\|\bm{\phi}_{i\cdot}\|_{1}+\|\bm{\delta}_{i\mathcal{S}_{i}^{c}}\|_{1}-\|\bm{\delta}_{i\mathcal{S}_{i}}\|_{1},$

where the last equality holds by decomposability of the $\ell_{1}$-norm. Thus
we obtain

$\displaystyle\|\bm{\phi}_{i\cdot}\|_{1}-\|\hat{\bm{\phi}}_{i\cdot}\|_{1}\leq\|\bm{\delta}_{i\mathcal{S}_{i}}\|_{1}-\|\bm{\delta}_{i\mathcal{S}_{i}^{c}}\|_{1}.$

From (A.11), we have

$\displaystyle 0\leq(2T)^{-1}\|\bm{\delta}_{i\cdot}\mathbf{X}\|_{2}^{2}$
$\displaystyle\leq
2^{-1}\lambda\|\bm{\delta}_{i\cdot}\|_{1}+\lambda\|\bm{\phi}_{i\cdot}\|_{1}-\lambda\|\hat{\bm{\phi}}_{i\cdot}\|_{1}$
$\displaystyle\leq
2^{-1}\lambda\|\bm{\delta}_{i\mathcal{S}_{i}}\|_{1}+2^{-1}\lambda\|\bm{\delta}_{i\mathcal{S}_{i}^{c}}\|_{1}+\lambda\|\bm{\delta}_{i\mathcal{S}_{i}}\|_{1}-\lambda\|\bm{\delta}_{i\mathcal{S}_{i}^{c}}\|_{1}$
$\displaystyle=(3/2)\lambda\|\bm{\delta}_{i\mathcal{S}_{i}}\|_{1}-2^{-1}\lambda\|\bm{\delta}_{i\mathcal{S}_{i}^{c}}\|_{1},$

giving the desired inequality, $\|\bm{\delta}_{i\mathcal{S}_{i}^{c}}\|_{1}\leq
3\|\bm{\delta}_{i\mathcal{S}_{i}}\|_{1}$. This completes the proof. ∎

######  Lemma 6.

Suppose the same conditions as Proposition 1. We have

$\displaystyle\min_{\bm{\delta}_{i\cdot}\in\mathcal{D}_{i}}\frac{T^{-1}\|\bm{\delta}_{i\cdot}\mathbf{X}\|_{\mathrm{F}}^{2}}{\|\bm{\delta}_{i\cdot}\|_{2}^{2}}\geq\gamma-16s_{i}\left\|T^{-1}\mathbf{X}\mathbf{X}^{\prime}-\bm{\Sigma}_{x}\right\|_{\max},$

where
$\mathcal{D}_{i}=\\{\bm{\delta}\in\mathbb{R}^{KN}:\|\bm{\delta}_{\mathcal{S}_{i}^{c}}\|_{1}\leq
3\|\bm{\delta}_{\mathcal{S}_{i}}\|_{1}\\}$, and $\gamma>0$ is some constant.

###### Proof of Lemma 6.

Let $\bm{\delta}_{i\cdot}^{\prime}\in\mathbb{R}^{KN}$ denote the $i$th column
vector of $\bm{\Delta}^{\prime}\in\mathbb{R}^{KN\times N}$. We see that

$\displaystyle\frac{T^{-1}\|\bm{\delta}_{i\cdot}\mathbf{X}\|_{2}^{2}}{\|\bm{\delta}_{i\cdot}\|_{2}^{2}}=\frac{\bm{\delta}_{i\cdot}\operatorname{\mathbb{E}}[T^{-1}\mathbf{X}\mathbf{X}^{\prime}]\bm{\delta}_{i\cdot}^{\prime}}{\bm{\delta}_{i\cdot}\bm{\delta}_{i\cdot}^{\prime}}-\frac{\bm{\delta}_{i\cdot}\left(\operatorname{\mathbb{E}}[T^{-1}\mathbf{X}\mathbf{X}^{\prime}]-T^{-1}\mathbf{X}\mathbf{X}^{\prime}\right)\bm{\delta}_{i\cdot}^{\prime}}{\bm{\delta}_{i\cdot}\bm{\delta}_{i\cdot}^{\prime}}.$
(A.15)

Then it is sufficient to show that uniformly in
$\bm{\delta}_{i\cdot}\in\mathcal{D}_{i}$ the first term is bounded from below
by some constant and the second term converges to zero with high probability.

We consider the first term of (A.15). In view of the Rayleigh quotient with
Condition 3, we obtain

$\displaystyle\min_{\bm{\delta}_{i\cdot}\in\mathcal{D}_{i}}\frac{\bm{\delta}_{i\cdot}\operatorname{\mathbb{E}}[T^{-1}\mathbf{X}\mathbf{X}^{\prime}]\bm{\delta}_{i\cdot}^{\prime}}{\bm{\delta}_{i\cdot}\bm{\delta}_{i\cdot}^{\prime}}\geq\gamma.$
(A.16)

We next consider the second term of (A.15). By applying Hölder’s inequality
twice, the numerator is bounded as

$\displaystyle\left|\bm{\delta}_{i\cdot}^{\prime}\left(\operatorname{\mathbb{E}}[T^{-1}\mathbf{X}\mathbf{X}^{\prime}]-T^{-1}\mathbf{X}\mathbf{X}^{\prime}\right)\bm{\delta}_{i\cdot}^{\prime}\right|\leq\left\|\bm{\delta}_{i\cdot}\left(\operatorname{\mathbb{E}}[T^{-1}\mathbf{X}\mathbf{X}^{\prime}]-T^{-1}\mathbf{X}\mathbf{X}^{\prime}\right)\right\|_{\max}\|\bm{\delta}_{i\cdot}\|_{1}$
$\displaystyle\quad\leq\left\|\operatorname{\mathbb{E}}[T^{-1}\mathbf{X}\mathbf{X}^{\prime}]-T^{-1}\mathbf{X}\mathbf{X}^{\prime}\right\|_{\max}\|\bm{\delta}_{i\cdot}\|_{1}^{2}.$
(A.17)

We further compute the upper bound of $\|\bm{\delta}_{i\cdot}\|_{1}^{2}$.
Lemma 5 yields

$\displaystyle\|\bm{\delta}_{i\cdot}\|_{1}^{2}=\left(\|\bm{\delta}_{\mathcal{S}_{i}}\|_{1}+\|\bm{\delta}_{\mathcal{S}_{i}^{c}}\|_{1}\right)^{2}\leq
16\|\bm{\delta}_{\mathcal{S}_{i}}\|_{1}^{2}\leq
16s_{i}\|\bm{\delta}_{\mathcal{S}_{i}}\|_{2}^{2}\leq
16s_{i}\|\bm{\delta}_{i\cdot}\|_{2}^{2}.$ (A.18)

Combining (A.17) and (A.18) and rearranging the terms yield

$\displaystyle\max_{\bm{\delta}_{i\cdot}\in\mathcal{D}}\frac{\bm{\delta}_{i\cdot}\left(\operatorname{\mathbb{E}}[T^{-1}\mathbf{X}\mathbf{X}^{\prime}]-T^{-1}\mathbf{X}\mathbf{X}^{\prime}\right)\bm{\delta}_{i\cdot}^{\prime}}{\|\bm{\delta}_{i\cdot}\|_{2}^{2}}\leq
16s_{i}\left\|\operatorname{\mathbb{E}}[T^{-1}\mathbf{X}\mathbf{X}^{\prime}]-T^{-1}\mathbf{X}\mathbf{X}^{\prime}\right\|_{\max}.$
(A.19)

From (A.15), (A.16), and (A.19), we obtain

$\displaystyle\min_{\bm{\delta}_{i\cdot}\in\mathcal{D}}\frac{T^{-1}\|\bm{\delta}_{i\cdot}\mathbf{X}\|_{2}^{2}}{\|\bm{\delta}_{i\cdot}\|_{2}^{2}}\geq\gamma-16s_{i}\left\|\operatorname{\mathbb{E}}[T^{-1}\mathbf{X}\mathbf{X}^{\prime}]-T^{-1}\mathbf{X}\mathbf{X}^{\prime}\right\|_{\max}.$

This completes the proof. ∎

######  Lemma 7.

If Conditions 1–5 are true, then the following inequality holds for all
$i\in[N]$ such that $\bar{r}_{i}=o(1)$ with probability at least $1-O((N\vee
T)^{-\nu})$:

$\displaystyle(a)~{}$
$\displaystyle~{}\max_{j\in[KN]}\left|\hat{\sigma}_{i}^{2}\hat{\bm{\omega}}_{j}^{\prime}\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}-\sigma_{i}^{2}\omega_{j}^{2}\right|\lesssim
M_{\omega}^{2}\lambda,$ $\displaystyle(b)~{}$
$\displaystyle~{}\max_{j\in[KN]}\left|\hat{\sigma}_{i}^{2}\hat{\omega}_{j}^{2}-\sigma_{i}^{2}\omega_{j}^{2}\right|\lesssim
M_{\omega}^{2}\lambda,$ $\displaystyle(c)~{}$
$\displaystyle~{}\max_{j\in[KN]}\left|\frac{\sigma_{i}\omega_{j}}{\hat{\sigma}_{i}\sqrt{\hat{\bm{\omega}}_{j}^{\prime}\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}}}-1\right|\lesssim\frac{M_{\omega}^{2}\lambda}{\left|\gamma^{2}-O(\sqrt{M_{\omega}^{2}\lambda})\right|},$
$\displaystyle(d)~{}$
$\displaystyle~{}\max_{j\in[KN]}\left|\frac{\sigma_{i}\omega_{j}}{\hat{\sigma}_{i}\hat{\omega}_{j}}-1\right|\lesssim\frac{M_{\omega}^{2}\lambda}{\left|\gamma^{2}-O(\sqrt{M_{\omega}^{2}\lambda})\right|}.$

In consequence, if $\bar{v}_{i}=\bar{r}_{i}+M_{\omega}^{2}\lambda$ is $o(1)$,
all the upper bounds are $o(1)$.

###### Proof of Lemma 7.

For any $(i,j)\in\mathcal{H}$, we consider two bounds:

$\displaystyle(a)$
$\displaystyle~{}~{}\left|\hat{\sigma}_{i}^{2}\hat{\bm{\omega}}_{j}^{\prime}\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}-\sigma_{i}^{2}\omega_{j}^{2}\right|\leq\left(\left|\hat{\sigma}_{i}^{2}-\sigma_{i}^{2}\right|+\sigma_{i}^{2}\right)\left|\hat{\bm{\omega}}_{j}^{\prime}\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}-\omega_{j}^{2}\right|+\left|\hat{\sigma}_{i}^{2}-\sigma_{i}^{2}\right|\omega_{j}^{2},$
$\displaystyle(b)$
$\displaystyle~{}~{}\left|\hat{\sigma}_{i}^{2}\hat{\omega}_{j}^{2}-\sigma_{i}^{2}\omega_{j}^{2}\right|\leq\left(\left|\hat{\sigma}_{i}^{2}-\sigma_{i}^{2}\right|+\sigma_{i}^{2}\right)\left|\hat{\omega}_{j}^{2}-\omega_{j}^{2}\right|+\left|\hat{\sigma}_{i}^{2}-\sigma_{i}^{2}\right|\omega_{j}^{2}.$

Thus to complete the proof, we derive upper bounds of

$\displaystyle(i):=\max_{j}|\hat{\bm{\omega}}_{j}^{\prime}\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}-\omega_{j}^{2}|,~{}~{}(i)^{\prime}:=\max_{j}\left|\hat{\omega}_{j}^{2}-\omega_{j}^{2}\right|,~{}~{}\text{and}~{}~{}(ii):=|\hat{\sigma}_{i}^{2}-\sigma_{i}^{2}|.$

where $(i)^{\prime}$ is bounded as

$\displaystyle(i)^{\prime}\lesssim\max_{j}\|{\bm{\omega}}_{j}\|\lambda_{1}\lesssim
M_{\omega}^{2}\lambda$

by Theorem 6 of Cai et al. (2011). Bound $(i)$. Since
$w_{jj}={\bm{\omega}}_{j}^{\prime}\mathbf{e}_{j}$, we have

$\displaystyle(i)$
$\displaystyle\leq\max_{j}|(\hat{\bm{\omega}}_{j}-{\bm{\omega}}_{j})^{\prime}\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}|+|{\bm{\omega}}_{j}^{\prime}(\mathbf{e}_{j}-\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j})|$
$\displaystyle\leq\max_{j}|(\hat{\bm{\omega}}_{j}-{\bm{\omega}}_{j})^{\prime}(\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}-\mathbf{e}_{j})|+\max_{j}|(\hat{\bm{\omega}}_{j}-{\bm{\omega}}_{j})^{\prime}\mathbf{e}_{j}|+\max_{j}|{\bm{\omega}}_{j}^{\prime}(\mathbf{e}_{j}-\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j})|$
$\displaystyle\leq\max_{j}\|\hat{\bm{\omega}}_{j}-{\bm{\omega}}_{j}\|_{1}\|\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}-\mathbf{e}_{j}\|_{\infty}+\max_{j}\|\hat{\bm{\omega}}_{j}-{\bm{\omega}}_{j}\|_{\infty}+\max_{j}\|{\bm{\omega}}_{j}\|_{1}\|\mathbf{e}_{j}-\hat{\bm{\Sigma}}_{x}\hat{\bm{\omega}}_{j}\|_{\infty}$
$\displaystyle\leq\max_{j}\|\hat{\bm{\omega}}_{j}-{\bm{\omega}}_{j}\|_{1}\lambda_{1}+2M_{\omega}\lambda_{1}.$

By the proof of Proposition 2, it holds with probability at least $1-O((N\vee
T)^{-\nu})$ that

$\displaystyle\max_{j}\|\hat{\bm{\omega}}_{j}-{\bm{\omega}}_{j}\|_{1}\lesssim(\max_{j}\|{\bm{\omega}}_{j}\|_{\ell_{1}}\lambda_{1})^{1-r}s_{\omega}\lesssim
M_{\omega}^{2-2r}\lambda^{1-r}s_{\omega}.$

where we have used $\lambda_{1}=b\max_{j}\|{\bm{\omega}}_{j}\|_{1}\lambda/2$.
Then Conditions 4 and 5 yield

$\displaystyle(i)\lesssim
M_{\omega}^{3-2r}\lambda^{2-r}s_{\omega}+M_{\omega}^{2}\lambda.$

Finally bound $(ii)$. We have

$\displaystyle(ii)\leq
T^{-1}\left|\sum_{t=1}^{T}(\hat{u}_{it}^{2}-u_{it}^{2})\right|+\left|T^{-1}\sum_{t=1}^{T}({u}_{it}^{2}-\operatorname{\mathbb{E}}u_{it}^{2})\right|.$

Consider the first term. Because
$\hat{\mathbf{u}}_{i\cdot}-\mathbf{u}_{i\cdot}=-(\hat{\bm{\phi}}_{i\cdot}^{L}-\bm{\phi}_{i\cdot})\mathbf{X}$,
we see that

$\displaystyle
T^{-1}\left|\sum_{t=1}^{T}(\hat{u}_{it}^{2}-u_{it}^{2})\right|\leq
T^{-1}\left|\sum_{t=1}^{T}(\hat{u}_{it}-u_{it})^{2}\right|+2T^{-1}\left|\sum_{t=1}^{T}(\hat{u}_{it}-u_{it})u_{it}\right|$