Matrix Autoregressive Model with Vector Time Series Covariates for Spatio-Temporal Data

We adopt the following notations throughout the article. We use calligraphic bold-face letters (e.g. $\mathbfcal{X},\mathbfcal{G}$) for tensors with at least three modes, uppercase bold-face letters (e.g. $\mathbf{X},\mathbf{G}$) for matrix and lowercase bold-face letters (e.g. $\mathbf{x},\mathbf{z}$) for vector and blackboard bold-faced letters for sets (e.g. $\mathbb{R},\mathbb{S}$). To subscript any tensor/matrix/vector, we use square brackets with subscripts such as $[\mathbfcal{G}]_{ijd},[\mathbf{z}_{t}]_{d},[\mathbf{X}_{t}]_{ij}$, and we keep the subscript $t$ inside the square bracket to index time. Any fibers and slices of tensor are subscript-ed with colons such as $[\mathbfcal{G}]_{ij:}$, $[\mathbfcal{G}]_{::d}$ and thus any row and column of a matrix is denoted as $[\mathbf{X}_{t}]_{i:}$ and $[\mathbf{X}_{t}]_{:j}$. If the slices of tensor/matrix are based on the last mode such as $[\mathbfcal{G}]_{::d}$ and $[\mathbf{X}_{t}]_{:j}$, we will often omit the colons and write as $[\mathbfcal{G}]_{d}$ and $[\mathbf{X}_{t}]_{j}$ for brevity. For any tensor $\mathbfcal{X}$, we use $\mathrm{\mathbf{vec}}\left(\mathbfcal{X}\right)$ to denote the vectorized tensor. For any two tensors $\mathbfcal{X},\mathbfcal{Y}$ with identical size, we define their inner product as: $\langle\mathbfcal{X},\mathbfcal{Y}\rangle=\mathrm{\mathbf{vec}}\left(\mathbfcal{X}\right)^{\top}\mathrm{\mathbf{vec}}\left(\mathbfcal{Y}\right)$, and we use $\|\mathbfcal{X}\|_{\mathrm{F}}$ to denote the Frobenius norm of a tensor and one has $\|\mathbfcal{X}\|_{\mathrm{F}}=\sqrt{\langle\mathbfcal{X},\mathbfcal{X}\rangle}$.

Following Li and Zhang (2017), the tensor-vector product between a tensor $\mathbfcal{G}$ of size $d_{1}\times\cdots\times d_{K+1}$ and a vector $\mathbf{z}\in\mathbb{R}^{d_{K+1}}$, denoted as $\mathbfcal{G}\bar{\times}_{(K+1)}\mathbf{z}$, or simply $\mathbfcal{G}\bar{\times}\mathbf{z}$, is a tensor of size $d_{1}\times\cdots\times d_{K}$ with $[\mathbfcal{G}\bar{\times}\mathbf{z}]_{i_{1}\ldots i_{K}}=\sum_{i_{K+1}}[\mathbfcal{G}]_{i_{1}\ldots i_{K}i_{K+1}}\cdot[\mathbf{z}]_{i_{K+1}}$. For tensor $\mathbfcal{X}\in\mathbb{R}^{d_{1}\times\cdots\times d_{K}}$, we use $\mathbf{X}_{(k)}\in\mathbb{R}^{d_{k}\times\prod_{m\neq k}d_{m}}$ to denote its $k$-mode matricization. The Kronecker product between matrices is denoted via $\mathbf{A}\otimes\mathbf{B}$ and the trace of a square matrix $\mathbf{A}$ is denoted as $\text{tr}\left(\mathbf{A}\right)$. We use $\bar{\rho}(\cdot),\underaccent{\bar}{\rho}(\cdot),\rho_{i}(\cdot)$ to denote the maximum, minimum and $i^{\rm{th}}$ largest eigenvalue of a matrix. We use $diag(\mathbf{C}_{1},\ldots,\mathbf{C}_{d})$ to denote a block diagonal matrix with $\mathbf{C}_{1},\ldots,\mathbf{C}_{d}$ along the diagonal. More details on the related tensor algebra can be found in Kolda and Bader (2009).

For the matrix time series $\mathbf{X}_{t}\in\mathbb{R}^{M\times N}$ in our modeling context, we assume that all $S=MN$ grid locations are points on an $M\times N$ grid within the domain $\bar{\mathbb{S}}=[0,1]^{2}$. The collection of all the spatial locations is denoted as $\mathbb{S}$ and any particular element of $\mathbb{S}$ corresponding to the $(i,j)^{\rm{th}}$ entry of the matrix is denoted as $s_{ij}$. We will often index the $(i,j)^{\rm{th}}$ entry of the matrix $\mathbf{X}_{t}$ with a single index $u=i+(j-1)M$ and thus $s_{ij}$ will be denoted as $s_{u}$. We use $[N]$ to denote index set, i.e. $[N]=\{1,2,\ldots,N\}$. Finally, we use $k(\cdot,\cdot):\bar{\mathbb{S}}\times\bar{\mathbb{S}}\mapsto\mathbb{R}$ to denote a spatial kernel function and $\mathbb{H}_{k}$ to denote the corresponding Reproducing Kernel Hilbert Space (RKHS).

Let $\{\mathbf{X}_{t},\mathbf{z}_{t}\}_{t=1}^{T}$ be a joint observation of the matrix and the auxiliary vector time series with $\mathbf{X}_{t}\in\mathbb{R}^{M\times N},\mathbf{z}_{t}\in\mathbb{R}^{D}$. To forecast $\mathbf{X}_{t}$, we propose our Matrix AutoRegression with Auxiliary Covariates, or MARAC, as:

$$\mathbf{X}_{t}=\sum_{p=1}^{P}\mathbf{A}_{p}\mathbf{X}_{t-p}\mathbf{B}_{p}^{\top}+\sum_{q=1}^{Q}\mathbfcal{G}_{q}\bar{\times}\mathbf{z}_{t-q}+\mathbf{E}_{t},$$

(1)

where $\mathbf{A}_{p}\in\mathbb{R}^{M\times M},\mathbf{B}_{p}\in\mathbb{R}^{N\times N}$ are the autoregressive coefficients for the lag-$p$ matrix predictor and $\mathbfcal{G}_{q}\in\mathbb{R}^{M\times N\times D}$ is the tensor coefficient for the lag-$q$ vector predictor, and $\mathbf{E}_{t}$ is a noise term whose distribution will be specified later. The lag parameters $P,Q$ are hyper-parameters of the model and we often refer to the model (1) as MARAC$(P,Q)$.

Based on model (1), for the $(i,j)^{\rm{th}}$ element of $\mathbf{X}_{t}$, the MARAC$(P,Q)$ specifies the following model:

$$[\mathbf{X}_{t}]_{ij}=\sum_{p=1}^{P}\left\langle[\mathbf{A}_{p}]_{i:}^{\top}[\mathbf{B}_{p}]_{j:},\mathbf{X}_{t-p}\right\rangle+\sum_{q=1}^{Q}[\mathbfcal{G}_{q}]_{ij:}^{\top}\mathbf{z}_{t-q}+[\mathbf{E}_{t}]_{ij},$$

(2)

where each autoregressive term is associated with a rank-$1$ coefficient matrix determined by the specific rows from $\mathbf{A}_{p},\mathbf{B}_{p}$ and each non-spatial auxiliary covariate is associated with a coefficient vector that is location-specific, i.e. $[\mathbfcal{G}_{q}]_{ij:}$. It now becomes more evident from (2) that the auxiliary vector covariates enter the model via an elementwise linear model. The autoregressive term utilizes $\mathbf{A}_{p},\mathbf{B}_{p}$ to transform each lag-$p$ predictor in a bi-linear form. Using such bi-linear transformation greatly reduces the total amount of parameters in that each lagged predictor that required $M^{2}N^{2}$ parameters previously now only requires $M^{2}+N^{2}$ parameters.

For the tensor coefficient $\mathbfcal{G}_{q}$, we assume that it is spatially smooth. More specifically, we assume that $[\mathbfcal{G}_{q}]_{ijd}$ and $[\mathbfcal{G}_{q}]_{uvd}$ are similar if $s_{ij},s_{uv}$ are spatially close. Formally, we assume that each $[\mathbfcal{G}_{q}]_{d}$, i.e. the coefficient matrix for the $d^{\rm{th}}$ covariate at lag-$q$, is a discrete evaluation of a function $g_{q,d}(\cdot):[0,1]^{2}\mapsto\mathbb{R}$ on $\mathbb{S}$. Furthermore, each $g_{q,d}(\cdot)$ comes from an RKHS $\mathbb{H}_{k}$ endowed with the spatial kernel function $k(\cdot,\cdot)$. The spatial kernel function specifies the spatial smoothness of the functional parameters $g_{q,d}(\cdot)$ and thus the tensor coefficient $\mathbfcal{G}_{q}$.

An alternative formulation for $\mathbfcal{G}_{q}$ would be a low-rank form (Li and Zhang, 2017). We choose locally smooth over low rank to explicitly model the spatial smoothness of the coefficients and avoid tuning the tensor rank of $\mathbfcal{G}_{q}$. We leave the low-rank model for future research and focus on the RKHS framework for the current paper.

Finally, for the additive noise term $\mathbf{E}_{t}$, we assume that it is i.i.d. from a multivariate normal distribution with a separable Kronecker-product covariance:

$$\mathrm{\mathbf{vec}}\left(\mathbf{E}_{t}\right)\mathrel{\overset{\mathrm{i.i.d.}}{\scalebox{1.5}[1.0]{$\sim$}}}\mathcal{N}\left(\mathbf{0},\boldsymbol{\Sigma}_{c}\otimes\boldsymbol{\Sigma}_{r}\right),\quad t\in[T]$$

(3)

where $\boldsymbol{\Sigma}_{r}\in\mathbb{R}^{M\times M},\boldsymbol{\Sigma}_{c}\in\mathbb{R}^{N\times N}$ are the row/column covariance components. Such a Kronecker-product covariance is commonly seen in the covariance models for multi-way data (Hoff, 2011; Fosdick and Hoff, 2014) with the merit of reducing the number of parameters significantly.

Compared to existing models that can only deal with either matrix or vector predictors, our model (1) can incorporate predictors with non-uniform modes. If one redefines $\mathbf{E}_{t}$ in our model as $\sum_{q=1}^{Q}\mathbfcal{G}_{q}\bar{\times}\mathbf{z}_{t-q}+\mathbf{E}_{t}$, i.e. all terms except the autoregressive term, then our model ends up specifying:

	$\displaystyle\mathrm{Cov}(\mathrm{\mathbf{vec}}\left(\mathbf{E}_{t}\right),\mathrm{\mathbf{vec}}\left(\mathbf{E}_{t^{\prime}}\right))$	$\displaystyle=\mathbbm{1}_{\{t=t^{\prime}\}}\cdot\boldsymbol{\Sigma}_{c}\otimes\boldsymbol{\Sigma}_{r}+\mathbf{F}\mathbf{M}\mathbf{F}^{\top}$
	$\displaystyle\mathbf{F}=[\left(\mathbfcal{G}_{1}\right)_{(3)}^{\top}:\cdots:\left(\mathbfcal{G}_{Q}\right)_{(3)}^{\top}],$	$\displaystyle\quad\mathbf{M}=\left[\mathrm{Cov}(\mathbf{z}_{t-q_{1}},\mathbf{z}_{t^{\prime}-q_{2}})\right]_{q_{1},q_{2}\in[Q]}$

where $(\mathbfcal{G}_{q})_{(3)}$ is the mode-3 matricization of $\mathbfcal{G}_{q}$ and we will use $\mathbf{G}_{q}$ to denote it for the rest of the paper. This new formulation reveals how our model differs from other autoregression models with matrix predictors. The covariance of $\mathbf{E}_{t},\mathbf{E}_{t^{\prime}}$ in our model contains a separable covariance matrix $\boldsymbol{\Sigma}_{c}\otimes\boldsymbol{\Sigma}_{r}$ that is based on the matrix grid geometry, a locally smooth coefficient matrix $\mathbf{F}$ that captures the local spatial dependency and an auto-covariance matrix $\mathbf{M}$ that captures the temporal dependency. Consequently, entries of $\mathbf{E}_{t}$ are more correlated if either they are spatially/temporally close or they share the same row/column index and are thus more flexible for spatial data distributed on a matrix grid.

As a comparison, in the kriging framework (Cressie, 1986), the covariance of $\mathbf{E}_{t},\mathbf{E}_{t^{\prime}}$ is characterized by a spatio-temporal kernel that captures the dependencies among spatial and temporal neighbors. Such kernel method can account for the local dependency but not the spatial dependency based on the matrix grid geometry. In the matrix autoregression model (Chen et al., 2021), the authors do not consider the local spatial dependencies among entries of $\mathbf{E}_{t}$ nor the temporal dependency across different $t$. In Hsu et al. (2021), the matrix autoregression model is generalized to adapt to spatial data via fixed-rank co-kriging (FRC) (Cressie and Johannesson, 2008) with $\mathrm{Cov}(\mathrm{\mathbf{vec}}\left(\mathbf{E}_{t}\right),\mathrm{\mathbf{vec}}\left(\mathbf{E}_{t^{\prime}}\right))=\mathbbm{1}_{\{t=t^{\prime}\}}\cdot\boldsymbol{\Sigma}_{c}\otimes\boldsymbol{\Sigma}_{r}+\mathbf{F}\mathbf{M}\mathbf{F}^{\top}$, where $\mathbf{M}$ is a $k\times k$ coefficient matrix and $\mathbf{F}$ is a pre-specified $MN\times k$ spatial basis matrix. Such a co-kriging framework does not account for the temporal dependency of the noises nor does it consider the auxiliary covariates. Our model generalizes these previous works to allow for temporally dependent noise with both local and grid spatial dependency.

The combination of (1) and (3) specifies the complete MARAC$(P,Q)$ model. Vectorizing both sides of (1) yields the vectorized MARAC$(P,Q)$ model:

$$\mathbf{x}_{t}=\sum_{p=1}^{P}\left(\mathbf{B}_{p}\otimes\mathbf{A}_{p}\right)\mathbf{x}_{t-p}+\sum_{q=1}^{Q}\mathbf{G}_{q}^{\top}\mathbf{z}_{t-q}+\mathbf{e}_{t},\quad\mathbf{e}_{t}\mathrel{\overset{\mathrm{i.i.d.}}{\scalebox{1.5}[1.0]{$\sim$}}}\mathcal{N}\left(\mathbf{0},\boldsymbol{\Sigma}_{c}\otimes\boldsymbol{\Sigma}_{r}\right)$$

(4)

where $\mathbf{x}_{t}=\mathrm{\mathbf{vec}}\left(\mathbf{X}_{t}\right),\mathbf{e}_{t}=\mathrm{\mathbf{vec}}\left(\mathbf{E}_{t}\right)$, and recall that $\mathbf{G}_{q}=(\mathbfcal{G}_{q})_{(3)}$. It is now more evident that the Kronecker-product structure of the autoregressive coefficient matrix and the noise covariance matrix greatly reduce the number of parameters, making the regression estimation feasible given finite samples. The spatially smooth structure of $\mathbf{G}_{q}$ leverages the spatial information of the spatial data. In the next section, we will discuss the estimating algorithm of the model parameters of MARAC.

To estimate the parameters of the MARAC$(P,Q)$ model, which we denote collectively as $\boldsymbol{\Theta}$, we propose a penalized maximum likelihood estimation (MLE) approach. Following the distribution assumption on $\mathbf{E}_{t}$ in (3), we can write the negative log-likelihood of $\{\mathbf{X}_{t}\}_{t=1}^{T}$ with a squared RKHS functional norm penalty, after dropping the constants, as:

$$\mathfrak{L}_{\lambda}(\boldsymbol{\Theta})=-\frac{1}{T}\sum_{t\in[T]}\ell\left(\mathbf{X}_{t}|\mathbf{X}_{t-1:P},\mathbf{z}_{t-1:Q};\boldsymbol{\Theta}\right)+\frac{\lambda}{2}\sum_{q\in[Q]}\sum_{d\in[D]}\|g_{q,d}\|_{\mathbb{H}_{k}}^{2},$$

(5)

where $\ell(\cdot)$ is the conditional log-likelihood of $\mathbf{X}_{t}$:

$$\ell\left(\mathbf{X}_{t}|\mathbf{X}_{t-1:P},\mathbf{z}_{t-1:Q};\boldsymbol{\Theta}\right)=-\frac{1}{2}\log|\boldsymbol{\Sigma}_{c}\otimes\boldsymbol{\Sigma}_{r}|-\frac{1}{2}\mathbf{r}_{t}^{\top}\left(\boldsymbol{\Sigma}_{c}^{-1}\otimes\boldsymbol{\Sigma}_{r}^{-1}\right)\mathbf{r}_{t},$$

(6)

and $\mathbf{r}_{t}=\mathbf{x}_{t}-\sum_{p}(\mathbf{B}_{p}\otimes\mathbf{A}_{p})\mathbf{x}_{t-p}-\sum_{q}\mathbf{G}_{q}^{\top}\mathbf{z}_{t-q}$ is the vectorized residual at $t$. To estimate the parameters, one needs to solve a constrained minimization problem:

$$\min_{\boldsymbol{\Theta}}\mathfrak{L}_{\lambda}(\boldsymbol{\Theta}),\quad\text{s.t. }g_{q,d}(s_{ij})=[\mathbfcal{G}_{q}]_{ijd},\text{ for all }s_{ij}\in\mathbb{S}.$$

(7)

We now define the functional norm penalty in (5) explicitly and derive a finite-dimensional equivalent of the optimization problem above. We assume that the spatial kernel function $k(\cdot,\cdot)$ is continuous and square-integrable, thus it has an eigen-decomposition following the Mercer’s Theorem (Williams and Rasmussen, 2006):

$$k(s_{ij},s_{uv})=\sum_{r=1}^{\infty}\lambda_{r}\psi_{r}(s_{ij})\psi_{r}(s_{uv}),\quad s_{ij},s_{uv}\in[0,1]^{2},$$

(8)

where $\lambda_{1}\geq\lambda_{2}\geq\ldots$ is a sequence of non-negative eigenvalues and $\psi_{1},\psi_{2},\ldots$ is a set of orthonormal eigen-functions on $[0,1]^{2}$. The functional norm of function $g$ from the RKHS $\mathbb{H}_{k}$ endowed with kernel $k(\cdot,\cdot)$ is defined as:

$$\|g\|_{\mathbb{H}_{k}}=\sqrt{\sum_{r=1}^{\infty}\frac{\beta_{r}^{2}}{\lambda_{r}}},\qquad\text{where }g(\cdot)=\sum_{r=1}^{\infty}\beta_{r}\psi_{r}(\cdot),$$

(9)

following van Zanten and van der Vaart (2008).

Given any $\lambda>0$ in (5), the generalized representer theorem (Schölkopf et al., 2001) suggests that the solution of the functional parameters, denoted as $\{\widetilde{g}_{q,d}\}_{q=1,d=1}^{Q,D}$, of the minimization problem (7), with all other parameters held fixed, is a linear combination of the representers $\{k(\cdot,s)\}_{s\in\mathbb{S}}$ plus a linear combination of the basis functions $\{\phi_{1},\ldots,\phi_{J}\}$ of the null space of $\mathbb{H}_{k}$, i.e.,

$$\widetilde{g}_{q,d}(\cdot)=\sum_{s\in\mathbb{S}}\gamma_{s}k(\cdot,s)+\sum_{j=1}^{J}\alpha_{j}\phi_{j}(\cdot),\quad\|\phi_{j}\|_{\mathbb{H}_{k}}=0,$$

(10)

where we omit the subscript $(q,d)$ for the coefficient $\gamma_{s},\alpha_{j}$ for brevity but they are essentially different for each $(q,d)$. We assume that the null space of $\mathbb{H}_{k}$ contains only the zero function for the remainder of the paper. As a consequence of (10), the minimization problem in (7) can be reduced to a finite-dimensional Kernel Ridge Regression (KRR) problem. We summarize the discussion above in the proposition below:

We give the proof in the supplemental material. After introducing the functional norm penalty, the original tensor coefficient is now converted to a linear combination of the representer functions with the relationship that $[\mathbfcal{G}_{q}]_{ijd}=\langle[\mathbf{K}]_{u:}^{\top},[\boldsymbol{\Gamma}_{q}]_{:d}\rangle$ where $u=i+(j-1)M$.

We attempt to solve the minimization problem in (11) with an alternating minimization algorithm (Attouch et al., 2013) where we update one block of parameters at a time while keeping the others fixed. We update the parameters following the order of: $\mathbf{A}_{1}\rightarrow\mathbf{B}_{1}\rightarrow\cdots\rightarrow\mathbf{A}_{P}\rightarrow\mathbf{B}_{P}\rightarrow\boldsymbol{\Gamma}_{1}\rightarrow\cdots\rightarrow\boldsymbol{\Gamma}_{Q}\rightarrow\boldsymbol{\Sigma}_{r}\rightarrow\boldsymbol{\Sigma}_{c}\rightarrow\mathbf{A}_{1}\rightarrow\cdots$. We choose the alternating minimization algorithm for its simplicity and efficiency. Each step of the algorithm conducts exact minimization over one block of the parameters, leading to a non-increasing sequence of the objective function, which guarantees the convergence of the algorithm towards a local stationary point.

To solve the optimization problem in (11) for $\mathbf{A}_{p}$ at the $(l+1)^{\rm{th}}$ iteration, it suffices to solve the following least-square problem:

$$\min_{\mathbf{A}_{p}}\left\{\sum_{t\in[T]}\text{tr}\left(\widetilde{\mathbf{X}}_{t}(\mathbf{A}_{p})^{\top}\left(\boldsymbol{\Sigma}_{r}^{(l)}\right)^{-1}\widetilde{\mathbf{X}}_{t}(\mathbf{A}_{p})\left(\boldsymbol{\Sigma}_{c}^{(l)}\right)^{-1}\right)\right\},$$

(12)

where $\widetilde{\mathbf{X}}_{t}(\mathbf{A}_{p})$ is the residual matrix when predicting $\mathbf{X}_{t}$:

	$\displaystyle\widetilde{\mathbf{X}}_{t}(\mathbf{A}_{p})$	$\displaystyle=\mathbf{X}_{t}-\sum_{p^{{}^{\prime}}<p}\mathbf{A}_{p^{{}^{\prime}}}^{(l+1)}\mathbf{X}_{t-p^{{}^{\prime}}}\left(\mathbf{B}_{p^{{}^{\prime}}}^{(l+1)}\right)^{\top}-\sum_{p^{{}^{\prime}}>p}\mathbf{A}_{p^{{}^{\prime}}}^{(l)}\mathbf{X}_{t-p^{{}^{\prime}}}\left(\mathbf{B}_{p^{{}^{\prime}}}^{(l)}\right)^{\top}$
		$\displaystyle-\sum_{q\in[Q]}\mathbfcal{G}_{q}^{(l)}\bar{\times}\mathbf{z}_{t-q}-\mathbf{A}_{p}\mathbf{X}_{t-p}\left(\mathbf{B}_{p}^{(l)}\right)^{\top}=\widetilde{\mathbf{X}}_{t,-p}-\mathbf{A}_{p}\mathbf{X}_{t-p}\left(\mathbf{B}_{p}^{(l)}\right)^{\top}$

and we use $\widetilde{\mathbf{X}}_{t,-p}$ to denote the partial residual excluding the term involving $\mathbf{X}_{t-p}$ and use $\mathbfcal{G}_{q}^{(l)}$ to denote the tensor coefficient satisfying $[\mathbfcal{G}_{q}^{(l)}]_{ijd}=\langle[\mathbf{K}]_{u:}^{\top},[\boldsymbol{\Gamma}_{q}^{(l)}]_{:d}\rangle$, with $u=i+(j-1)M$. The superscript $l$ represents the value at the $l^{\rm{th}}$ iteration. To simplify the notation, we define $\boldsymbol{\Phi}(\mathbf{A}_{t},\mathbf{B}_{t},\boldsymbol{\Sigma})=\sum_{t}\mathbf{A}_{t}^{\top}\boldsymbol{\Sigma}^{-1}\mathbf{B}_{t}$, where $\boldsymbol{\Sigma},\mathbf{A}_{t},\mathbf{B}_{t}$ are arbitrary matrices/vectors with conformal matrix sizes and we simply write $\boldsymbol{\Phi}(\mathbf{A}_{t},\boldsymbol{\Sigma})$ if $\mathbf{A}_{t}=\mathbf{B}_{t}$. Solving (12) yields the following updating formula for $\mathbf{A}_{p}^{(l+1)}$:

$$\mathbf{A}_{p}^{(l+1)}\leftarrow\boldsymbol{\Phi}\left(\widetilde{\mathbf{X}}_{t,-p}^{\top},\mathbf{B}_{p}^{(l)}\mathbf{X}_{t-p}^{\top},\boldsymbol{\Sigma}_{c}^{(l)}\right)\boldsymbol{\Phi}\left(\mathbf{B}_{p}^{(l)}\mathbf{X}_{t-p}^{\top},\boldsymbol{\Sigma}_{c}^{(l)}\right)^{-1}$$

(13)

Similarly, we have the following updating formula for $\mathbf{B}_{p}^{(l+1)}$:

$$\mathbf{B}_{p}^{(l+1)}\leftarrow\boldsymbol{\Phi}\left(\widetilde{\mathbf{X}}_{t,-p},\mathbf{A}_{p}^{(l+1)}\mathbf{X}_{t-p},\boldsymbol{\Sigma}_{r}^{(l)}\right)\boldsymbol{\Phi}\left(\mathbf{A}_{p}^{(l+1)}\mathbf{X}_{t-p},\boldsymbol{\Sigma}_{r}^{(l)}\right)^{-1}$$

(14)

For updating $\boldsymbol{\Gamma}_{q}$, or its vectorized version $\boldsymbol{\gamma}_{q}=\mathrm{\mathbf{vec}}\left(\boldsymbol{\Gamma}_{q}\right)$, it is required to solve the following kernel ridge regression problem:

$$\min_{\boldsymbol{\gamma}_{q}}\left\{\frac{1}{2T}\boldsymbol{\Phi}\left(\widetilde{\mathbf{x}}_{t,-q}-\left(\mathbf{z}_{t-q}^{\top}\otimes\mathbf{K}\right)\boldsymbol{\gamma}_{q},\boldsymbol{\Sigma}^{(l)}\right)+\frac{\lambda}{2}\boldsymbol{\gamma}_{q}^{\top}\left(\mathbf{I}_{D}\otimes\mathbf{K}\right)\boldsymbol{\gamma}_{q}\right\},$$

where $\boldsymbol{\Sigma}^{(l)}=\boldsymbol{\Sigma}_{c}^{(l)}\otimes\boldsymbol{\Sigma}_{r}^{(l)}$ and $\widetilde{\mathbf{x}}_{t,-q}$ is the vectorized partial residual of $\mathbf{X}_{t}$ by leaving out the lag-$q$ auxiliary predictor, defined in a similar way as $\widetilde{\mathbf{X}}_{t,-p}$. Solving the kernel ridge regression leads to the following updating formula for $\boldsymbol{\gamma}_{q}^{(l+1)}$:

$$\boldsymbol{\gamma}_{q}^{(l+1)}\leftarrow\left[\left(\sum_{t\in[T]}\mathbf{z}_{t-q}\mathbf{z}_{t-q}^{\top}\right)\otimes\mathbf{K}+\lambda T\left(\mathbf{I}_{D}\otimes\boldsymbol{\Sigma}^{(l)}\right)\right]^{-1}\left[\sum_{t\in[T]}\left(\mathbf{z}_{t-q}\otimes\widetilde{\mathbf{x}}_{t,-q}\right)\right].$$

(15)

The step in (15) can be slow since one needs to invert a square matrix of size $MND\times MND$. In Section 3.2, we propose an approximation to (15) to speed up the computation under high dimensionality.

The updating rule of $\boldsymbol{\Sigma}_{r}^{(l+1)}$ and $\boldsymbol{\Sigma}_{c}^{(l+1)}$ can be easily derived by taking their derivative in (11) and setting it to zero. Specifically, we have:

	$\displaystyle\boldsymbol{\Sigma}_{r}^{(l+1)}$	$\displaystyle\leftarrow\frac{1}{NT}\boldsymbol{\Phi}\left(\widetilde{\mathbf{X}}_{t}^{\top},\boldsymbol{\Sigma}_{c}^{(l)}\right)$		(16)
	$\displaystyle\boldsymbol{\Sigma}_{c}^{(l+1)}$	$\displaystyle\leftarrow\frac{1}{MT}\boldsymbol{\Phi}\left(\widetilde{\mathbf{X}}_{t},\boldsymbol{\Sigma}_{r}^{(l+1)}\right).$		(17)

where $\widetilde{\mathbf{X}}_{t}$ is the full residual when predicting $\mathbf{X}_{t}$.

The algorithm cycles through (13), (14), (15), (16) and (17) and terminates when $\mathbf{B}_{p}^{(l)}\otimes\mathbf{A}_{p}^{(l)}$, $\mathbfcal{G}_{q}^{(l)}$, $\boldsymbol{\Sigma}_{c}^{(l)}\otimes\boldsymbol{\Sigma}_{r}^{(l)}$ have their relative changes between iterations fall under a pre-specified threshold. We make two additional remarks on the algorithm:

The iterative algorithm in Section 3.1 requires inverting an $MND\times MND$ matrix in (15) when updating $\boldsymbol{\gamma}_{q}$, i.e., the coefficients of the representer functions $k(\cdot,s)$. One way to reduce the computational complexity without any approximation is to divide the step of updating $\boldsymbol{\gamma}_{q}=[\boldsymbol{\gamma}_{q,1}^{\top}:\cdots:\boldsymbol{\gamma}_{q,D}^{\top}]^{\top}$ to updating one block of parameters at a time following the order of $\boldsymbol{\gamma}_{q,1}\rightarrow\cdots\rightarrow\boldsymbol{\gamma}_{q,D}$. However, such a procedure requires inverting a matrix of size $MN\times MN$, which could still be high-dimensional.

To circumvent the issue of inverting large matrices, we can approximate the linear combination of all $MN$ representers using a set of $R<<MN$ basis functions, i.e., $\mathbf{K}\boldsymbol{\gamma}_{q,d}\approx\mathbf{K}_{R}\boldsymbol{\theta}_{q,d}$, where $\mathbf{K}_{R}\in\mathbb{R}^{MN\times R},\boldsymbol{\theta}_{q,d}\in\mathbb{R}^{R}$. For example, one can reduce the spatial resolution by subsampling a fraction of the rows and columns of the matrix and only use the representers at the subsampled “knots” as the basis functions. In this subsection, we consider an alternative approach by truncating the Mercer decomposition in (8). A similar technique can be found in Kang et al. (2018).

Given the eigen-decomposition of $k(\cdot,\cdot)$ in (8), one can truncate the decomposition at the $R^{\rm{th}}$ largest eigenvalue $\lambda_{R}$ and get an approximation: $k(\cdot,\cdot)\approx\sum_{r\leq R}\lambda_{r}\psi_{r}(\cdot)\psi_{r}(\cdot)$. We will use the set of eigen-functions $\{\psi_{1}(\cdot),\ldots,\psi_{R}(\cdot)\}$ for faster computation. The choice of $R$ depends on the decaying rate of the eigenvalue sequence $\{\lambda_{r}\}_{r=1}^{\infty}$ (thus the smoothness of the underlying functional parameters). Our simulation result shows that the estimation and prediction errors shrink monotonically as $R\rightarrow\infty$. Therefore, $R$ can be chosen based on the computational resources available. The kernel truncation speeds up the computation at the cost of providing an overly-smoothed estimator, as we demonstrate empirically in the supplemental material.

Given the kernel truncation, any functional parameter $g_{q,d}(\cdot)$ is now approximated as: $g_{q,d}(\cdot)\approx\sum_{r\in[R]}[\boldsymbol{\theta}_{q,d}]_{r}\psi_{r}(\cdot)$. The parameter to be estimated now is $\boldsymbol{\Theta}_{q}=[\boldsymbol{\theta}_{q,1};\cdots;\boldsymbol{\theta}_{q,D}]\in\mathbb{R}^{R\times D}$, whose dimension is much lower than before ($\boldsymbol{\Gamma}_{q}\in\mathbb{R}^{MN\times D}$). Estimating $\boldsymbol{\Theta}_{q}$ requires solving a ridge regression problem, and the updating formula for generating $\mathrm{\mathbf{vec}}\left(\boldsymbol{\Theta}_{q}\right)$ at the $(l+1)^{\rm{th}}$ iteration can be written as:

$$\left[\boldsymbol{\Phi}\left(\mathbf{z}_{t-q}^{\top}\otimes\mathbf{K}_{R},\boldsymbol{\Sigma}^{(l)}\right)+\lambda T\left(\mathbf{I}_{D}\otimes\boldsymbol{\Lambda}_{R}^{-1}\right)\right]^{-1}\boldsymbol{\Phi}\left(\mathbf{z}_{t-q}^{\top}\otimes\mathbf{K}_{R},\widetilde{\mathbf{x}}_{t,-q},\boldsymbol{\Sigma}^{(l)}\right),$$

where $\mathbf{K}_{R}\in\mathbb{R}^{MN\times R}$ satisfies $[\mathbf{K}_{R}]_{ur}=\psi_{r}(s_{ij}),u=i+(j-1)M$, and $\boldsymbol{\Lambda}_{r}=\text{diag}(\lambda_{1},\ldots,\lambda_{R})$, with $\lambda_{r}$ being the $r^{\rm{th}}$ largest eigenvalue of the Mercer decomposition of $k(\cdot,\cdot)$. Now we only need to invert a matrix of size $RD\times RD$, which speeds up the computation.

The MARAC$(P,Q)$ model (1) has three hyper-parameters: the autoregressive lag $P$, the auxiliary covariate lag $Q$, and the RKHS norm penalty weight $\lambda$. In practice, $\lambda$ can be chosen by cross-validation while choosing $P$ and $Q$ requires a more formal model selection criterion. We propose to select $P$ and $Q$ by using information criterion, including the Akaike Information Criterion (AIC) (Akaike, 1998) and the Bayesian Information Criterion (BIC) (Schwarz, 1978). We formally define the AIC and BIC for the MARAC$(P,Q)$ model here and empirically validate their consistency via simulation experiments in Section 5.

Let $\widehat{\boldsymbol{\Theta}}$ be the set of the estimated parameters of the MARAC$(P,Q)$ model, and $\mathbf{df}_{P,Q,\lambda}$ be the effective degrees of the freedom of the model. We can then define the AIC and the BIC as follows:

	$\displaystyle\mathrm{AIC}(P,Q,\lambda)$	$\displaystyle=-2\sum_{t\in[T]}\ell(\mathbf{X}_{t}|\mathbf{X}_{t-1:P},\mathbf{z}_{t-1:Q},\widehat{\boldsymbol{\Theta}})+2\cdot\mathbf{df}_{P,Q,\lambda},$		(19)
	$\displaystyle\mathrm{BIC}(P,Q,\lambda)$	$\displaystyle=-2\sum_{t\in[T]}\ell(\mathbf{X}_{t}|\mathbf{X}_{t-1:P},\mathbf{z}_{t-1:Q},\widehat{\boldsymbol{\Theta}})+\log(T)\cdot\mathbf{df}_{P,Q,\lambda}.$		(20)

To calculate $\mathbf{df}_{P,Q,\lambda}$, we decompose it into the sum of three components: 1) for each pair of the autoregressive coefficient $\widehat{\mathbf{A}}_{p},\widehat{\mathbf{B}}_{p}$, the model has $(M^{2}+N^{2}-1)$ degrees of freedom; 2) for the noise covariance $\widehat{\boldsymbol{\Sigma}}_{r},\widehat{\boldsymbol{\Sigma}}_{c}$, the model has $(M^{2}+N^{2})$ degrees of freedom; and 3) for the auxiliary covariate functional parameters $\widehat{g}_{q,1},\ldots,\widehat{g}_{q,D}$, inspired by the kernel ridge regression estimator in (15), we define the sum of their degrees of freedom as:

$$\mathbf{df}_{q}(\widehat{g})=\text{tr}\left\{\left[\widetilde{\mathbf{K}}+\lambda\left(\mathbf{I}_{D}\otimes\widehat{\boldsymbol{\Sigma}}_{c}\otimes\widehat{\boldsymbol{\Sigma}}_{r}\right)\right]^{-1}\widetilde{\mathbf{K}}\right\},$$

where $\widetilde{\mathbf{K}}=\left(T^{-1}\sum_{t\in[T]}\mathbf{z}_{t-q}\mathbf{z}_{t-q}^{\top}\right)\otimes\mathbf{K}$. As $\lambda\rightarrow 0$, we have $\mathbf{df}_{q}(\widehat{g})\rightarrow MND$; namely each covariate has $MN$ free parameters, which then reduces to the element-wise linear regression model. Empirically, we find that the BIC is a consistent lag selection criterion for our model.

To facilitate the theoretical analysis, we make an assumption for the vector time series $\mathbf{z}_{t}$, which greatly simplifies the theoretical analysis.

Given the joint vector autoregressive model in (22), we derive the conditions for $\mathbf{x}_{t}$ and $\mathbf{z}_{t}$ to be jointly stationary in Theorem 6.

The proof of Theorem 6 follows the proof of the stationarity of the joint autoregressive model in (22). As a special case where $P=\widetilde{Q}=1$, the stationarity condition in (23) is equivalent to $\bar{\rho}(\mathbf{A}_{1})\cdot\bar{\rho}(\mathbf{B}_{1})<1$ and $\bar{\rho}(\mathbf{C}_{1})<1$, where $\bar{\rho}(\cdot)$ is the spectral radius of a square matrix.

Based on Theorem 6, the stationarity of the matrix and vector time series relies on the stationarity of the autoregressive coefficients of the MARAC$(P,Q)$ and VAR$(\widetilde{Q})$ models and the tensor coefficients $\mathbfcal{G}_{1},\ldots,\mathbfcal{G}_{Q}$ do not affect the stationarity. The MARAC model can be extended to the joint autoregressive process (22) where the dynamics of $\mathbf{X}_{t}$ and $\mathbf{z}_{t}$ are modeled jointly. We stick to the simpler case here and leave the joint autoregression to future work.

In this subsection, we establish the consistency and asymptotic normality of the MARAC model estimators under the scenario that $M,N$ are fixed. Given a fixed matrix dimensionality, the functional parameters $g_{q,d}\in\mathbb{H}_{k}$ can only be estimated at $S=MN$ fixed locations, and thus the asymptotic normality result is established for the corresponding tensor coefficient $\widehat{\mathbfcal{G}}_{q}$. In Section 4.3, we will discuss the double asymptotics when both $S,T\rightarrow\infty$. For the remainder of the paper, we denote the true model coefficient with an asterisk superscript, such as $\mathbf{A}_{1}^{*},\mathbf{B}_{1}^{*},\mathbfcal{G}_{1}^{*}$ and $\boldsymbol{\Sigma}^{*}$.

To start with, we make an assumption on the Gram matrix $\mathbf{K}$:

As a result of Assumption 7, every $\mathbfcal{G}_{q}^{*}$ has a unique kernel decomposition: $\mathrm{\mathbf{vec}}\left(\mathbfcal{G}_{q}^{*}\right)=(\mathbf{I}_{D}\otimes\mathbf{K})\boldsymbol{\gamma}_{q}^{*}$. With this additional assumption, the first theoretical result we establish is the consistency of the covariance matrix estimator $\widehat{\boldsymbol{\Sigma}}=\widehat{\boldsymbol{\Sigma}}_{c}\otimes\widehat{\boldsymbol{\Sigma}}_{r}$, which we summarize in Proposition 8.

Given this result, we can further establish the asymptotic normality of the other model estimators: