# Sparse Generalized Yule–Walker Estimation for Large Spatio-temporal
Autoregressions with an Application to NO2 Satellite Data

Hanno Reuvers  Department of Econometrics, Erasmus University Rotterdam, 3062
PA Rotterdam, The Netherlands  Etienne Wijler  Department of Quantitative
Economics, Maastricht University, 6200 MD Maastricht, The Netherlands

###### Abstract

We consider sparse estimation of a class of high-dimensional spatio-temporal
models. Unlike classical spatial autoregressive models, we do not rely on a
predetermined spatial interaction matrix. Instead, under the assumption of
sparsity, we estimate the relationships governing both the spatial and
temporal dependence in a fully data-driven way by penalizing a set of Yule-
Walker equations. While this regularization can be left unstructured, we also
propose a customized form of shrinkage to further exploit diagonally
structured forms of sparsity that follow intuitively when observations
originate from spatial grids such as satellite images. We derive finite sample
error bounds for this estimator, as well estimation consistency in an
asymptotic framework wherein the sample size and the number of spatial units
diverge jointly. A simulation exercise shows strong finite sample performance
compared to competing procedures. As an empirical application, we model
satellite measured $\text{NO}_{2}$ concentrations in London. Our approach
delivers forecast improvements over a competitive benchmark and we discover
evidence for strong spatial interactions between sub-regions.

##  1 Introduction

Spatial econometrics provides models to describe the cross-sectional
dependence among entities observed at different locations. A key quantity is
the so-called spatial weight matrix. Parametric choices for this spatial
weight matrix, e.g. adjacency matrix, inverse physical distance and spatial
contiguity weights, are very popular and give rise to interesting spatial
models and novel estimation techniques. The formulation of these spatial
models typically involves two steps. First, the product of the spatial weight
matrix and the vector of the contemporaneous observations is constructed to
create a spatial lag. Second, spatial dependencies are modeled by augmenting
an existing econometric model with the spatial lag. The necessity for
dedicated estimation techniques stems from the endogeneity problem caused by
the contemporaneous variable appearing simultaneously in the dependent
variable and in the spatial lag. Examples of this two-step approach include:
the spatial autoregressive model (a linear regression model with spatial lag)
with a Gaussian quasi-maximum likelihood estimator (QMLE), see e.g. Lee,
(2004); the QMLE estimation of stationary spatial panel with fixed effects
detailed in Yu et al., (2008); the extension of these spatial panels to
include spatially autoregressive disturbances as in Lee and Yu, (2010); Yu et
al., (2012) cover a non-stationary case where units can be spatially
cointegrated; and the computationally beneficial generalized method of moments
(GMM) estimator by Lee and Yu, (2014). Possible application areas for spatial
econometric models are wide-spread. A non-exhaustive list of examples include
product demand across U.S. states, pollution concentrations across measurement
sites, housing prices with location information, co-movement between stock
indices, and interactions between social network users have all been studied
using spatial econometric methods.

The choice of the spatial weight matrix is a key element of the model
specification. Nevertheless, it also involves a selection process that can
feel somewhat arbitrary and/or tedious. The arbitrariness might prevail when
practical considerations fail to suggest a particular mechanism for the
spatial interactions. Alternatively, in large data sets with many spatial
sites, there are simply too many spatial interactions to be modeled. Two
solutions have been proposed in the literature. First, there is the option to
use multiple weight matrices. These weight matrices can be combined through
either a convex combination (see, e.g. Debarsy and LeSage, (2018)) or a model
averaging approach as in Zhang and Yu, (2018). Second, at the expense of
estimating many parameters, we might decide to directly infer the spatial
interactions from the data. In both cases, the estimation method has to
account for the endogeneity issue.

We propose the SPatial LAsso-type SHrinkage (SPLASH) estimator. Building upon
previous works by Dou et al., (2016) and Gao et al., (2019), we study spatio-
temporal regressions. Apart from a more generous bandwidth upper bound, this
SPLASH estimator leaves the spatial weight matrix and autoregressive matrix
unspecified while employing a lasso approach to recover sparse solutions. The
endogeneity problem is resolved by using an estimator which solves the
generalized Yule-Walker equations. We make five contributions. First, assuming
sparsity in the coefficient matrices and general mixing conditions on the
innovations, we derive finite-sample performance bounds for the estimation and
prediction error of our estimator. We subsequently utilize these bounds to
derive asymptotic consistency in a variety of settings. For example, in the
special case of a finite bandwidth and unstructured sparsity, it is
demonstrated that the number of spatial units $N$ may grow at any polynomial
rate of the number of temporal observations $T$. Second, we adopt a banded
estimation procedure for the autocovariance matrices that underlie the
generalized Yule-Walker equations. The faster convergence rates of these
banded autocovariance matrix estimators are shown to translate into better
convergence rates of our SPLASH estimator. Third, we show that dependence
between neighbouring units that are ordered on a spatial grid translates to
diagonally structured forms of sparsity in the spatial weight matrix, for
which we develop a tailored regularization procedure reminiscent of the sparse
group lasso. Fourth and final, we consider an empirical application in which
we employ SPLASH to predict $\text{NO}_{2}$ concentrations in London based on
satellite data.

Elaborating on the empirical application, we collect daily $\text{NO}_{2}$
column densities from August 2018 to October 2020, recorded by the
TROPOspheric Monitoring Instrument (TROPOMI) on board of the Corpernicus
Sentinel-5 Precursor satellite. Each spatial unit in this application
corresponds to an aggegration of a small number of pixels on the satellite
image. We find that SPLASH constructs more accurate one-step ahead predictions
for all spatial units compared to the procedure in Gao et al., (2019), while
outperforming a competitive penalized VAR benchmark for the majority of
spatial units. In addition, we find evidence for spatial interactions between
first-order neighbours and second-order neighbours (i.e. neighbours of
neighbours).

There are two strands of literature that are closely linked to this work: the
literature on the estimation of (nonparametric) spatial weight matrices and
the literature on spatio-temporal vector autoregressions. Some similarities
and differences are as follows. Lam and Souza, (2014) consider a model
specification where the spatial units depend linearly on a spatial lag and a
set of exogenous regressors. The adaptive lasso is proven to select the
correct sparsity pattern. To solve the endogeneity issue, they require the
error variance to decay to zero as the time dimension grows large. Ahrens and
Bhattacharjee, (2015) solve the endogeneity problem using external
instruments. Their two-step lasso estimation procedure selects the relevant
instruments in the first step and the relevant spatial interactions in the
second step. The theoretical properties of this estimator are derived using
moderate deviation theory as in Jing et al., (2003). This approach requires
the instruments and the idiosyncratic component to be serially independent.
Clearly, a serial independence assumption is unrealistic for the spatio-
temporal models we consider here. Finally, Lam and Souza, (2019) augment a
spatial lag model with a set of potentially endogenous variables (the
augmenting set). They decompose the spatial weight matrix into a pre-
determined component based on expert knowledge and a sparse adjustment matrix
that represents specification errors. The adjustment matrix is sparsely
estimated based on a penalized version of instrumental variables (IV)
regression. If these instrumental variables are selected as temporal lags of
the dependence variable, then their IV regressions are similar to generalized
Yule-Walker estimation. In contrast to our approach, Lam and Souza, (2019) do
not regularize the interactions between the dependent variables and the
variables in the augmenting set, and they assume the number of such
interactions to be fixed. A fixed number of interactions is inappropriate in
high-dimensional settings in which the number of spatial units is allowed to
diverge.

Most closely related to the work presented in this paper are Gao et al.,
(2019), and Ma et al., (2021). Both of these paper consider the same model and
an estimation procedure that relies on the generalized Yule-Walker equations.
The key difference with our paper lies in the method by which the model
complexity is controlled during estimation. Gao et al., (2019) assume the
coefficient matrices to be banded with a bandwidth that is small in relation
to the number of spatial units. The bandwidth is determined from the data and
all parameters within the selected bandwidth are left unregularized. Our
SPLASH estimator, however, has the ability to exploit (structured) sparsity
within the bandwidth and thereby improve estimation and forecasting
performance. In addition, after imposing a general upper bound on the
bandwidth to ensure identification, no a priori choice regarding the bandwidth
is required when using SPLASH. The recently developed bagging approach in Ma
et al., (2021) does allow for sparsity within the bands, yet it also requires
the calculation of so-called solution paths. That is, a forward addition and
backward deletion stage are needed to determine the variables that enter the
final model specification. In contrast, the SPLASH estimator provides this
solution at once. Furthermore, their approach is not designed to detect
diagonally structured forms of sparsity, while the ability to do so results in
clear performance improvements of SPLASH in both the simulations and the
empirical application considered below.

This paper is organized as follows. Section 2 introduces the spatio-temporal
vector autoregression and the banded autocovariance estimator that underlies
the generalized Yule-Walker estimation approach. The SPLASH estimator and its
theoretical properties are discussed in Section 3. The simulation results in
Section 4 and the empirical application in Section 5 demonstrate the benefits
of the SPLASH estimator. Section 6 concludes.

### Notation

The indicator function $\mathbbm{1}_{\\{A\\}}$ equals 1 if $A$ is true and
zero otherwise. For a vector $\bm{x}\in\mathbb{R}^{N}$, the $L_{p}$-norm of
$\bm{x}$ is denoted
$\left\lVert\bm{x}\right\rVert_{p}=\big{(}\sum_{i=1}^{N}|x_{i}|^{p}\big{)}^{1/p}$,
with $\left\lVert\bm{x}\right\rVert_{\infty}=\max_{i}\left\lvert
x_{i}\right\rvert$ as an important special case. The total number elements in
$\bm{x}$ is denoted by $\left\lvert\bm{x}\right\rvert$, while the number of
non-zero elements in $\bm{x}$ is denoted by
$\mathscr{M}(\bm{x})=\sum_{i=1}^{N}\mathbbm{1}\\{x_{i}\neq 0\\}$. The Orlicz
norm is defined as
$\left\lVert\cdot\right\rVert_{\psi}=\inf\left\\{c>0:\operatorname{\mathbb{E}}\big{[}\psi\left(|\cdot|/c\right)\big{]}\leq
1\right\\}$ for any $\psi(\cdot):\mathbb{R}^{+}\to\mathbb{R}^{+}$ being a
convex, increasing function with $\psi(0)=0$ and $\psi(x)\to\infty$ as
$x\to\infty$. In addition, we rely on several types of matrix norms. For a
matrix $\bm{A}\in\mathbb{R}^{M\times N}$, the matrix norms induced by the
vector $L_{p}$-norms are given by
$\left\lVert\bm{A}\right\rVert_{p}=\sup_{\bm{x}\in\mathbb{R}^{M}}\big{(}\left\lVert\bm{A}\bm{x}\right\rVert_{p}/\left\lVert\bm{x}\right\rVert_{p}\big{)}$.
Noteworthy examples are: $\left\lVert\bm{A}\right\rVert_{1}=\max_{1\leq j\leq
N}\sum_{i=1}^{M}\left\lvert a_{ij}\right\rvert$, the spectral norm
$\left\lVert\bm{A}\right\rVert_{2}=\big{[}\lambda_{\max}(\bm{A}^{\prime}\bm{A})\big{]}^{1/2}$
where $\lambda_{\max}(\cdot)$ stands for the maximum eigenvalue, and
$\left\lVert\bm{A}\right\rVert_{\infty}=\max_{1\leq i\leq
M}\sum_{j=1}^{N}\left\lvert a_{ij}\right\rvert$. The Frobenius norm of
$\bm{A}$ is
$\left\lVert\bm{A}\right\rVert_{F}=\big{(}\sum_{i=1}^{M}\sum_{j=1}^{N}|a_{ij}|^{2}\big{)}^{1/2}$.
Finally, we define $\left\lVert\bm{A}\right\rVert_{\max}=\max_{i,j}\left\lvert
a_{ij}\right\rvert$ and
$\left\lVert\bm{A}\right\rVert_{\vdash}=\max\left\\{\left\lVert\bm{A}\right\rVert_{1},\left\lVert\bm{A}\right\rVert_{\infty}\right\\}$.
Let $S\subseteq\\{1,\ldots,N\\}$ denote an index set with cardinality denoted
by $\left\lvert S\right\rvert$. Then, $\bm{x}_{S}$ denotes the $\left\lvert
S\right\rvert$-dimensional vector with the elements of $\bm{x}$ indexed by
$S$, whereas $\bm{A}_{S}$ denotes the $(M\times\left\lvert
S\right\rvert)$-dimensional matrix containing the columns of $\bm{A}$ indexed
by $S$. In addition, we define $\mathscr{D}_{A}(k)=\left\\{a_{ij}\ |\
\left\lvert i-j\right\rvert=k\right\\}$ as the collections of elements lying
on (pairs of) the diagonals in the matrix $\bm{A}$. Finally, $C$ is a generic
constant that can change value from line-to-line.

##  2 The Spatio-temporal Vector Autoregression

As in the recent paper by Gao et al., (2019), we consider the spatio-temporal
vector autoregression

$\bm{y}_{t}=\bm{A}\bm{y}_{t}+\bm{B}\bm{y}_{t-1}+\bm{\epsilon}_{t},\qquad\qquad\qquad\qquad
t=1,\ldots,T,$ (1)

where $\bm{y}_{t}=(y_{1t},\ldots,y_{Nt})^{\prime}$ stacks the observations at
time $t$ over a collection of $N$ spatial units. The contemporaneous spatial
dependence between these spatial units is governed by the matrix
$\bm{A}=(a_{ij})_{i,j=1}^{N}$ with $a_{ii}=0$ for $i=1,\ldots,N$. The matrix
$\bm{B}=(b_{ij})_{i,j=1}^{N}$ incorporates dependence on past realizations.
Finally, we have the innovation vector $\bm{\epsilon}_{t}$. We impose the
following assumptions on the DGP in (1).

######  Assumption 1 (Stability).

  1. (a)

$\left\lVert\bm{A}\right\rVert_{\vdash}=\max\left\\{\left\lVert\bm{A}\right\rVert_{1},\left\lVert\bm{A}\right\rVert_{\infty}\right\\}\leq\delta_{A}<1$.

  2. (b)

$\left\lVert\bm{B}\right\rVert_{\vdash}\leq C_{B}$ and
$\frac{C_{B}}{1-\delta_{A}}<1$.

######  Remark 1.

Assumption 1 is defined in terms of $\left\lVert\cdot\right\rVert_{\vdash}$.
Since
$\left\lVert\bm{A}\right\rVert_{1}=\left\lVert\bm{A}^{\prime}\right\rVert_{\infty}\leq\left\lVert\bm{A}\right\rVert_{\vdash}$
for any matrix $\bm{A}$, the norm $\left\lVert\cdot\right\rVert_{\vdash}$ is
convenient when bounding products of matrices containing transposes.

######  Assumption 2 (Innovations).

  1. (a)

The sequence $\\{\bm{\epsilon}_{t}\\}$ is a covariance stationary, martingale
difference process with respect to the filtration
$\mathscr{F}_{t-1}=\sigma\left(\bm{\epsilon}_{t-1},\bm{\epsilon}_{t-2},\ldots\right)$,
and geometrically strong mixing ($\alpha$-mixing). That is, the mixing
coefficients $\\{\alpha_{m}\\}$ satisfy $\alpha_{m}\leq
c_{2}e^{-\gamma_{\alpha}m}$ for all $m$ and some constants
$c_{2},\gamma_{\alpha}>0$. The largest and smallest eigenvalues of
$\bm{\varSigma}_{\epsilon}=\operatorname{\mathbb{E}}(\bm{\epsilon}_{1}\bm{\epsilon}_{1}^{\prime})=(\sigma_{ij})_{i,j=1}^{N}$
are bounded away from $0$ and $\infty$, uniformly in $T\in\mathbb{N}$.

  2. (b)

Either one of the following assumptions holds:

    1. (b1)

For $\psi(x)=x^{d}$, we require
$\sup_{i,t}\left\lVert\epsilon_{it}\right\rVert_{\psi}=\left(\operatorname{\mathbb{E}}|\epsilon_{it}|^{d}\right)^{1/d}\leq\mu_{d}<\infty$
for $d\geq 4$.

    2. (b2)

For $\psi(x)=\exp(x)-1$, we have
$\sup_{i,t}\left\lVert\epsilon_{it}\right\rVert_{\psi}\leq\mu_{\infty}<\infty$.

Assumption 1 ensures that
$\bm{y}_{t}=\bm{A}\bm{y}_{t}+\bm{B}\bm{y}_{t-1}+\bm{\epsilon}_{t}$ has a
stable reduced form VAR(1) specification. This follows from the following two
observations. First, Assumption 2(a) bounds the maximum row and column sums of
$\bm{A}$ and thereby constraints the contemporaneous dependence between the
time series. This assumption reminds of the spatial econometrics literature in
which the spatial parameter $\lambda$ is bounded from above and the
prespecified spatial weight matrix $\bm{W}_{N}$ is standardized (see, e.g.
Lee, (2004) and Lee and Yu, (2010)). Typically, the product
$\lambda\bm{W}_{N}$, the natural counter part of the matrix $\bm{A}$, is
required to fulfil conditions similar to
$\left\lVert\bm{A}\right\rVert_{\vdash}\leq\delta_{A}<1$.111For instance, it
is not uncommon to row-normalize $\bm{W}_{N}$ (each absolute row sum equal to
1) and restrict $\lambda<1$, see pages 1903-1904 of Lee, (2004). If
$\bm{W}_{N}$ is symmetric, then also
$\left\lVert\lambda\bm{W}_{N}\right\rVert_{\vdash}<1$. Invertibility of
$\bm{I}_{N}-\bm{A}$ is also guaranteed because
$\left\lVert\bm{A}\right\rVert_{2}\leq\sqrt{\left\lVert\bm{A}\right\rVert_{1}\,\left\lVert\bm{A}\right\rVert_{\infty}}\leq\delta_{A}\leq
1$ and we have the reduced-form representation
$\bm{y}_{t}=\bm{C}\bm{y}_{t-1}+\bm{D}\bm{\epsilon}_{t}$ with
$\bm{C}=(\bm{I}_{N}-\bm{A})^{-1}\bm{B}$ and $\bm{D}=(\bm{I}_{N}-\bm{A})^{-1}$.
From
$\left\lVert\bm{D}\right\rVert_{\vdash}\leq\sum_{j=0}^{\infty}\left\lVert\bm{A}\right\rVert_{\vdash}^{j}=\frac{1}{1-\delta_{A}}$,
we infer that the absolute row and column sum of $\bm{I}_{N}-\bm{A}$ are
bounded. The latter is the logical counterpart of assumption B2 in Dou et al.,
(2016). Second, Assumption 2(b) controls serial dependence. Indeed, we
conclude from
$\left\lVert\bm{C}\right\rVert_{2}\leq\left\lVert\bm{C}\right\rVert_{\vdash}\leq\frac{C_{B}}{1-\delta_{A}}<1$
that both unit root and explosive behaviour of the reduced form specification
are ruled out. The resulting stable VAR(1) representation is convenient to
study the theoretical properties of our penalized estimator.

The assumptions on the innovation process $\\{\bm{\epsilon}_{t}\\}$,
Assumption 2, are closely related to those in Masini et al., (2019).
Assumption 2(a) places restrictions on the time series properties of the error
term through martingale difference (m.d.) and mixing assumptions. The m.d.
assumption implies that
$\operatorname{\mathbb{E}}(\bm{\epsilon}_{t}\bm{y}_{t-j}^{\prime})=\bm{0}$,
while the mixing assumption controls the serial correlation in the data.
Polynomial or exponential tail decay of the distribution of the innovations
are imposed through Assumptions 2(b1) and 2(b2), respectively. The type of
tail decay will directly influence the growth rates we can allow for $N$ and
$T$. The discussions on pages 4–8 of Masini et al., (2019) demonstrate that
Assumption 2 allows for a wide range of innovation models.

Any further structure being absent, there are $(2N-1)N$ unknown parameters in
$\bm{A}$ and $\bm{B}$ to estimate. Three complications are encountered when
estimating these parameters. First, if $\bm{A}\neq\mathbf{O}$, then
$\bm{y}_{t}$ occurs on both sides of the equation, and we face an endogeneity
problem which renders the ordinary least squares estimator inconsistent.
Second, the number of unknown parameters grows quadratically in the cross-
sectional dimension $N$. The model thus quickly becomes too large to estimate
accurately without regularization. Finally, the multitude of parameters raises
concerns about identifiability. These three complications are addressed by:
(1) imposing structure on the matrices $\bm{A}$ and $\bm{B}$, and (2)
estimating the unknown coefficients using the Yule-Walker equations (e.g.
Brockwell and Davis,, 1991, p. 420).

There are several possibilities to introduce structure into $\bm{A}$ and/or
$\bm{B}$. Early spatial econometrics models, e.g. the spatial autoregressive
(SAR) model or spatial Durbin model (SDM), incorporate spatial effects through
the product $\lambda\bm{W}_{N}$ (with $\bm{W}_{N}$ pre-specified). The
specification $\bm{A}=\lambda\bm{W}_{N}$ imposes substantial structure on
$\bm{A}$ and leaves only the single parameter $\lambda$ to estimate. Dou et
al., (2016) consider a more general setting in which each row of $\bm{W}_{N}$
receives its own scale parameter. Specifically, they set
$\bm{A}=\operatorname{\text{diag}}(\bm{\lambda}_{0})\bm{W}_{N}$ and
$\bm{B}=\operatorname{\text{diag}}(\bm{\lambda}_{1})+\operatorname{\text{diag}}(\bm{\lambda}_{2})\bm{W}_{N}$,
and estimate the $3N$ coefficients in
$(\bm{\lambda}_{0}^{\prime},\bm{\lambda}_{1}^{\prime},\bm{\lambda}_{2}^{\prime})^{\prime}$.
Gao et al., (2019) require $\bm{A}$ and $\bm{B}$ to be banded matrices.222We
say that the matrix $\bm{A}$ has bandwidth $k$ if the total number of nonzero
entries in the any row or column is at most $k$. In other words, if the matrix
$\bm{A}$ satisfies $a_{ij}=0$ for all $|i-j|>k_{0}$, then $\bm{A}$ has a
bandwidth of $2k_{0}+1$ (at most). We employ a similar assumption.

######  Assumption 3 (Banded matrices).

Recall $\bm{A}=(a_{ij})_{i,j=1}^{N}$, $\bm{B}=(b_{ij})_{i,j=1}^{N}$, and
$\bm{\varSigma}_{\epsilon}=(\sigma_{ij})_{i,j=1}^{N}$. We have: (a)
$a_{ij}=b_{ij}=0$ for all $\left\lvert i-j\right\rvert>k$ with $k<\lfloor
N/4\rfloor$, and (b) $\sigma_{ij}=0$ for all $\left\lvert
i-j\right\rvert>l_{0}$.

Assumption 3 serves two purposes. First, for each spatial unit $i=1,\ldots,N$,
the matrices $\bm{A}$ and $\bm{B}$ are banded to have no more than $N$ unknown
parameters per equation. With $N$ moment conditions for each $i$, Assumption
3(a) is key in identifying the parameters. Our discussions in Section 3
illustrate that this assumption is realistic when the data is observed on a
regular grid. The combination of Assumptions 3(a)–(b) is exploited in the
Yule-Walker estimation approach. This approach requires estimation of the
$(N\times N)$ autocovariance matrices
$\bm{\varSigma}_{j}=\operatorname{\mathbb{E}}(\bm{y}_{t}\bm{y}_{t-j}^{\prime})$.
Especially in our large $N$ settings, it is crucial to rely on covariance
matrix estimators that converge at a fast rate. If $\bm{A}$, $\bm{B}$, and
$\bm{\varSigma}_{\epsilon}$ are banded, then the following result applies.

######  Theorem 1 (Convergence rates for banded sample autocovariance
matrices).

For any matrix $\bm{M}=(m_{ij})$, its $h$-banded counterpart is defined as
$\mathscr{B}_{h}\big{(}\bm{M}\big{)}=(m_{ij}\mathbbm{1}_{\\{|i-j|\leq h\\}})$.
Define the $(N\times 2N)$ matrix
$\hat{\bm{V}}_{h}=\big{[}\mathscr{B}_{h}\big{(}\hat{\bm{\varSigma}}_{1}\big{)}^{\prime}\;\mathscr{B}_{h}\big{(}\hat{\bm{\varSigma}}_{0}\big{)}\big{]}$
with
$\hat{\bm{\varSigma}}_{1}=\frac{1}{T}\sum_{t=2}^{T}\bm{y}_{t}\bm{y}_{t-1}^{\prime}$
and
$\hat{\bm{\varSigma}}_{0}=\frac{1}{T}\sum_{t=2}^{T}\bm{y}_{t}\bm{y}_{t}^{\prime}$,
and choose

$\displaystyle
h=h(\epsilon,N)=\left(\frac{\log(C_{4}N/\epsilon)}{|\log(\delta_{C})|}+1\right)$
$\displaystyle\left(3\max\Big{\\{}s^{*},\frac{\log\left(C_{4}N/(1-\delta_{C})\epsilon\right)}{|\log(\delta_{A})|}\Big{\\}}+1\right)k_{0}$
(2)
$\displaystyle\qquad\qquad+6\max\Big{\\{}s^{*},\frac{\log\left(C_{4}N/(1-\delta_{C})\epsilon\right)}{|\log(\delta_{A})|}\Big{\\}}k_{0}+l_{0},$

then $\left\lVert\widehat{\bm{V}}_{h}-\bm{V}\right\rVert_{\vdash}\leq
6\epsilon$ with a probability of at least

  1. (a)

$1-2\mathscr{P}_{1}(\epsilon,N,T)$ under Assumptions 1–3 using Assumption
2(b1) (polynomial tails),

  2. (b)

$1-2\mathscr{P}_{2}(\epsilon,N,T)$ under Assumptions 1–3 using Assumption
2(b2) (exponential tails),

where

$\mathscr{P}_{1}(\epsilon,N,T)=N^{2}\left[\left(b_{1}T^{(1-\delta)/3}+\frac{[2h(\epsilon,N)+1]b_{3}}{\epsilon}\right)\exp\left(-\frac{T^{(1-\delta)/3}}{2b_{1}^{2}}\right)+\frac{b_{2}[2h(\epsilon,N)+1]^{d}}{\epsilon^{d}T^{\frac{\delta}{2}(d-1)}}\right],$

for some $0<\delta<1$, and

$\mathscr{P}_{2}(\epsilon,N,T)=N^{2}\left[\frac{\beta_{1}[2h(\epsilon,N)+1]}{\epsilon}+\frac{2}{\beta_{2}}\left(\frac{T\epsilon^{2}}{[2h(\epsilon,N)+1]^{2}}\right)^{\frac{1}{7}}\right]\exp\left(-\frac{1}{\beta_{3}}\left(\frac{T\epsilon^{2}}{[2h(\epsilon,N)+1]^{2}}\right)^{\frac{1}{7}}\right).$

All constants ($C_{4}$, $\delta_{C}$, $s^{*}$, etc.) are positive and
independent of $N$ and $T$. The proof (see the Appendix) shows how these
constants are related to quantities in Assumptions 1–3.

Theorem 1 shows that banded estimators for $\bm{\varSigma}_{0}$ and
$\bm{\varSigma}_{1}$ can provide an accurate approximation to
$\bm{V}=\big{[}\bm{\varSigma}_{1}^{\prime}\;\bm{\varSigma}_{0}\big{]}^{\prime}$.
Each of these banded matrices has at most $2h(\epsilon,N)+1$ nonzero elements
in their columns/rows. For fixed $\epsilon$, the term
$2h(\epsilon,N)+1=O\left([\log(N/\epsilon)]^{2}\right)$ increases at a rate
much slower than $N^{2}$. This will eventually lead to better convergence
rates for our SPLASH estimator.

##  3 Sparse Estimation

###  3.1 The SPLASH($\alpha$,$\lambda$) Estimator

The matrices $\bm{A}$ and $\bm{B}$ in (1) contain a total of $2N(N-1)$
parameters. For large $N$, the accurate estimation of all these parameters
becomes infeasible rather quickly. To alleviate this curse of dimensionality,
we rely on sparsity. Some sparsity has been imposed already through Assumption
3 which was made for identifiability purposes. But this only partially solves
the issue, since the number of parameters continues to grow quadratically in
$N$. Further dimension reductions are possible whenever the spatial units are
ordered in a structured way.

Figure 1: The left figures show $(5\times 5)$ grids of spatial units with
arrows depicting horizontal (blue), vertical (red), and/or diagonal
interactions (grey). The right figures illustrate the sparsity pattern
associated with these spatial layouts. For example, spatial unit $y_{1}$
interacts horizontally with $y_{2}$, vertically with $y_{6}$, and diagonally
with $y_{7}$. If these are the only possible interactions for spatial unit
$y_{1}$, then only $a_{12}$, $a_{16}$, and $a_{17}$ in the first row of
$\bm{A}$ are possibly nonzero (the interactions for $y_{14}$ are indicated as
well). The right figures show the implied corresponding sparsity pattern for
all spatial units.

As an illustrative example, let us consider repeated measurements on the
$(5\times 5)$ spatial grids shown in the left column of Figure 1. The $N=25$
spatial units are labelled $y_{1}$ up to $y_{25}$ and enumerated row-wise.
This ordering of the spatial entities creates an implicit notion of proximity
and we intuitively expect economic and/or physical interactions to be most
pronounced at short length scales. In Figure 1(a) we depart from the situation
in which the spatial units are restricted to communicate horizontally. Blue
arrows indicate explicitly that $y_{1}$ interacts with $y_{2}$, and $y_{14}$
interacts with both $y_{13}$ and $y_{15}$. Of course, such interactions will
occur among all elements in the grid. More importantly, if only these
horizontal interaction exist, then the $(25\times 25)$ matrices $\bm{A}$ and
$\bm{B}$ feature a sparsity pattern as shown in Figure 1(b). The blue elements
are potentially nonzero whereas uncolored elements are zero. The nonzero
elements in $\bm{A}$ and $\bm{B}$ are thus clustered in specific, dense
diagonals with the occasional zero when horizontal neighbours are absent (i.e.
on the boundary of the grid). The diagonal nature of the sparsity is not an
artifact of allowing horizontal interactions only. Figures 1(c) adds the
vertical interactions and the accompanying sparsity pattern again manifests
itself along diagonals (Figure 1(d)). Finally, with diagonal nearest
neighbours being horizontal neighbors of vertical elements, we observe a
“thickening” of the diagonals in Figure 1(f). Guided by these considerations
we combine generalized Yule-Walker estimation with a sparse group penalty
(e.g. Simon et al.,, 2013). The Yule-Walker estimator will control for
endogeneity. The sparse group penalty will shrink towards diagonal structures
by including/omitting complete diagonals and thus selecting the required
interactions. Compared to Gao et al., (2019), we hereby gain the ability to
exploit sparsity within banded matrices.

A formal definition of our estimator requires further notation. Part of this
notation comes naturally if we briefly review the generalized Yule-Walker
estimator. After post-multiplying by $\bm{y}_{t-1}^{\prime}$ and taking
expectations, we find
$\bm{\varSigma}_{1}=\bm{A}\bm{\varSigma}_{1}+\bm{B}\bm{\varSigma}_{0}$ or,
equivalently,

$\bm{\varSigma}_{1}^{\prime}=\begin{bmatrix}\bm{\varSigma}_{1}^{\prime}&\bm{\varSigma}_{0}\end{bmatrix}\begin{bmatrix}\bm{A}&\bm{B}\end{bmatrix}^{\prime}=:\bm{V}\bm{C}^{\prime}.$
(3)

The $i$th column of $\bm{C}^{\prime}$ contains all coefficients that belong to
the $i$th equation in (1). Assumption 3 requires several of these coefficients
to be zero so we exclude these from the outset. We collect all remaining
(possibly) nonzero coefficients in the $i$th equation in the vector
$\bm{c}_{i}$, and define $\bm{V}_{i}$ as the matrix containing the
corresponding columns from $\bm{V}$. In the population, we have
$\bm{V}_{i}\bm{c}_{i}=\bm{\varSigma}_{1}^{\prime}\bm{e}_{i}=:\bm{\sigma}_{i}$
for $i=1,\ldots,N$. Sample counterparts of $\bm{V}_{i}$ and $\bm{\sigma}_{i}$
are readily available from the sample autocovariance matrices. More
explicitly, Gao et al., (2019) set
$\hat{\bm{\sigma}}_{i}=\frac{1}{T}\sum_{t=2}^{T}\bm{y}_{t-1}y_{it}$ and
construct $\hat{\bm{V}}_{i}$ from the appropriate columns of
$\hat{\bm{V}}=\big{[}\hat{\bm{\varSigma}}_{1}^{\prime}\;\hat{\bm{\varSigma}}_{0}\big{]}$.
Motivated by $\bm{\sigma}_{i}-\bm{V}_{i}\bm{c}_{i}=\bm{0}$, they define the
estimator $\hat{\bm{c}}_{i}$ as the following minimizer:

$\hat{\bm{c}}_{i}=\operatorname*{arg~{}min}_{\bm{c}}\left\lVert\hat{\bm{\sigma}}_{i}-\hat{\bm{V}}_{i}\bm{c}\right\rVert_{2}^{2}.$
(4)

$3$$1$$\bm{A}$786$\bm{B}$(a)13678$\hat{\bm{V}}_{h}=\left[\mathscr{B}_{h}\big{(}\hat{\bm{\varSigma}}_{1}\big{)}^{\prime}\;\mathscr{B}_{h}\big{(}\hat{\bm{\varSigma}}_{0}\big{)}\right]$(b)13678$\hat{\bm{V}}_{2,h}$(c)
Figure 2: A visualization on the construction of $\hat{\bm{V}}_{2,h}$ for
$N=5$. (a) If $h=1$, then grey elements in $\bm{A}$ and $\bm{B}$ are
(potentially) nonzero whereas white elements are zero by construction.
Enumerating along the second row, the active elements are in the set
$\\{1,3,6,7,8\\}$. (b) We select the columns from $\hat{\bm{V}}_{h}$
corresponding to the active set. (c) The matrix $\hat{\bm{V}}_{2,h}$ is the
submatrix of $\hat{\bm{V}}_{h}$ with only active columns.

We will adjust this objective function in three ways. First, we define our
estimator in terms of banded estimated covariance matrices, which allows us to
exploit the results in Theorem 1. Second, our group penalty penalizes
parameters across equations so we can no longer estimate the parameters
equation-by-equations. We therefore define
$\hat{\bm{\sigma}}_{h}=\operatorname{\mathrm{vec}}\big{(}\mathscr{B}_{h}\big{(}\hat{\bm{\varSigma}}_{1}\big{)}^{\prime}\big{)}$
and
$\hat{\bm{V}}_{h}^{(d)}=\operatorname{\text{diag}}\big{(}\hat{\bm{V}}_{i,h},\ldots,\hat{\bm{V}}_{N,h}\big{)}$,
with $\hat{\bm{V}}_{i,h}$ being constructed similarly to $\bm{V}_{i}$ (see
Figure 2 for an illustration).333In general, all quantities derived from the
autocovariance matrices come in three versions: (1) the population quantity,
(2) the estimated counterpart without banding denoted with an additional
“hat”, and (3) the estimated counterpart with banding featuring both a “hat”
and the subscript $h$. In this notation, the expression
$\big{\|}\hat{\bm{\sigma}}_{h}-\hat{\bm{V}}_{h}^{(d)}\bm{c}\big{\|}_{2}^{2}$
defines the joint objective function that sums the individual contributions in
(4) over all equations. Finally, we construct the penalty function. We define
an index set that partitions the vector $\bm{c}$ into sub-vectors, denoted
$\\{\bm{c}_{g}\\}$, that contain the non-zero diagonals of $\bm{A}$ and
$\bm{B}$ that are admissible under Assumption 3 as

$\begin{split}\mathscr{G}_{A}&:=\left\\{g\subset\mathbb{N}:\bm{c}_{g}=\mathscr{D}_{A}(k),k\in\\{1,\ldots,\lfloor
N/4\rfloor\right\\},\\\
\mathscr{G}_{B}&:=\left\\{g\subset\mathbb{N}:\bm{c}_{g}=\mathscr{D}_{B}(k),k\in\\{0,\ldots,\lfloor
N/4\rfloor\right\\},\end{split}$ (5)

respectively, and let $\mathscr{G}=\mathscr{G}_{A}\cup\mathscr{G}_{B}$. Based
on this notation, we define our objective function as

$\mathscr{L}_{\alpha}(\bm{c};\lambda)=\left\lVert\hat{\bm{\sigma}}_{h}-\hat{\bm{V}}_{h}^{(d)}\bm{c}\right\rVert_{2}^{2}+\lambda\Bigg{(}\underbrace{(1-\alpha)\sum_{g\in\mathscr{G}}\sqrt{\left\lvert
g\right\rvert}\left\lVert\bm{c}_{g}\right\rVert_{2}+\alpha\left\lVert\bm{c}\right\rVert_{1}}_{=:P_{\alpha}(\bm{c})}\Bigg{)}=\left\lVert\hat{\bm{\sigma}}_{h}-\hat{\bm{V}}_{h}^{(d)}\bm{c}\right\rVert_{2}^{2}+\lambda
P_{\alpha}(\bm{c}).$ (6)

The spatial lasso-type shrinkage estimator, abbreviated
SPLASH($\alpha$,$\lambda$) or SPLASH in short, is defined as the minimizer of
(6), i.e.
$\hat{\bm{c}}=\operatorname*{arg~{}min}_{\bm{c}}\mathscr{L}_{\alpha}(\bm{c};\lambda)$.
The importance of the penalty function $P_{\alpha}(\bm{c})$ is governed by the
penalty parameter $\lambda$ and a second hyperparameter $\alpha$ balances
group-structured sparsity versus individual sparsity. At the extremities of
$\alpha\in[0,1]$ we find the group lasso ($\alpha=0$) and the lasso
($\alpha=1$). Intermediate values of $\alpha$ will shrink both groups of
diagonal coefficients in $\bm{A}$ and $\bm{B}$ and individual parameters. The
SPLASH solution promotes completely sparse diagonals and sparse elements
within nonzero diagonals, and thus shrinks towards sparsity patterns of the
type displayed in Figure 1(b). As the structure of our estimator is similar to
that of the sparse group lasso (SGL), efficient algorithms are available to
compute its solution (see, e.g. Simon et al.,, 2013).444An R/C++
implementation of the SPLASH estimator based on this algorithm is available on
https://sites.google.com/view/etiennewijler/code.

###  3.2 Theoretical Properties of the SPLASH($\alpha$,$\lambda$) Estimator

In this section we derive the theoretical properties of the SPLASH estimator.
First, however, we require an additional assumption on the DGP in order to
ensure that $\bm{A}$ and $\bm{B}$ in (1) are uniquely identified. To this end,
we leverage the bandedness assumption in Assumption 3, which enables unique
identification of $\bm{A}$ and $\bm{B}$ via a straightforward full-rank
condition on sub-matrices of the autocovariance matrices that appear in the
generalized Yule-Walker equations:

######  Assumption 4 (Restricted minimum eigenvalue).

Assume that

$\phi_{\min}(\bm{x}):=\underset{\bm{x}\in\mathbb{R}^{2N}:\mathscr{M}(\bm{x})\leq
N}{\text{min}}\frac{\left\lVert\bm{V}\bm{x}\right\rVert_{2}}{\left\lVert\bm{x}\right\rVert_{2}}\geq\phi_{0}.$

Essentially, Assumption 4 states that every sub-matrix containing $N$ columns
from $\bm{V}$ has full column-rank and a minimum singular value bounded away
from zero. Related assumptions appear in Bickel et al., (2009, Section 4), who
refer to $\phi_{\min}(\bm{x})$ as a restricted eigenvalue and use this
quantity to construct sufficient conditions for their restricted eigenvalue
assumptions. Assumption 4 fits our framework particularly well, as the assumed
maximum bandwidth of the matrices $\bm{A}$ and $\bm{B}$ in Assumption 3 imply
that the diagonal blocks of the matrix $\bm{V}^{(d)}$ never contain more than
$N$ unique columns of $\bm{V}$. Using this property, we show in Lemma 1 in
Appendix A that the SGL compatibility condition is implied by Assumption 4.

Equipped with Assumption 4, we are now able to derive the following finite-
sample performance bounds on the prediction and estimation error of SPLASH.

######  Theorem 2.

Under Assumptions 1-4, it holds that

$\left\lVert\hat{\bm{V}}_{h}^{(d)}(\hat{\bm{c}}-\bm{c})\right\rVert_{2}^{2}+\lambda\left((1-\alpha)\sum_{g\in\mathscr{G}}\sqrt{\left\lvert
g\right\rvert}\left\lVert\hat{\bm{c}}_{g}-\bm{c}_{g}\right\rVert_{2}+\alpha\left\lVert\hat{\bm{c}}-\bm{c}\right\rVert_{1}\right)\leq\frac{64\bar{\omega}_{\alpha}^{2}\lambda^{2}}{\phi_{0}^{2}}$

with a probability of at least

  1. (a)

$1-10\mathscr{P}_{1}\left(f(\lambda,\phi_{0}),N,T\right)$ when Assumption
2(b1) (polynomial tail decay) is valid, or

  2. (b)

$1-10\mathscr{P}_{2}\left(f(\lambda,\phi_{0}),N,T\right)$ when Assumption
2(b2) (exponential tail decay) is valid,

where $\mathscr{P}_{1}\left(x,N,T\right)$ and
$\mathscr{P}_{2}\left(x,N,T\right)$ are defined in Theorem 1 and
$f(\lambda,\phi_{0})=\min\left(\frac{\lambda^{1/2}}{24},\frac{\lambda}{96C_{v}},\frac{\phi_{0}}{12}\right)$.

Theorem 2 contains a finite-sample performance bound for the prediction and
estimation error for the SPLASH($\alpha$,$\lambda$) estimator. The result
demonstrates that both errors are increasing in $\bar{\omega}_{\alpha}$, which
in turn is increasing in the bandwidths $k_{0}$ and $l_{0}$ and the group
sizes when $\alpha<1$ and increasing in the number of relevant interactions
$\left\lvert S\right\rvert$ when $\alpha>0$. Furthermore, the prediction and
estimation error increase along with the imposed amount of penalization
$\lambda$. While this seems to suggest minimizing the amount of penalization
to obtain better performance bounds, we emphasize that the effect of
regularization in Theorem 2 is two-fold: increasing the regularization
($\lambda$) worsens the performance bound, but increases the probability of
the set on which the performance bound holds. Intuitively, the imposed
shrinkage induces finite-sample bias in the estimates which worsen the
prediction and estimation accuracy. At the same time, however, shrinkage
reduces sensitivity to noise, such that performance guarantees can be provided
with more certainty.

The aforementioned interactions can be intuitively demonstrated by means of an
asymptotic analysis. Accordingly, based on Theorem 2, we derive the conditions
for convergence of the prediction and estimation errors, as well as their
convergence rates, in the following corollary.

######  Corollary 1.

Let $\lambda\in O\left(T^{-q_{\lambda}}\right)$, $N\in
O\left(T^{q_{N}}\right)$, $\left\lvert\mathscr{G}_{S}\right\rvert\in
O\left(T^{q_{g}}\right)$, $\left\lvert S\right\rvert\in
O\left(T^{q_{s}}\right)$, $k_{0},l_{0}\in O\left(T^{q_{k}}\right)$, where
$q_{\lambda}$, $q_{N}$, $q_{s}$, and $q_{k}$ are fixed and positive constants.
Maintain Assumptions 1-4 and assume that either (i)
$q_{\lambda}<-\frac{2q_{N}}{d}-q_{k}+\frac{\delta(d-1)}{2d}$ for some
$0<\delta<1$ and Assumption 2(b1) holds, or (ii)
$q_{\lambda}<\frac{1}{2}-q_{k}$ and Assumption 2(b2) holds. Then,

  1. (a)

$\left\lVert\hat{\bm{V}}_{h}^{(d)}\left(\hat{\bm{c}}-\bm{c}\right)\right\rVert_{2}^{2}=O_{p}\left((1-\alpha)T^{2q_{g}+q_{N}-2q_{\lambda}}+\alpha
T^{q_{s}-2q_{\lambda}}\right)$,

  2. (b)

$(1-\alpha)\sum_{g\in\mathscr{G}}\sqrt{\left\lvert
g\right\rvert}\left\lVert\hat{\bm{c}}_{g}-\bm{c}_{g}\right\rVert_{2}+\alpha\left\lVert\hat{\bm{c}}-\bm{c}\right\rVert_{1}=O_{p}\left((1-\alpha)T^{2q_{g}+q_{N}-q_{\lambda}}+\alpha
T^{q_{s}-q_{\lambda}}\right)$.

Corollary 1 provides detailed insights into the determinants of the
convergence rate. In particular, the result confirms that the convergence rate
decreases in the bandwidths $k_{0}$ and $l_{0}$, the number of spatial units
$N$, the number of interactions $\left\lvert S\right\rvert$ and the degree of
penalization $\lambda$.555Recall that $\lambda\in
O\left(T^{-q_{\lambda}}\right)$, such that a higher $q_{\lambda}$ implies a
faster decay of the penalty term. To ensure that the set on which the
performance bound in Theorem 2 holds occurs with probability converging to
one, conditions (i) and (ii) impose that the degree of penalization does not
decay too fast. The optimal convergence rate is obtained by choosing
$q_{\lambda}$ as large as possible without violating these conditions. Several
examples are provided in Remark 2.

######  Remark 2.

Several insightful cases can be examined based on Corollary 1. For the sake of
brevity, we consider the estimation error, as measured by
$P_{\alpha}\left(\hat{\bm{c}}-\bm{c}\right)$, for the case in which the errors
possess at least $d$ finite moments (Assumption 2(b1)) in some special cases.
First, in the absence of within-group shrinkage (i.e. $\alpha=0$), Corollary 1
demonstrates that
$P_{0}\left(\hat{\bm{c}}-\bm{c}\right)=O_{p}\left(T^{q_{N}-q_{\lambda}}\right)$,
with $q_{\lambda}<\frac{1}{2}-\frac{2q_{N}}{d}-\frac{1}{2d}$ by choosing
$\delta$ arbitrarily close to 1. Hence, the estimator converges almost at the
rate $\frac{T^{1/2-1/2d}}{N^{1+2/d}\left\lvert\mathscr{G}_{S}\right\rvert}$.
For fixed $N$ and large $d$, this would indeed be close to the optimal
$\sqrt{T}$ convergence familiar from the fixed-dimensional case without
regularization. Alternatively, when the shrinkage is imposed at the individual
interaction level only ($\alpha=1)$, then
$P_{1}\left(\hat{\bm{c}}-\bm{c}\right)=O_{p}\left(T^{q_{s}-q_{\lambda}}\right)$
such that the estimation error converges almost at the rate
$\frac{T^{1/2-1/2d}}{\left\lvert S\right\rvert N^{2/d}}$, instead.

##  4 Simulations

###  4.1 Simulation setting

In this section, we explore the finite sample performance of our estimator by
Monte Carlo simulation. The data generating process underlying the simulations
is the spatio-temporal VAR in (1). We study $T\in\\{500,1000,2000\\}$ and draw
all errors $\epsilon_{it}$ independently and $N(0,1)$ distributed. The
matrices $\bm{A}$ and $\bm{B}$ and the cross-sectional dimension $N$ are
specified in the two designs below. All simulation results are based on
$N_{sim}=500$ Monte Carlo replications.

Design A (Banded specification): We revisit simulation Case 1 in Gao et al.,
(2019). The matrices $\bm{A}$ and $\bm{B}$ are banded with a bandwidth of
$k_{0}=3$. Specifically, the elements in the matrices $(\bm{A})_{i,j=1}^{N}$
and $(\bm{B})_{i,j}^{N}$ are generated according to the following two steps:

  1. Step 1:

If $|i-j|=k_{0}$, then $a_{ij}$ and $b_{ij}$ are drawn independently from a
uniform distributions on the two points $\\{-2,2\\}$. All remaining elements
within the bandwidth are drawn from the mixture distribution $\omega
I_{\\{0\\}}+(1-\omega)N(0,1)$ with $\operatorname{\mathbb{P}}(\omega=1)=0.4$
and $\operatorname{\mathbb{P}}(\omega=0)=0.6$.

  2. Step 2:

Rescale the matrices $\bm{A}$ and $\bm{B}$ from Step 1 to
$\eta_{1}\times\bm{A}/\left\lVert\bm{A}\right\rVert_{2}$ and
$\eta_{2}\times\bm{B}/\left\lVert\bm{B}\right\rVert_{2}$, where $\eta_{1}$ and
$\eta_{2}$ are drawn independently from $U[0.4,0.8]$.666This rescaling does
not necessarily imply
$\left\lVert(\bm{I}_{N}-\bm{A})^{-1}\bm{B}\right\rVert_{2}<1$ (stability).
During the simulations we redraw the matrices $\bm{A}$ and $\bm{B}$ whenever
$\left\lVert(\bm{I}_{N}-\bm{A})^{-1}\bm{B}\right\rVert_{2}>0.95$.

We vary the cross-sectional dimension over $N\in\\{25,50\\}$.

Design B (Spatial grid with neighbor interactions): As in Figure 1, we
consider an $(m\times m)$ grid of spatial units. For $m=5$ ($m=7$),this
results in a cross-sectional dimension of $N=25$ ($N=49$). The matrix $\bm{A}$
contains interactions between first horizontal and first vertical neighbours
while all other coefficients are zero. The magnitude of these nonzero
interactions are 0.2. For $m=5$ ($m=7$), the temporal matrix $\bm{B}$ is a
diagonal matrix with elements 0.25 (0.23) on the diagonal. The reduced form
VAR matrix $\bm{C}=(\bm{I}_{N}-\bm{A})^{-1}\bm{B}$ has a maximum eigenvalue of
0.814 (0.882).

For each design, we report simulation results for three sets of estimators.
The first set includes the estimators developed in this paper: (1) the
SPLASH(0,$\lambda$) estimator promotes non-sparse groups only, (2) the
SPLASH(1,$\lambda$) encourages unstructured sparsity, and (3) the
SPLASH($\alpha$,$\lambda$) balances group formation and sparse solutions by
selecting $\alpha$ in a data-driven way. In congruence with Theorems 1 and 2,
we rely on banded autocovariance matrices
$\mathscr{B}_{h}\big{(}\hat{\bm{\varSigma}}_{0}\big{)}$ and
$\mathscr{B}_{h}\big{(}\hat{\bm{\varSigma}}_{1}\big{)}$. The bandwidth choice
is determined by the bootstrap procedure described in Guo et al., (2016, p.
7). Second, we include two unpenalized estimators in the spirit of Gao et al.,
(2019). The estimator GMWY$(k_{0})$ implements generalized Yule-Walker
estimation for banded $\bm{A}$ and $\bm{B}$ with the true bandwidth $k_{0}$.
GMWY($\bm{c}_{S}$) is even more informed about the data generating process and
estimates true nonzero coefficients only. To allow for easy comparison with
the simulation results by the aforementioned authors, we implement these GMWY
estimators without banding the covariance matrix estimators
$\hat{\bm{\varSigma}}_{0}=\frac{1}{T}\sum_{t=2}^{T}\bm{y}_{t}\bm{y}_{t}^{\prime}$
and
$\hat{\bm{\varSigma}}_{1}=\frac{1}{T}\sum_{t=2}^{T}\bm{y}_{t}\bm{y}_{t-1}^{\prime}$.777In
unreported simulation results, that are available upon request, we find that
the results are insensitive to this choice. Note that GMWY$(k_{0})$ and
GMWY($\bm{c}_{S}$) are infeasible in practice and are thus given a comparative
advantage.888A (consistent) ratio-based estimator for $k_{0}$ is described on
page 214 of Gao et al., (2019). A feasible two-step estimation procedures is:
(1) estimate $k_{0}$ using the ratio-based estimator (i.e. compute $\hat{k}$),
and (2) estimate the spatio-temporal VAR assuming $\bm{A}$ and $\bm{B}$ have
bandwidth $\hat{k}$. As in Gao et al., (2019), we also observed a (small)
tendency to overestimate $k_{0}$. In high-dimensional settings, this causes
highly volatile parameter estimates and performance issues. We thus decided to
provide either $k_{0}$ or the true sparsity structure. Gao et al., (2019) also
propose estimators that use multiple Yule-Walker equations or pre-select the
estimating equations that are most informative about the parameters. In their
simulation section, see Gao et al., (2019, p. 219), they state that these
estimators do not provide clear benefits. As such, we decided to not include
these estimators in our simulation study. The third set contains the
$L_{1}$-penalized reduced form VAR(1) estimator (abbreviated PVAR) only. That
is, we consider the reduced form VAR($1$) specification
$\bm{y}_{t}=\bm{C}\bm{y}_{t-1}+\bm{u}_{t}$ and estimate $\bm{C}$ by minimizing
$\mathscr{L}_{pvar}(\bm{C})=\sum_{t=2}^{T}\left\lVert\bm{y}_{t}-\bm{C}\bm{y}_{t-1}\right\rVert_{2}^{2}+\lambda\sum_{i,j=1}^{N}\left\lvert
c_{ij}\right\rvert$. This estimator is well-researched in the literature (see,
e.g. Kock and Callot,, 2015; Gelper et al.,, 2016; Masini et al.,, 2019),
albeit in different settings. It will serve as a competitive benchmark for the
forecasting performance of our proposed estimation procedure.

The forecasting performance of each estimator is assessed by the relative
mean-squared forecast error (RMSFE). Using a superscript $j$ to index the
Monte Carlo replications, the RMSFE is calculated as

$\text{RMSFE}=\frac{\sum_{j=1}^{N_{sim}}\left\lVert\bm{y}_{T+1}^{j}-\hat{\bm{C}}^{j}\bm{y}_{T}^{j}\right\rVert_{2}^{2}}{\sum_{j=1}^{N_{sim}}\left\lVert\bm{y}_{T+1}^{j}-\bm{C}\bm{y}_{T}^{j}\right\rVert_{2}^{2}}.$

As the SPLASH and GMWY procedures estimate $\bm{A}$ and $\bm{B}$, we can also
compare the estimation accuracy. Using the superscript $j$ as before, the
_Estimation Error (EE)_ of the coefficient matrices are:

$\text{EE}_{A}=\frac{1}{N_{sim}}\sum_{j=1}^{N_{sim}}\left\lVert\hat{\bm{A}}^{j}-\bm{A}\right\rVert_{2},\text{
and
}\text{EE}_{B}=\frac{1}{N_{sim}}\sum_{j=1}^{N_{sim}}\left\lVert\hat{\bm{B}}^{j}-\bm{B}\right\rVert_{2}.$
(7)

Finally, the hyperparameter, $\lambda$ (and also $\alpha$ for the
SPLASH($\alpha$,$\lambda$) estimator), for the SPLASH and PVAR estimators are
selected based on time series cross-validation. The grid for $\lambda$
consists of 20 evenly spaced values on a logarithmic scale. The maximum value
for $\lambda$ is the smallest value that gives the zero solution and the
minimum value is set to $10^{-4}$ times the maximum. For $\alpha$, the
admissible values are $\alpha\in\\{0,0.25,0.5,0.75,1\\}$. For each choice of
hyperparameters, we estimate the model on the first 80% of the data and
evaluate the out-of-sample fit on the last 20% of the data. The
hyperparameters with the lowest out-of-sample model error are chosen.

###  4.2 Simulation results

Table 1: Simulation results for Design A (Banded specification).

$N$ | $T$ | SPLASH($0$,$\lambda$) | SPLASH($\alpha$,$\lambda$) | SPLASH($1$,$\lambda$) | GMWY($k_{0}$) | GMWY($c_{S}$) | PVAR
---|---|---|---|---|---|---|---
Relative Mean-Squared Forecast Error (RMSFE)
25 | 500 | 1.013 | 1.014 | 1.023 | 1.039 | 1.018 | 1.106
| 1,000 | 1.010 | 1.011 | 1.019 | 1.016 | 1.008 | 1.108
| 2,000 | 1.004 | 1.005 | 1.010 | 1.004 | 1.002 | 1.122
50 | 500 | 1.019 | 1.019 | 1.035 | 1.033 | 1.018 | 1.108
| 1,000 | 1.010 | 1.010 | 1.020 | 1.018 | 1.010 | 1.101
| 2,000 | 1.009 | 1.009 | 1.013 | 1.008 | 1.006 | 1.101
100 | 500 | 1.022 | 1.023 | 1.040 | 1.034 | 1.021 | 1.115
| 1,000 | 1.014 | 1.014 | 1.025 | 1.018 | 1.013 | 1.106
| 2,000 | 1.009 | 1.009 | 1.015 | 1.010 | 1.007 | 1.103
Estimation Error in A (EEA)
25 | 500 | 0.549 | 0.552 | 0.836 | 1.208 | 0.689 |
| 1,000 | 0.527 | 0.527 | 0.854 | 1.095 | 0.588 |
| 2,000 | 0.489 | 0.491 | 0.854 | 1.018 | 0.491 |
50 | 500 | 0.558 | 0.557 | 0.954 | 0.899 | 0.634 |
| 1,000 | 0.529 | 0.529 | 0.968 | 0.783 | 0.537 |
| 2,000 | 0.506 | 0.504 | 1.010 | 0.702 | 0.450 |
100 | 500 | 0.555 | 0.554 | 0.928 | 0.764 | 0.593 |
| 1,000 | 0.521 | 0.519 | 0.917 | 0.689 | 0.519 |
| 2,000 | 0.502 | 0.499 | 0.976 | 0.598 | 0.440 |
Estimation Error in B (EEB)
25 | 500 | 0.282 | 0.284 | 0.386 | 0.417 | 0.244 |
| 1,000 | 0.243 | 0.244 | 0.356 | 0.342 | 0.187 |
| 2,000 | 0.207 | 0.208 | 0.330 | 0.295 | 0.146 |
50 | 500 | 0.305 | 0.306 | 0.436 | 0.368 | 0.250 |
| 1,000 | 0.257 | 0.257 | 0.394 | 0.277 | 0.186 |
| 2,000 | 0.222 | 0.220 | 0.380 | 0.218 | 0.139 |
100 | 500 | 0.323 | 0.323 | 0.449 | 0.359 | 0.263 |
| 1,000 | 0.273 | 0.272 | 0.404 | 0.267 | 0.196 |
| 2,000 | 0.234 | 0.231 | 0.383 | 0.207 | 0.147 |

  * •

Note: The relative mean-squared forecast error (RMSFE) and estimation errors
($\text{EE}_{A}$ and $\text{EE}_{B}$) of design A. Lower numbers indicate
better performance. As PVAR estimates a reduced form VAR, there are no model
errors for $\bm{A}$ and $\bm{B}$ to report for this method.

The results for Design A are reported in Table 1. First, we consider the
predictive performance. Across all specifications, the penalized VAR estimator
is outperformed by all other estimators, indicating that direct estimation of
the contemporaneous spatial interactions can deliver forecast improvements
over regularized reduced form estimation. A plausible explanation for this
finding is that the sparsity patterns in the spatial representation are less
prevalent, and thus more difficult to exploit, in the reduced form
representation. Overall, and in line with expectations, the best predictor is
the infeasible GMWY($c_{s}$) estimator, which requires knowledge of the exact
sparsity pattern. Alternatively, the best feasible estimators are
SPLASH(0,$\lambda$) and SPLASH($\alpha$,$\lambda$), who obtain nearly
equivalent performance, since the cross-validation approach decides on
$\alpha=0$ in the overwhelming majority of simulation trials. Moreover, for
smaller samples ($T=500$), the SPLASH estimators obtain a forecast accuracy is
notable higher than that obtained by GMWY($k_{0}$). However, as the latter
improves rapidly in $T$, the performance gains of SPLASH are less notable for
$T=1000$ and absent for $T=2000$.

Next, we explore the estimation accuracy. For SPLASH(0,$\lambda$) and
SPLASH($\alpha$,$\lambda$), the estimation error decreases in $T$. This is to
be expected since an increasing sample size will result in more accurate
estimates for $\bm{\varSigma}_{0}$ and $\bm{\varSigma}_{1}$ and thus (through
the Yule-Walker equations) better parameter estimates and forecasts. Somewhat
counter-intuitively, however, SPLASH(1,$\lambda$) seems to attain worse
estimation accuracy on larger samples. Since the penalty is chosen by cross-
validation, which solely focuses on out-of-sample predictive performance, we
conjecture that the most accurate estimates may not always the one that
delivers the most accurate predictive performance within a given sample.
Regarding the performance differentials across methods, the implied
performance ranking largely coincides with the ranking based on the relative
mean-squared forecast error. SPLASH(0,$\lambda$) and
SPLASH($\alpha$,$\lambda$) obtain the lowest estimation error, narrowly
followed by the infeasible GMWY($c_{s}$). The GMWY($k_{0}$) performs notably
worse than SPLASH in small samples, but closes the gap somewhat when $T$
increases.

Table 2: Simulation results for Design B (Spatial grid with neighbor
interactions).

$N$ | $T$ | SPLASH($0$,$\lambda$) | SPLASH($\alpha$,$\lambda$) | SPLASH($1$,$\lambda$) | GMWY($k_{0}$) | GMWY($c_{S}$) | PVAR
---|---|---|---|---|---|---|---
Relative Mean-Squared Forecast Error (RMSFE)
25 | 500 | 1.011 | 1.025 | 1.046 | 1,573.644 | 1.011 | 1.109
| 1,000 | 1.004 | 1.016 | 1.032 | 23.176 | 1.004 | 1.068
| 2,000 | 1.008 | 1.025 | 1.019 | 7.441 | 1.004 | 1.043
100 | 500 | 1.023 | 1.029 | 1.086 | 131.346 | 1.015 | 1.168
| 1,000 | 1.011 | 1.023 | 1.037 | 37.423 | 1.011 | 1.116
| 2,000 | 1.007 | 1.017 | 1.015 | 1.170 | 1.004 | 1.070
Estimation Error in A (EEA)
25 | 500 | 0.322 | 0.351 | 0.832 | 4.187 | 0.435 |
| 1,000 | 0.277 | 0.327 | 0.718 | 3.968 | 0.313 |
| 2,000 | 0.232 | 0.318 | 0.648 | 3.773 | 0.224 |
100 | 500 | 0.351 | 0.349 | 1.066 | 2.078 | 0.549 |
| 1,000 | 0.328 | 0.339 | 0.913 | 2.030 | 0.473 |
| 2,000 | 0.291 | 0.316 | 0.769 | 2.010 | 0.388 |
Estimation Error in B
25 | 500 | 0.102 | 0.141 | 0.338 | 1.269 | 0.137 |
| 1,000 | 0.087 | 0.137 | 0.298 | 1.098 | 0.095 |
| 2,000 | 0.072 | 0.134 | 0.248 | 0.966 | 0.069 |
100 | 500 | 0.105 | 0.114 | 0.361 | 0.906 | 0.218 |
| 1,000 | 0.095 | 0.123 | 0.324 | 0.687 | 0.160 |
| 2,000 | 0.080 | 0.121 | 0.282 | 0.553 | 0.120 |

  * •

Note: The relative mean-squared forecast error (RMSFE) and estimation errors
($\text{EE}_{A}$ and $\text{EE}_{B}$) of design B. Lower numbers indicate
better performance. As PVAR estimates a reduced form VAR, there are no model
errors for $\bm{A}$ and $\bm{B}$ to report for this method.

Simulation results for Design B are shown in Table 2.The high RMSFEs for
GMWY($k_{0}$) in comparison to all competing models are most striking.
However, do recall that this approach cannot exploit any sparsity as present
in the sparse representations for $\bm{A}$ and $\bm{B}$ in this design. For
instance, looking at Figure 3(a) we see that the entries $a_{ij}$ of $\bm{A}$
with $|i-j|\leq 5$ contain 165 zeros and 80 nonzeros for $N=25$. The matrix
$\bm{A}$ within the correct bandwidth is thus 51.52% sparse. Any unpenalized
approach is likely overparameterized which explains the bad forecasting
performance and high estimation errors that is especially pronounced in small
samples. The other results tell the same story as in Design A;
SPLASH(0,$\lambda$) and SPLASH($\alpha$,$\lambda$) are performing very close
to the optimal, though infeasible, GMWY($\bm{c}_{S})$ estimator, and forecast
notably better than the PVAR. SPLASH(1,$\lambda$) ranks in-between these
approaches, as it benefits from the direct estimation of spatial interactions,
but does so less efficiently by not taking into account the diagonally
structured form of sparsity.

A small visual analysis provides further evidence on the good estimation
accuracy obtained by SPLASH($\alpha$,$\lambda$). That is, we can visualize the
capability of recovering the correct sparsity pattern by displaying the
absolute value of the coefficients as averaged across all $N_{sim}$ simulation
runs. Figure 3 illustrates the similarity between the true matrix $\bm{A}$ and
the average magnitude of the estimated coefficients.

Figure 3: Visualizations of the true and estimated spatial weight matrix
$\bm{A}$ for Design B. (a) The true spatial weight matrix $\bm{A}$ implied by
the $(5\times 5)$ spatial grid design (Design B with $m=5$). (b) The average
absolute values of the entries in $\hat{\bm{A}}$ as computed by
SPLASH($\alpha$,$\lambda$) for $N=25$ and $T=1000$. That is, the $(i,j)$th
entry in the matrix on the right equals
$\frac{1}{N_{sim}}\sum_{k=1}^{N_{sim}}|\hat{a}_{ij}^{k}|$ with
$\hat{a}_{ij}^{k}$ being the estimated $(i,j)$th entry of $\bm{A}$ in the
$k$th Monte Carlo replication.

######  Remark 3.

In elaborate, though unreported, visual analysis, we find that most zero
diagonals are actually not estimated as exactly zero by the (sparse) group
lasso. When tuning the penalty parameter by the BIC criterion, in which the
number of estimated non-zeros is used as a proxy for the degrees of freedom,
the true zero diagonals are typically estimated as exactly zeros. However, the
increased amount of regularization that is required to effectuate this has a
negative impact on the forecast performance.

##  5 Empirical Application

We examine the empirical performance of the SPLASH estimator when predicting
daily $\text{NO}_{2}$ concentrations in (Greater) London. The data is publicly
available via the Copernicus Open Access Hub and we consider the time span
from 1 August 2018 to 18 October 2020.999See http://www.tropomi.eu/ for more
info on the TROPOMI data products, and use https://scihub.copernicus.eu/ to
access the database. The original $\text{NO}_{2}$ concentrations are reported
in mol/m2 but we convert to mol/cm2 to avoid numerical instabilities caused by
the use of small-scale numbers. The far majority of measurements are captured
sometime between 11:00 and 14:00 UTC. The area of interest is divided into a
($5\times 9$) grid, implying that longitudes and latitudes increment by
approximately 0.2 from cell to cell. All available $\text{NO}_{2}$
measurements are averaged within each cell and within the same day. The
resulting data set has 0.8% missing observations. These missing values are
imputed using the Multivariate Time Series Data Imputation (mtsdi) R
package.101010This imputation method is proposed by Junger and Ponce De Leon,
(2015) to impute missing values in time series for air pollutants. The package
is written by the same authors and currently maintained by W. L. Junger.

A rolling-window approach is used to assess the predictive power of the SPLASH
estimator. Each window contains 80% of the data, or 641 days, such that a
total of 160 one-step ahead forecasts can be produced. For each window, we
proceed along the following four steps: (i) de-mean the data, (ii) determine
the hyperparameters and estimate each model, (iii) produce a forecast for the
de-meaned data, and (iv) add the means back to the forecast. In addition to
the estimators described in the simulation section (Section 4, see page 4), we
add another forecast: the window’s mean. This new forecast is abbreviated
CONST and all other forecasts follow the notational conventions from the
simulation section. Similarly, all tuning parameters are chosen via the time-
series cross-validation procedure outlined in Section 4. The spatial grid
contains $N=5\times 9=45$ spatial units, such that the
SPLASH($\alpha$,$\lambda$) models contain $2N^{2}-N=4,005$ parameters. For the
purpose of identifiability, we band the spatial matrix $\bm{A}$ and
autoregressive matrix $\bm{B}$ such that $a_{ij}=b_{ij}=0$ for $\left\lvert
i-j\right\rvert>\lfloor N/4\rfloor=11$. By ordering the spatial units
vertically, this banding puts no restrictions on the degree of vertical
interaction but allows no more than second-order interaction between
horizontal neighbours (see Figure 5 in Appendix C.1 for details).

The forecast performance is measured along three metrics and always relative
to the $L_{1}$-penalized reduced form VAR($1$) (PVAR) benchmark. That is, we
report: (i) the number of spatial units that are predicted more accurately
than the PVAR method (#wins), (ii) the number of spatial units that are
predicted significantly more accurately based on a Diebold-Mariano test at the
the 5% significance level (#sign.), and (iii) the average loss relative to the
penalized VAR over all spatial units. These three metrics are calculated based
on two loss functions for the forecast errors, namely the mean squared
forecast error (MSFE) and the mean absolute forecast error (MAFE). We report
the MAFE because the $\text{NO}_{2}$ column densities display several abrupt
spikes which may carry to much weight when relying on a squared loss function.
The results are reported in Table 3.

Table 3: Forecast performance of various methods for $\text{NO}_{2}$ satellite
data on a $(5\times 9)$ grid of observations.

| MSFE |  | MAFE
---|---|---|---
| #wins | #sign. wins | RMSFE |  | #wins | #sign. wins | RMAFE
CONST | 0 | 0 | 1.185 |  | 0 | 0 | 1.139
GMWY | 0 | 0 | 153.122 |  | 0 | 0 | 4.851
SPLASH($0$,$\lambda$) | 45 | 42 | 0.917 |  | 45 | 43 | 0.939
SPLASH($\alpha$,$\lambda$) | 44 | 30 | 0.941 |  | 44 | 32 | 0.961
SPLASH($1$,$\lambda$) | 34 | 14 | 0.969 |  | 36 | 20 | 0.975

  * •

Note: Number of grid points (out of $N=45$) with lower prediction errors
(#wins) and significantly lower prediction errors (#sign.) compared to the
$L_{1}$-penalized reduced form VAR($1$) estimator (PVAR). Relative MAFE and
MSFE are abbreviated by RMAFE and RMSFE, respectively. Values below (above) 1
indicate superior (inferior) performance compared to PVAR.

We first look at the mean squared forecast errors (MSFEs). The window-mean
forecast (CONST) clearly does not improve the benchmark PVAR forecast for any
spatial unit. However, this forecast still attains a RMSFE of 1.185,
potentially indicating a low predictability of $\text{NO}_{2}$ column
densities. The GMWY approach obtains the worst forecast performance, possibly
because relatively large bandwidths are needed to allow for second-order
horizontal interaction, thereby paralleling our findings from the simulations.
SPLASH(0,$\lambda$) improves upon the benchmark for all 45 spatial units, 42
of which are found to be significant by a Diebold-Mariano test (based on the
absolute forecast errors). Allowing for sparsity within groups does not
deliver additional forecast improvements as SPLASH($\alpha$,$\lambda$) attains
a slightly worse forecast performance than the SPLASH specification with group
penalty only ($\alpha=0$). Completely omitting regularization at the group
level further deteriorate the forecast performance, with SPLASH(1,$\lambda$)
only beating the benchmark for 34 out 45 spatial units, 14 of which are
statistically significant. We take this as evidence that the ability to
promote diagonally structured sparsity is indeed beneficial in real-life
spatial applications.

Next, we focus on the mean absolute forecast error (MAFE). The results are
qualitatively similar to those obtained based on the MSFE. In particular, GMWY
obtains the worst forecast accuracy for GMWY, whereas SPLASH(0,$\lambda$))
forecasts best. The GMWY method, while still standing out, does not score as
poorly anymore based on the RMAFE. We conjecture that the absence of
regularization may increase sensitivity to noise, thereby resulting in
particularly high squared forecast errors at periods of atypical NO2
concentrations.

Finally, we illustrate the second key benefit of SPLASH-type estimators:
interpretability. Recall that we convoluted satellite images to a ($5\times
9$) grid of spatial units. To examine the relevant interactions between these
spatial units, we provide several visualizations off the spatial weight
matrices estimated by the SPLASH(0.5,$\lambda$) estimator in Figure 4. First,
in panel (a), we visualize the absolute magnitude of the spatial interactions.
A clear diagonal pattern emerges, with the two diagonals closest to the
principal diagonal and the two outer diagonals containing the largest
interactions. These four diagonals correspond to first-order vertical and
second-order horizontal interactions, respectively. The additional two
diagonals, that are sandwiched in between the former, contain the first-order
horizontal interactions between spatial units, which surprisingly seem to be
smaller in magnitude. In panel (b), each cell indicates the proportion of
rolling windows the corresponding spatial interaction is estimated as being
non-zero. It turns out, that these proportions are either one (yellow) or zero
(purple), indicating very stable selection across samples. It becomes apparent
that in addition to the six diagonals that were clear from panel (a), two
additional diagonals are always selected, which contain the first order
diagonal interactions between spatial units. To facilitate interpretation of
this sparsity pattern, we provide a spatial plot of our region of interest in
panel (c), with the spatial grid overlaid. For two pixels on this vertically
labelled spatial grid, we visualize the interactions implied by panel (b),
namely pixel 1 (left-top) and 23 (center) with the use of arrows whose
thickness is determined by the average absolute magnitudes estimated in panel
(a). The pattern of spatial interactions that emerges clearly shows the
interactions between NO2 concentrations of neighbouring districts in London.
The wider horizontal interactions, as well as the diagonal interactions, may
be explainable by the “prevailing winds”, which come from the West or South-
West and are the most commonly occurring winds in London. Overall, the
intuitive sparsity patterns that arise, in combination with the improvement in
forecast performance, are encouraging and provide empirical validation for the
use of SPLASH on spatial data, especially when the spatial units follow a
natural ordering on a spatial grid.

(a) Average absolute magnitudes (b) Selection proportions

(c) Spatial plot of London

Figure 4: Sparsity patterns for estimates of $\bm{A}$ based on rolling window
samples.

##  6 Conclusion

In this paper, we develop the Spatial Lasso-type Shrinkage (SPLASH) estimator,
a novel estimation procedure for high-dimensional spatio-temporal models. The
SPLASH estimator is designed to promote the recovery of diagonally structured
forms of sparsity, which are argued to occur in numerous realistic setting,
whilst not imposing said structure a priori. We derive consistency of our
estimator in a joint asymptotic framework in which both the number of spatial
units and the temporal observations diverge. To solve the identifiability
issue, we rely on a relatively non-restrictive assumption that the coefficient
matrices in the spatio-temporal model are sufficiently banded. Based on this
assumption, we consider banded estimation of high-dimensional spatio-temporal
autocovariance matrices, for which we derive novel convergence rates that are
likely to be of independent interest. As an application, we use SPLASH to
model NO2 concentrations in London based on satellite data and we find
convincing evidence for the presence spatial interactions between horizontally
and vertically neighbouring regions. In addition, our estimator obtains
superior forecast accuracy compared to a number of competitive benchmarks,
including the recently introduced spatio-temporal estimator by Gao et al.,
(2019) that formed the inspiration for the development of SPLASH.

##  7 Acknowledgements

This paper (or earlier versions hereof) has been presented at the seminar of
the Quantitative Economics department of Maastricht University, the
Econometrics Internal Seminar (EIS) at Erasmus University Rotterdam, and the
online Bernoulli-IMS One World Symposium. We gratefully acknowledge comments
and feedback from the participants. Suggestions by Stephan Smeekes and Ines
Wilms have been particularly helpful so we thank them explicitly. All
remaining errors are our own.

## References

  * Ahrens and Bhattacharjee, (2015) Ahrens, A. and Bhattacharjee, A. (2015).  Two-step lasso estimation of the spatial weights matrix.  Econometrics, 3:128–155.
  * Bickel et al., (2009) Bickel, P. J., Ritov, Y., and Tsybakov, A. B. (2009).  Simultaneous analysis of lasso and Dantzig selector.  The Annals of Statistics, 37:1705–1732.
  * Brockwell and Davis, (1991) Brockwell, P. J. and Davis, R. A. (1991).  Time Series: Theory and Methods.  Springer-Verlag, New York, 2nd edition.
  * Debarsy and LeSage, (2018) Debarsy, N. and LeSage, J. (2018).  Flexible dependence modeling using convex combinations of different types of connectivity structures.  Regional Science and Urban Economics, 69:48–68.
  * Dou et al., (2016) Dou, B., Parrell, M. L., and Yao, Q. (2016).  Generalized Yule-Walker estimation for spatio-temporal models with unknown diagonal coefficients.  Journal of Econometrics, 194:369–382.
  * Gao et al., (2019) Gao, Z., Ma, Y., Wang, H., and Yao, Q. (2019).  Banded spatio-temporal autoregressions.  Journal of Econometrics, 208:211–230.
  * Gelper et al., (2016) Gelper, S., Wilms, I., and Croux, C. (2016).  Identifying demand effects in a large network of product categories.  Journal of Retailing, 92:25–39.
  * Guo et al., (2016) Guo, S., Wang, Y., and Yao, Q. (2016).  High-dimensional and banded vector autoregressions.  Biometrika, 103:889–903.
  * Jing et al., (2003) Jing, B.-Y., Shao, Q.-M., and Wang, Q. (2003).  Self-normalized Cramér-type large deviations for independent random variables.  Annals of Probability, 31:2167–2215.
  * Junger and Ponce De Leon, (2015) Junger, W. L. and Ponce De Leon, A. (2015).  Imputation of missing data in time series for air pollutants.  Atmospheric Environment, 102:96–104.
  * Kock and Callot, (2015) Kock, A. B. and Callot, L. (2015).  Oracle inequalities for high dimensional vector autoregressions.  Journal of Econometrics, 186:325–344.
  * Lam and Souza, (2014) Lam, C. and Souza, P. C. L. (2014).  Regularization for spatial panel time series using the adaptive lasso.  Working paper, London School of Economics and Political Science.
  * Lam and Souza, (2019) Lam, C. and Souza, P. C. L. (2019).  Estimation and selection of spatial weight matrix in a spatial lag model.  Journal of Business & Economic Statistics, 38:1–41.
  * Lee, (2004) Lee, L.-F. (2004).  Asymptotic distributions of quasi-maximum likelihood estimators for spatial autoregressive models.  Econometrica, 72:1899–1925.
  * Lee and Yu, (2010) Lee, L.-F. and Yu, J. (2010).  Estimation of spatial autoregressive panel data models with fixed effects.  Journal of Econometrics, 154:165–185.
  * Lee and Yu, (2014) Lee, L.-f. and Yu, J. (2014).  Efficient GMM estimation of spatial dynamic panel data models with fixed effects.  Journal of Econometrics, 180:174–197.
  * Ma et al., (2021) Ma, Y., Guo, S., and Wang, H. (2021).  Sparse spatio-temporal autoregressions by profiling and bagging.  Journal of Econometrics, X:XX–XX.
  * Masini et al., (2019) Masini, R. P., Medeiros, M. C., and Mendes, E. F. (2019).  Regularized estimation of high-dimensional vector autoregressions with weakly dependent innovations.  arXiv e-prints 1912.09002, arXiv.
  * Simon et al., (2013) Simon, N., Friedman, J., Hastie, T., and Tibshirani, R. (2013).  A sparse-group lasso.  Journal of Computational and Graphical Statistics, 22:231–245.
  * van der Vaart and Wellner, (1996) van der Vaart, A. W. and Wellner, J. A. (1996).  Weak Convergence and Empirical Processes.  Springer.
  * Wainwright, (2019) Wainwright, M. J. (2019).  High-dimensional Statistics: A Non-asymptotic Viewpoint.  Cambridge University Press.
  * Yu et al., (2012) Yu, J., de Jong, R., and Lee, L.-f. (2012).  Estimation for spatial dynamic panel data with fixed effects: The case of spatial cointegration.  Journal of Econometrics, 167:16–37.
  * Yu et al., (2008) Yu, J., de Jong, R. M., and Lee, L.-F. (2008).  Quasi-maximum likelihood estimators for spatial dynamic panel data with fixed effects when both n and T are large.  Journal of Econometrics, 146:118–134.
  * Zhang and Yu, (2018) Zhang, X. and Yu, J. (2018).  Spatial weight matrix selection and model averaging for spatial autoregressive models.  Journal of Econometrics, 203:1–18.

##  Appendix A Lemmas

######  Lemma 1.

Define the quantities $N_{c}=\left\lvert\bm{c}\right\rvert$, $S=\\{j:c_{j}\neq
0\\}$,
$\mathscr{G}_{S}=\left\\{g\in\mathscr{G}:\bm{c}_{g}\neq\bm{0}\right\\}$,
$\mathscr{G}_{S}^{c}=\left\\{g\in\mathscr{G}:\bm{c}_{g}=\bm{0}\right\\}$,
$\bar{\omega}_{\alpha}=\max\left\\{(1-\alpha)\sum_{g\in\mathscr{G}_{S}}\sqrt{\left\lvert
g\right\rvert},\alpha\sqrt{\left\lvert S\right\rvert}\right\\}$ and consider

$\bm{\varDelta}\in\mathscr{C}_{N_{c}}(\mathscr{G},S):=\left\\{\bm{\varDelta}\in\mathbb{R}^{N_{c}}:P_{\alpha,S^{c}}(\bm{\varDelta})\leq
3P_{\alpha,S}(\bm{\varDelta})\right\\},$

where

$\begin{split}P_{\alpha,S^{c}}(\bm{\varDelta})&=(1-\alpha)\sum_{g\in\mathscr{G}_{S}^{c}}\sqrt{\left\lvert
g\right\rvert}\left\lVert\Delta_{g}\right\rVert_{2}+\alpha\left\lVert\bm{\varDelta}_{S^{c}}\right\rVert_{1},\text{
and }\\\
P_{\alpha,S}(\bm{\varDelta})&=(1-\alpha)\sum_{g\in\mathscr{G}_{S}}\sqrt{\left\lvert
g\right\rvert}\left\lVert\Delta_{g}\right\rVert_{2}+\alpha\left\lVert\bm{\varDelta}_{S}\right\rVert_{1}.\end{split}$

Then, under Assumption 4, it holds that

$\min_{\bm{\varDelta}\in\mathscr{C}_{N_{c}}(\mathscr{G},S)}\frac{\bar{\omega}_{\alpha}\left\lVert\bm{V}^{(d)}\bm{\varDelta}\right\rVert_{2}}{P_{\alpha,S}(\bm{\varDelta})}\geq\frac{\phi_{0}}{2}.$
(8)

###### Proof.

First, we show that (8) is bounded from below by $\tfrac{1}{2}$ times the
smallest singular value of $\bm{V}^{(d)}$. For $0\leq\alpha<1$, it holds that

$\frac{(1-\alpha)\sum_{g\in\mathscr{G}_{S}}\sqrt{\left\lvert
g\right\rvert}\left\lVert\bm{\varDelta}_{g}\right\rVert_{2}}{\bar{\omega}_{\alpha}}\leq\frac{\left\lVert\bm{\varDelta}\right\rVert_{2}\sum_{g\in\mathscr{G}_{S}}\sqrt{\left\lvert
g\right\rvert}}{\sum_{g\in\mathscr{G}_{S}}\sqrt{\left\lvert
g\right\rvert}}=\left\lVert\bm{\varDelta}\right\rVert_{2}$

Alternatively, for $\alpha=1$, we have
$\frac{\alpha\left\lVert\bm{\varDelta}_{S}\right\rVert_{1}}{\bar{\omega}_{\alpha}}=\frac{\left\lVert\bm{\varDelta}_{S}\right\rVert_{1}}{\sqrt{\left\lvert
S\right\rvert}}\leq\left\lVert\bm{\varDelta}\right\rVert_{2}$. Combining the
two previous results yields

$\frac{P_{\alpha,S}(\bm{\varDelta})}{\bar{\omega}_{\alpha}}=\frac{(1-\alpha)\sum_{g\in\mathscr{G}_{S}}\sqrt{\left\lvert
g\right\rvert}\left\lVert\bm{\varDelta}_{g}\right\rVert_{2}+\alpha\left\lVert\bm{\varDelta}_{S}\right\rVert_{1}}{\bar{\omega}_{\alpha}}\leq
2\left\lVert\bm{\varDelta}\right\rVert_{2},$

for any $0\leq\alpha\leq 1$.

Next, the result in Lemma 1 follows by noting that
$\bm{V}^{(d)}=\operatorname{\text{diag}}\left(\bm{V}_{1},\ldots,\bm{V}_{N}\right)$
is a block-diagonal matrix whose singular values correspond to those of its
sub-blocks. Let $N_{i}$ denote the number of columns of $\bm{V}_{i}$ and note
that $N_{i}<N$ by construction. Then,

$\displaystyle\min_{\bm{\varDelta}\in\mathscr{C}_{N_{c}}(\mathscr{G},S)}\frac{\omega_{\alpha}\left\lVert\bm{V}^{(d)}\bm{\varDelta}\right\rVert_{2}}{P_{\alpha,S}(\bm{\varDelta})}$
$\displaystyle\geq\min_{\bm{\varDelta}\in\mathbb{R}^{N_{c}}}\frac{\left\lVert\bm{V}^{(d)}\bm{\varDelta}\right\rVert_{2}}{2\left\lVert\bm{\varDelta}\right\rVert_{2}}=\min_{1\leq
i\leq
N}\min_{\bm{\varDelta}\in\mathbb{R}^{N_{i}}}\frac{\left\lVert\bm{V}_{i}\bm{\varDelta}\right\rVert_{2}}{2\left\lVert\bm{\varDelta}\right\rVert_{2}}$
$\displaystyle\overset{(i)}{\geq}\min_{\bm{\varDelta}\in\mathbb{R}^{2N}:\left\lVert\bm{\varDelta}\right\rVert_{0}\leq
N}\frac{\left\lVert\bm{V}\bm{\varDelta}\right\rVert_{2}}{2\left\lVert\bm{\varDelta}\right\rVert_{2}}\overset{(ii)}{\geq}\frac{\phi_{0}}{2},$

where (i) follows since $N_{i}<N$ for all $i=1,\ldots,N$ and (ii) holds by
Assumption 4. ∎

######  Lemma 2.

Define the set
$\mathscr{V}(x):=\left\\{\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{2}\leq
x\right\\}$. Then, under Assumption 4, it holds on
$\mathscr{V}\left(\frac{\phi_{0}}{4}\right)$ that

$\min_{\bm{x}\in\mathscr{C}_{N_{c}}(\mathscr{G},S)}\frac{\bar{\omega}_{\alpha}\left\lVert\hat{\bm{V}}_{h}^{(d)}\bm{x}\right\rVert_{2}}{P_{\alpha,S}(\bm{x})}\geq\frac{\phi_{0}}{4}.$

###### Proof.

First, recall the construction of
$\hat{\bm{V}}_{h}^{(d)}=\operatorname{\text{diag}}\left(\hat{\bm{V}}_{1,h},\ldots,\hat{\bm{V}}_{N,h}\right)$
with $\hat{\bm{V}}_{i,h}$ containing at most $N-1$ columns of the matrix
$\hat{\bm{V}}_{h}=\left(\mathscr{B}_{h}\left(\hat{\bm{\varSigma}}_{1}^{\prime}\right),\mathscr{B}_{h}\left(\hat{\bm{\varSigma}}_{0}\right)\right)$.
From the block-diagonal construction, it follows that

$\left\lVert\hat{\bm{V}}_{h}^{(d)}-\bm{V}^{(d)}\right\rVert_{2}=\max_{1\leq
i\leq
N}\left\lVert\hat{\bm{V}}_{i,h}-\bm{V}_{i}\right\rVert_{2}\leq\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{2}.$

Then,

$\left\lVert\hat{\bm{V}}_{h}^{(d)}\bm{x}\right\rVert_{2}\geq\left\lVert\bm{V}^{(d)}\bm{x}\right\rVert_{2}-\left\lVert\left(\hat{\bm{V}}_{h}^{(d)}-\bm{V}^{(d)}\right)\bm{x}\right\rVert_{2}\geq\left\lVert\bm{V}^{(d)}\bm{x}\right\rVert_{2}-\left\lVert\hat{\bm{V}}-\bm{V}\right\rVert_{2}\left\lVert\bm{x}\right\rVert_{2}\geq\left\lVert\bm{V}^{(d)}\bm{x}\right\rVert_{2}-\frac{\phi_{0}}{2}\left\lVert\bm{x}\right\rVert_{2},$
(9)

where the last inequality follows holds on the set
$\mathscr{V}\left(\frac{\phi_{0}}{4}\right)$. Consequently,

$\min_{\bm{x}\in\mathscr{C}_{N_{c}}(\mathscr{G},S)}\frac{\omega_{\alpha}\left\lVert\hat{\bm{V}}_{h}^{(d)}\bm{x}\right\rVert_{2}}{P_{\alpha,S}(\bm{\varDelta})}\overset{(i)}{\geq}\min_{\bm{x}\in\mathbb{R}^{N_{c}}}\frac{\left\lVert\hat{\bm{V}}_{h}^{(d)}\bm{x}\right\rVert_{2}}{2\left\lVert\bm{x}\right\rVert_{2}}\overset{(ii)}{\geq}\min_{\bm{x}\in\mathbb{R}^{N_{c}}}\frac{\left\lVert\bm{V}^{(d)}\bm{x}\right\rVert_{2}}{2\left\lVert\bm{x}\right\rVert_{2}}-\frac{\phi_{0}}{4}\overset{(iii)}{\geq}\frac{\phi_{0}}{4},$

where (i) follows from the proof of Lemma 1, (ii) by (9) and (iii) from Lemma
1. ∎

######  Lemma 3.

Let Assumptions 1 and 2 hold. Moreover, let $i\in\\{1,2,\ldots,N\\}$ be an
arbitrary index and let $x$ be non-negative.

  1. (a)

Under polynomial tail decay of the innovations, Assumption 2(b1), we have

$\operatorname{\mathbb{P}}\left(\left\lVert\frac{1}{T}\sum_{t=2}^{T}\bm{y}_{t-1}\epsilon_{it}\right\rVert_{max}>x\right)\leq
C_{5}\frac{N}{x^{d}T^{d/2}}.$

  2. (b)

Under exponential tail decay of the innovations, Assumption 2(b2), we have

$\operatorname{\mathbb{P}}\Bigg{(}\left\lVert\frac{1}{T}\sum_{t=2}^{T}\bm{y}_{t-1}\epsilon_{it}\right\rVert_{max}>x\Bigg{)}\leq
C_{6}N\left(1+\frac{1}{x}\right)\exp\left(-C_{7}(xT)^{1/3}\right).$

The constants $C_{5}$, $C_{6}$, and $C_{7}$ are positive and independent of
$N$ and $T$. The exact expression can be found in the proof below.

###### Proof.

(a) The required upper bound follows from

$\displaystyle\operatorname{\mathbb{P}}$
$\displaystyle\Bigg{(}\left\lVert\frac{1}{T}\sum_{t=2}^{T}\bm{y}_{t-1}\epsilon_{it}\right\rVert_{max}>x\Bigg{)}=\operatorname{\mathbb{P}}\left(\max_{1\leq
k\leq N}\left|\frac{1}{T}\sum_{t=2}^{T}y_{k,t-1}\epsilon_{it}\right|\geq
x\right)\leq\sum_{k=1}^{N}\operatorname{\mathbb{P}}\left(\left|\frac{1}{T}\sum_{t=2}^{T}y_{k,t-1}\epsilon_{it}\right|\geq
x\right)$
$\displaystyle\overset{(i)}{\leq}\frac{1}{x^{d}T^{d}}\sum_{k=1}^{N}\operatorname{\mathbb{E}}\left|\sum_{t=2}^{T}y_{k,t-1}\epsilon_{it}\right|^{d}\overset{(ii)}{\leq}\frac{(d-1)^{d}}{x^{d}T^{d}}\sum_{k=1}^{N}\operatorname{\mathbb{E}}\left(\sum_{t=2}^{T}y_{k,t-1}^{2}\epsilon_{it}^{2}\right)^{d/2}\overset{(iii)}{\leq}\frac{(d-1)^{d}}{x^{d}T^{d/2+1}}\sum_{k=1}^{N}\sum_{t=2}^{T}\left\lVert
y_{k,t-1}\epsilon_{it}\right\rVert_{d}^{d}$
$\displaystyle\overset{(iv)}{\leq}\frac{(d-1)^{d}}{x^{d}}\frac{N}{T^{d/2}}(\mu_{2d}c_{D}c_{1})^{d}\mu_{2d}^{d},$

by (i) the Markov inequality, (ii) Burkhölder-Davis-Gundy inequality, (iii)
the $C_{r}$-inequality, and (iv) a combination of the Cauchy-Schwartz
inequality and the moment inequalities in Assumption 2(b1). Set
$C_{5}=(d-1)^{d}(\mu_{2d}c_{D}c_{1})^{d}\mu_{2d}^{d}$. (b) For some $C>0$ to
be decided later, we have

$\displaystyle\operatorname{\mathbb{P}}$
$\displaystyle\Bigg{(}\left\lVert\frac{1}{T}\sum_{t=2}^{T}\bm{y}_{t-1}\epsilon_{it}\right\rVert_{max}>x\Bigg{)}\leq\sum_{k=1}^{N}\operatorname{\mathbb{P}}\left(\left|\sum_{t=2}^{T}y_{k,t-1}\epsilon_{it}\right|\geq
xT\right)$ (10)
$\displaystyle\leq\sum_{k=1}^{N}\operatorname{\mathbb{P}}\left(\left|\sum_{t=2}^{T}y_{k,t-1}\epsilon_{it}\mathbbm{1}_{\\{|y_{k,t-1}\epsilon_{it}|\leq
C\\}}\right|\geq\frac{xT}{2}\right)+\sum_{k=1}^{N}\frac{2}{xT}\operatorname{\mathbb{E}}\left(\left|\sum_{t=2}^{T}y_{k,t-1}\epsilon_{it}\mathbbm{1}_{\\{|y_{k,t-1}\epsilon_{it}|>C\\}}\right|\right)$
$\displaystyle\leq
2N\exp\left(-\frac{x^{2}T}{8C^{2}}\right)+\frac{2}{xT}\sum_{k=1}^{N}\sum_{t=2}^{T}\operatorname{\mathbb{E}}\left|y_{k,t-1}\epsilon_{it}\mathbbm{1}_{\\{|y_{k,t-1}\epsilon_{it}|>C\\}}\right|,$

where the first term was bounded using the Azuma-Hoeffding inequality.
Subsequently, using Cauchy-Schwartz on the expectation in the RHS,

$\displaystyle\operatorname{\mathbb{E}}\left|y_{k,t-1}\epsilon_{it}\mathbbm{1}_{\\{|y_{k,t-1}\epsilon_{it}|>C\\}}\right|$
$\displaystyle\leq\left\lVert
y_{k,t-1}\epsilon_{it}\right\rVert_{2}\Big{[}\operatorname{\mathbb{P}}\left(|y_{k,t-1}\epsilon_{it}|>C\right)\Big{]}^{1/2}$
$\displaystyle\leq\left\lVert
y_{k,t-1}\right\rVert_{4}\,\left\lVert\epsilon_{it}\right\rVert_{4}\left[\operatorname{\mathbb{P}}\left(|y_{k,t-1}|>\sqrt{C}\right)+\operatorname{\mathbb{P}}\left(|\epsilon_{it}|>\sqrt{C}\right)\right]^{1/2}.$

With exponential tail decay, we can bound moments and tail probabilities (see,
e.g. page 95 in van der Vaart and Wellner, (1996) and Exercise 2.18 in
Wainwright, (2019), and combine with Lemma LABEL:lemma:boundsonmoments). The
resulting bound is

$\displaystyle\operatorname{\mathbb{E}}\left|y_{k,t-1}\epsilon_{it}\mathbbm{1}_{\\{|y_{k,t-1}\epsilon_{it}|>C\\}}\right|$
$\displaystyle\leq(4!)^{2}\mu_{\infty}^{2}c_{D}c_{1}\left[2\exp(-C/\mu_{\infty}c_{D}c_{1})+2\exp(-C/\mu_{\infty})\right]^{1/2}$
$\displaystyle\leq
2(4!)^{2}\mu_{\infty}^{2}c_{D}c_{1}\exp\left(-C/2\mu_{\infty}(1+c_{D}c_{1})\right),$

and (10) can be continued as

$\operatorname{\mathbb{P}}\Bigg{(}\left\lVert\frac{1}{T}\sum_{t=2}^{T}\bm{y}_{t-1}\epsilon_{it}\right\rVert_{max}>x\Bigg{)}\leq
2N\exp\left(-\frac{x^{2}T}{8C^{2}}\right)+4(4!)^{2}\mu_{\infty}^{2}c_{D}c_{1}\frac{N}{x}\exp\left(-\frac{C}{2\mu_{\infty}(1+c_{D}c_{1})}\right).$

The choice $C=(\frac{1}{4}\mu_{\infty}(1+c_{D}c+1)xT)^{1/3}$ sets the
exponents equal, and for this $C$ we find

$\operatorname{\mathbb{P}}\Bigg{(}\left\lVert\frac{1}{T}\sum_{t=2}^{T}\bm{y}_{t-1}\epsilon_{it}\right\rVert_{max}>x\Bigg{)}\leq\left(2+4(4!)^{2}\mu_{\infty}^{2}c_{D}c_{1}\right)N\left(1+\frac{1}{x}\right)\exp\left(-\left(\frac{xT}{128[\mu_{\infty}(1+c_{D}c_{1})]^{2}}\right)^{\frac{1}{3}}\right).$

This establishes the desired result with
$C_{6}=\left(2+4(4!)^{2}\mu_{\infty}^{2}c_{D}c_{1}\right)$ and
$C_{7}=\left(128[\mu_{\infty}(1+c_{D}c_{1})]^{2}\right)^{-1/3}$. ∎

##  Appendix B Proofs of Main Results

###### Proof of Theorem 1.

We first prove various intermediate results, see (a)–(d) below. We afterwards
combine these results and recover Theorem 1.

  1. (a)

The matrix
$\widetilde{\bm{C}}_{s}=\left(\sum_{j=0}^{s-1}\bm{A}^{j}\right)\bm{B}=:\widetilde{\bm{D}}_{s}\bm{B}$
has a maximum bandwidth of $(s+1)(k_{0}-1)+1$ and satisfies

$\left\lVert\widetilde{\bm{C}}_{s}-\bm{C}\right\rVert_{\vdash}\leq
C_{1}\delta_{A}^{s}.$

  2. (b)

Define
$\bm{\varSigma}_{0}^{r,s}=\sum_{j=0}^{r}\widetilde{\bm{C}}_{s}^{j}\widetilde{\bm{D}}_{s}\bm{\varSigma}_{\epsilon}\widetilde{\bm{D}}_{s}^{\prime}(\widetilde{\bm{C}}_{s}^{\prime})^{j}$
with $\widetilde{\bm{C}}_{s}$ as in Theorem 1(a). The matrix
$\bm{\varSigma}_{0}^{r,s}$ is a banded matrix with bandwidth no larger than
$2(rs+r+s)(k_{0}-1)+2l_{0}+1$. Moreover,

$\left\lVert\bm{\varSigma}_{0}^{r,s}-\bm{\varSigma}_{0}\right\rVert_{\vdash}\leq
C_{2}\left(\frac{\delta_{A}^{s}}{\big{[}1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}\big{]}^{3}}+\delta_{C}^{2(r+1)}\right),$

whenever $s$ is large enough such that $C_{1}\delta_{A}^{s}+\delta_{C}<1$.

  3. (c)

Define
$\bm{\varSigma}_{1}^{r,s}=\widetilde{\bm{C}}_{s}\bm{\varSigma}_{0}^{r,s}$ with
$\widetilde{\bm{C}}_{s}$ and $\bm{\varSigma}_{0}^{r,s}$ as in Theorems 1(a)
and 1(b), respectively. The matrix $\bm{\varSigma}_{1}^{r,s}$ is a banded
matrix with bandwidth no larger than $(2rs+2r+3s+1)(k_{0}-1)+2l_{0}+1$.
Moreover,

$\left\lVert\bm{\varSigma}_{1}^{r,s}-\bm{\varSigma}_{1}\right\rVert_{\vdash}\leq
C_{3}\left(\frac{\delta_{A}^{s}}{\big{[}1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}\big{]}^{3}}+\delta_{C}^{2(r+1)}\right),$

whenever $s$ is large enough such that $C_{1}\delta_{A}^{s}+\delta_{C}<1$.

  4. (d)

Take any $h_{1}\geq 2(rs+r+s)(k_{0}-1)+2l_{0}+1$ , then

$\big{\|}\mathscr{B}_{h_{1}}\big{(}\widehat{\bm{\varSigma}}_{0}\big{)}-\bm{\varSigma}_{0}\big{\|}_{2}\leq\epsilon+2C_{2}(h_{1}+1)\left(\frac{\delta_{A}^{s}}{\big{[}1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}\big{]}^{3}}+\delta_{C}^{2(r+1)}\right),$

with a probability of at least $1-\mathscr{P}_{1}(\epsilon,N,T)$ (for
polynomial tail decay) or $1-\mathscr{P}_{2}(\epsilon,N,T)$ (for exponential
tail decay).

  5. (e)

Take any $h_{2}\geq(2rs+2r+3s+1)(k_{0}-1)+2l_{0}+1$ , then

$\big{\|}\mathscr{B}_{h_{2}}\big{(}\widehat{\bm{\varSigma}}_{1}\big{)}-\bm{\varSigma}_{1}\big{\|}_{2}\leq\epsilon+2C_{3}(h_{2}+1)\left(\frac{\delta_{A}^{s}}{\big{[}1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}\big{]}^{3}}+\delta_{C}^{2(r+1)}\right),$

with a probability of at least $1-\mathscr{P}_{1}(\epsilon,N,T)$ (for
polynomial tail decay) or $1-\mathscr{P}_{2}(\epsilon,N,T)$ (for exponential
tail decay).

Explicit expressions for $C_{1}$, $C_{2}$, $C_{3}$, and $0\leq\delta_{C}<1$
are provided in the proofs below.

(a) The proof builds upon results from Guo et al., (2016) on banded vector
autoregressions. Recall that $\bm{C}=\bm{D}\bm{B}$ with
$\bm{D}=(\bm{I}_{N}-\bm{A})^{-1}$.
$\widetilde{\bm{D}}_{s}=\sum_{j=0}^{s-1}\bm{A}^{j}$ has a bandwidth of at most
$s(k_{0}-1)+1$ and satisfies111111If matrices $\bm{F}_{1}$ and $\bm{F}_{2}$
are banded matrices with bandwidths $k_{1}$ and $k_{2}$, respectively, then
the product $\bm{F}_{1}\bm{F}_{2}$ is again a banded matrix with a bandwidth
of at most $k_{1}+k_{2}-1$.

$\left\lVert\widetilde{\bm{D}}_{s}-\bm{D}\right\rVert_{\vdash}=\left\lVert\sum_{j=0}^{s-1}\bm{A}^{j}-(\bm{I}_{N}-\bm{A})^{-1}\right\rVert_{\vdash}=\left\lVert-\sum_{j=s}^{\infty}\bm{A}^{j}\right\rVert_{\vdash}\leq\sum_{j=s}^{\infty}\left\lVert\bm{A}\right\rVert_{\vdash}^{j}\leq\frac{\delta_{A}^{s}}{1-\delta_{A}}.$

The product $\widetilde{\bm{C}}_{s}=\widetilde{\bm{D}}_{s}\bm{B}$ has a
maximal bandwidth of $(s+1)(k_{0}-1)+1$. Since
$\left\lVert\bm{B}\right\rVert_{\vdash}\leq C_{B}$ (Assumption 3(c)), we also
have

$\left\lVert\widetilde{\bm{C}}_{s}-\bm{C}\right\rVert_{\vdash}=\left\lVert\left(\widetilde{\bm{D}}_{s}-\bm{D}\right)\bm{B}\right\rVert_{\vdash}\leq\left\lVert\widetilde{\bm{D}}_{s}-\bm{D}\right\rVert_{\vdash}\left\lVert\bm{B}\right\rVert_{\vdash}\leq
C_{B}\frac{\delta_{A}^{s}}{1-\delta_{A}}.$

and the claim follows with $C_{1}=\frac{C_{B}}{1-\delta_{A}}$. (b) Iterating
on the observation in footnote 11, we conclude that the bandwidth of
$\widetilde{\bm{C}}_{s}^{r}$ is at most
$r\Big{[}(s+1)(k_{0}-1)+1\Big{]}-(r-1)=r(s+1)(k_{0}-1)+1$. The bandwidth of
$\bm{\varSigma}_{0}^{r,s}$ therefore does not exceed

$2\Big{[}r(s+1)(k_{0}-1)+1\Big{]}+2(s(k_{0}-1)+1)+(2l_{0}+1)-4=2(rs+r+s)(k_{0}-1)+2l_{0}+1.$

We now bound
$\left\lVert\bm{\varSigma}_{0}^{r,s}-\bm{\varSigma}_{0}\right\rVert_{1}$.
Assumption 1(b), imposes
$\left\lVert\bm{C}\right\rVert_{\vdash}\leq\delta_{C}$ for some
$0\leq\delta_{C}<1$. Because
$\bm{\varSigma}_{0}=\sum_{j=0}^{\infty}\bm{C}^{j}\bm{D}\bm{\varSigma}_{\epsilon}\bm{D}^{\prime}(\bm{C}^{\prime})^{j}$,
it holds that

$\displaystyle\left\lVert\bm{\varSigma}_{0}^{r,s}-\bm{\varSigma}_{0}\right\rVert_{\vdash}=\left\lVert\sum_{j=0}^{r}\left[\widetilde{\bm{C}}_{s}^{j}\widetilde{\bm{D}}_{s}\bm{\varSigma}_{\epsilon}\widetilde{\bm{D}}_{s}^{\prime}(\widetilde{\bm{C}}_{s}^{\prime})^{j}-\bm{C}^{j}\bm{D}\bm{\varSigma}_{\epsilon}\bm{D}^{\prime}(\bm{C}^{\prime})^{j}\right]-\sum_{j=r+1}^{\infty}\bm{C}^{j}\bm{D}\bm{\varSigma}_{\epsilon}\bm{D}^{\prime}(\bm{C}^{\prime})^{j}\right\rVert_{\vdash}$
(11)
$\displaystyle\leq\sum_{j=0}^{r}\left\lVert\bm{\varSigma}_{\epsilon}\right\rVert_{\vdash}\left\lVert\widetilde{\bm{C}}_{s}^{j}\right\rVert_{\vdash}^{2}\left\lVert\widetilde{\bm{D}}_{s}-\bm{D}\right\rVert_{\vdash}^{2}+2\sum_{j=0}^{r}\left\lVert\bm{\varSigma}_{\epsilon}\right\rVert_{\vdash}\left\lVert\widetilde{\bm{C}}_{s}^{j}\right\rVert_{\vdash}^{2}\left\lVert\bm{D}\right\rVert_{\vdash}\left\lVert\widetilde{\bm{D}}_{s}-\bm{D}\right\rVert_{\vdash}^{2}$
$\displaystyle\qquad+\sum_{j=1}^{r}\left\lVert\bm{\varSigma}_{\epsilon}\right\rVert_{\vdash}\left\lVert\widetilde{\bm{C}}_{s}^{j}-\bm{C}^{j}\right\rVert_{\vdash}^{2}\left\lVert\bm{D}\right\rVert_{\vdash}^{2}+2\sum_{j=1}^{r}\left\lVert\bm{\varSigma}_{\epsilon}\right\rVert_{\vdash}\left\lVert\widetilde{\bm{C}}_{s}^{j}-\bm{C}^{j}\right\rVert_{\vdash}\left\lVert\bm{C}^{j}\right\rVert_{\vdash}\left\lVert\bm{D}\right\rVert_{\vdash}^{2}$
$\displaystyle\qquad+\left\lVert\sum_{j=r+1}^{\infty}\bm{C}^{j}\bm{D}\bm{\varSigma}_{\epsilon}\bm{D}^{\prime}(\bm{C}^{\prime})^{j}\right\rVert_{\vdash}.$

An inspection of (11) shows that additional upper bounds are required on
$\left\lVert\widetilde{\bm{C}}_{s}^{j}\right\rVert_{\vdash}$,
$\left\lVert\bm{C}^{j}-\widetilde{\bm{C}}_{s}^{j}\right\rVert_{\vdash}$, and
$\left\lVert\sum_{j=r+1}^{\infty}\bm{C}^{j}\bm{D}\bm{\varSigma}_{\epsilon}\bm{D}^{\prime}(\bm{C}^{\prime})^{j}\right\rVert_{\vdash}$.
For $\left\lVert\widetilde{\bm{C}}_{s}^{j}\right\rVert_{\vdash}$, using
properties of matrix norms and Theorem 1(a), we obtain

$\left\lVert\widetilde{\bm{C}}_{s}^{j}\right\rVert_{\vdash}=\left\lVert\left(\widetilde{\bm{C}}_{s}-\bm{C}+\bm{C}\right)^{j}\right\rVert_{\vdash}\leq\left(\left\lVert\widetilde{\bm{C}}_{s}-\bm{C}\right\rVert_{\vdash}+\left\lVert\bm{C}\right\rVert_{\vdash}\right)^{j}\leq\Big{[}C_{1}\delta_{A}^{s}+\delta_{C}\Big{]}^{j}.$
(12)

Furthermore, expanding the matrix powers provides

$\displaystyle\left\lVert\bm{C}^{j}-\widetilde{\bm{C}}_{s}^{j}\right\rVert_{\vdash}=\left\lVert\bm{C}^{j}-\left(\widetilde{\bm{C}}_{s}-\bm{C}+\bm{C}\right)^{j}\right\rVert_{\vdash}$
(13)
$\displaystyle=\left\lVert\bm{C}^{j}-\left((\widetilde{\bm{C}}_{s}-\bm{C})^{j}+(\widetilde{\bm{C}}_{s}-\bm{C})^{j-1}\bm{C}+(\widetilde{\bm{C}}_{s}-\bm{C})^{j-2}\bm{C}(\widetilde{\bm{C}}_{s}-\bm{C})+\ldots+\bm{C}(\widetilde{\bm{C}}_{s}-\bm{C})^{j-1}\ldots+\bm{C}^{j}\right)\right\rVert_{\vdash}$
$\displaystyle\leq\sum_{k=1}^{j}\binom{j}{k}\left\lVert\bm{C}\right\rVert_{\vdash}^{j-k}\left\lVert\widetilde{\bm{C}}_{s}-\bm{C}\right\rVert_{\vdash}^{k}=\sum_{k=0}^{j-1}\binom{j}{k+1}\left\lVert\bm{C}\right\rVert_{\vdash}^{(j-1)-k}\left\lVert\widetilde{\bm{C}}_{s}-\bm{C}\right\rVert_{\vdash}^{k+1}$
$\displaystyle=\left\lVert\widetilde{\bm{C}}_{s}-\bm{C}\right\rVert_{\vdash}\sum_{k=0}^{j-1}\frac{j}{k+1}\binom{j-1}{k}\left\lVert\bm{C}\right\rVert_{\vdash}^{(j-1)-k}\left\lVert\widetilde{\bm{C}}_{s}-\bm{C}\right\rVert_{\vdash}^{k}\leq
C_{1}j\delta_{A}^{s}\left[\left\lVert\widetilde{\bm{C}}_{s}-\bm{C}\right\rVert_{\vdash}+\left\lVert\bm{C}\right\rVert_{\vdash}\right]^{j-1}$
$\displaystyle\leq
C_{1}j\delta_{A}^{s}\Big{[}C_{1}\delta_{A}^{s}+\delta_{C}\Big{]}^{j-1}$

Finally, since $\left\lVert\bm{\varSigma}_{\epsilon}\right\rVert_{\vdash}\leq
C_{\epsilon}$ and
$\left\lVert\bm{D}\right\rVert_{\vdash}=\left\lVert\sum_{j=0}^{\infty}\bm{A}^{j}\right\rVert_{\vdash}=\sum_{j=0}^{\infty}\left\lVert\bm{A}\right\rVert_{\vdash}^{j}\leq
c_{D}$ (Assumption 1(a)), we have

$\left\lVert\sum_{j=r+1}^{\infty}\bm{C}^{j}\bm{D}\bm{\varSigma}_{\epsilon}\bm{D}^{\prime}(\bm{C}^{\prime})^{j}\right\rVert_{\vdash}\leq\sum_{j=r+1}^{\infty}\left\lVert\bm{\varSigma}_{\epsilon}\right\rVert_{\vdash}\left\lVert\bm{D}\right\rVert_{\vdash}^{2}\left\lVert\bm{C}\right\rVert_{\vdash}^{2j}\leq
C_{\epsilon}c_{D}^{2}\sum_{j=r+1}^{\infty}\delta_{C}^{2j}=C_{\epsilon}c_{D}^{2}\frac{\delta_{C}^{2(r+1)}}{1-\delta_{C}^{2}}.$
(14)

Returning to (11) and inserting the upper bounds in (12)–(14), we find

$\displaystyle\left\lVert\bm{\varSigma}_{0}^{r,s}-\bm{\varSigma}_{0}\right\rVert_{\vdash}\leq\frac{C_{\epsilon}}{(1-\delta_{A})^{2}}\delta_{A}^{2s}\sum_{j=0}^{r}\left(C_{1}\delta_{A}^{s}+\delta_{C}\right)^{2j}+2\frac{C_{\epsilon}c_{D}}{(1-\delta_{A})^{2}}\delta_{A}^{2s}\sum_{j=0}^{r}\left(C_{1}\delta_{A}^{s}+\delta_{C}\right)^{2j}$
$\displaystyle+C_{\epsilon}C_{1}^{2}c_{D}^{2}\delta_{A}^{2s}\sum_{j=1}^{r}j^{2}\left(C_{1}\delta_{A}^{s}+\delta_{C}\right)^{2(j-1)}+2C_{\epsilon}C_{1}c_{D}^{2}\delta_{A}^{s}\sum_{j=1}^{r}j\delta_{C}^{j}\left(C_{1}\delta_{A}^{s}+\delta_{C}\right)^{j-1}+C_{\epsilon}c_{D}^{2}\frac{\delta_{C}^{2(r+1)}}{1-\delta_{C}^{2}}.$

Assuming $s$ is sufficient large, that is assuming
$C_{1}\delta_{A}^{s}+\delta_{C}<1$, we subsequently use result on geometric
series and conclude121212Specifically,
$\sum_{j=1}^{\infty}jz^{j}=\frac{z}{(1-z)^{2}}$ and
$\sum_{j=1}^{\infty}j^{2}z^{j}=\frac{z}{(1-z)^{3}}$ for $|z|<1$.

$\displaystyle\left\lVert\bm{\varSigma}_{0}^{r,s}-\bm{\varSigma}_{0}\right\rVert_{\vdash}\leq\frac{C_{\epsilon}}{(1-\delta_{A})^{2}}\delta_{A}^{2s}\frac{1}{1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}}+2\frac{C_{\epsilon}c_{D}}{(1-\delta_{A})^{2}}\delta_{A}^{2s}\frac{1}{1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}}$
$\displaystyle\qquad+\frac{C_{\epsilon}C_{1}^{2}c_{D}^{2}\delta_{A}^{2s}}{(C_{1}\delta_{A}^{s}+\delta_{C})^{2}}\frac{(C_{1}\delta_{A}^{s}+\delta_{C})^{2}}{\big{[}1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}\big{]}^{3}}+2\frac{C_{\epsilon}C_{1}c_{D}^{2}\delta_{A}^{s}}{C_{1}\delta_{A}^{s}+\delta_{C}}\frac{\delta_{C}(C_{1}\delta_{A}^{s}+\delta_{C})}{\big{[}1-\delta_{C}(C_{1}\delta_{A}^{s}+\delta_{C})\big{]}^{2}}+C_{\epsilon}c_{D}^{2}\frac{\delta_{C}^{2(r+1)}}{1-\delta_{C}^{2}}$
$\displaystyle=\frac{C_{\epsilon}}{(1-\delta_{A})^{2}}\frac{1}{1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}}\delta_{A}^{2s}+2\frac{C_{\epsilon}c_{D}}{(1-\delta_{A})^{2}}\frac{1}{1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}}\delta_{A}^{2s}+\frac{C_{\epsilon}C_{1}^{2}c_{D}^{2}}{\big{[}1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}\big{]}^{3}}\delta_{A}^{2s}$
$\displaystyle\qquad+2\frac{C_{\epsilon}C_{1}c_{D}^{2}\delta_{C}}{\big{[}1-\delta_{C}(C_{1}\delta_{A}^{s}+\delta_{C})\big{]}^{2}}\delta_{A}^{s}+C_{\epsilon}c_{D}^{2}\frac{\delta_{C}^{2(r+1)}}{1-\delta_{C}^{2}}$
$\displaystyle\leq\frac{1}{(1-\delta_{A})^{2}}\frac{C_{\epsilon}+2C_{\epsilon}c_{D}}{1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}}\delta_{A}^{2s}+\frac{C_{\epsilon}C_{1}^{2}c_{D}^{2}+2C_{\epsilon}C_{1}c_{D}^{2}\delta_{C}}{\big{[}1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}\big{]}^{3}}\delta_{A}^{s}+\frac{C_{\epsilon}c_{D}^{2}}{1-\delta_{C}^{2}}\delta_{C}^{2(r+1)}$
$\displaystyle\leq\left[\frac{C_{\epsilon}+2C_{\epsilon}c_{D}}{(1-\delta_{A})^{2}}+C_{\epsilon}C_{1}^{2}c_{D}^{2}+2C_{\epsilon}C_{1}c_{D}^{2}\delta_{C}\right]\frac{\delta_{A}^{s}}{\big{[}1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}\big{]}^{3}}+\frac{C_{\epsilon}c_{D}^{2}}{1-\delta_{C}^{2}}\delta_{C}^{2(r+1)}.$

The claim is thus indeed valid with
$C_{2}=\max\left\\{\frac{C_{\epsilon}+2C_{\epsilon}c_{D}}{(1-\delta_{A})^{2}}+C_{\epsilon}C_{1}^{2}c_{D}^{2}+2C_{\epsilon}C_{1}c_{D}^{2}\delta_{C},\frac{C_{\epsilon}c_{D}^{2}}{1-\delta_{C}^{2}}\right\\}$.
(c) We have
$\bm{\varSigma}_{1}=(\bm{I}_{N}-\bm{A})^{-1}\bm{B}\bm{\varSigma}_{0}=\bm{C}\bm{\varSigma}_{0}$,
and hence

$\left\lVert\bm{\varSigma}_{1}^{r,s}-\bm{\varSigma}_{1}\right\rVert_{\vdash}=\left\lVert\widetilde{\bm{C}}_{s}(\bm{\varSigma}_{0}^{r,s}-\bm{\varSigma}_{0})+(\widetilde{\bm{C}}_{s}-\bm{C})\bm{\varSigma}_{0}\right\rVert_{\vdash}\leq\left\lVert\widetilde{\bm{C}}_{s}\right\rVert_{\vdash}\left\lVert\bm{\varSigma}_{0}^{r,s}-\bm{\varSigma}_{0}\right\rVert_{\vdash}+\left\lVert\widetilde{\bm{C}}_{s}-\bm{C}\right\rVert_{\vdash}\left\lVert\bm{\varSigma}_{0}\right\rVert_{\vdash}$
(15)

Now combine $\left\lVert\widetilde{\bm{C}}_{s}\right\rVert_{\vdash}\leq
C_{1}\delta_{\kappa}^{s}+\delta_{C}\leq 1$ (since $s$ is taken sufficiently
large), Theorem 1(b), (13) with $j=1$, and
$\left\lVert\bm{\varSigma}_{0}\right\rVert_{\vdash}\leq\sum_{j=0}^{\infty}\left\lVert\bm{C}^{j}\bm{D}\bm{\varSigma}_{\epsilon}\bm{D}^{\prime}(\bm{C}^{\prime})^{j}\right\rVert_{\vdash}\leq\frac{C_{\epsilon}c_{D}^{2}}{1-\delta_{C}^{2}}$
to obtain the stated result with
$C_{3}=C_{2}+\frac{C_{1}C_{\epsilon}c_{D}^{2}}{1-\delta_{C}^{2}}$. (d) We have

$\displaystyle\left\lVert\mathscr{B}_{h_{1}}\big{(}\widehat{\bm{\varSigma}}_{0}\big{)}-\bm{\varSigma}_{0}\right\rVert_{\vdash}$
$\displaystyle\leq\left\lVert\mathscr{B}_{h_{1}}\big{(}\widehat{\bm{\varSigma}}_{0}\big{)}-\mathscr{B}_{h_{1}}\big{(}\bm{\varSigma}_{0}\big{)}\right\rVert_{\vdash}+\left\lVert\mathscr{B}_{h_{1}}\big{(}\bm{\varSigma}_{0}\big{)}-\mathscr{B}_{h_{1}}\big{(}\bm{\varSigma}_{0}^{r,s}\big{)}\right\rVert_{\vdash}+\left\lVert\mathscr{B}_{h_{1}}\big{(}\bm{\varSigma}_{0}^{r,s}\big{)}-\bm{\varSigma}_{0}\right\rVert_{\vdash}$
(16)
$\displaystyle=\left\lVert\mathscr{B}_{h_{1}}\big{(}\widehat{\bm{\varSigma}}_{0}-\bm{\varSigma}_{0}\big{)}\right\rVert_{\vdash}+\left\lVert\mathscr{B}_{h_{1}}\big{(}\bm{\varSigma}_{0}-\bm{\varSigma}_{0}^{r,s}\big{)}\right\rVert_{\vdash}+\left\lVert\bm{\varSigma}_{0}^{r,s}-\bm{\varSigma}_{0}\right\rVert_{\vdash},$

because $\bm{\varSigma}_{0}^{r,s}$ is a banded matrix already. We consider the
three terms in the RHS of (16) separately. First, there are at most $2h_{1}+1$
nonzero elements in any column/row of
$\mathscr{B}_{h_{1}}\big{(}\widehat{\bm{\varSigma}}_{0}-\bm{\varSigma}_{0}\big{)}$
and thus

$\left\lVert\mathscr{B}_{h_{1}}\big{(}\widehat{\bm{\varSigma}}_{0}-\bm{\varSigma}_{0}\big{)}\right\rVert_{\vdash}\leq(2h_{1}+1)\left\lVert\widehat{\bm{\varSigma}}_{0}-\bm{\varSigma}_{0}\right\rVert_{max}$
(17)

and

$\displaystyle\operatorname{\mathbb{P}}\left(\left\lVert\mathscr{B}_{h_{1}}\big{(}\widehat{\bm{\varSigma}}_{0}-\bm{\varSigma}_{0}\big{)}\right\rVert_{\vdash}\leq
x\right)\geq\operatorname{\mathbb{P}}\left((2h_{1}+1)\left\lVert\widehat{\bm{\varSigma}}_{0}-\bm{\varSigma}_{0}\right\rVert_{max}\leq
x\right)=1-\operatorname{\mathbb{P}}\left(\left\lVert\widehat{\bm{\varSigma}}_{0}-\bm{\varSigma}_{0}\right\rVert_{max}>\frac{x}{2h_{1}+1}\right)$
(18) $\displaystyle\geq 1-\operatorname{\mathbb{P}}\left(\bigcup_{1\leq
i,j\leq
N}\left[\widehat{\bm{\varSigma}}_{0}-\bm{\varSigma}_{0}\right]_{ij}>\frac{x}{2h_{1}+1}\right)\geq
1-\sum_{i=1}^{N}\sum_{j=1}^{N}\operatorname{\mathbb{P}}\left(\left|\sum_{t=2}^{T}x_{it}x_{jt}-\operatorname{\mathbb{E}}(x_{it}x_{jt})\right|>\frac{Tx}{2h_{1}+1}\right)$
$\displaystyle\geq\begin{cases}1-N^{2}\left[\left(b_{1}T^{(1-\delta)/3}+\frac{(2h_{1}+1)b_{3}}{x}\right)\exp\left(-\frac{T^{(1-\delta)/3}}{2b_{1}^{2}}\right)+\frac{b_{2}(2h_{1}+1)^{d}}{x^{d}T^{\frac{\delta}{2}(d-1)}}\right]&\qquad(\text{polynomial
tails}),\\\
1-N^{2}\left[\frac{\beta_{1}(2h_{1}+1)}{x}+\frac{2}{\beta_{2}}\left(\frac{Tx^{2}}{(2h_{1}+1)^{2}}\right)^{1/7}\right]\exp\left(-\frac{1}{\beta_{3}}\left(\frac{Tx^{2}}{(2h_{1}+1)^{2}}\right)^{1/7}\right)&\qquad(\text{exponential
tails}),\end{cases}$

where the last inequality exploits Lemma LABEL:lemma:covmatrixelements (see
Supplement). Note that the probabilities in (18) coincide with the
probabilities defined as $1-\mathscr{P}_{1}(x,N,T)$ and
$1-\mathscr{P}_{2}(x,N,T)$ in Theorem 1. Second, using
$\left\lVert\cdot\right\rVert_{max}\leq\left\lVert\cdot\right\rVert_{\vdash}$,
we have
$\left\lVert\mathscr{B}_{h_{1}}\big{(}\bm{\varSigma}_{0}-\bm{\varSigma}_{0}^{r,s}\big{)}\right\rVert_{\vdash}\leq(2h_{1}+1)\left\lVert\bm{\varSigma}_{0}-\bm{\varSigma}_{0}^{r,s}\right\rVert_{max}\leq(2h_{1}+1)\left\lVert\bm{\varSigma}_{0}^{r,s}-\bm{\varSigma}_{0}\right\rVert_{\vdash}$.
Overall, if
$\left\lVert\mathscr{B}_{h_{1}}\big{(}\widehat{\bm{\varSigma}}_{0}-\bm{\varSigma}_{0}\big{)}\right\rVert_{\vdash}\leq\epsilon$
holds, then continuing from (16),

$\displaystyle\left\lVert\mathscr{B}_{h_{1}}\big{(}\widehat{\bm{\varSigma}}_{0}\big{)}-\bm{\varSigma}_{0}\right\rVert_{\vdash}\leq\epsilon+2(h_{1}+1)\left\lVert\bm{\varSigma}_{0}^{r,s}-\bm{\varSigma}_{0}\right\rVert_{\vdash}\leq\epsilon+2C_{2}N\left(\frac{\delta_{A}^{s}}{\big{[}1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}\big{]}^{3}}+\delta_{C}^{2(r+1)}\right),$

where the last inequality follows from part (b) of this proof and the fact
that $h_{1}+1\leq N$. A lower bound on the probability of
$\left\\{\left\lVert\mathscr{B}_{h_{1}}\big{(}\widehat{\bm{\varSigma}}_{0}-\bm{\varSigma}_{0}\big{)}\right\rVert_{\vdash}\leq\epsilon\right\\}$
occurring is immediately available from (18). (e) Mimicking the steps from
(d), we find

$\left\lVert\mathscr{B}_{h_{2}}\big{(}\widehat{\bm{\varSigma}}_{1}\big{)}-\bm{\varSigma}_{1}\right\rVert_{\vdash}\leq\epsilon+2(h_{2}+1)\left\lVert\bm{\varSigma}_{1}^{r,s}-\bm{\varSigma}_{1}\right\rVert_{\vdash}\leq\epsilon+2C_{3}N\left(\frac{\delta_{A}^{s}}{\big{[}1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}\big{]}^{3}}+\delta_{C}^{2(r+1)}\right)$

if we use part (c) and if
$\big{\\{}\left\lVert\mathscr{B}_{h_{2}}\big{(}\widehat{\bm{\varSigma}}_{1}-\bm{\varSigma}_{1}\big{)}\right\rVert_{\vdash}\leq\epsilon\big{\\}}$
holds. The probability of the latter event is
$1-\mathscr{P}_{1}(\epsilon,N,T)$ (polynomial tails) or
$1-\mathscr{P}_{2}(\epsilon,N,T)$ (exponential tails).

We now combine all these intermediate results to recover the result from the
theorem. Parts (d) and (e) are both applicable since $h=\max\\{h_{1},h_{2}\\}$
. Lemma LABEL:lemma:L2normineqs implies that

$\displaystyle\left\lVert\widehat{\bm{V}}_{h}-\bm{V}\right\rVert_{\vdash}$
$\displaystyle\leq\left\lVert\mathscr{B}_{h_{2}}\big{(}\widehat{\bm{\varSigma}}_{1}\big{)}-\bm{\varSigma}_{1}\right\rVert_{\vdash}+\left\lVert\mathscr{B}_{h_{1}}\big{(}\widehat{\bm{\varSigma}}_{0}\big{)}-\bm{\varSigma}_{0}\right\rVert_{\vdash}$
$\displaystyle\leq
2\epsilon+2C_{4}N\left(\frac{\delta_{\kappa}^{s}}{\big{[}1-(C_{1}\delta_{\kappa}^{s}+\delta_{C})^{2}\big{]}^{3}}+\delta_{C}^{2(r+1)}\right),$

when
$\left\\{\left\lVert\mathscr{B}_{h_{2}}\big{(}\widehat{\bm{\varSigma}}_{1}-\bm{\varSigma}_{1}\big{)}\right\rVert_{\vdash}\leq\epsilon\right\\}\cap\left\\{\left\lVert\mathscr{B}_{h_{1}}\big{(}\widehat{\bm{\varSigma}}_{0}-\bm{\varSigma}_{0}\big{)}\right\rVert_{\vdash}\leq\epsilon\right\\}$
takes place and $C_{4}=C_{2}+C_{3}$. From $\operatorname{\mathbb{P}}(A\cap
B)=1-\operatorname{\mathbb{P}}(A^{c}\cup B^{c})\geq
1-\operatorname{\mathbb{P}}(A^{c})-\operatorname{\mathbb{P}}(B^{c})$ it
follows that this joint event takes place with probabilities
$1-2\mathscr{P}_{1}(\epsilon,N,T)$ and $1-2\mathscr{P}_{2}(\epsilon,N,T)$ in
the cases of polynomial and exponential tail decay, respectively.

It remains to determine $s$ and $r$ such that
$C_{4}N\delta_{A}^{s}/\big{[}1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}\big{]}^{3}\leq\epsilon$
and $C_{4}N\delta_{C}^{2(r+1)}\leq\epsilon$ hold. For any $s\geq s^{*}$ with
$s^{*}=\log\left(\frac{1}{C_{1}}\big{[}(1-(1-\delta_{C})^{1/3})^{1/2}-\delta_{C}\big{]}\right)/\log(\delta_{A})$,
we have $1/[1-(C_{1}\delta_{A}^{s}+\delta_{C})^{2}]^{3}\leq 1/[1-\delta_{C}]$.
Under the latter assumption, we determine the $s$ such that

$\displaystyle\frac{C_{4}N\delta_{A}^{s}}{1-\delta_{C}}\leq\epsilon\iff
s\log(\delta_{A})\leq\log\left(\frac{(1-\delta_{C})\epsilon}{C_{4}N}\right)\iff
s\geq\frac{\log\left(C_{4}N/(1-\delta_{C})\epsilon\right)}{|\log(\delta_{A})|}.$

Similarly, for $r$, the choice $r=\log(C_{4}N/\epsilon)/2|\log(\delta_{C})|$
suffices. ∎

###### Proof of Theorem 2.

The proof of the theorem relies on the properties of dual norms. Recall
$P_{\alpha}(\bm{c})=\alpha\sum_{g\in\mathscr{G}}\sqrt{|g|}\left\lVert\bm{c}_{g}\right\rVert_{2}+(1-\alpha)\left\lVert\bm{c}\right\rVert_{1}$.
Exploiting the properties of $\left\lVert\cdot\right\rVert_{1}$ and
$\left\lVert\cdot\right\rVert_{2}$, it is straightforward to verify that
$P_{\alpha}(\cdot)$ is a norm for any $0\leq\alpha\leq 1$. For any norm
$\left\lVert\cdot\right\rVert$, we define its dual norm
$\left\lVert\cdot\right\rVert^{*}$ through
$\left\lVert\bm{c}\right\rVert^{*}=\sup_{\bm{x}\neq\bm{0}}\frac{\left\lvert\bm{c}^{\prime}\bm{x}\right\rvert}{\left\lVert\bm{x}\right\rVert}$.
The dual-norm inequality states that

$\bm{c}^{\prime}\bm{x}\leq\left\lVert\bm{c}\right\rVert^{*}\,\left\lVert\bm{x}\right\rVert\qquad\qquad\text{for
all conformable vectors $\bm{c}$ and $\bm{x}$}.$ (19)

For the norm $P_{\alpha}(\bm{c})$, its dual norm $P_{\alpha}^{*}(\bm{c})$ is
bounded by

$\displaystyle P_{\alpha}^{*}(\bm{c})$
$\displaystyle=\sup_{\bm{x}\neq\bm{0}}\frac{|\bm{c}^{\prime}\bm{x}|}{P_{\alpha}(\bm{x})}=\sup_{\bm{x}\neq\bm{0}}\frac{|\bm{c}^{\prime}\bm{x}|}{\alpha\sum_{g\in\mathscr{G}}\sqrt{|g|}\left\lVert\bm{x}_{g}\right\rVert_{2}+(1-\alpha)\left\lVert\bm{x}\right\rVert_{1}}$
(20)
$\displaystyle\stackrel{{\scriptstyle(i)}}{{\leq}}\alpha\sup_{\bm{x}\neq\bm{0}}\frac{|\bm{c}^{\prime}\bm{x}|}{\sum_{g\in\mathscr{G}}\sqrt{|g|}\left\lVert\bm{x}_{g}\right\rVert_{2}}+(1-\alpha)\sup_{\bm{x}\neq\bm{0}}\frac{|\bm{c}^{\prime}\bm{x}|}{\left\lVert\bm{x}\right\rVert_{1}}\stackrel{{\scriptstyle(ii)}}{{\leq}}\alpha\max_{g\in\mathscr{G}}\frac{\left\lVert\bm{c}_{g}\right\rVert_{2}}{\sqrt{\left\lvert
g\right\rvert}}+(1-\alpha)\left\lVert\bm{c}\right\rVert_{\infty},$

by convexity of the function $f(x)=x^{-1}$ in step (i), and using for step
(ii) both
$\left\lVert\bm{c}\right\rVert_{1}^{*}=\left\lVert\bm{c}\right\rVert_{\infty}$
and

$\displaystyle\sup_{\bm{x}\neq\bm{0}}\frac{\left\lvert\bm{c}^{\prime}\bm{x}\right\rvert}{\sum_{g\in\mathscr{G}}\sqrt{\left\lvert
g\right\rvert}\left\lVert\bm{x}_{g}\right\rVert_{2}}$
$\displaystyle=\sup_{\bm{x}\neq\bm{0},\sum_{g\in\mathscr{G}}\sqrt{\left\lvert
g\right\rvert}\left\lVert\bm{x}_{g}\right\rVert_{2}=1}\left\lvert\sum_{g\in\mathscr{G}}\bm{c}_{g}^{\prime}\bm{x}_{g}\right\rvert$
$\displaystyle\leq\sup_{\bm{x}\neq\bm{0},\sum_{g\in\mathscr{G}}\sqrt{\left\lvert
g\right\rvert}\left\lVert\bm{x}_{g}\right\rVert_{2}=1}\sum_{g\in\mathscr{G}}\left\lVert\frac{\bm{c}_{g}}{\sqrt{\left\lvert
g\right\rvert}}\right\rVert_{2}\left\lVert\sqrt{\left\lvert
g\right\rvert}\bm{x}_{g}\right\rVert_{2}\leq\max_{g\in\mathscr{G}}\frac{\left\lVert\bm{c}_{g}\right\rVert_{2}}{\sqrt{\left\lvert
g\right\rvert}}.$

We now start the actual proof. Recall
$\hat{\bm{\sigma}}_{h}=\operatorname{\mathrm{vec}}\left(\mathscr{B}_{h}\big{(}\hat{\bm{\varSigma}}_{1}\big{)}^{\prime}\right)$,
$\hat{\bm{V}}_{h}=\begin{bmatrix}\mathscr{B}_{h}\big{(}\hat{\bm{\varSigma}}_{1}\big{)}^{\prime}&\mathscr{B}_{h}\big{(}\hat{\bm{\varSigma}}_{0}\big{)}\end{bmatrix}$,
$\bm{V}=\begin{bmatrix}\bm{\varSigma}_{1}^{\prime}&\bm{\varSigma}_{0}\end{bmatrix}$,
and $\bm{C}=\begin{bmatrix}\bm{A}&\bm{B}\end{bmatrix}$. Then,

$\displaystyle\hat{\bm{\sigma}}_{h}-\hat{\bm{V}}^{(d)}_{h}\bm{c}$
$\displaystyle=\operatorname{\mathrm{vec}}\left(\mathscr{B}_{h}\big{(}\hat{\bm{\varSigma}}_{1}\big{)}^{\prime}-\hat{\bm{V}}_{h}\bm{C}^{\prime}\right)=\operatorname{\mathrm{vec}}\left(\left[\mathscr{B}_{h}\big{(}\hat{\bm{\varSigma}}_{1}\big{)}-\bm{\varSigma}_{1}\right]^{\prime}-\left[\hat{\bm{V}}_{h}-\bm{V}\right]\bm{C}^{\prime}+\underbrace{\left[\bm{\varSigma}_{1}^{\prime}-\bm{V}\bm{C}^{\prime}\right]}_{=\mathbf{O},\text{
see }\eqref{eq:YWpopulation}}\right)$ (21)
$\displaystyle=\underbrace{\operatorname{\mathrm{vec}}\left(\mathscr{B}_{h}\big{(}\hat{\bm{\varSigma}}_{1}\big{)}^{\prime}-\bm{\varSigma}_{1}^{\prime}\right)}_{\hat{\bm{\varDelta}}_{\Sigma}}-\underbrace{\left[\hat{\bm{V}}_{h}^{(d)}-\bm{V}^{(d)}\right]}_{\hat{\bm{\varDelta}}_{V}}\bm{c}=\hat{\bm{\varDelta}}_{\Sigma}-\hat{\bm{\varDelta}}_{V}\bm{c},$

Using (21), we rewrite

$\displaystyle\left\lVert\hat{\bm{\sigma}}_{h}-\hat{\bm{V}}_{h}^{(d)}\hat{\bm{c}}\right\rVert_{2}^{2}$
$\displaystyle=\left\lVert\Big{[}\hat{\bm{\sigma}}_{h}-\hat{\bm{V}}_{h}^{(d)}\bm{c}\Big{]}-\Big{[}\hat{\bm{V}}_{h}^{(d)}(\hat{\bm{c}}-\bm{c})\Big{]}\right\rVert_{2}^{2}$
$\displaystyle=\left\lVert\hat{\bm{\sigma}}_{h}-\hat{\bm{V}}_{h}^{(d)}\bm{c}\right\rVert_{2}^{2}+\left\lVert\hat{\bm{V}}_{h}^{(d)}(\hat{\bm{c}}-\bm{c})\right\rVert_{2}^{2}-2\left(\hat{\bm{c}}-\bm{c}\right)^{\prime}\hat{\bm{V}}_{h}^{(d)\prime}\left(\hat{\bm{\varDelta}}_{\Sigma}-\hat{\bm{\varDelta}}_{V}\bm{c}\right).$

Recalling the objective function
$\mathscr{L}_{\alpha}(\bm{c};\lambda)=\left\lVert\hat{\bm{\sigma}}-\hat{\bm{V}}_{h}^{(d)}\bm{c}\right\rVert_{2}^{2}+\lambda
P_{\alpha}(\bm{c})$ and noting that
$\mathscr{L}_{\alpha}(\hat{\bm{c}};\lambda)\leq\mathscr{L}_{\alpha}(\bm{c};\lambda)$
by construction, it follows that

$\displaystyle\left\lVert\hat{\bm{V}}_{h}^{(d)}(\hat{\bm{c}}-\bm{c})\right\rVert_{2}^{2}+\lambda
P_{\alpha}\left(\hat{\bm{c}}\right)$ $\displaystyle\leq
2\left(\hat{\bm{c}}-\bm{c}\right)^{\prime}\hat{\bm{V}}_{h}^{(d)\prime}\left(\hat{\bm{\varDelta}}_{\Sigma}-\hat{\bm{\varDelta}}_{V}\bm{c}\right)+\lambda
P_{\alpha}\left(\bm{c}\right)$ (22) $\displaystyle\leq
2P_{\alpha}\left(\hat{\bm{c}}-\bm{c}\right)P_{\alpha}^{*}\left(\hat{\bm{V}}_{h}^{(d)\prime}\left(\hat{\bm{\varDelta}}_{\Sigma}-\hat{\bm{\varDelta}}_{V}\bm{c}\right)\right)+\lambda
P_{\alpha}\left(\bm{c}\right)$ $\displaystyle\leq
P_{\alpha}\left(\hat{\bm{c}}-\bm{c}\right)\left[2P_{\alpha}^{*}\left(\hat{\bm{V}}_{h}^{(d)\prime}\hat{\bm{\varDelta}}_{\Sigma}\right)+2P_{\alpha}^{*}\left(\hat{\bm{V}}_{h}^{(d)\prime}\hat{\bm{\varDelta}}_{V}\bm{c}\right)\right]+\lambda
P_{\alpha}\left(\bm{c}\right),$

where we used the dual-norm inequality (see (19)) and the triangle property of
(dual) norms in the second and third inequality, respectively. Define the sets

$\mathscr{H}_{1}(x)=\left\\{2P_{\alpha}^{*}\left(\hat{\bm{V}}_{h}^{(d)\prime}\hat{\bm{\varDelta}}_{\Sigma}\right)\leq
x\right\\}\qquad\text{ and
}\qquad\mathscr{H}_{2}(x)=\left\\{2P_{\alpha}^{*}\left(\hat{\bm{V}}_{h}^{(d)\prime}\hat{\bm{\varDelta}}_{V}\bm{c}\right)\leq
x\right\\}.$ (23)

On the set
$\mathscr{H}_{1}(\frac{\lambda}{4})\cap\mathscr{H}_{2}(\frac{\lambda}{4})$, we
can scale (22) by a factor 2 to obtain

$2\left\lVert\hat{\bm{V}}_{h}^{(d)}(\hat{\bm{c}}-\bm{c})\right\rVert_{2}^{2}+2\lambda
P_{\alpha}\left(\hat{\bm{c}}\right)\leq\lambda
P_{\alpha}(\hat{\bm{c}}-\bm{c})+2\lambda P_{\alpha}(\bm{c}).$ (24)

We subsequently manipulate $P_{\alpha}(\hat{\bm{c}})$ and
$P_{\alpha}\left(\hat{\bm{c}}-\bm{c}\right)$. Using the reverse triangle
inequality, we have

$\displaystyle
P_{\alpha}(\hat{\bm{c}})=\alpha\sum_{g\in\mathscr{G}}\sqrt{\left\lvert
g\right\rvert}\left\lVert\hat{\bm{c}}_{g}\right\rVert_{2}+(1-\alpha)\left\lVert\hat{\bm{c}}\right\rVert_{1}$
(25) $\displaystyle\geq\alpha\sum_{g\in\mathscr{G}_{S}}\sqrt{\left\lvert
g\right\rvert}\Big{[}\left\lVert\bm{c}_{g}\right\rVert_{2}-\left\lVert\hat{\bm{c}}_{g}-\bm{c}_{g}\right\rVert_{2}\Big{]}+\alpha\sum_{g\in\mathscr{G}_{S^{c}}}\sqrt{\left\lvert
g\right\rvert}\left\lVert\hat{\bm{c}}_{g}\right\rVert_{2}+(1-\alpha)\Big{[}\left\lVert\bm{c}_{S}\right\rVert_{1}-\left\lVert\hat{\bm{c}}_{S}-\bm{c}_{S}\right\rVert_{1}+\left\lVert\hat{\bm{c}}_{S^{c}}\right\rVert_{1}\Big{]}$
$\displaystyle=P_{\alpha,S}(\bm{c})+P_{\alpha,S^{c}}\left(\hat{\bm{c}}-\bm{c}\right)-P_{\alpha,S}\left(\hat{\bm{c}}-\bm{c}\right),$

where $\mathscr{G}_{S}$ and $\mathscr{G}_{S^{c}}$ are defined in Lemma 1.
Simple rewriting provides

$\displaystyle P_{\alpha}\left(\hat{\bm{c}}-\bm{c}\right)$
$\displaystyle=\alpha\sum_{g\in\mathscr{G}}\sqrt{\left\lvert
g\right\rvert}\left\lVert\hat{\bm{c}}_{g}-\bm{c}_{g}\right\rVert_{2}+(1-\alpha)\left\lVert\hat{\bm{c}}_{S}-\bm{c}_{S}\right\rVert_{1}+\alpha\sum_{g\in\mathscr{G}_{S^{c}}}\sqrt{\left\lvert
g\right\rvert}\left\lVert\hat{\bm{c}}_{g}\right\rVert_{2}+(1-\alpha)\left\lVert\hat{\bm{c}}_{S^{c}}\right\rVert_{1}$
(26)
$\displaystyle=P_{\alpha,S}\left(\hat{\bm{c}}-\bm{c}\right)+P_{\alpha,S^{c}}\left(\hat{\bm{c}}-\bm{c}\right).$

Combining results (24)–(26) yields

$2\left\lVert\hat{\bm{V}}_{h}^{(d)}(\hat{\bm{c}}-\bm{c})\right\rVert_{2}^{2}+\lambda
P_{\alpha,S^{c}}\left(\hat{\bm{c}}-\bm{c}\right)\leq 3\lambda
P_{\alpha,S}(\hat{\bm{c}}-\bm{c}),$ (27)

and $\hat{\bm{c}}-\bm{c}$ is thus a member of the set
$\mathscr{C}_{N_{c}}(\mathscr{G},S)$ as defined in Lemma 1. In combination
with Lemma 2, thus requiring
$\mathscr{H}_{1}(\frac{\lambda}{4})\cap\mathscr{H}_{2}(\frac{\lambda}{4})\cap\mathscr{V}\left(\frac{\phi_{0}}{2}\right)$
to hold, we conclude

$\displaystyle
2\left\lVert\hat{\bm{V}}_{h}^{(d)}(\hat{\bm{c}}-\bm{c})\right\rVert_{2}^{2}$
$\displaystyle+\lambda
P_{\alpha}(\hat{\bm{c}}-\bm{c})=2\left\lVert\hat{\bm{V}}_{h}^{(d)}(\hat{\bm{c}}-\bm{c})\right\rVert_{2}^{2}+\lambda
P_{\alpha,S}(\hat{\bm{c}}-\bm{c})+\lambda
P_{\alpha,S^{c}}(\hat{\bm{c}}-\bm{c})$ (28)
$\displaystyle\overset{(i)}{\leq}4\lambda
P_{\alpha,S}(\hat{\bm{c}}-\bm{c})\overset{(ii)}{\leq}16\left\lVert\hat{\bm{V}}_{h}^{(d)}(\hat{\bm{c}}-\bm{c})\right\rVert_{2}\left(\frac{\bar{\omega}_{\alpha}\lambda}{\phi_{0}}\right)\overset{(iii)}{\leq}\left\lVert\hat{\bm{V}}_{h}^{(d)}(\hat{\bm{c}}-\bm{c})\right\rVert_{2}^{2}+\frac{64\bar{\omega}_{\alpha}^{2}\lambda^{2}}{\phi_{0}^{2}},$

where step (i) follows from (27), step (ii) is implied by
$P_{\alpha,S}(\hat{\bm{c}}-\bm{c})\leq\frac{4\bar{\omega}_{0}\left\lVert\hat{\bm{V}}_{h}^{(d)}(\hat{\bm{c}}-\bm{c})\right\rVert_{2}}{\phi_{0}}$
for $\hat{\bm{c}}-\bm{c}\in\mathscr{C}_{N_{c}}(\mathscr{G},S)$ (Lemma 2), and
step (iii) uses the elementary inequality $16uv\leq u^{2}+64v^{2}$ (i.e.
manipulating $(u-8v)^{2}\geq 0$). A straightforward rearrangement of (28)
provides the inequality of Theorem 2.

It remains to determine a lower bound on the probability of
$\mathscr{H}_{1}(\frac{\lambda}{4})\cap\mathscr{H}_{2}(\frac{\lambda}{4})\cap\mathscr{V}\left(\frac{\phi_{0}}{2}\right)$.
We rely on the elementary inequality
$\operatorname{\mathbb{P}}\left(\mathscr{H}_{1}(\frac{\lambda}{4})\cap\mathscr{H}_{2}(\frac{\lambda}{4})\cap\mathscr{V}\left(\frac{\phi_{0}}{2}\right)\right)\geq
1-\operatorname{\mathbb{P}}\left(\mathscr{H}_{1}(\frac{\lambda}{4})^{c}\right)-\operatorname{\mathbb{P}}\left(\mathscr{H}_{2}(\frac{\lambda}{4})^{c}\right)-\operatorname{\mathbb{P}}\left(\mathscr{V}\left(\frac{\phi_{0}}{2}\right)^{c}\right)$
to bound the individual probabilities.

We start with
$\operatorname{\mathbb{P}}\left(\mathscr{H}_{1}(\frac{\lambda}{4})^{c}\right)=\operatorname{\mathbb{P}}\left(2P_{\alpha}^{*}\left(\hat{\bm{V}}_{h}^{(d)\prime}\hat{\Delta}_{\Sigma}\right)>\frac{\lambda}{4}\right)\leq\operatorname{\mathbb{P}}\left(\left\lVert\hat{\bm{V}}_{h}^{(d)\prime}\hat{\Delta}_{\Sigma}\right\rVert_{\infty}>\frac{\lambda}{8}\right)$.
The last inequality is true because continuing from (20), we have

$P_{\alpha}^{*}(\bm{c})\leq\alpha\max_{g\in\mathscr{G}}\frac{\left\lVert\bm{c}_{g}\right\rVert_{2}}{\sqrt{\left\lvert
g\right\rvert}}+(1-\alpha)\left\lVert\bm{c}\right\rVert_{\infty}\leq\alpha\max_{g\in\mathscr{G}}\frac{\sqrt{|g|}\left\lVert\bm{c}_{g}\right\rVert_{\infty}}{\sqrt{\left\lvert
g\right\rvert}}+(1-\alpha)\left\lVert\bm{c}\right\rVert_{\infty}=\left\lVert\bm{c}\right\rVert_{\infty}$
(29)

for any vector $\bm{c}$. Subsequently, we have

$\displaystyle\operatorname{\mathbb{P}}$
$\displaystyle\left(\left\lVert\hat{\bm{V}}_{h}^{(d)\prime}\hat{\Delta}_{\Sigma}\right\rVert_{\infty}>\frac{\lambda}{8}\right)=\operatorname{\mathbb{P}}\left(\left\lVert\left[(\hat{\bm{V}}_{h}^{(d)}-\bm{V}^{(d)})+\bm{V}^{(d)}\right]^{\prime}\hat{\bm{\varDelta}}_{\Sigma}\right\rVert_{\infty}>\frac{\lambda}{8}\right)$
(30)
$\displaystyle\leq\operatorname{\mathbb{P}}\left(\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{1}\left\lVert\hat{\bm{\varDelta}}_{\Sigma}\right\rVert_{\infty}+C_{V}\left\lVert\hat{\bm{\varDelta}}_{\Sigma}\right\rVert_{\infty}>\frac{\lambda}{8}\right)\leq\operatorname{\mathbb{P}}\left(\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{\vdash}^{2}+C_{V}\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{\vdash}>\frac{\lambda}{8}\right)$
$\displaystyle\leq\operatorname{\mathbb{P}}\left(\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{\vdash}>\frac{\lambda^{1/2}}{4}\right)+\operatorname{\mathbb{P}}\left(\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{\vdash}>\frac{\lambda}{16C_{V}}\right),$

exploiting block-diagonality of $\hat{\bm{V}}_{h}^{(d)}-\bm{V}^{(d)}$ and
$\bm{V}^{(d)}$ such that
$\left\lVert\hat{\bm{V}}_{h}^{(d)}-\bm{V}^{(d)}\right\rVert_{1}=\max_{1\leq
i\leq
N}\left\lVert\hat{\bm{V}}_{i,h}-\bm{V}_{i}\right\rVert_{1}\leq\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{\vdash}$
and
$\left\lVert\bm{V}^{(d)}\right\rVert_{1}\leq\left\lVert\bm{V}\right\rVert_{\vdash}\leq
C_{V}$ (Assumption X). Bounds for the final RHS terms in (30) are available
from Theorem 1.

Second, we have

$\displaystyle\operatorname{\mathbb{P}}\left(\mathscr{H}_{2}\left(\frac{\lambda}{4}\right)^{c}\right)$
$\displaystyle\leq\operatorname{\mathbb{P}}\left(\left\lVert\hat{\bm{V}}_{h}^{(d)\prime}\hat{\bm{\varDelta}}_{V}\bm{c}\right\rVert_{\infty}>\frac{\lambda}{8}\right)\leq\operatorname{\mathbb{P}}\left(\left\lVert\hat{\bm{V}}_{h}^{(d)\prime}\left[\hat{\bm{V}}_{h}^{(d)}-\bm{V}^{(d)}\right]\right\rVert_{\infty}>\frac{\lambda}{8\left\lVert\bm{c}\right\rVert_{\infty}}\right)$
(31)
$\displaystyle\leq\operatorname{\mathbb{P}}\left(\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{1}\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{\infty}+C_{V}\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{\infty}>\frac{\lambda}{8\left\lVert\bm{c}\right\rVert_{\infty}}\right)$
$\displaystyle\leq\operatorname{\mathbb{P}}\left(\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{\vdash}>\frac{\lambda^{1/2}}{4}\right)+\operatorname{\mathbb{P}}\left(\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{\vdash}>\frac{\lambda}{16C_{V}}\right),$

where the last line relies on the union bound and the fact that
$\left\lVert\bm{c}\right\rVert_{\infty}<1$ (implied by Assumption 1). Hence,
the sets $\mathscr{H}_{1}\left(\frac{\lambda}{4}\right)^{c}$ and
$\mathscr{H}_{2}\left(\frac{\lambda}{4}\right)^{c}$ admit the same probability
bound.

Finally,
$\operatorname{\mathbb{P}}\left(\mathscr{V}\left(\frac{\phi_{0}}{2}\right)^{c}\right)=\operatorname{\mathbb{P}}\left(\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{1}>\frac{\phi_{0}}{2}\right)\leq\operatorname{\mathbb{P}}\left(\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{\vdash}>\frac{\phi_{0}}{2}\right)$.
Combining all previous results, we conclude

$\displaystyle\operatorname{\mathbb{P}}$
$\displaystyle\left(\mathscr{H}_{1}\left(\frac{\lambda}{4}\right)\cap\mathscr{H}_{2}\left(\frac{\lambda}{4}\right)\cap\mathscr{V}\left(\frac{\phi_{0}}{2}\right)\right)$
(32) $\displaystyle\geq
1-2\operatorname{\mathbb{P}}\left(\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{\vdash}>\frac{\lambda^{1/2}}{4}\right)-2\operatorname{\mathbb{P}}\left(\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{\vdash}>\frac{\lambda}{16C_{v}}\right)-\operatorname{\mathbb{P}}\left(\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{\vdash}>\frac{\phi_{0}}{2}\right)$
$\displaystyle\geq
1-5\operatorname{\mathbb{P}}\left(\left\lVert\hat{\bm{V}}_{h}-\bm{V}\right\rVert_{\vdash}>6f(\lambda,\phi_{0})\right),$

where
$f(\lambda,\phi_{0})=\min\left(\frac{\lambda^{1/2}}{24},\frac{\lambda}{96C_{v}},\frac{\phi_{0}}{12}\right)$.
The proof is completed by evaluating the final probability in (32) with the
use of Theorem 1.

∎

###### Proof of Corollary 1.

First, we derive the conditions under which the set on which the performance
bound in Theorem 2 holds occurs with probability converging to one. Under
Assumption 2(b1), along with the remaining assumptions in Theorem 2, this
probability is given by $1-\mathscr{P}_{1}(f(\lambda,\phi_{0}),N,T)$, where we
recall from Theorem 1 that

$\begin{split}\mathscr{P}_{1}(f(\lambda,\phi_{0}),N,T)&=N^{2}\left(b_{1}T^{(1-\delta)/3}+\frac{[2h(f(\lambda,\phi_{0}),N)+1]b_{3}}{f(\lambda,\phi_{0})}\right)\exp\left(-\frac{T^{(1-\delta)/3}}{2b_{1}^{2}}\right)\\\
&\qquad+\frac{b_{2}N^{2}[2h(f(\lambda,\phi_{0}),N)+1]^{d}}{f(\lambda,\phi_{0})^{d}T^{\frac{\delta}{2}(d-1)}}\\\
\end{split}$ (33)

Given that $\lambda\in O\left(T^{-q_{\lambda}}\right)$ with $q_{\lambda}>0$,
it follows immediately that

$f(\lambda,\phi_{0})=\min\left(\frac{\lambda^{1/2}}{24},\frac{\lambda}{96C_{v}},\frac{\phi_{0}}{12}\right)=O\left(T^{-q_{\lambda}}\right).$
(34)

In addition, following the remark below Theorem 1, it holds that

$h\left(\frac{\lambda}{96C_{v}},N\right)=O\left(\log\left(\frac{N}{\lambda}\right)^{2}k_{0}\right)=O\left(\log(T)^{2}T^{q_{k}}\right).$
(35)

Based on (34) and (35), it follows that the first RHS-term in (33) converges
to zero exponentially in $T$ for any $\delta<1$. The second RHS-term, however,
converges to zero at most at a polynomial rate. Accordingly,

$\begin{split}\mathscr{P}_{1}(f(\lambda,\phi_{0}),N,T)=O\left(\frac{N^{2}[h\left(\frac{\lambda}{96C_{v}},N\right)+1]^{d}}{\lambda^{d}T^{\frac{\delta}{2}(d-1)}}\right)=O\left(\log(T)^{2d}T^{2q_{N}+dq_{k}+dq_{\lambda}-\frac{\delta(d-1)}{2}}\right),\end{split}$
(36)

where the second equality holds from the observation that for any $\delta<1$
the first two RHS terms in the first equality converge to zero at an
exponential rate, whereas the third term may converge to zero at most at a
polynomial rate. From (36), it follows that
$\mathscr{P}_{1}(f(\lambda,\phi_{0}),N,T)\to 0$ if

$2q_{n}+dq_{k}+dq_{\lambda}-\frac{\delta(d-1)}{2}<0\Rightarrow
q_{\lambda}<-\frac{2q_{N}}{d}-q_{k}+\frac{\delta(d-1)}{2d}.$

In a similar fashion, we derive the conditions under which
$\mathscr{P}_{2}(f(\lambda,\phi_{0})\to 0$, by noting that

$\begin{split}\mathscr{P}_{2}(f(\lambda,\phi_{0}),N,T)=&N^{2}\left[\frac{\beta_{1}[2h(f(\lambda,\phi_{0}),N)+1]}{f(\lambda,\phi_{0})}+\frac{2}{\beta_{2}}\left(\frac{Tf(\lambda,\phi_{0})^{2}}{[2h(x,N)+1]^{2}}\right)^{\frac{1}{7}}\right]\times\\\
&\exp\left(-\frac{1}{\beta_{3}}\left(\frac{Tf(\lambda,\phi_{0})^{2}}{[2h(f(\lambda,\phi_{0}),N)+1]^{2}}\right)^{\frac{1}{7}}\right)\end{split}$
(37)

converges to zero exponentially fast in $T$ if
$\frac{Tf(\lambda,\phi_{0})^{2}}{[2h(f(\lambda,\phi_{0}),N)+1]^{2}}$ diverges
at a polynomial rate in $T$. Making use of (34) and (35), it follows that

$\frac{Tf(\lambda,\phi_{0})^{2}}{[2h(f(\lambda,\phi_{0}),N)+1]^{2}}=O\left(\log(T)^{-4}T^{1-2q_{\lambda}-2q_{k}}\right),$

which translates to the condition

$1-2q_{\lambda}-2q_{k}>0\Rightarrow q_{\lambda}<\frac{1}{2}-q_{k}.$

This establishes conditions (i) and (ii) in Corollary 1.

We proceed by deriving the order of the performance bound in Theorem 2. Noting
that

$\sum_{g\in\mathscr{G}_{S}}\sqrt{\left\lvert
g\right\rvert}\leq\left\lvert\mathscr{G}_{S}\right\rvert\max_{g\in\mathscr{G}_{S}}\sqrt{\left\lvert
g\right\rvert}=O\left(T^{q_{g}+q_{N}/2}\right),$

it follows that

$\bar{\omega}=O\left((1-\alpha)T^{q_{g}+q_{N}/2}+\alpha T^{q_{S}/2}\right).$

Then, by Theorem 2,

$\left\lVert\hat{\bm{V}}_{h}^{(d)}(\hat{\bm{c}}-\bm{c})\right\rVert_{2}^{2}\leq\frac{4\bar{\omega}_{\alpha}^{2}\lambda^{2}}{\phi_{0}^{2}}=O\left((1-\alpha)T^{2q_{g}+q_{N}-2q_{\lambda}}+\alpha
T^{q_{S}-2q_{\lambda}}\right)$

and

$(1-\alpha)\sum_{g\in\mathscr{G}}\sqrt{\left\lvert
g\right\rvert}\left\lVert\hat{\bm{c}}_{g}-\bm{c}_{g}\right\rVert_{2}+\alpha\left\lVert\hat{\bm{c}}-\bm{c}\right\rVert_{1}\leq\frac{4\bar{\omega}_{\alpha}^{2}\lambda}{\phi_{0}^{2}}=O\left((1-\alpha)T^{2q_{g}+q_{N}-q_{\lambda}}+\alpha
T^{q_{S}-q_{\lambda}}\right),$

on a set with probability $1-\mathscr{P}_{1}(f(\lambda,\phi_{0}),N,T)$ or
$1-\mathscr{P}_{2}(f(\lambda,\phi_{0}),N,T)$, depending on whether Assumption
2(b1) or 2(b2) applies, respectively. Since we have shown that both
$\mathscr{P}_{1}(f(\lambda,\phi_{0}),N,T)\to 0$ and
$\mathscr{P}_{2}(f(\lambda,\phi_{0}),N,T)\to 0$ under the conditions imposed
in Corollary 1, the proof is complete.

∎

##  Appendix C Empirical Application: Additional Material

###  C.1 Satellite data

$y_{1}$$y_{2}$$y_{3}$$y_{4}$$y_{5}$$y_{6}$$y_{7}$$y_{8}$$y_{9}$$y_{10}$$y_{11}$$y_{12}$$y_{13}$$y_{14}$$y_{15}$$y_{16}$$y_{17}$$y_{18}$$y_{19}$$y_{20}$$y_{21}$$y_{22}$$y_{23}$$y_{24}$$y_{25}$$y_{26}$$y_{27}$$y_{28}$$y_{29}$$y_{30}$$y_{31}$$y_{32}$$y_{33}$$y_{34}$$y_{35}$$y_{36}$$y_{37}$$y_{38}$$y_{39}$$y_{40}$$y_{41}$$y_{42}$$y_{43}$$y_{44}$$y_{45}$
Figure 5: Vertically ordered spatial grid. The arrows display part of the
admissible interactions for spatial units $y_{1}$ and $y_{23}$ under the
identification constraint that spatial units $y_{i},y_{j}$ may only interact
if $\left\lvert i-j\right\rvert<\lfloor N/4\rfloor=11$. Figure 6: The
coefficient sizes estimated by SPLASH(0,$\lambda_{\text{CV}}$). While there is
no sparsity within the admissible bandwidth, the diagonal that were estimated
zero by $\lambda_{\text{AIC}}$ are estimated close to zero.