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

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

Introduction

Spatio-temporal models are powerful tools to explain and exploit dependencies between variables that are observed over both time and space, but they come with a number of challenges. In particular, endogeneity issues arise because the contemporaneous observations occur on both sides of the model equation. Furthermore, the inclusion of both spatial and temporal lags quickly results in heavily parameterized models. To circumvent these issues, a large part of the literature incorporates predetermined spatial weight matrices that govern the contemporaneous interactions between spatial units. Examples of this modelling strategies are: the spatial autoregressive model with a Gaussian quasi-maximum likelihood estimator (QMLE) by lee2004; the QMLE estimation of stationary spatial panels with fixed effects detailed in yudejonglee2008; the extension of these spatial panels to include spatially autoregressive disturbances as in leeyu2010; a further extension to a non-stationary setting in which units can be spatially cointegrated in Yu2012; and the computationally beneficial generalized method of moments (GMM) estimator by Lee2014. While the choice of the spatial weight matrix is a key element of the model specification, its selection process can feel somewhat arbitrary and/or tedious. The arbitrariness might prevail when practical considerations fail to suggest a particular mechanism for the spatial interactions. Accordingly, more recent literature focuses on either incorporating multiple weight matrices DebarsyLeSage2018,zhangyu2018 or, at the expense of estimating many parameters, directly inferring all spatial interactions from the data lamsouza2019,Gao2019,MaGuoWang2021. We contribute to the latter strand of literature with the development of a new estimator that provides several unique benefits.

In this paper, we propose the SPatial LAsso-type SHrinkage (SPLASH) estimator as a fully data-driven estimator of spatio-temporal interactions. Apart from a 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. Building upon previous works by DouParrellaYao2016 and Gao2019, we resolve the endogeneity problem in our spatio-temporal regressions by estimating the generalized Yule-Walker equations. Our contributions are five-fold. 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 finitely bounded bandwidth and unstructured sparsity, it follows 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. A tailored regularization procedure reminiscent of the sparse group lasso is proposed. Fourth, we generalize the model to include exogenous variable and demonstrate that an extended system of Yule-Walker equations continues to provide consistent estimators. Our simulation study confirms these points. Fifth and final, we employ SPLASH to predict $\text{NO}_2$ concentrations in London from 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 is an aggregation 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 Gao2019, 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. lamsouza2014 consider a model specification where the spatial units depend linearly on a spatial lag and 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. Ahrens2015 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 Jing2003. 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, lamsouza2019 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, lamsouza2019 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.

Closest related to the our work are Gao2019, and MaGuoWang2021. Both papers consider the same model and an estimation procedure that relies on generalized Yule-Walker equations. The key difference with our paper lies in the method by which the model complexity is controlled during estimation. Gao2019 assume the coefficient matrices to be banded with a bandwidth that is small compared 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, apart from a generous upper bound on the bandwidth to ensure identification, SPLASH does not require an a priori choice regarding the bandwidth. The recently developed bagging approach in MaGuoWang2021 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 (ref) 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 (ref). The simulation results in Section (ref) and the empirical application in Section (ref) demonstrate the benefits of the SPLASH estimator. Section (ref) 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\lvertx_i\right\rvert$ as an important special case. The total number elements in $\bm x$ is denoted by $\left\lvert\bm x\right\rvert$ and the number of non-zero elements in $\bm x$ is denoted by $\mathcal{M}(\bm x) = \sum_{i=1}^N \mathbbm{1}\lbrace x_i \neq 0\rbrace$. 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 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\lverta_{ij}\right\rvert$, the spectral norm $\left\lVert\bm A\right\rVert_2 = \big[\lambda_{\max}(\bm A'\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\lverta_{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\lverta_{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 \lbrace 1, \ldots, N \rbrace$ denote an index set with cardinality $\left\lvertS\right\rvert$. Then, $\bm x_S$ denotes the $\left\lvertS\right\rvert$-dimensional vector with the elements of $\bm x$ indexed by $S$, whereas $\bm A_S$ denotes the $(M \times \left\lvertS\right\rvert)$-dimensional matrix containing the columns of $\bm A$ indexed by $S$. In addition, we define $\mathcal{D}_A(k) = \left\lbrace a_{ij} \ \vert \ \left\lverti-j\right\rvert = k\right\rbrace$ 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.

The Spatio-temporal Vector Autoregression

As in the recent paper by Gao2019, we consider the spatio-temporal vector autoregression

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 (ref).

Assumption (ref) 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 (ref)(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. lee2004 and leeyu2010). Typically, the product $\lambda \bm W_N$ -- the natural counterpart of the matrix $\bm A$ -- is required to fulfil conditions similar to $\left\lVert\bm A\right\rVert_{\vdash}\leq \delta_A <1$.\footnote{For 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 lee2004. 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 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 $\normoneinf \bm D \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 DouParrellaYao2016. Second, Assumption (ref)(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 (ref), are closely related to those in Masini2019. Assumption (ref)(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}')=\boldsymbol{0}$ while the mixing assumption controls the serial correlation in the data. Polynomial or exponential tail decay of the distribution of the innovations is imposed through either Assumption (ref)(b1) or Assumption (ref)(b2), respectively. The type of tail decay will directly influence the growth rates we can allow for $N$ and $T$. The discussions in Masini2019 demonstrate that Assumption (ref) 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 OLS estimation 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 Brockwell1991.

There are several possibilities to introduce structure into $\bm A$ and $\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. DouParrellaYao2016 consider a more general setting in which each row of $\bm W_N$ receives its own spatial autoregressive 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',\bm \lambda_1',\bm \lambda_2')'$. Gao2019 require $\bm A$ and $\bm B$ to be banded matrices.\footnote{The matrix $\bm A$ has bandwidth $k$ if the total number of nonzero entries in any row or column is at most $k$.} We employ a similar assumption.

Assumption (ref) 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 (ref)(a) is key in identifying the parameters. Our discussions in Section (ref) illustrate that this assumption is realistic when the data is observed on a regular grid. The combination of Assumptions (ref)(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}')$. 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 (ref) shows that banded estimators for $\bm \varSigma_0$ and $\bm \varSigma_1$ provide an accurate approximation to $\bm V = \big[ \bm \varSigma_1' \; \bm \varSigma_0 \big]'$. Each of these banded matrices has at most $2h(\epsilon)+1$ nonzero elements in their columns/rows. In other words, given $\epsilon$, $l_0$ and $k_0$, Assumptions (ref)--(ref) guarantee that $\bm \varSigma_0$ and $\bm \varSigma_1$ can be well-approximated by matrices with bandwidths smaller than $N$. This improves the convergence rate of our estimator.

Sparse Estimation

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

Even under Assumption (ref), the number of unknown parameters in $\bm A$ and $\bm B$ continues to grow quadratically in $N$. For large $N$, the accurate estimation of all these parameters becomes infeasible rather quickly. To alleviate this curse of dimensionality, we rely on sparsity. Sparsity naturally occurs when two spatial units do not interact with each other. We demonstrate, however, that a special, and exploitable, sparsity pattern arises whenever the spatial units are ordered in a structured way.

As an illustrative example, let us consider repeated measurements on the $(5\times5)$ spatial grids shown in the left column of Figure (ref). 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/physical interactions to be most pronounced at short length scales. In Figure (ref)(a) we start 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}$. Such interactions 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 (ref)(b). The blue elements are potentially nonzero whereas uncolored elements are zero. The nonzero elements in $\bm A$ and $\bm B$ are seen to cluster in specific, dense diagonals with the occasional zero when horizontal neighbours are absent (at the boundary of the grid). This diagonal sparsity pattern is not an artifact of allowing horizontal interactions only. Figures (ref)(c) adds the vertical interactions and the accompanying sparsity pattern again manifests itself along diagonals (Figure (ref)(d)). Finally, with diagonal nearest neighbours being horizontal neighbors of vertical elements, we observe a “thickening” of the diagonals in Figure (ref)(f). Guided by these considerations we combine generalized Yule-Walker estimation with a sparse group penalty Simon2013. The Yule-Walker estimator will control for endogeneity, while the sparse group penalty will shrink towards diagonal structures by including/omitting complete diagonals and thus selecting the required interactions. Compared to Gao2019, 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}'$ and taking expectations, we find $\bm \varSigma_1=\bm A \bm \varSigma_1 + \bm B \bm \varSigma_0$ or, equivalently,

The $i$\textsuperscript{th} column of $\bm C'$ contains all coefficients that belong to the $i$\textsuperscript{th} equation in (ref). Assumption (ref) 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$\textsuperscript{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' \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, Gao2019 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'\;\hat\bm \varSigma_0 \big]$. Motivated by $\bm \sigma_i - \bm V_i \bm c_i=\boldsymbol{0}$, they define their estimator $\hat{\bm c}^{GMWY}_i$ as the following minimizer:

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 (ref). 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(\mathcal{B}_{h_{}}\big(\hat\bm \varSigma_1\big)' \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 (ref) for an illustration).\footnote{In general, all quantities derived from 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 (ref) over all equations. Finally, we construct the penalty function. We define an index set that partitions the vector $\bm c$ into sub-vectors, denoted $\lbrace\bm c_g\rbrace$, that contain the non-zero diagonals of $\bm A$ and $\bm B$ that are admissible under Assumption (ref) as

respectively, and let $\mathcal{G} = \mathcal{G}_A \cup \mathcal{G}_B$. Based on this notation, we define our objective function as

The spatial lasso-type shrinkage estimator, abbreviated SPLASH($\alpha$,$\lambda$) or SPLASH in short, is defined as the minimizer of (ref), i.e. $\hat{\bm c}=\operatorname*{arg~min}_\bm c \mathcal{L}_\alpha(\bm c;\lambda)$. The importance of the penalty function $P_\alpha(\bm c)$ is governed by the penalty parameter $\lambda$ and the 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 (ref)(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 Simon2013. An R/C++ implementation of the SPLASH estimator based on this algorithm is available on one of the author's website.\footnote{\url{https://sites.google.com/view/etiennewijler/code}}.

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 (ref) are uniquely identified. To this end, we leverage the bandedness assumption in Assumption (ref), 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 (ref) 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 Bickel2009, 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 (ref) fits our framework particularly well, as the assumed maximum bandwidth of the matrices $\bm A$ and $\bm B$ in Assumption (ref) 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 (ref) of Appendix (ref) that a Sparse Group Lasso compatibility condition is implied by Assumption (ref).

Equipped with Assumption (ref), we find the following finite-sample performance bounds on the prediction and estimation error of SPLASH.

Theorem (ref) contains a finite-sample performance bound on the prediction and estimation error for the SPLASH($\alpha$,$\lambda$) estimator. It offers some interesting insights. First, we focus on the probability with which inequality (ref) holds. For VAR estimation with a penalized least-squares objective function, such probabilities are governed by tail probabilities of the process $\{\frac{1}{T} \sum_{t=1}^T y_{it} \epsilon_{jt}\}$ (see, e.g. lemma 4 in KockCallot2015, or lemmas 5--6 in Medeiros2016). Because Yule-Walker estimation relies primarily on autocovariance matrix estimation, our probability depends on the tail decay of the distribution of $\{\big\|\widehat{\bm V}_h - \bm V\big\|_\vdash \}$. Overall, the probability of (ref) improves through faster tail decay of the innovation distribution (compare cases (a) and (b)) and banded autocovariance matrix estimation (Theorem (ref)). Second, we look closer at the performance upper bound itself. The right-hand side of (ref) demonstrates that the upper bound of the prediction and estimation error is increasing in $\bar{\omega}_\alpha$, which in turn is increasing in the bandwidths $k_0$ and $l_0$, increasing in the group sizes ($\alpha<1$), and increasing in the number of relevant interactions $\left\lvertS\right\rvert$ ($\alpha>0$). Furthermore, the prediction and estimation error increases in the degree of penalization. Whereas this seemingly suggests to minimize $\lambda$ as to improve performance bounds, we emphasize that the effect of regularization in Theorem (ref) is two-fold: increasing regularization deteriorates the performance bound, but increases the probability of the set on which the performance bound holds. Intuitively, shrinkage induces finite-sample bias which worsens accuracy, but simultaneously reduces sensitivity to noise, thereby enabling performance guarantees at higher degrees of certainty.

The aforementioned effects can also be demonstrated by means of an asymptotic analysis. Based on Theorem (ref), we derive the conditions for convergence of the prediction and estimation errors in the following corollary. The exact convergence rates are also provided.

Corollary (ref) provides 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\lvertS\right\rvert$ and the degree of penalization $\lambda$.\footnote{Recall 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 (ref) 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. Some concrete examples are provided in Remark (ref).

Exogenous variables

We generalize model specification (ref) by accommodating $K$ exogenous variables, i.e.

Each vector $\bm x_{t,k}=(x_{1t,k},\ldots,x_{Nt,k})'$ augments the spatio-temporal vector autoregression with an extra regressor. This regressor may vary over time and it is exogenous, i.e. we have $ \operatorname{\mathbb{E}}(\bm x_{t,k} \bm \epsilon_t')=\mathbf{O}$ for $k=1,\ldots,K$. For notational brevity, we consider the situation in which the exogenous regressors $x_{it,1}\ldots,x_{it,K}$ can only directly influence spatial unit $i$. This explains the diagonal structure in $\operatorname{\text{diag}}(\bm \beta_k)$. In Remark (ref) we argue that this simplification does not greatly hinder generality. In contrast to MaGuoWang2021, we allow $\bm \beta_k=(\beta_{1k},\ldots,\beta_{Nk})'$ to vary with location. We keep $K$ fixed.

To account for the exogenous variables, we modify the generalized Yule-Walker estimator of Section (ref). We recall $\bm \varSigma_j = \operatorname{\mathbb{E}}(\bm y_t \bm y_{t-j}')$, and define the matrices $\bm \varSigma_{j}^{x_k y}=\operatorname{\mathbb{E}}(\bm x_{t,k} \bm y_{t-j}')$ and $\bm \varSigma_j^{x_k x_\ell}=\operatorname{\mathbb{E}}(\bm x_{t,k} \bm x_{t-j,\ell}')$. Two sets of Yule-Walker equations, namely

are derived by post-multiplying the model by respectively $\bm y_{t-1}'$ and $\bm x_{t,k}'$, and taking expectations. Compared to (ref), the Yule-Walker equations in (ref) contain the additional term $\sum_{k=1}^K \operatorname{\text{diag}}(\bm \beta_k) \bm \varSigma_1^{x_k y}$ to provide information on $\bm \beta_1,\ldots,\bm \beta_K$. However, if $\bm \varSigma_1^{x_k y}=\mathbf{O}$ (e.g. when $\{\bm x_{t,k}\}$ and $\{\bm y_t\}$ are independent and $\bm \beta_k=\boldsymbol{0}$), then (ref) alone will not identify $\bm \beta_k$. We therefore add the additional Yule-Walker equations in (ref). To develop the estimator, we combine (ref) and (ref) into

From this point onward, the development of the SPLASHX($\alpha$,$\lambda$) estimator mimics the reasoning of page (ref) closely. First, we focus on the $i$\textsuperscript{th} spatial unit and collect all the nonzero coefficients of $\bm A$ and $\bm B$ (as stipulated by Assumption (ref)) in $\bm c_i$. Letting $\bm V_i^*$ denote the columns in $\bm V^*$ related to $\bm c_i$ and defining both $\bm \sigma_i^* =

' \bm e_i$ and $\bm w_{ik}^*=\bm W_k^* \bm e_i$, result \eqref{eq:EXOsystem} implies $$ \bm V_i^* \bm c_i + \sum_{k=1}^K \bm w_{ik}^* \beta_{ik} = \bm \sigma_i^*. $$ Second, we define (a) the sample counterparts of $\bm \varSigma_j$, $\bm \varSigma_j^{x_k y}$ and $\bm \varSigma_j^{x_k x_l}$ as respectively $\hat \bm \varSigma_j = \frac{1}{T} \sum_{t=j+1}^T \bm y_t \bm y_{t-j}'$, $\hat \bm \varSigma_j^{x_k y }=\frac{1}{T} \sum_{t=j+1}^T \bm x_{t,k} \bm y_{t-j}'$ and $\hat\bm \varSigma_j^{x_k x_\ell}= \frac{1}{T} \sum_{t=j+1}^T \bm x_{t,k} \bm x_{t-j,\ell}'$, and (b) define the quantities $\hat\bm \sigma_i^*$, $\hat\bm w_{ik}^*$ and $\hat\bm V_i^*$ based on their underlying sample covariance matrix estimators. Finally, set $\hat\bm \sigma^*=(\hat\bm \sigma_1^{*\prime},\ldots,\hat\bm \sigma_N^{*\prime})'$, $\hat\bm V^{*(d)}=\operatorname{diag}(\hat\bm V_1^*,\ldots,\hat\bm V_N^*)$, and $\hat\bm W_k^{*(d)}=\operatorname{diag}(\hat\bm w_{1k}^*,\ldots,\hat\bm w_{Nk}^*)$. The SPLASHX($\alpha$,$\lambda$) objective function is

This objective function allows for the estimation of $\bm \beta_1,\ldots,\bm \beta_K$, sparse coefficients, completely sparse vectors $\bm \beta_k$, and completely sparse diagonals in the coefficient matrices $\bm A$ and $\bm B$. There is a clear mathematical resemblance between the SPLASH and SPLASHX estimators. Accordingly, under appropriate modifications to Assumptions (ref)--(ref), a finding similar to Theorem is attainable. We provide this result as Theorem (ref) and refer the reader to Supplement (ref) for detailed assumptions and proofs.

Simulations

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 (ref). 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 Gao2019. 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:

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

Design B (Spatial grid with neighbor interactions): As in Figure (ref), we consider an $(m\times m)$ grid of spatial units. For $m=5$ ($m=10$), this results in a cross-sectional dimension of $N=25$ ($N=100$). 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=10$), the temporal matrix $\bm B$ is a diagonal matrix with elements 0.25 (0.21) 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.904).

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) SPLASH($0.5$,$\lambda$) provides equal weight to sparsity at the group and individual level, and (3) SPLASH(1,$\lambda$) encourages unstructured sparsity only.\footnote{The choice for $\alpha=0.5$ is solely made to illustrate the effect of combining both group and individual penalties. For different designs, this choice may or may not be optimal.} In congruence with Theorems (ref) and (ref), we rely on banded autocovariance matrices $\mathcal{B}_{h_{}}\big(\hat\bm \varSigma_0\big)$ and $\mathcal{B}_{h_{}}\big(\hat\bm \varSigma_1\big)$. The bandwidth choice is determined by the bootstrap procedure described in guowangyao2016. Second, we include two unpenalized estimators in the spirit of Gao2019: GMWY and GMWY($k_0$). The GMWY estimator implements generalized Yule-Walker estimation for banded $\bm A$ and $\bm B$ with the bandwidth being chosen by the selection rule proposed by Gao2019, whereas GMWY$(k_0)$ is based on the true bandwidth $k_0$. 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 '$ and $\hat \bm \varSigma_1 = \frac{1}{T} \sum_{t=2}^T \bm y_t \bm y_{t-1}'$.\footnote{In unreported simulation results (available upon request), we find that the results are insensitive to this choice.} As GMWY$(k_0)$ is infeasible in practice, it is given a comparative advantage. The third set solely contains the $L_1$-penalized reduced form VAR(1) estimator (abbreviated PVAR). In detail, 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 $\mathcal{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\lvertc_{ij}\right\rvert$. This estimator is well-researched in the literature KockCallot2015,Gelper2016,Masini2019, 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 will be assessed using the Relative Mean-Squared Forecast Error (RMSFE). Using a superscript $j$ to index a specific Monte Carlo replication, the RMSFE is calculated as

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

Finally, a word on the selection of the the penalty parameter. For the SPLASH estimator, we calculate the maximum penalty, $\lambda_\max$, as the smallest value producing the zero solution for all values of $\alpha$, i.e.

Given $\lambda_\max$, we define the smallest penalty $\lambda_\min$ as $10^{-4}\lambda_\max$ ($10^{-6}\lambda_\max$) for Design A (B) and construct an ordered grid of 20 equidistant values on a log-scale, say $\lambda_\max = \lambda_1 > \lambda_2 > \ldots > \lambda_{20} = \lambda_\min$. Estimating SPLASH solutions for each $\lambda_i$, a grid of $\alpha$-values, and each individual simulation trial is computationally expensive (especially for large $N$). We instead perform a small-scale preliminary analysis in which we draw a small set of simulations from Designs A and B on which we estimate all solutions for a given value of $T$. Then, we choose the order $i_T \in \lbrace 1,\ldots,20\rbrace$ that minimizes the RMSFE in this preliminary set of simulations. This process of choosing the order $i_T$ on a log-equidistant grid for each value of $T$, is equivalent to setting $\lambda = m_T\lambda_\max$ with $m_T = 10^{-4(i_T-1)/20}$ or $m_T = 10^{-6(i_T-1)/20}$ for designs A and B, respectively. For Design A (B), our selected orders for $T=\{500,1000,2000\}$ are $i_T = \{9,10,11\}$ ($i_T =\{10,11,12\}$), corresponding to $m_T \approx 0.025,0.015,0.01$ ($m_T \approx 0.002,0.001,0.0005$), respectively. Having fixed the preferred order or multiplier, it remains to estimate a single solution per $\alpha$-value, thus resulting in substantial reductions in computation time. The penalized VAR is computationally less expensive. Accordingly, we choose its penalty parameter based on a time series cross-validation (TSCV) scheme Hyndman2018. In our implementation of TSCV, the first 80% of the data is used to fit multiple solutions on, which are then evaluated based on the MSFE obtained on the latter 20% of the data. The preferred penalty is chosen as the solution that attains the smallest MSFE.\footnote{We also tried to select the penalty for the PVAR as the sparsest solution whose prediction error lies within one standard error of the minimum prediction error. This selection rule, however, did not lead to an improvement in forecast or estimation accuracy.}

Simulation results

The results for Design A are reported in Table (ref). First, we consider the predictive performance in Panel 1. For all methods, we observe a monotonic decrease in RMSFE when $T$ increases. The SPLASH estimators and GMWY($k_0$) exhibit the best overall forecast performance, with SPLASH outperforming for smaller sample sizes ($T=500$). Among the SPLASH estimators, SPLASH($0$,$\lambda$) attains the lowest RMSFE in the majority of specifications but differences are generally marginal. The penalized VAR forecasts are less accurate than the aforementioned methods. An explanation is that sparsity patterns in the reduced form representation are less prevalent and thus more difficult to exploit. Direct estimation of the contemporaneous spatial interactions thus delivers forecast improvements over regularized reduced form estimation. The GMWY estimator is highly competitive when $T=2000$ but performs notably worse for small $N$ and $T$. As GMWY has a tendency to select a too large bandwidth (as in Gao2019, table 1), this is probably caused by the estimation of redundant parameters. Given that the majority of sparsity in this design comes from the small bandwidth of $\bm A$ and $\bm B$, which is fully exploited by the infeasible GMWY($k_0$) estimator, we consider it reassuring that the SPLASH estimators attains comparable, and occasionally better, forecast performance without necessitating an a priori specification of the bandwidth.

Next, we explore the estimation accuracy for $\bm A$ and $\bm B$ in Panels 2 and 3, respectively. As before, all estimators display an improvement in estimation accuracy when $T$ increases. The SPLASH($0$,$\lambda$) attains a lower estimation error than the SPLASH($0.5$,$\lambda$) estimator, which in turn performs better than the unstructured sparsity variant SPLASH($1$,$\lambda$). The tight bandwidth in this design implies that many diagonals ought to be set to zero, which seems to be best effectuated by means of the group penalty. The GMWY($k_0$) estimator appears to deliver somewhat more accurate estimates than SPLASH for larger values of $T$. This apparently slower convergence of the SPLASH estimator might, at least partly, be considered the price of not knowing the true sparsity pattern, as represented by the term $\bar{\omega}_\alpha$ in Theorem (ref). It is worth mentioning, however, that the choice of penalty parameter is motivated based on the predictive performance, which may not be optimal from the perspective of estimation accuracy. Indeed, in an unreported analysis we find that the penalty that minimizes the estimation error is typically higher and delivers sparser solutions. Regarding the GMWY estimator, we note that the detrimental effect of overestimating the bandwidth in smaller sample sizes is again visible, with the estimation error being substantially larger for the $T=500$ setting.

Simulation results for Design B are shown in Table (ref). The high RMSFEs for the GMWY estimators are most striking. In the setting $N=25$ and $T=500$, the GMWY estimator frequently selects a bandwidth equal to 1, translating to inferior performance across all metrics. The GMWY($k_0$) estimator, on the other hand, is based on the correct bandwidth. This method, however, forecasts far worse, while its estimation accuracy instead is competitive to SPLASH. Upon closer inspection, we find that the high RMSFE in this case is driven by a few extreme prediction errors. These prediction outliers in turn correspond to simulation trials in which the smallest absolute eigenvalues of the estimated matrix $\bm I - \hat{\bm A}$ are close to zero (see Fig (ref) in the Supplementary Appendix). This implies that the GMWY estimator may be prone to stability issues when the bandwidth is large relative to the dimension.\footnote{Recall that converting the spatial representation to the reduced form representation requires inverting $\bm I-\hat{\bm A}$.} Apparently, owing to the implementation of sparsity, the SPLASH estimator does not suffer from such stability issues. For $N=25$ and $T=2,000$, the bandwidth selection in GMWY improves, while its forecast performance ironically worsens as a result of the increasing stability issues. The remaining results tell the same story as in Design A; SPLASH(0,$\lambda$) and SPLASH(0.5,$\lambda$) are forecasting very close to the optimal forecast, and forecast notably better than the PVAR. While the forecast performance of SPLASH(1,$\lambda$) comes across as equivalent to the SPLASH implementation with group-penalization, the estimation accuracy is superior for the latter. Hence, the group penalty seems especially valuable for the purpose of model interpretation.

A small visual analysis provides further evidence on the favourable estimation accuracy obtained by SPLASH with shrinkage towards diagonally structured sparsity. We 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 (ref) illustrates the similarity between the true matrix $\bm A$ and the average magnitude of the estimated coefficients.

Simulations with exogenous regressors

In this section, we examine the estimation performance of our estimator in the presence of exogenous regressors. Simulated data is drawn from

where $\bm A$ and $\bm B$ are generated analogously to designs A and B in Section (ref), and the coefficients of the exogenous regressors are given by $\bm \beta_1 = \bm \iota_N$ and $\bm \beta_2 = \bm{0}_N$. Hence, only $\bm x_{t,1}$ contributes to the variation in $\bm y_t$. Accordingly, we henceforth refer to $\bm x_{t,1}$ and $\bm x_{t,2}$ as the relevant and irrelevant exogenous regressor, respectively. All elements of exogenous variables and innovations are drawn i.i.d. from $N(0,1)$. At each simulation trial, we implement the same estimators as considered in Section (ref), with the exception of the penalized VAR which is omitted here. The selection of $\lambda$ for the SPLASH estimator is again done via the construction $\lambda_T = m_T\lambda_\max$, where the sequence of multipliers $m_T$ are the same as those in Section (ref). Forecasts under (ref) require predictions of the exogenous variables. This leads to two complications. First, in the reduced-form VAR(1) representation, $\bm y_t = \bm B \bm y_{t-1} + \big[\bm D\operatorname{\text{diag}}(\bm \beta_1) \big]\bm x_{t,1} + \big[\bm D\operatorname{\text{diag}}(\bm \beta_2) \big] \bm x_{t,2} + \bm D\bm \epsilon_t$, the coefficient matrices in front of $\bm x_{t,1}$ and $\bm x_{t,2}$ are no longer diagonal. This would cause an unfair comparison with PVAR so we decided to omit the penalized VAR approach from the comparison. Second, to avoid results that depend on the prediction method employed, we focus solely on the estimation accuracy.

The estimation accuracy for the estimates of $\bm A$ and $\bm B$ is compared on the basis of the metrics $EE_A$ and $EE_B$, as given in (ref). In addition, we also report the average estimation errors of the coefficients for the relevant and irrelevant exogenous regressors, calculated as

The results for Design A and B are reported in Tables (ref) and (ref), respectively.

First, we consider the results for Design A. The first two panels display the average estimation errors in $\bm A$ and $\bm B$, respectively. Reassuringly, all estimators display a clear monotonic decrease in estimation accuracy with growing sample size. Comparing the SPLASH estimators among each other, we observe that shrinkage towards group sparsity is most beneficial in high-dimensional settings ($T=500$ or $N=100$). In these instances, the SPLASH($0$,$\lambda$) and SPLASH($0.5$,$\lambda$) estimators obtain the lowest estimation error across all methods. Conversely, when the dimension is small ($N=25$) and sample size is large ($T=2000$), we find little gain in penalizing towards structured sparsity and the SPLASH($1$,$\lambda$) outperforms all other methods. Furthermore, the SPLASH estimators attain a lower estimation error than the GMWY estimators for almost all settings, with the performance gains attained by SPLASH being most pronounced in the case where the sample size is small, i.e. when the exploitation of sparsity matters most. Comparing the GMWY estimators, we find that, in lower-dimensional settings, using a data-driven selection of the bandwidth performs comparable to relying on the true bandwidth. However, when $N=100$ and $T=500$, we find that the bandwidth selection procedure over-estimates the true bandwidth in roughly 40% of the simulation trials. Accordingly, the GMWY estimator attains inferior estimation accuracy in this particular setting. Regarding the exogenous regressors, we observe a similar monotonic decrease in the estimation error for the coefficients of both the relevant and irrelevant exogenous regressor. The SPLASH estimator outperforms the GMWY estimators across all dimensions and sample sizes, with the performance gain again being most prominent in the higher-dimensional settings. The estimation error obtained by SPLASH for the irrelevant exogenous regressor is remarkably small, further demonstrating the benefits of the incorporated shrinkage.

The results for Design B depict a similar, if not more compelling, story. The SPLASH estimators again outperform across all settings, with the performance differentials between SPLASH and GMWY being more pronounced compared to Design A. Again, we observe that the SPLASH(1,$\lambda$) estimator seems to outperform based on $EE_A$ and $EE_B$ for $N=25$, whereas the SPLASH(0,$\lambda$) and SPLASH($0.5$,$\lambda$) estimators do better when $N=100$. We conjecture that shrinkage towards structured sparsity only becomes beneficial when the group sizes are sizable enough, at which point the accumulation of selection errors by SPLASH(1,$\lambda$) starts to deteriorate the overall estimation accuracy. Contrasting the performance of SPLASH to the GMWY estimators, we observe that the exploitation of sparsity within the bandwidth results in substantial performance gains across all specifications and coefficient matrices. Moreover, the bandwidth selection procedure of the GMWY estimator now frequently selects very small bandwidths. This negatively impacts the estimation accuracy when $N=25$, whilst having a positive impact when $N=100$. In the latter case, the number of parameters to estimate is simply too large without further regularization, such that one might be better off by forcing most diagonals to zero, even if some of those are relevant. Interestingly, the inability to exploit sparsity also affects the estimation accuracy for the relevant exogenous regressors, as the third panel reveals a sizeable difference in the $EE_R$ between the SPLASH and GMWY estimators. Regarding the irrelevant exogenous regressor, we find that the $EE_{IR}$ is substantially larger for GMWY, but comparable across the SPLASH and GMWY($k_0$) estimators.

Overall, SPLASH unambiguously attains more accurate estimates of all coefficient matrices in the spatial VAR with exogenous regressors. In line with expectations, the performance gain of SPLASH is most notable in high-dimensional designs with substantial degrees of sparsity. However, even in lower-dimensional designs in which the degree of sparsity is less, SPLASH remains competitive to the GMWY estimators.

Empirical Application

Nitrogen dioxide (NO\textsubscript{2}) is emitted during combustion of fossil fuels (e.g. by motor vehicles) and it has been associated with adverse effects on the respiratory system.\footnote{The direct health effect of nitrogen dioxide is difficult to determine because its emission process is typically accompanied with the emission of other air pollutants (see, e.g. brunekreefholgate2002).} The Air Quality Standards Regulations 2010 requires a regular monitoring of NO\textsubscript{2} concentration levels in the UK.\footnote{Source: https://www.legislation.gov.uk/uksi/2010/1001/contents/made.} Using satellite data, we examine the empirical performance of the SPLASH estimator when predicting daily NO\textsubscript{2} concentrations in Greater London. This satellite data is publicly available via the Copernicus Open Access Hub and we consider the time span from 1 August 2018 to 18 October 2020.\footnote{See \url{http://www.tropomi.eu/} for more info on TROPOMI data products and use \url{https://scihub.copernicus.eu/} to access the database.} The original $\text{NO}_2$ concentrations are reported in mol/m$^2$, which we convert to mol/cm$^2$ to avoid numerical instabilities caused by small-scale numbers. The far majority of measurements are captured 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 (see Figure (ref), part c). All available $\text{NO}_2$ measurements are averaged within each cell and within the same day. The resulting data set contains 0.8% missing observations, which we impute using the Multivariate Time Series Data Imputation (mtsdi) R package.\footnote{This imputation method is proposed by JungerDeLeon2015 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 (641 days) allowing 160 one-step ahead forecasts to be made. 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 (ref), see page (ref)), 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. For SPLASH, we follow the procedure described in Section (ref) and set $\lambda = 1.8\times 10^{-4}\lambda_\max$. 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\lverti-j\right\rvert > \lfloor N/4 \rfloor = 11$. By ordering the spatial units vertically, this banding puts no restrictions on the vertical interactions but allows no more than second-order interaction between horizontal neighbours (see Figure (ref) in the Supplementary Appendix (ref) for details).

The forecast performance is measured along three metrics and is always expressed 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 a 5% significance level (\#sign. wins), 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 additionally report the MAFE because the $\text{NO}_2$ column densities display several abrupt spikes which may carry too much weight when relying on a squared loss function. The results are reported in Table (ref).

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 a large bandwidth is needed to allow for second-order horizontal interaction and, consequently, a large number of parameters to estimate. With an RMSFE of 0.919, SPLASH(0,$\lambda$) does manage to improve upon the benchmark. In fact, the MSFEs for all 45 spatial units are smaller than that of the benchmark, 42 of which are found to be significant by a Diebold-Mariano test based on the squared forecast errors. Allowing for sparsity within groups does not seem to deliver additional forecast improvements, as SPLASH(0.5,$\lambda$) attains a slightly worse forecast performance and significantly outperforms the benchmark for only 39 locations. Completely omitting regularization at the group level results in a further deterioration of the forecast performance, with SPLASH(1,$\lambda$) attaining an RMSFE of 0.94 and significantly beating the benchmark for 30 out 45 spatial units. We take this as evidence that the ability to promote diagonally structured sparsity is indeed beneficial in real-life spatial applications, although even estimating the spatial VAR with unstructured sparsity manages to improve upon regularized reduced form VAR estimation.

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 still has the worst forecast accuracy for GMWY and SPLASH(0,$\lambda$) continues to perform 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 NO$_2$ 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. First, in Figure (ref)(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 Figure (ref)(b), each cell indicates the proportion of rolling windows the corresponding spatial interaction is estimated as being non-zero. 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 with the spatial grid overlaid (Figure (ref)(c)). We explicitly visualize the interactions implied by (ref)(b) for two pixels -- pixel 1 (left-top) and pixel 23 (center) -- using arrows whose thickness is determined by the average absolute magnitudes estimated in Figure (ref)(a). The emerging pattern of spatial interactions shows clearly the interactions between NO$_2$ 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.

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 structured forms of sparsity without imposing such structure a priori. We derive consistency of our estimator in a joint asymptotic framework in which the number of both spatial units and 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. The SPLASHX extension explains how to include exogenous variables. As an application, we use SPLASH to predict satellite-measured NO\textsubscript{2} concentrations in London. We find evidence for spatial interactions between 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 Gao2019 (the inspiration for the development of SPLASH).

Acknowledgements

This paper (or earlier versions hereof) has been presented during the internal seminar of the Quantitative Economics department of Maastricht University, the Econometrics Internal Seminar (EIS) at Erasmus University Rotterdam, the Bernoulli-IMS One World Symposium, the workshop on Dimensionality Reduction and Inference in High-dimensional Time Series at Maastricht University, the 2021 Annual Conference of the International Association for Applied Econometrics (IAAE), the 5\textsuperscript{th} Conference on Econometric Models of Climate Change, the internal seminar at Tor Vergata University of Rome, and the 2021 $(\text{EC})^2$ Conference. We gratefully acknowledge comments and feedback from the participants. Suggestions by Stephan Smeekes and Ines Wilms were particularly helpful so we thank them explicitly. All remaining errors are our own.