A Two-Way Transformed Factor Model for Matrix-Variate Time Series

A Two-Way Transformed Factor Model for Matrix-Variate Time Series

{\sl Keywords}: Structured factor, Eigen-analysis, Projected PCA, Kronecker product, Diverging eigenvalues, High-dimensional white noise test.

Introduction

Modern scientific studies often collect data under combinations of multiple factors. For example, neuroimaging experiments record brain activities at multiple spatial locations and multiple time points under a variety of experimental stimuli. Studies of social networks record social links for a variety of settings from multiple initiators of social activity to multiple receivers of the activity. Data such as these are naturally represented not as lists or tables of numbers, but as multi-indexed arrays, or tensors. As many types of such data are collected over time, it is natural to view them as tensor-valued time series. The matrix-variate time series is a sequence of second-order random tensors. For example, financial and economic studies often collect data from different countries with a number of economic indicators (e.g., growth rate of the gross domestic product, unemployment rate, etc.) every quarter. Therefore, it is important and interesting to develop appropriate statistical methods to analyze such data. The most commonly used approach to modeling such data is to stack the matrix into a long vector and to apply the standard multivariate methods. However, such an approach ignores the matrix structure of the data and often overlooks some important patterns embedded in the data. For example, werner2008 pointed out that after vectorizing the matrices the resulting vectors have a Kronecker structure, and ignoring this structure means that a much larger number of parameters need to be estimated. Furthermore, the dimension of a matrix-variate time series itself can become large in the current era of big data. Therefore, it is important to make use of the matrix structure and to find an effective way to reduce the number of parameters, especially when the dimension is high. There are some works on tensor time series, e.g., rogers2013 and surana2016, but these articles focus on data processing rather than on statistical properties or the high dimensional case.

In modeling vector time series, the available methods to reduce the number of parameters can be classified in two categories: regularization and dimension reduction. The former imposes some conditions on the structure of a vector autoregressive moving-average (VARMA) model, and the latter assumes there is a lower dimensional representation for the high-dimensional process. For the regularization methods, some special structures are often imposed on the VARMA model. For example, Chapter 4 of Tsay_2014 and the references therein discussed two different canonical structures. Davis2012 studied the VAR model with sparse coefficient matrices based on partial spectral coherence. The Lasso regularization has also been applied to VAR models, see ShojaieMichailidis_2010, SongBickel_2011, and HanTsay2020, among others. For dimension reduction, popular methods include the canonical correlation analysis (CCA) of BoxTiao_1977, the principle component analysis (PCA) of StockWatson_2002, the scalar component analysis of TiaoTsay_1989. The factor model approach can be found in BaiNg_Econometrica_2002, StockWatson_2005, forni2000,forni2005, panyao2008, LamYaoBathia_Biometrika_2011, lamyao2012, gaotsay2018a,gaotsay2020a,gaotsay2018b, among others. However, none of the methods mentioned above can directly be used to model matrix-variate time series without vectorization. The matrix-variate time series has not been well studied in the literature; walden2002 handled this type of data in signal and image processing, wang2018 proposed a factor model for matrix-variate time series, which maintains and utilizes the matrix structure to achieve the dimension reduction, and chentsaychen2018 studied the constrained matrix-variate factor models by incorporating domain or prior knowledge in the model through linear constraints. However, the mechanism of the proposed matrix factor model deserves a further study and the bounded eigenvalue assumption of the covariance matrix of the vectorized idiosyncratic term is often violated in the high-dimensional setting, especially for the notable case of low signal-to-noise ratio commonly seen in finance and economics. See, for example, black1986.

The goal of this paper is to study the common dynamic dependence of matrix-variate time series from a new perspective. We first illustrate our primitive idea below and propose our approach in Section 2. Let ${\mathbf Y}_t\in\mathbb{R}^{p_1\times p_2}$ be an observable matrix-variate time series. For simplicity, we assume that ${\mathbf Y}_t$ is weakly stationary with $E({\mathbf Y}_t)$ = ${\mathbf 0}$. We postulate that there exist two full-rank matrices ${\mathbf T}_1\in R^{p_1\times p_1}$ and ${\mathbf T}_2\in R^{p_2\times p_2}$ such that ${\mathbf T}_1{\mathbf Y}_t{\mathbf T}_2'$ is of the form

where ${\mathbf F}_t\in\mathbb{R}^{r_1\times r_2}$ is a matrix-variate factor that captures the dynamic dependence of ${\mathbf Y}_t$, and ${\mathbf Z}_{12,t}$, ${\mathbf Z}_{21,t}$ and ${\mathbf Z}_{22,t}$ are matrix-variate idiosyncratic components, which are white noise processes. Equivalently, model ((ref)) is to seek two nonsingular transformation matrices ${\mathbf L}=({\mathbf L}_1,{\mathbf L}_2):={\mathbf T}_1^{-1}$ and ${\mathbf R}=({\mathbf R}_1,{\mathbf R}_2):={\mathbf T}_2^{-1}$ with ${\mathbf L}_1\in \mathbb{R}^{p_1\times r_1}$ and ${\mathbf R}_1\in \mathbb{R}^{p_2\times r_2}$ such that ${\mathbf L}_1$ and ${\mathbf R}_1$ are the front and back loading matrices associated with the common factors. To see the rationale of model ((ref)), let $\textnormal{vec}(\cdot)$ be the conventional vectorization operator that converts a matrix to a vector by stacking columns of the matrix on top of each other. By the basic properties of Kronecker product, we rewrite the model in the following vector form:

where ${\mathbf A}=[{\mathbf R}_1\otimes{\mathbf L}_1,{\mathbf R}_1\otimes{\mathbf L}_2,{\mathbf R}_2\otimes{\mathbf L}_1,{\mathbf R}_2\otimes{\mathbf L}_2]\in\mathbb{R}^{p_1p_2\times p_1p_2}$, ${\mathbf f}_t=\textnormal{vec}({\mathbf F}_t)\in\mathbb{R}^{r_1r_2}$ and ${\mathbf z}_t=[\textnormal{vec}({\mathbf Z}_{21,t})',\textnormal{vec}({\mathbf Z}_{12,t})',\textnormal{vec}({\mathbf Z}_{22,t})']'\in\mathbb{R}^{p_1p_2-r_1r_2}$. For identifiability, we assume that both ${\mathbf f}_t$ and ${\mathbf z}_t$ have zero mean and identity covariance matrices. This is a special case of the model considered in gaotsay2018b for vector time series by assuming that the covariance of the vectorized data has a Kronecker structure. That is, we expect there exists a transformation matrix ${\mathbf A}^{-1}$ with a Kronecker structure such that ${\mathbf A}^{-1}{\mathbf y}_t=({\mathbf f}_t',{\mathbf z}_t')'$, and this can be done via canonical correlation analysis between ${\mathbf y}_t$ and its past lagged variables, and the resulting vector $({\mathbf f}_t',{\mathbf z}_t')'$ are contemporaneously uncorrelated with an identity covariance matrix. See the discussions in gaotsay2018b and TiaoTsay_1989. The structure of ${\mathbf A}$ is different from that in gaotsay2018b in order to preserve the structure of the matrix-valued data. Consequently, the main task of the proposed method is to estimate ${\mathbf L}_1$, ${\mathbf R}_1$ and to recover the matrix factor ${\mathbf F}_t$.

To summarize, we propose in this paper a new framework for statistical modeling of matrix-variate time series based on the aforementioned motivation and the concepts of factor models. We reparametrize the model by compressing the strengths of the linear transformation matrices to the corresponding factors and idiosyncratic components, and the resulting front and back loading matrices associated with the common factor and the idiosyncratic terms are all half-orthonormal. Our first step is to find common orthonormal projections for the row and column vectors respectively based on an eigen-analysis of certain matrices, and the top few projected coordinates form a matrix-variate common factor process. The rest of the projected coordinates form a matrix-variate white noise process. When recovering the factor matrix, we introduce a two-way projected principal component analysis (PCA) to estimate the loading matrices associated with the idiosyncratic matrix; see Section 2 for details. In the presence of diverging noise components, the projected PCA helps to mitigate the effect of the idiosyncratic component in estimating the common factor matrix. Furthermore, we propose a diagonal-path selection method to estimate the order (dimension) of the factor matrix based on a white noise testing procedure. The testing procedure is more reasonable and statistically interpretable than the ratio-based method in wang2018, which essentially follows the method in LamYaoBathia_Biometrika_2011. Consequently, the extracted matrix-variate factors capture most of the dynamic dependence of the data and is useful if one is interested in out-of-sample forecasting of matrix-variate time series. An autoregressive type of model can be used to model the low-dimensional common factor process. See, for example, the model in chenxiaoyang2020. Asymptotic properties of the proposed method are established for both fixed and diverging dimensions as the sample size $n$ tends to infinity. We use simulated and real examples to assess the performance of the proposed method.

The rest of the paper is organized as follows. We introduce the proposed model and estimation methodology in Section 2. In Section 3, we study the theoretical properties of the proposed model and its associated estimates. Numerical illustrations with both simulated and real data sets are reported in Section 4. Section 5 provides concluding remarks. All technical proofs are given in an Appendix. Throughout the article, we use the following notation. For a $p\times 1$ vector ${\mathbf u}=(u_1,..., u_p)',$ $||{\mathbf u}||_2 =\|{\mathbf u}'\|_2= (\sum_{i=1}^{p} u_i^2)^{1/2} $ is the Euclidean norm, and ${\mathbf I}_p$ denotes a $p\times p$ identity matrix. For a matrix ${\mathbf H}=(h_{ij})$, $|{\mathbf H}|_\infty=\max_{i,j}|h_{ij}|$, $\|{\mathbf H} \|_F=\sqrt{\sum_{i,j}h_{ij}^2}$ is the Frobenius norm, $\|{\mathbf H} \|_2=\sqrt{\lambda_{\max} ({\mathbf H}' {\mathbf H} ) }$ is the operator norm, where $\lambda_{\max} (\cdot) $ denotes for the largest eigenvalue of a matrix, and $\|{\mathbf H}\|_{\min}$ is the square root of the minimum non-zero eigenvalue of ${\mathbf H}'{\mathbf H}$. The superscript ${'}$ denotes the transpose of a vector or matrix. We also use the notation $a\asymp b$ to denote $a=O(b)$ and $b=O(a)$.

Models and Methodology

Setting

Let ${\mathbf Y}_t= [y_{ij,t}] = ({\mathbf y}_{1,t},...,{\mathbf y}_{p_2, t})$ be an observable $p_1\times p_2$ matrix-variate time series with ${\mathbf y}_{jt}=(y_{1j,t},...,y_{p_1 j, t})'\in \mathbb{R}^{p_1}$ and $E({\mathbf y}_{jt})={\bf 0}$ for $1\leq j\leq p_2$. We assume ${\mathbf Y}_t$ admits a latent structure:

where ${\mathbf F}_t\in \mathbb{R}^{r_1\times r_2}$ is a matrix-variate common factor process, ${\mathbf Z}_{12,t}\in \mathbb{R}^{r_1\times v_2}$, ${\mathbf Z}_{21,t}\in\mathbb{R}^{v_1\times r_2}$, and ${\mathbf Z}_{22,t}\in \mathbb{R}^{v_1\times v_2}$ are matrix-variate idiosyncratic noise processes with $r_1+v_1=p_1$ and $r_2+v_2=p_2$. ${\mathbf L}=({\mathbf L}_1,{\mathbf L}_2)\in \mathbb{R}^{p_1\times p_1}$ is the front loading matrix with ${\mathbf L}_1\in \mathbb{R}^{p_1\times r_1}$ and ${\mathbf L}_2\in \mathbb{R}^{p_1\times v_1}$, and ${\mathbf R}=({\mathbf R}_1,{\mathbf R}_2)\in \mathbb{R}^{p_2\times p_2}$ is the back loading matrix with ${\mathbf R}_1\in \mathbb{R}^{p_2\times r_2}$ and ${\mathbf R}_2\in \mathbb{R}^{p_2\times v_2}$. We assume ${\mathbf L}$ and ${\mathbf R}$ are full-rank so that ${\mathbf F}_t$, ${\mathbf Z}_{12,t}$, ${\mathbf Z}_{21,t}$, and ${\mathbf Z}_{22,t}$ can be viewed as transformed processes by applying the inverses of ${\mathbf L}$ and ${\mathbf R}'$, respectively, to the left and right of the data matrix ${\mathbf Y}_t$ as discussed in Section 1. Furthermore, letting ${\mathbf f}_t$ and ${\mathbf z}_t$ be the vectorized factor and idiosyncratic terms, we assume that $\textnormal{Cov}({\mathbf f}_t)={\mathbf I}_{r_1r_2}$ and $\textnormal{Cov}({\mathbf z}_t)={\mathbf I}_{p_1p_2-r_1r_2}$. This assumption holds because one can adjust the scales of ${\mathbf L}$ and ${\mathbf R}$ accordingly. Therefore, the three noise terms are uncorrelated with each other and individually identified. Model ((ref)) is general if one allows $r_1$, $r_2$, $v_1$ and $v_2$ to be zero, but for effective dimension reduction, $r_1$ and $r_2$ should be small and fixed positive integers. In addition, we assume ${\mathbf f}_t$ and ${\mathbf z}_s$ are uncorrelated for any $t$ and $s$. This is only for the simplicity in illustration, and it can be relaxed by imposing some dynamic dependence between ${\mathbf f}_t$ and ${\mathbf z}_s$. See gaotsay2018b for details. We do not pursue it here. Note that ${\mathbf L}$ and ${\mathbf R}$ are not uniquely identified because $c{\mathbf L}$ and ${\mathbf R}/c$, where $c\neq 0$, also holds for Equation ((ref)).

To proceed, we further decompose ${\mathbf L}$ and ${\mathbf R}$ as follows: \[{\mathbf L}_1={\mathbf A}_1{\mathbf W}_1,\quad{\mathbf L}_2={\mathbf A}_2{\mathbf W}_2,\quad{\mathbf R}_1={\mathbf P}_1{\mathbf G}_1,\,\, \text{and}\,\,{\mathbf R}_2={\mathbf P}_2{\mathbf G}_2,\] where ${\mathbf A}_i$ and ${\mathbf P}_i$ ($i=1,2$) are half orthonormal matrices, i.e., ${\mathbf A}_i'{\mathbf A}_i={\mathbf I}_{r_i}$ and ${\mathbf P}_i'{\mathbf P}_i={\mathbf I}_{v_i}$. This can be done via QR or singular value decomposition. Furthermore, let ${\mathbf X}_t={\mathbf W}_1{\mathbf F}_t{\mathbf G}_1'$, ${\mathbf E}_{21,t}={\mathbf W}_2{\mathbf Z}_{21,t}{\mathbf G}_1'$, ${\mathbf E}_{12,t}={\mathbf W}_1{\mathbf Z}_{12,t}{\mathbf G}_2'$, and ${\mathbf E}_{22,t}={\mathbf W}_2{\mathbf Z}_{22,t}{\mathbf G}_2'$, then model ((ref)) can be rewritten as

Even though ${\mathbf L}$ and ${\mathbf R}$ are full rank, ${\mathbf A}_1$ (${\mathbf P}_1$) is not orthogonal to ${\mathbf A}_2$ (${\mathbf P}_2$) in general. Note that model ((ref)) is still not identified since we can replace the triplets (${\mathbf A}_1$,${\mathbf X}_t$,${\mathbf P}_1$) by (${\mathbf A}_1{\mathbf H}_1$,${\mathbf H}_1'{\mathbf X}_t{\mathbf H}_2$, ${\mathbf P}_1{\mathbf H}_2$) for any orthonormal matrices ${\mathbf H}_1\in\mathbb{R}^{r_1\times r_1}$ and ${\mathbf H}_2\in\mathbb{R}^{r_2\times r_2}$ without altering the data generating process. The same issue exists for the idiosyncratic terms. Nevertheless the linear spaces spanned by the columns of ${\mathbf A}_i$ and ${\mathbf P}_i$, denoted respectively by $\mathcal{M}({\mathbf A}_i)$ and $\mathcal{M}({\mathbf P}_i)$, are uniquely defined and $\mathcal{M}({\mathbf A}_i)=\mathcal{M}({\mathbf L}_i)$ and $\mathcal{M}({\mathbf P}_i)=\mathcal{M}({\mathbf R}_i)$.

Common Orthonormal Projections

To illustrate our estimation method, we first introduce some notation. For $i=1,2$, let ${\mathbf B}_i$ and ${\mathbf Q}_i$ be the orthonormal complements of ${\mathbf A}_i$ and ${\mathbf P}_i$, respectively, i.e., ${\mathbf B}_i\in\mathbb{R}^{p_i\times r_i}$ and ${\mathbf Q}_i=\mathbb{R}^{p_i\times v_i}$ are half orthonormal matrices with ${\mathbf B}_i'{\mathbf A}_i=\bf 0$ and ${\mathbf Q}_i'{\mathbf P}_i=\bf{0}$. Furthermore, denote $\boldsymbol{\ell}_{i,j}$, ${\mathbf r}_{i,j}$, ${\mathbf a}_{i,j}$, ${\mathbf b}_{i,j}$, ${\mathbf p}_{i,j}$ and ${\mathbf q}_{i,j}$ the $j$-th columns of ${\mathbf L}_i$, ${\mathbf R}_i$, ${\mathbf A}_i$, ${\mathbf B}_i$, ${\mathbf P}_i$ and ${\mathbf Q}_i$, respectively, where the range of $j$ depends on the dimension of the corresponding matrix.

Let $\boldsymbol{\eta}_{t}=[\textnormal{vec}({\mathbf Y}_{t-1})',...,\textnormal{vec}({\mathbf Y}_{t-k_0})']'$ be the vector of past $k_0$ lagged values of ${\mathbf Y}_t$, where $\textnormal{vec}({\mathbf Y}_t)=({\mathbf y}_{1,t}',...,{\mathbf y}_{p_2,t}')'$ and $k_0$ is a prescribed positive integer. Define $\boldsymbol{\Sigma}_{y,ij}(k)=\mbox{Cov}({\mathbf y}_{i,t},{\mathbf y}_{j,t-k})$. We seek the direction ${\mathbf a}\in \mathbb{R}^{p_1}$ that solves the following optimization problem:

That is, we look for a common direction ${\mathbf a}$ with ${\mathbf a}'{\mathbf a}=1$ such that it maximizes the sum of the covariance between ${\mathbf a}'{\mathbf y}_{i,t}$ and the past lagged variables, which characterize the dynamic dependence of the columns. Note that \[\sum_{i=1}^{p_2}\|\mbox{Cov}({\mathbf a}'{\mathbf y}_{i,t},\boldsymbol{\eta}_t)\|_2^2={\mathbf a}'\left[\sum_{k=1}^{k_0}\sum_{i=1}^{p_2}\sum_{j=1}^{p_2}\boldsymbol{\Sigma}_{y,ij}(k)\boldsymbol{\Sigma}_{y,ij}(k)'\right]{\mathbf a}.\] Then, ${\mathbf a}$ is an eigenvector of the matrix

On the other hand, under model ((ref)), let ${\mathbf p}_{1,i\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}}$ be the $i$-th row vectors of ${\mathbf P}_1$ and define $\boldsymbol{\Sigma}_{xp,ij}(k)=\mbox{Cov}({\mathbf X}_t{\mathbf p}_{1,i\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}}',{\mathbf X}_{t-k}{\mathbf p}_{1,j\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}}')$. Then

where we assume ${\mathbf f}_t$ and ${\mathbf z}_s$ are uncorrelated for any $t$ and $s$. Therefore,

We observe that ${\mathbf M}_1{\mathbf B}_1=\bf 0$, that is, the columns of ${\mathbf B}_1$ are the eigenvectors associated with the zero eigenvalues of ${\mathbf M}_1$, and the front factor loading space $\mathcal{M}({\mathbf A}_1)$ is spanned by the eigenvectors corresponding to the $r_1$ non-zero eigenvalues of ${\mathbf M}_1$. Equivalently, the space spanned by the first $r_1$ solutions to the problem ((ref)) are just the front factor loading space $\mathcal{M}({\mathbf A}_1)$.

The $r_2$ orthonormal directions of the columns of ${\mathbf P}_1$ can be obtained by performing the same procedure on the the transpose of ${\mathbf Y}_t$'s. We can similarly construct ${\mathbf M}_2$ as ${\mathbf M}_1$ in ((ref)) such that ${\mathbf M}_2{\mathbf Q}_1=\bf 0$, and therefore, $\mathcal{M}({\mathbf P}_1)$ is the space spanned by the first $r_2$ non-zero eigenvectors of ${\mathbf M}_2$. We omit the details here.

Two-Way Projected Principal Component Analysis

In this section, we introduce the idea of a 2-way projected PCA in order to recover the true factor matrix ${\mathbf X}_t$ and, hence, ${\mathbf F}_t$. Using the notation in Section 2.2, it follows from model ((ref)) that

which implies that ${\mathbf B}_1'{\mathbf Y}_t{\mathbf Q}_1$ is a matrix-variate white noise process, and hence $\{{\mathbf b}_{1,i}'{\mathbf Y}_t{\mathbf q}_{1,j}|t=0,\pm 1,...\}$ is a univariate white noise process for all $1\leq i\leq v_1$ and $1\leq j\leq v_2$. Furthermore,

Therefore, ${\mathbf B}_2'{\mathbf Y}_t$ and ${\mathbf Y}_t{\mathbf Q}_2$ are uncorrelated with ${\mathbf B}_1'{\mathbf Y}_t{\mathbf Q}_1$ defined in ((ref)). Let $\boldsymbol{\Omega}_{y_i}=\mbox{Cov}({\mathbf y}_{i,t},\textnormal{vec}({\mathbf Y}_t))$ and $\boldsymbol{\Omega}_{e_{22,i}p}=\mbox{Cov}({\mathbf E}_{22,t}{\mathbf p}_{2,i\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}}',\textnormal{vec}({\mathbf E}_{22,t}))$, where ${\mathbf p}_{2,i\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}}$ is the $i$-th row vector of ${\mathbf P}_2$. It follows from ((ref)) and ((ref)) that

For each $1\leq i\leq p_2$, ${\mathbf B}_2'{\mathbf y}_{i,t}$ is uncorrelated with $\textnormal{vec}({\mathbf B}_1'{\mathbf Y}_t{\mathbf Q}_1)$ and we define

from which we can see, via ((ref)), that ${\mathbf S}_1{\mathbf B}_2={\mathbf 0}$. In addition, the rank of ${\mathbf S}_1\in \mathbb{R}^{p_1\times p_1}$ is $v_1$, therefore, ${\mathbf B}_2$ contains all the eigenvectors corresponding to the zero eigenvalues of ${\mathbf S}_1$. From the form of ${\mathbf S}_1$, we can see that the middle component contains the information of the noise, and we seek the direction ${\mathbf b}_{2,j}\in \mathbb{R}^{p_2}$ such that ${\mathbf b}_{2,j}'{\mathbf y}_{i,t}$ minimizes the covariance between the projected direction and the noise, and therefore, ${\mathbf b}_{2,j}'{\mathbf y}_{i,t}$ contains the information of the signal ${\mathbf F}_t$.

Similarly, we can construct ${\mathbf S}_2$ such that ${\mathbf S}_2{\mathbf Q}_2=0$, and ${\mathbf Q}_2$ contains all the eigenvectors associated with the zero eigenvalues of ${\mathbf S}_2$ at the population level. Furthermore, if ${\mathbf A}_1$, ${\mathbf P}_1$, ${\mathbf B}_2$ and ${\mathbf Q}_2$ are known, it follows from ((ref)) that

and consequently,

where ${\mathbf B}_2'{\mathbf A}_1\in\mathbb{R}^{r_1\times r_1}$ and ${\mathbf P}_1'{\mathbf Q}_2\in\mathbb{R}^{r_2\times r_2}$ are two invertible matrices. To see this, note that ${\mathbf L}$ is a matrix of full rank, thus there exist matrices ${\mathbf H}_1\in \mathbb{R}^{r_1\times r_1}$ and ${\mathbf H}_2\in\mathbb{R}^{v_1\times r_1}$ such that \[{\mathbf B}_2={\mathbf L}_1{\mathbf H}_1+{\mathbf L}_2{\mathbf H}_2={\mathbf A}_1{\mathbf W}_1{\mathbf H}_1+{\mathbf A}_2{\mathbf W}_2{\mathbf H}_2.\] Then, \[{\mathbf I}_{r_1}={\mathbf B}_2'{\mathbf B}_2={\mathbf B}_2'{\mathbf A}_1{\mathbf W}_1{\mathbf H}_1,\] implying that rank$({\mathbf B}_2'{\mathbf A}_1)=r_1$ which is of full rank. The invertibility of ${\mathbf P}_1'{\mathbf Q}_2$ follows from a similar argument.

Estimation

In practice, given a sample $\{{\mathbf Y}_t:t=1,...,n\}$, the goal is to estimate ${\mathbf A}_1$ and ${\mathbf P}_1$ or equivalently $\mathcal{M}({\mathbf A}_1)$ and $\mathcal{M}({\mathbf P}_1)$, the dimension $(r_1,r_2)$ of the matrix factor, and to recover the latent factor matrix process ${\mathbf X}_t$. To illustrate the main idea, we first assume $(r_1,r_2)$ is known, but propose a way to estimate them in the next subsection.

For the estimation of ${\mathbf A}_1$ and ${\mathbf P}_1$, we construct the sample version of ${\mathbf M}_1$ defined in ((ref)) as follows:

where

and $\bar{{\mathbf y}}_i=n^{-1}\sum_{t=1}^n{\mathbf y}_{i,t}$ which is essentially ${\bf 0}$ if the data are centered. Then, $\mathcal{M}({\mathbf A}_1)$ can be estimated by $\mathcal{M}(\widehat{\mathbf A}_1)$, where $\widehat{\mathbf A}_1=(\widehat{\mathbf a}_{1,1},...,\widehat{\mathbf a}_{1,r_1})$ with $\widehat{\mathbf a}_{1,1},...,\widehat{\mathbf a}_{1,r_1}$ being the eigenvectors corresponding to the $r_1$ largest eigenvalues of $\widehat{\mathbf M}_1$. Consequently, the orthogonal space $\mathcal{M}(\widehat{\mathbf B}_1)$ can be similarly obtained by $\widehat{\mathbf B}_1=(\widehat{\mathbf b}_{1,1},...,\widehat{\mathbf b}_{1,v_1})$, where $\widehat{\mathbf b}_{1,1},...,\widehat{\mathbf b}_{1,v_1}$ are the eigenvectors corresponding to the $v_1$ smallest eigenvalues of $\widehat{\mathbf M}_1$.

By a similar procedure on $\{{\mathbf Y}_t',t=1,...,n\}$, we can construct $\widehat{\mathbf M}_2$, and the estimator $\widehat{\mathbf P}_1$ for ${\mathbf P}_1$ is then obtained. Once we have the estimators $\widehat{\mathbf A}_1$ and $\widehat{\mathbf P}_1$, we consider methods for obtaining the estimators of ${\mathbf B}_2$ and ${\mathbf Q}_2$. The choices of the estimators $\widehat{\mathbf B}_2$ and $\widehat{\mathbf Q}_2$ are different for small and large dimensions. We only discuss the cases when $p_1$ and $p_2$ are both small or large, and the case when one of them is small can be solved by applying both methods jointly. Let

where $\widehat\boldsymbol{\Omega}_{y_i}$ is the sample estimator of $\boldsymbol{\Omega}_{y_i}$ defined in Section 2.3. When $p_1$ and $p_2$ are small, we perform an eigen-analysis on $\widehat{\mathbf S}_1$, and $\widehat{\mathbf B}_2=(\widehat{\mathbf b}_{2,1},...,\widehat{\mathbf b}_{2,r_1})$, where $\widehat{\mathbf b}_{2,1},...,\widehat{\mathbf b}_{2,r_1}$ are the eigenvectors of ${\mathbf S}_1$ corresponding to its $r_1$ smallest eigenvalues. We can similarly obtain $\widehat{\mathbf Q}_1$ based on the eigen-analysis on $\widehat{\mathbf S}_2$, which is calculated based on the transposed data ${\mathbf Y}_t'$.

When the dimensions $p_1$ and $p_2$ are relatively large, the choices of $\widehat{\mathbf B}_2$ and $\widehat{\mathbf Q}_2$ by selecting the eigenvectors associated with the smallest $r_1$ and $r_2$ eigenvalues of $\widehat{\mathbf S}_1$ and $\widehat{\mathbf S}_2$ respectively may not fare well because the linear spaces spanned by the chosen vectors are not consistent to the true ones in the high-dimensional case. Suppose that the elements $z_{ij,t}$ of ${\mathbf Z}_{22, t}$ are independent of each other for $1\leq i\leq v_1$ and $1\leq j\leq v_2$. A reasonable assumption is that the top eigenvalues of the covariance matrix of the idiosyncratic component $\textnormal{vec}({\mathbf L}_2{\mathbf Z}_{22,t}{\mathbf R}_2')$ or equivalently $\textnormal{vec}({\mathbf A}_2{\mathbf E}_{22,t}{\mathbf P}_2')$ are diverging. Thus, we assume the top singular values of ${\mathbf L}_2$ and ${\mathbf R}_2$ are diverging. See also Assumption 4 in Section 3. We can partition the singular vectors ${\mathbf A}_2=({\mathbf A}_{21},{\mathbf A}_{22})$ and ${\mathbf P}_2=({\mathbf P}_{21},{\mathbf P}_{22})$ with ${\mathbf A}_{21}\in \mathbb{R}^{p_1\times k_1}$ and ${\mathbf P}_{21}\in \mathbb{R}^{p_2\times k_2}$ which correspond to the $k_1$ and $k_2$ diverging singular values of ${\mathbf L}_2$ and ${\mathbf R}_2$, respectively. Let ${\mathbf B}_2^*=({\mathbf A}_{22},{\mathbf B}_{2})\in \mathbb{R}^{p_1\times (p_1-k_1)}$ and ${\mathbf Q}_2^*=({\mathbf P}_{22},{\mathbf Q}_2)\in \mathbb{R}^{p_2\times(p_2-k_2)}$. Under the assumption that the top $k_1$ singular values of ${\mathbf L}_2$ and $k_2$ of ${\mathbf R}_2$ are diverging, we can consistently estimate the spaces $\mathcal{M}({\mathbf A}_{21})$ and $\mathcal{M}({\mathbf Q}_{21})$ and hence their orthogonal parts $\mathcal{M}({\mathbf B}_2^*)$ and $\mathcal{M}({\mathbf Q}_2^*)$. From the above discussion, $\mathcal{M}({\mathbf B}_2)$ and $\mathcal{M}({\mathbf Q}_2)$ are subspaces of $\mathcal{M}({\mathbf B}_2^*)$ and $\mathcal{M}({\mathbf Q}_2^*)$, respectively. Once we have the consistent estimators for ${\mathbf B}_2^*$ and ${\mathbf Q}_2^*$, denoted by $\widehat{{\mathbf B}}_2^*$ and $\widehat{{\mathbf Q}}_2^*$, respectively, there are half orthonormal matrices $\boldsymbol{\Xi}_1\in\mathbb{R}^{(p_1-k_1)\times r_1 }$ and $\boldsymbol{\Xi}_2\in\mathbb{R}^{(p_2-k_2)\times r_2}$ such that $\widehat{\mathbf B}_2=\widehat{{\mathbf B}}_2^*\boldsymbol{\Xi}_1$ and $\widehat{\mathbf Q}_2=\widehat{{\mathbf Q}}_2^*\boldsymbol{\Xi}_2$. In practice, it is not easy to find $\boldsymbol{\Xi}_1$ and $\boldsymbol{\Xi}_2$ such that $\widehat{\mathbf B}_2$ and $\widehat{\mathbf Q}_2$ are consistent to ${\mathbf B}_2$ and ${\mathbf Q}_2$. Nevertheless, any choices of $\boldsymbol{\Xi}_1$ and $\boldsymbol{\Xi}_2$ can mitigate the diverging effect of the top eigenvalues since they are all orthogonal to $\widehat{\mathbf A}_{21}$ and $\widehat{\mathbf P}_{21}$, respectively. Thus, we only need to guarantee the invertiblities of the the matrices $\widehat{\mathbf B}_2'\widehat{\mathbf A}_1$ and $\widehat{\mathbf P}_1'\widehat{\mathbf Q}_2$ in order to recover the latent factors.

In practice, with the estimators $\widehat{{\mathbf B}}_2^*$ and $\widehat{{\mathbf Q}}_2^*$, the columns of $\boldsymbol{\Xi}_1$ are chosen as the $r_1$ eigenvectors of $\widehat{{\mathbf B}}_2^*{'}\widehat{\mathbf A}_1\widehat{\mathbf A}_1'\widehat{{\mathbf B}}_2^*$ corresponding the $r_1$ largest eigenvalues, and the columns of $\boldsymbol{\Xi}_2$ are the $r_2$ eigenvectors of $\widehat{{\mathbf Q}}_2^*{'}\widehat{\mathbf P}_1\widehat{\mathbf P}_1'\widehat{{\mathbf Q}}_2^*$ corresponding to the largest $r_2$ eigenvalues. These choices guarantee that both $\widehat{\mathbf B}_2'\widehat{\mathbf A}_1$ and $\widehat{\mathbf P}_1'\widehat{\mathbf Q}_2$ behave well in practical calculations. Finally, we recover the latent factor matrix as

With $\widehat{\mathbf A}_1$, $\widehat{\mathbf P}_1$ and the estimated factor process $\widehat{\mathbf X}_t$, we can make an $h$-step ahead prediction for the ${\mathbf Y}_t$ series using the formula $\widehat{\mathbf Y}_{n+h}=\widehat{\mathbf A}_1\widehat{\mathbf X}_{n+h}\widehat{\mathbf P}_1'$, where $\widehat{\mathbf X}_{t+h}$ is an $h$-step ahead forecast for ${\mathbf X}_t$ based on the estimated past values $\widehat{\mathbf X}_1,...,\widehat{\mathbf X}_n$. This can be done, for example, by fitting a matrix-autoregressive model to $\{\widehat{\mathbf X}_1,...,\widehat{\mathbf X}_n\}$ as, for example, the one introduced in chenxiaoyang2020.

Diagonal-Path Selections of the Order of Factor Matrix

The estimation of ${\mathbf A}_1$, ${\mathbf P}_1$, and ${\mathbf X}_t$ of the prior sections are based on given $r_1$ and $r_2$, which are unknown in practice. To the best of our knowledge, there is no efficient method available to estimate them in the literature. The most relevant one is the ratio-based method of wang2018, but it can be shown that the method is not appropriate when the top eigenvalues of the covariance of the idiosyncratic term are diverging. See the simulation results in Section 4. For the vector factor models, there are some methods available. See, for example, the information criterion in Bai and Ng (2002) and Bai (2003), the random matrix theory method in onatski2010, the ratio-based method in Lam and Yao (2012), the canonical correlation analysis in gaotsay2018a, and the white noise testing approach in gaotsay2018b, among others. However, those methods cannot apply to the matrix-factor models directly.

In this section, we propose a diagonal-path method to search the dimension $(r_1,r_2)$ by modifying the approach of gaotsay2018b. The idea of our method follows from equation ((ref)) that ${\mathbf B}_1'{\mathbf Y}_t{\mathbf Q}_1$ is a matrix-variate white noise process. Let $\widehat \boldsymbol{\Gamma}_1$ and $\widehat\boldsymbol{\Gamma}_2$ be the matrices of eigenvectors (in the decreasing order of corresponding eigenvalues) of the sample matrix $\widehat{\mathbf M}_1$ in ((ref)) and $\widehat {\mathbf M}_2$, respectively. Define $\widehat{\mathbf W}_t=\widehat\boldsymbol{\Gamma}_1{'}{\mathbf Y}_t\widehat\boldsymbol{\Gamma}_2$ and let $\widehat{\mathbf W}_t(i,j)\in\mathbb{R}^{(p_1-i+1)\times (p_2-j+1)}$ be the lower-right submatrix consisting of the $i$-th to the $p_1$-th rows and the $j$-th to the $p_2$-th columns of $\widehat{\mathbf W}_t$, and $\widehat{\mathbf W}_t^*(i,j)\in \mathbb{R}^{(i-1)\times (j-1)}$ be the upper-left submatrix of $\widehat{\mathbf W}_t$. Our test procedure searches the order $(i,j)$ such that $\widehat{\mathbf W}_t^*(i,j)$ consists of all the factors and the remaining elements of $\widehat{\mathbf W}_t$ are white noises. The estimate of $(r_1,r_2)$ is then $(i-1,j-1)$. The testing procedure is discussed below, and the test statistic used depends on the dimension $p_1p_2$.

If the dimension $p_1p_2$ is small, implying that ${\mathbf Y}_t$ is a low dimensional matrix, we recommend using the well-known Ljung-Box statistic $Q_s(m)$ for multivariate time series, where $s$ and $m$ denote the dimension of the vector and the number of lags used. See, for example, hosking1980 and Tsay_2014. Specifically, we first search the minimum of $r_1$ and $r_2$ along the diagonal of $\widehat{\mathbf W}_t$. Consider the null hypothesis \[H_0(l): \textnormal{vec}(\widehat{\mathbf W}_t(l,l))\,\, \text{is a vector white noise}, \] with type-I error $\alpha$. $H_0(l)$ is rejected if $Q_{d_l}(m)\geq \chi_{d_l^2 m,1-\alpha}^2$, where $d_l=(p_1-l+1)(p_2-l+1)$ is the dimension of $\textnormal{vec}(\widehat{\mathbf W}_t(l,l))$ and $\chi_{d_l^2 m,1-\alpha}^2$ is the $(1-\alpha)$-th quantile of a chi-squared distribution with $d_l^2m$ degrees of freedom. We start with $l=1$. If $H_0(1)$ is rejected, we increase $l$ by 1 and repeat the test until we cannot reject $H_0(l)$, and denote the resulting order as $l^*$. Two situations can happen. If $l^*=\min(p_1,p_2)$ and we still reject $H_0(l^*)$, we fix one dimension (say $p_1$ when $p_1=l^*$), and test whether $\textnormal{vec}(\widehat{\mathbf W}_t(p_1,p_1+j))$ is white noise or not by starting with $j=1$ until we cannot reject $H_0$. If $l^*<\min(p_1,p_2)$, then we perform a back testing to determine the maximum order of the factor matrix. That is, we first test whether $\textnormal{vec}(\widehat{\mathbf W}_t(l^*-1+i,l^*-1))$ is a vector white noise starting with $i=1$. Increase $i$ by 1 and repeat the test until we cannot reject $H_0$ at $i=i^*$. Second, we test whether $\textnormal{vec}(\widehat{\mathbf W}_t(l^*+i^*-2,l^*-1+j))$ is a vector white noise starting with $j=1$. Increase $j$ by 1 and repeat the test until we reject $H_0$ at $j=j^*$. Then, we have $\widehat r_1=l^*+i^*-2$ and $\widehat r_2=l^*+j^*-2$. Finally, $\widehat \boldsymbol{\Gamma}_1=[\widehat{\mathbf A}_1,\widehat{\mathbf B}_1]$ and $\widehat\boldsymbol{\Gamma}_2=[\widehat{\mathbf P}_1,\widehat{\mathbf Q}_1]$, where $\widehat{\mathbf A}_1\in\mathbb{R}^{p_1\times \widehat r_1}$ and $\widehat{\mathbf P}_1\in\mathbb{R}^{p_2\times \widehat r_2}$.

For large $p_1$ and/or $p_2$, we use the same testing procedure, but the multivariate white noise test statistics are no longer adequate. Instead, some methods have been developed in recent years to test high-dimensional white noise series. We consider two such methods in this paper. The first method is introduced by changyaozhou2017 and makes use of the maximum absolute auto- and cross-correlations of the component series. The second method of high-dimensional white noise test is by Tsay_2018 and uses rank correlations and the extreme value theory. The test is simple and easy to use with a close-form limiting distribution under some weak assumptions. Details of the two test statistics can be found in changyaozhou2017 and Tsay_2018, respectively. See also the formulation and a brief discussion of the two test statistics $T_n$ and $T(m)$ in Section 2.3 of gaotsay2018b.

Theoretical Properties

In this section, we first present the asymptotic theory for the estimation method described in Section 2 assuming $r_1$ and $r_2$ are fixed. The consistency of the white noise test to determine $r_1$ and $r_2$ of the matrix factor is shown thereafter. The conventional asymptotic properties are established under the setting that the sample size $n$ tends to $\infty$ and everything else is fixed. Modern time series analysis encounters the situation that the number of time series $p_1p_2$ is as large as, or even larger than, the sample size $n$. We deal with these two settings separately in Sections 3.1 and 3.2 below.

Asymptotics When $n\rightarrow\infty$ But $p_1$ and $p_2$ Are Fixed

We first consider asymptotic properties under the assumption that $n\rightarrow\infty$ with $p_1$ and $p_2$ being fixed. These properties reflect the behavior of our estimation method when $n$ is large and the dimensions $p_1$ and $p_2$ are relatively small. We begin with some assumptions.

Assumption 1 is standard for dependent random processes. See gaoetal2017 for a theoretical justification for VAR models. The conditions in Assumption 2 imply that $E|y_{ij,t}|^{2\gamma}<C$ under the setting that $p_1$ and $p_2$ are fixed. To this end, we adopt the discrepancy measure used by panyao2008: for two $p\times r$ half orthogonal matrices ${\bf H}_1$ and ${\bf H}_2$ satisfying the condition ${\bf H}_1'{\bf H}_1={\bf H}_2'{\bf H}_2={\mathbf I}_{r}$, the difference between the two linear spaces $\mathcal{M}({\bf H}_1)$ and $\mathcal{M}({\bf H}_2)$ is measured by

Note that $D(\mathcal{M}({\bf H}_1),\mathcal{M}{\bf H}_2) \in [0,1].$ It is equal to $0$ if and only if $\mathcal{M}({\bf H}_1)=\mathcal{M}({\bf H}_2)$, and to $1$ if and only if $\mathcal{M}({\bf H}_1)\perp \mathcal{M}({\bf H}_2)$. The following theorem establishes the consistency of the estimated loading matrices $\widehat{\mathbf A}_1$ and $\widehat{\mathbf P}_1$, their orthonormal complements $\widehat{\mathbf B}_1$ and $\widehat{\mathbf Q}_1$, the matrices $\widehat{\mathbf B}_2$ and $\widehat{\mathbf Q}_2$, and the extracted common factor $\widehat{\mathbf A}_1\widehat{\mathbf X}_t\widehat{\mathbf P}_1'$.

From Theorem 1 and as expected, the convergence rates of all estimates are standard at $\sqrt{n}$, which is commonly seen in the traditional statistical theory. If the largest $r_1$ and $r_2$ eigenvalues of ${\mathbf M}_1$ and ${\mathbf M}_2$ are distinct, then ${\mathbf A}_1$ and ${\mathbf P}_1$ are uniquely defined up to a change of signs in columns. In fact, the consistency of the linear spaces of $\mathcal{M}({\mathbf B}_1)$ and $\mathcal{M}({\mathbf B}_2)$ is more meaningful since their columns correspond to the zero eigenvalues of ${\mathbf M}_1$ and ${\mathbf S}_1$, respectively, and they cannot be uniquely characterized.

Asymptotics When $n\rightarrow\infty$ and $p_1,p_2\rightarrow\infty$

Turn to the case of high-dimensional matrices. For vectorized variables, it is well known that if the dimension $p_1p_2$ diverges faster than $n^{1/2}$, the sample covariance matrix is no longer a consistent estimate of the population covariance matrix. When $p_1p_2=o(n^{1/2})$, it is still possible to consistently estimate the factor loading matrix and the number of common factors. See gaotsay2018a for details. Therefore, without any additional assumptions on the underlying structure of time series, $p_1p_2$ can only be as large as $o(n^{1/2})$. To deal with the case of large $p_1p_2$, we impose some conditions on the transformation matrices ${\mathbf L}$ and ${\mathbf R}$ of Equation ((ref)) and the cross dependence of ${\mathbf Y}_t$.

The quantity $\delta_1$ of Assumption 3 is used to quantify the strength of the factors. If $\delta_1=0$, the corresponding factors are called strong factors, since it includes the case where each element of $\boldsymbol{\ell}_i$ and ${\mathbf r}_{j}$ is $O(1)$. If $\delta_1>0$, the corresponding factors are weak factors and the smaller the $\delta_1$ is, the stronger the factors are. One advantage of using index $\delta_1$ is to link the convergence rates of the estimated factors explicitly to the strength of the factors. This assumption is slightly different from Condition 4 in wang2018, which actually impose two different strengths $\varsigma_1$ and $\varsigma_2$ on the front and back loading matrices, respectively. Due to the non-uniqueness of the loading matrices, we can always choose $\delta_1$ such that $(p_1p_2)^{(1-\delta_1)/2}\asymp p_1^{(1-\varsigma_1)/2}p_2^{(1-\varsigma_2)/2}$. Hence Assumption 4 ensures that all common factor components in ${\mathbf F}_t$ are of equal strength $\delta_1$. There are many sufficient conditions for Assumption 4 to hold. See the discussion of Assumption 5 in gaotsay2018b. Assumption 5(i) is mild and includes the standard normal distribution as a special case. Assumption 5(ii) is reasonable since ${\mathbf B}_2$ is a subspace of ${\mathbf B}_2^*$, $\widehat{\mathbf Q}_2$ is a subspace of $\widehat{\mathbf Q}_2^*$, and the discussion in Section 2.3 implies that that $\boldsymbol{\Xi}_1'{\mathbf B}_2^{*}{'}{\mathbf A}_1$ and $\boldsymbol{\Xi}_2'{\mathbf Q}_2^{*}{'}{\mathbf P}_1$ are invertible. The choices of $\widehat\boldsymbol{\Xi}_1$ and $\widehat\boldsymbol{\Xi}_2$, and hence $\widehat{\mathbf B}_2=\widehat{{\mathbf B}}_2^*\widehat\boldsymbol{\Xi}_1$ and $\widehat{\mathbf Q}_2=\widehat{{\mathbf Q}}_2^*\widehat\boldsymbol{\Xi}_2$ will be discussed later.

If $p_1$ and $p_2$ are large, it is not possible to consistently estimate ${\mathbf B}_2$ (also ${\mathbf Q}_2$) or even $\mathcal{M}({\mathbf B}_2)$ (also $\mathcal{M}({\mathbf Q}_2)$). Instead, we will estimate ${\mathbf B}_2^*=({\mathbf A}_{22},{\mathbf B}_2)$ or equivalently $\mathcal{M}({\mathbf B}_2^*)$, which is the subspace spanned by the eigenvectors associated with the $p_1-k_1$ smallest eigenvalues of ${\mathbf S}_1$. Assume $\widehat{{\mathbf B}}_2^*$ consists of the eigenvectors corresponding to the smallest $p-k_1$ eigenvalues of $\widehat{\mathbf S}_1$. Under some conditions, we can show that $\mathcal{M}(\widehat{{\mathbf B}}_2^*)$ is consistent to $\mathcal{M}({\mathbf B}_2^*)$. This is also the case in the literature on high-dimensional PCA with i.i.d. data. See, for example, shenetal2016 and the references therein. Therefore, the choice of $\widehat{\mathbf B}_2$ should be a subspace of $\widehat{{\mathbf B}}_2^*$, and we discuss it before Theorem 3 below.

Once we have $\widehat{{\mathbf B}}_2^*$ and $\widehat{\mathbf Q}_2^*$, we suggest to choose $\widehat{\mathbf B}_2$ and $\widehat{\mathbf Q}_2$ as $\widehat{\mathbf B}_2=\widehat{{\mathbf B}}_2^*\widehat\boldsymbol{\Xi}_1$ and $\widehat{\mathbf Q}_2=\widehat{{\mathbf Q}}_2^*\widehat\boldsymbol{\Xi}_2$, where $\widehat\boldsymbol{\Xi}_1=(\widehat\boldsymbol{\xi}_{1,1},..,\widehat\boldsymbol{\xi}_{1,r_1})\in \mathbb{R}^{(p_1-k_1)\times r_1}$ and $\widehat\boldsymbol{\Xi}_2=(\widehat\boldsymbol{\xi}_{2,1},..,\widehat\boldsymbol{\xi}_{2,r_2})\in \mathbb{R}^{(p_2-k_2)\times r_2}$, where $\widehat\boldsymbol{\xi}_{1,i}$ is the vector associated with the $i$-th largest eigenvalues of $\widehat{{\mathbf B}}_2^*{'}\widehat{\mathbf A}_1\widehat{\mathbf A}_1'\widehat{{\mathbf B}}_2^*$ and $\widehat\boldsymbol{\xi}_{2,j}$ is the vector associated with the $j$-th largest eigenvalues of $\widehat{{\mathbf Q}}_2^*{'}\widehat{\mathbf P}_1\widehat{\mathbf P}_1'\widehat{{\mathbf Q}}_2^*$. These choices can guarantee that the matrices $(\widehat{\mathbf B}_2'\widehat{\mathbf A}_1)^{-1}$ and $(\widehat{\mathbf Q}_2'\widehat{\mathbf P}_1)^{-1}$ behave well when recovering the factor $\widehat{\mathbf X}_t$. On the other hand, they could still eliminate the diverging part of the noise covariance matrix and give prominent convergence rate, as shown in Theorem 3. There are many ways to choose the numbers of components $k_1$ and $k_2$ in Assumption 4 so long as $p_1-k_1>r_1$ and $p_2-k_2>r_2$. We discuss the choices of $k_1$ and $k_2$ in Remark 2 below. The following theorem states the convergence rate of the extracted common factors.

Next, we study the consistency of the white noise tests described in Section 2. In fact, the consistency conditions depend on which test statistic we use. We only consider the two test statistics $T_n$ and $T(m)$ discussed in Section 2.3 of gaotsay2018b and present the consistency when $p_1$ and $p_2$ are large since the case of small $p_1$ and $p_2$ is trivial. For any random vector ${\mathbf x}_t$ to be sub-Gaussian we mean there exists a constant $C>0$ such that $P(|{\mathbf v}'({\mathbf x}_t-E{\mathbf x}_t)|>x)\leq C\exp(-Cx^2)$ for any constant vector $\|{\mathbf v}\|_2=1$. We need an additional assumption.

With the estimator $\widehat r_1$, we may define the estimator for ${\mathbf A}_1$ as $\widehat{\mathbf A}_1=(\widehat{\mathbf a}_1,...,\widehat{\mathbf a}_{\widehat r_1})$, where $\widehat{\mathbf a}_1,...,\widehat{\mathbf a}_{\widehat r_1}$ are the orthonormal eigenvectors of $\widehat {\mathbf M}_1$, defined in ((ref)), corresponding to the $\widehat r_1$ largest eigenvalues. In addition, we may also replace $r_1$ by $\widehat r_1$ in the whole methodology described in Section 2. We can define $\widehat{\mathbf P}_1$ in a similar way.

Numerical Properties

Simulation

In this section, we illustrate the finite-sample properties of the proposed methodology under different choices of $p_1$ and $p_2$. Because the actual dimension is $p_1p_2$ which can easily go to hundreds for even relatively small $p_1$ and $p_2$, we focus on the case of high dimension, which is of more interest. As the dimensions of $\widehat{\mathbf A}_1$ and ${\mathbf A}_1$ are not necessarily the same, and ${\mathbf L}_1$ is not an orthogonal matrix in general, we first extend the discrepancy measure in Equation ((ref)) to a more general form below. Let ${\mathbf H}_i$ be a $p\times h_i$ matrix with rank$({\mathbf H}_i) = h_i$, and ${\mathbf P}_i = {\mathbf H}_i({\mathbf H}_i'{\mathbf H}_i)^{-1} {\mathbf H}_i'$, $i=1,2$. Define

Then $\bar{D} \in [0,1]$. Furthermore, $\bar{D}(\mathcal{M}({\mathbf H}_1),\mathcal{M}({\mathbf H}_2))=0$ if and only if either $\mathcal{M}({\mathbf H}_1)\subset \mathcal{M}({\mathbf H}_2)$ or $\mathcal{M}({\mathbf H}_2)\subset \mathcal{M}({\mathbf H}_1)$, and it is 1 if and only if $\mathcal{M}({\mathbf H}_1) \perp \mathcal{M}({\mathbf H}_2)$. When $h_1 = h_2=h$ and ${\mathbf H}_i'{\mathbf H}_i= {\mathbf I}_r$, $\bar{D}(\mathcal{M}({\mathbf H}_1),\mathcal{M}({\mathbf H}_2)) $ reduces to that in Equation ((ref)). We only present the simulation results for $k_0=2$ in Equation ((ref)) to save space since other choices of $k_0$ produce similar patterns.\\

{ \bf Example 1.} Consider model ((ref)) with common factors satisfying \[{\mathbf F}_t=\boldsymbol{\Phi}{\mathbf F}_{t-1}\boldsymbol{\Psi}'+{\mathbf N}_t,\] where ${\mathbf N}_t$ is a matrix-variate white noise process with independent entries, $\boldsymbol{\Phi}\in \mathbb{R}^{r_1\times r_1}$ and $\boldsymbol{\Psi}\in \mathbb{R}^{r_2\times r_2}$ are two diagonal coefficient matrices. We set the true dimension of the matrix factors $(r_1,r_2)=(2,3)$, the orders of the diverging noise components $(k_1,k_2)=(1,2)$ as defined in Assumption 4, the dimensions $(p_1,p_2)=(7,7)$, $(10,15)$, $(20,20)$ and $(20,30)$, and the sample sizes are $n=300$, $500$, $1000$, $1500$, $3000$. We consider three scenarios for $\delta_1$ and $\delta_2$: $(\delta_1,\delta_2)=(0,0.9)$, $(0.2,0.8)$ and $(0.5,0.5)$. We can also obtain similar results for other settings but omit the details to save space. For each scenario mentioned above, the elements of ${\mathbf L}$ and ${\mathbf R}$ are drawn independently from $U(-2,2)$, and then we divide ${\mathbf L}_1$ (also ${\mathbf R}_1$) by $p_1^{\delta_1/2}$ (also $p_2^{\delta_1/2}$), the first $k_1$ (also $k_2$) columns of ${\mathbf L}_2$ (also ${\mathbf R}_2$) by $p_1^{\delta_2/2}$ (also $p_2^{\delta_2/2}$) and the rest $v_1-k_1$ (also $v_2-k_2$) columns by $p_1$ (also $p_2$) to satisfy Assumptions 3 and 4. $\boldsymbol{\Phi}$ and $\boldsymbol{\Psi}$ are diagonal matrices with their diagonal elements drawn independently from $U(0.5,0.9)$, $\textnormal{vec}({\mathbf Z}_{12,t})\sim N(0,{\mathbf I}_{r_1v_2})$, $\textnormal{vec}({\mathbf Z}_{21,t})\sim N(0,{\mathbf I}_{v_1r_2})$, $\textnormal{vec}({\mathbf Z}_{22,t})\sim N(0,{\mathbf I}_{v_1v_2})$, $\textnormal{vec}({\mathbf N}_t)\sim N(0,{\mathbf I}_{r_1r_2})$. We use $500$ replications in each experiment.

We first study the performance of estimating the dimension of the matrix-variate factors. For simplicity, we only report the results of the test statistic $T(m)$ with $m=10$ defined in gaotsay2018b, and the results for the other test are similar. When $p_1p_2>n$, we only keep the upper ${{\varepsilon}}\sqrt{n}$ row- and column-transformed series of $\widehat\boldsymbol{\Gamma}_1'{\mathbf Y}_t\widehat\boldsymbol{\Gamma}_2$ with ${{\varepsilon}}=0.9$ in the testing. Similar results can be obtained for other choices of ${{\varepsilon}}$ and we do not report them here. The testing results are given in Table (ref). From the table, we see that for each setting of $(\delta_1,\delta_2)$ and fixed $(p_1,p_2)$, the performance of the white noise test improves as the sample size increases. The performance is also quite satisfactory for moderately large $p_1p_2$ when the factor strength is stronger than that of the noises. When $(\delta_1,\delta_2)=(0.5,0.5)$, we see that the test does not perform well for small sample sizes, which is understandable since the factors and the noises have the same level of strength but the diverging noise effect is much more prominent by Equation ((ref)), yet the performance improves significantly when the sample size increases.

Next, we study the accuracy of the estimated loading matrices. The boxplots of $\bar{D}(\mathcal{M}(\widehat{\mathbf A}_1),\mathcal{M}({\mathbf L}_1))$ and $\bar{D}(\mathcal{M}(\widehat{\mathbf P}_1),\mathcal{M}({\mathbf R}_1))$ are shown in Figure (ref)(a) and (b), respectively. From Figure (ref), we see that the estimation accuracy of the loading matrix improves as the sample size increases even for moderately large $p_1p_2$, which is in line with our asymptotic theory. Furthermore, we study the estimation accuracy of the estimated factor process by

The results are shown in Figure (ref), from which we see that, for fixed $(p_1,p_2)$, the estimation accuracy also improves as the sample size increases. This result is consistent with our Theorem 3 of Section 3.

To see the advantages of the proposed method, we compare it with that of wang2018 (denoted by WLC) in selecting the order of the matrix-variate factors. For the ratio-based method in WLC, let $\widehat\lambda_{i,1},...,\widehat\lambda_{i,p_i}$ be the $p_i$ eigenvalues of $\widehat{\mathbf M}_i$ for $i=1,2$, define

Figure (ref)(a)-(b) present the boxplots of $\widehat r_1$ and $\widehat r_2$, respectively. We see from Figure (ref) that the estimated number of factors $\widehat r_i$ tend to be the sum of the number of common factors $r_i$ and the number of spiked components of the noises $k_i$ in most of the scenarios. The result indicates that the ratio-based method of wang2018 may fail to identify the correct dimension of the matrix-variate factor process with dynamic dependence if the covariance of the noise has diverging eigenvalues, while the proposed white noise test continues to work well, as shown in Table (ref).

Finally, we compare our method with the one of wang2018 in recovering the common factors since a key difference between the two methods is that we allow some of the eigenvalues of the noise covariance to diverge. We denote our method by GT and the results are reported in Table (ref) for $(r_1,r_2)=(2,3)$, $(k_1,k_2)=(1,2)$, and $(\delta_1,\delta_2)=(0.5,0.5)$. From the table, we see that, because the ratio-based method tends to overestimate the dimension of the common factors, the estimation error of our method is much smaller than that obtained by WLC. In addition, for a given $(p_1,p_2)$, the estimation error by our method tends to decrease as the sample size increases, which is in agreement with our asymptotic theory. Overall, under the assumption that the noise effect is prominent, the proposed method outperforms the existing one in the literature.

Real Data Analysis

{\bf Example 2.} In this example, we use the Fama-French return series to illustrate application of the proposed method. The data contain monthly returns of 100 portfolios, structured in a $10$ by $10$ matrix according to ten levels of market capitalization (Size, in rows from small to large) and ten levels of investment (Inv, in columns from low to high) both of which are factors for average stock returns considered in famafrench2015. The return series spans from July 1963 to December 2019 and consists of 678 monthly observations for each individual process. Therefore, the series forms a $10\times 10\times 678$ tensor-valued data set. The data and relevant information are available at

\url{http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html}.

Following sharpe1964 and famafrench2015, we adjust each of the return series by subtracting the corresponding risk-free asset returns, which are also available from the above website. The missing values were imputed by a simple exponential smoothing method. Time plots of the adjusted $10\times 10$ series are shown in Figure (ref) with $p_1=p_2=10$ and $n=678$.

We first applied the method of Section 2.4 with $k_0=2$ and found that the test statistic $T(m)$ with $m=10$ for testing the null hypothesis $H_0(2)$ defined in Section 2.5 is 4.85, which exceeds the critical value 4.81 based on the limiting Gumbel distribution in Tsay_2018 with $\alpha=0.05$. But the null hypothesis is not rejected if we increase the order in either the column or row direction. Therefore, $\widehat r_1=2$ and $\widehat r_2=2$ implying that a $2\times 2$ matrix-variate latent factor process is detected. The estimated front and back loading matrices after being multiplied by $30$ are reported in Table (ref), which has several implications. First, for Size, it seems that the 10 rows of the portfolios can be divided into two or three groups. The one with the smallest size (corresponding to S1) depends on both the first and the second factors heavier than the others, the second smallest size portfolio depends more on the first row factor and less on the second one, and the 3rd to the 10th size portfolios have similar dependence on both the first and the second rows of the matrix-variate factors. Second, for Investment, all the portfolios have similar dependence on the first column of the factor matrix, and the dependence on the second columns seems to have three groups; the lowest investment portfolio (corresponding to Inv1) seems to depend heavily on the second row of the factors, the 5th to the 8th and the 10-th investment portfolios have similar dependence on the first and the second rows of the factors, whereas the 2nd to the 4th and the 9th investment portfolios depend more on the first column of the factors. In addition, the signs of the first coefficients of the size loading and the investment loading are the same, which implies that each return series have a co-movement with respect to the $[1,1]$-factor series. This is understandable since we can treat this common factor as representing the market factor of the capital asset pricing model (CAPM) of sharpe1964. The product of the first coefficients of the size loading and the investment loading can be treated as a market beta, even though the starting point of our approach is different from that of CAPM. The usefulness of the detected market factor, however, deserves a further investigation.

To obtain the extracted factors, by the two-way projected PCA of Sections 2.3 and 2.4, we first examine the eigenvalues of the sample covariance matrices $\widehat {\mathbf S}_1$ and $\widehat{\mathbf S}_2$. From Figure (ref), we see that the first eigenvalues of $\widehat{\mathbf S}_1$ and $\widehat{\mathbf S}_2$ are much larger than the others. Therefore, we choose $\widehat k_1=\widehat k_2=1$, and the recovered matrix-variate factors are shown in Figure (ref)(a) as well as their corresponding spectrum in Figure (ref)(b). From Figure (ref), we see that there are three series with non-trivial spectra, which are all dynamically dependent and they capture most of the dynamic information of the data, and the $[1,2]$-factor appears to be not serially correlated by itself because its sample spectrum is flat. However, this does not imply that the $[1,2]$-factor captures no dynamic information in the detected matrix-variate common factors. For example, the lag-1 cross-correlation between the $[1,2]$-factor and the $[2,2]$-factor is 0.08. If we test for the zero lag-1 corss correlation between these two series using the long-run covariance matrix calculated by the method in andrews1991, the $p$-value is 0.038 implying that the two factors are lag-1 cross-correlated. Therefore, the detected 2-by-2 matrix-variate common factor process does not violate the assumptions of the proposed model.

Next we examine and compare the forecasting performance of the extracted factors via the proposed method (denoted by GT) and those by wang2018 (denoted by WLC). We estimate the models using the data in the time span $[1,\tau]$ with $\tau=558,...,678-h$ for the $h$-step ahead forecasts, i.e., we use returns of the last ten years for out-of-sample forecasting. For the method of wang2018, the estimated dimension of the matrix-variate factor is $(\widehat r_1,\widehat r_2)=(1,1)$. For simplicity, we employ a simple AR(1) model for each detected common factor to produce forecasts. We also fit a scalar AR(1) (denoted by SAR) model to each individual return series as a benchmark approach in out-of-sample forecasting. The following two criteria are used to measure the forecast errors:

and

where $p_1=p_2=10$. Table (ref) reports the 1-step to 4-step ahead forecast errors of Equations ((ref)) and ((ref)) for the methods GT, WLC, and SAR. The smallest forecast error of each step is shown in boldface. From the table, we see that our proposed method is capable of producing accurate forecasts and the associated forecast errors based on the extracted factors by our method are smaller than that based on the factor extracted by WLC or the benchmark approach SAR. Although the difference in forecasting errors between the three methods used in Table (ref) is small, it is generally not easy to produce accurate forecasts in asset returns and the improvements by our proposed method could have substantial implications to practitioners, especially over the ten-year horizon.

In conclusion, for the monthly excess return series considered, our method not only produces interpretable factors, but also improves the out-of-sample forecasting. We like to emphasize that the proposed method is different from the traditional factor model analysis, especially those based on the conventional principal component analysis. The proposed model explores a different aspect of the data via a two-way transformation. Finally, the out-of-sample forecasting can be improved if we adopt some regularization method, but we do not pursue it here. The proposed method is intended as another tool for modeling high-dimensional and possibly highly-correlated matrix-variate time series.

Concluding Remarks

This paper proposed a new approach to analyze high-dimensional, dynamically dependent matrix-variate data in the presence of prominent noise effect. The proposed approach is an extension of that for vector time series in TiaoTsay_1989 and gaotsay2018b. The approach not only can reduce the dimensionality of the matrix-variate data, but also preserves the structure of the matrix to mitigate loss in information. The proposed approach is easy to implement for high-dimensional matrix-variate time series data and empirical results show that it can effectively extract the number of common factors from complex data. In addition, the extracted common factors could be useful in out of sample predictions.