Vector Autoregressive Models with Spatially Structured Coefficients for Time Series on a Spatial Grid

Motivated by the need to analyze readily available data collected in space and time, especially in environmental sciences, we propose a parsimonious spatiotemporal model for time series data on a spatial grid. In essence, our model is a vector autoregressive model that utilizes the spatial structure to achieve parsimony of autoregressive matrices at two levels. The first level ensures the sparsity of the autoregressive matrices using a lagged-neighborhood scheme. The second level performs a spatial clustering of the nonzero autoregressive coefficients such that within some subregions, nearby locations share the same autoregressive coefficients while across different subregions the coefficients may have distinct values. The model parameters are estimated using the penalized maximum likelihood with an adaptive fused Lasso penalty. The estimation procedure can be tailored to accommodate the need and prior knowledge of a modeler. Performance of the proposed estimation algorithm is examined in a simulation study. Our method gives reliable estimation results that are interpretable and especially useful to identify geographical subregions, within each of which, the time series have similar dynamical behavior with homogeneous autoregressive coefficients. We apply our model to a wind speed time series dataset generated from a climate model over Saudi Arabia to illustrate its power in explaining the dynamics by the spatially structured coefficients. Moreover, the estimated model can be useful for building stochastic weather generators as an approximation of the computationally expensive climate model.

Authors and Affiliations

1. Department of Mathematics & Statistics, Dalhousie University, Halifax, NS, B3H 4R2, Canada

   Yuan Yan

2. Institute of Statistical Science Academia Sinica, Taipei 115, Taiwan

   Hsin-Cheng Huang

3. Statistics Program, King Abdullah University of Science and Technology (KAUST), Thuwal, 23955-6900, Saudi Arabia

   Marc G. Genton

Corresponding author

Correspondence to Yuan Yan.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This publication is based upon work supported by the King Abdullah University of Science and Technology (KAUST) Office of Sponsored Research (OSR) under Award No: OSR-2018-CRG7-3742 and Academia Sinica Investigator Award AS-IA-109-M05.

A Appendix

Derivation of (5):

For the two terms involved in the above equation,

Then, it is easy to see that the log-likelihood is a quadratic form of \(\varvec{\alpha }\), i.e.,

where \(\varvec{X}\) and \(\varvec{y}\) satisfy (3) and (4). \(\square \)

Asymptotic Properties:

Some notations are needed to state the asymptotic properties. We consider an undirected graph \({\mathcal {G}}=({\mathcal {V}},{\mathcal {E}})\), where \({\mathcal {V}}=\{1,\ldots ,m\}\) is a set of vertices corresponding to the \(m=K n_I + n_B\) parameters in \(\varvec{\alpha }\) to be estimated, and \({\mathcal {E}}\subset {\mathcal {V}}\times {\mathcal {V}}\) is a set of undirected edges that depends on the structure of the spatial grid. Here \({\mathcal {V}}\) can be partitioned into \(n_B+K\) disjoint components, including \(n_B\) isolated nodes (corresponding to the non-regularized coefficients in \(\varvec{\alpha }_\text {B}\)) and K identical disjoint components each with \(n_I\) vertices and the same inner grid structure (corresponding to the penalized coefficients in \(\varvec{\alpha }_k\), for \(k=1,\ldots ,K\)). With these notations, the penalty term in (6) can be rewritten as

where \(w_{{i},{j}}=|{\tilde{\alpha }}_i-{\tilde{\alpha }}_j|^{-\gamma }\).

We let \(\varvec{\alpha }^0\) be the true parameter vector and define \(\varvec{\xi }=\varvec{y}-\varvec{X}\varvec{\alpha }^0\). We introduce \({\mathcal {A}}=\left\{ (i,j) \in {\mathcal {E}}: \alpha ^0_i=\alpha ^0_j, i,j=1,\ldots ,n_I \right\} \) and consider the subgraph \({\mathcal {G}}_{{\mathcal {A}}}=({\mathcal {V}},{\mathcal {A}})\) of \({\mathcal {G}}\). We denote by \(m_0\) the number of its connected components, that is, the number of distinct values in \(\varvec{\alpha }^0\) supported by \({\mathcal {G}}\). We further denote by \({\mathcal {V}}_1,\ldots , {\mathcal {V}}_{m_0}\) the sets of nodes of each connected components of \({\mathcal {G}}_{{\mathcal {A}}}\) and set \(l_i=\min {{\mathcal {V}}_i}\) for \(i=1,\ldots ,m_0\). We define \(\varvec{\alpha }_{{\mathcal {A}}}^0=(\alpha ^0_{l_1},\ldots ,\alpha ^0_{l_{m_0}})'\) and \(\varvec{X}_{\mathcal {A}}\) a matrix whose i-th column, for \(i=1,\ldots ,m_0\), is \(\varvec{X}_{{\mathcal {A}}_i}=\sum _{j\in {\mathcal {V}}_i}\varvec{X}_j\), where \(\varvec{X}_j\) is the j-th column of \(\varvec{X}\). We assume that \(\varvec{X}'\varvec{X}/(T-1)\rightarrow \mathbf {C}\) as \(T\rightarrow \infty \) for some positive-definite \(m\times m\) matrix \(\mathbf {C}\), which depends on \(\varvec{\alpha }^0\) and \(\varvec{\Psi }\). We denote by \(\mathbf {C}_{{\mathcal {A}}}\) the limiting \(m_0\times m_0\) matrix of \(\varvec{X}_{{\mathcal {A}}}'\varvec{X}_{{\mathcal {A}}}/(T-1)\) as \(T\rightarrow \infty \).

We derive the asymptotic properties of the adaptive Lasso estimator \(\hat{\varvec{\alpha }}\) from solving (7). We let \({\mathcal {A}}_n=\left\{ (i,j) \in {\mathcal {E}}: {\hat{\alpha }}_i={\hat{\alpha }}_j,i,j=1,\ldots ,n_I \right\} \) and \(\hat{\varvec{\alpha }}_{{\mathcal {A}}}=({\hat{\alpha }}_{l_1},\ldots ,{\hat{\alpha }}_{l_{m_0}})'\).

Theorem 1

Suppose that \(\lambda /\sqrt{T}\rightarrow 0,\, \lambda T^{(\gamma -1)/2}\rightarrow \infty \), and \(\mathbf {C}=\lim _{T\rightarrow \infty }\) \(\varvec{X}'\varvec{X}/(T-1)\) is positive-definite. Then as \(T\rightarrow \infty \), the following are satisfied by the adaptive Lasso estimator \(\hat{\varvec{\alpha }}\):

1. 1.

     Consistency in grouping: \(\Pr ({\mathcal {A}}_n={\mathcal {A}})\rightarrow 1\).

2. 2.

     Asymptotic normality: \(\sqrt{T}\big (\hat{\varvec{\alpha }}_{{\mathcal {A}}}-\varvec{\alpha }^0_{{\mathcal {A}}}\big )\xrightarrow {d}{\mathcal {N}}_{m_0}(\varvec{0},\mathbf {C}_{{\mathcal {A}}}^{-1})\).

Proof

The following proof is obtained by adapting the Proof of Theorem 2 in Zou (2006) and Theorem 3 in Viallon et al. (2013).

First, it is known that \(\hat{\varvec{\Psi }}\) in (8) is a consistent estimator of \(\varvec{\Psi }\), as T goes to infinity; see Chapter 5 of Lütkepohl (2005). We define \(V_T(\mathbf {u})=F(\varvec{\alpha }^0)-F(\varvec{\alpha }^0+\mathbf {u}/\sqrt{T})\), with F defined in (6). It is obvious that \(V_T(\mathbf {u})\) is minimized at \(\sqrt{T}(\hat{\varvec{\alpha }}-\varvec{\alpha }^0)\) and

We have

We denote \(\varvec{\mathbf {u}}_{{\mathcal {A}}}=(\mathbf {u}_{l_1},\ldots ,\mathbf {u}_{l_{m_0}})'\) and then, with an application of the Martingale difference central limit theory to \(\varvec{\xi }'\varvec{X}\), we obtain

for \(\mathbf {u}\in {\mathbb {R}}^m\), where \({\mathbf{W }}_{\mathcal {A}}\sim {\mathcal {N}}_m(\varvec{0},{\mathbf{C }}_{{\mathcal {A}}})\); V is convex and has a unique minimum satisfying \(u_i=u_j\) for all \((i,j)\in {\mathcal {A}}\) and \(\mathbf {u}_{\mathcal {A}}=\mathbf {C}_{{\mathcal {A}}}^{-1}\mathbf {W}_{\mathcal {A}}\). The asymptotic normality part can be derived as in Zou (2006) by using the epi-convergence results.

For consistency in grouping, we need to show that for all \((i,j)\not \in {\mathcal {A}}, \Pr ((i,j)\in {\mathcal {A}}_n^c)\rightarrow 1\) and for all \((i,j) \in {\mathcal {A}}, \Pr ((i,j)\in {\mathcal {A}}_n^c)\rightarrow 0\). The first part is implied by the asymptotic normality result. To prove the second part, we apply the subgradient equations for the optimality condition, for \(i=1,\ldots ,m\):

where \(t_{ij}={{\,\mathrm{sign}\,}}({\hat{\alpha }}_i-{\hat{\alpha }}_j)\) if \({\hat{\alpha }}_i\ne {\hat{\alpha }}_j\) and \(t_{ij}\) is some real number in \([-1,1]\) if \({\hat{\alpha }}_i= {\hat{\alpha }}_j.\) We prove by contradiction. Suppose that for \({\mathcal {V}}_k\) that contains at least two vertices, there exist \(i,j \in {\mathcal {V}}_k\) such that \({\hat{\alpha }}_i\ne {\hat{\alpha }}_j\). We define \(\displaystyle a^{\min }=\min _{i \in {\mathcal {V}}_k} {\hat{\alpha }}_i\) and \({\mathcal {V}}^{\min }=\left\{ i: i \in {\mathcal {V}}_k \text { and } {\hat{\alpha }}_i=a^{\min }\right\} \). Summing up the optimality conditions over the indices in \({\mathcal {V}}^{\min }\), we get:

where in the right-hand side, the first sum converges to 0 in probability, while the second sum tends to \(-\infty \). However, the left-hand side is \(O_{p}(1)\), since it can be decomposed as the sum of two asymptotically normal variables as in Zou (2006) with an application of the martingale central limit theorem. Therefore, \(\Pr ((i,j)\in {\mathcal {A}}_n^c)\rightarrow 0\). \(\square \)

Reprints and permissions