The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
83,836 characters
A Two-Way Transformed Factor Model for Matrix-Variate Time Series
\title{\bf A Two-Way Transformed Factor Model for Matrix-Variate Time Series}
\author{
Zhaoxing Gao$^1$ and Ruey S. Tsay$^2$ \\
$^1$Department of Mathematics, Lehigh University\\
$^2$Booth School of Business, University of Chicago
}
\date{}
\maketitle
\begin{abstract}
We propose a new framework for modeling high-dimensional matrix-variate time series by a
two-way transformation, where the transformed data consist of a matrix-variate factor process, which is dynamically dependent, and three other blocks of white noises. Specifically, for a given $p_1\times p_2$ matrix-variate time series, we seek common nonsingular transformations to project the rows and columns onto another $p_1$ and $p_2$ directions according to the strength of the dynamic dependence of the series on the past values. Consequently, we treat the data as nonsingular linear row and column transformations of dynamically dependent common factors and white noise idiosyncratic components. We propose a common orthonormal projection method to estimate the front and back loading matrices of the matrix-variate factors. Under the setting that the largest eigenvalues of the covariance of the vectorized idiosyncratic term diverge for large $p_1$ and $p_2$,
we introduce a two-way projected Principal Component Analysis (PCA) to estimate the associated loading matrices of the idiosyncratic terms to mitigate such diverging noise effects. A diagonal-path white noise testing procedure is proposed to estimate the order of the factor matrix.
Asymptotic properties of the proposed method are established for both fixed and diverging dimensions as the sample size increases to infinity. We use simulated and real examples to assess the performance of the proposed method. We also compare our method with some existing ones in the literature and find that the proposed approach not only provides interpretable results but also performs well in out-of-sample forecasting.
\end{abstract}
\noindent {\sl Keywords}: Structured factor, Eigen-analysis, Projected PCA, Kronecker product, Diverging eigenvalues, High-dimensional white noise test.
\newpage
\section{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, \cite{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., \cite{rogers2013} and \cite{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 \cite{Tsay_2014} and the references therein discussed two different canonical structures. \cite{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 \cite{ShojaieMichailidis_2010}, \cite{SongBickel_2011}, and \cite{HanTsay2020}, among others.
For dimension reduction, popular methods include the canonical correlation analysis (CCA) of \cite{BoxTiao_1977}, the principle component analysis (PCA) of \cite{StockWatson_2002}, the scalar component analysis of \cite{TiaoTsay_1989}. The factor model approach can be found in \cite{BaiNg_Econometrica_2002}, \cite{StockWatson_2005}, \cite{forni2000,forni2005}, \cite{panyao2008}, \cite{LamYaoBathia_Biometrika_2011}, \cite{lamyao2012}, \cite{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; \cite{walden2002} handled this type of data in signal and image processing, \cite{wang2018} proposed a factor model for matrix-variate time series, which maintains and utilizes the matrix structure to achieve the dimension reduction, and \cite{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, \cite{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
\begin{equation}\label{int:m}
{\mathbf T}_1{\mathbf Y}_t{\mathbf T}_2'=\left[\begin{array}{cc}
{\mathbf F}_t&{\mathbf Z}_{12,t}\\
{\mathbf Z}_{21,t}&{\mathbf Z}_{22,t}
\end{array}\right],
\end{equation}
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{int:m}) 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{int:m}), 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:
\begin{equation}\label{vec:m}
{\mathbf y}_t:=\textnormal{vec}({\mathbf Y}_t)={\mathbf A}\left[\begin{array}{c}
{\mathbf f}_t\\
{\mathbf z}_t
\end{array}\right],
\end{equation}
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 \cite{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 \cite{gaotsay2018b} and \cite{TiaoTsay_1989}. The structure of ${\mathbf A}$ is different from that in \cite{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 \cite{wang2018}, which essentially follows the method in \cite{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 \cite{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)$.
\section{Models and Methodology}
\subsection{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:
\begin{equation}\label{m-factor}
{\mathbf Y}_t={\mathbf L}\left[\begin{array}{cc}
{\mathbf F}_t&{\mathbf Z}_{12,t}\\
{\mathbf Z}_{21,t}&{\mathbf Z}_{22,t}
\end{array}\right]{\mathbf R}'={\mathbf L}_1{\mathbf F}_t{\mathbf R}_1'+{\mathbf L}_2{\mathbf Z}_{21,t}{\mathbf R}_1'+{\mathbf L}_1{\mathbf Z}_{12,t}{\mathbf R}_2'+{\mathbf L}_2{\mathbf Z}_{22,t}{\mathbf R}_2',
\end{equation}
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{m-factor}) 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 \cite{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{m-factor}).
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{m-factor}) can be rewritten as
\begin{equation}\label{mf}
{\mathbf Y}_t={\mathbf A}_1{\mathbf X}_t{\mathbf P}_1'+{\mathbf A}_2{\mathbf E}_{21,t}{\mathbf P}_1'+{\mathbf A}_1{\mathbf E}_{12,t}{\mathbf P}_2'+{\mathbf A}_2{\mathbf E}_{22,t}{\mathbf P}_2'.
\end{equation}
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{mf}) 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)$.
\subsection{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:
\begin{equation}\label{opm}
\max_{{\mathbf a}\in\mathbb{R}^{p_1}}\sum_{i=1}^{p_2}\|\mbox{Cov}({\mathbf a}'{\mathbf y}_{i,t},\boldsymbol{\eta}_t)\|_2^2,\quad\text{subject to}\quad {\mathbf a}'{\mathbf a}=1.
\end{equation}
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
\begin{equation}\label{m1}
{\mathbf M}_1=\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)'.
\end{equation}
On the other hand, under model (\ref{mf}), 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
\begin{equation}\label{sigma:k}
\boldsymbol{\Sigma}_{y,ij}(k)={\mathbf A}_1\boldsymbol{\Sigma}_{xp,ij}(k){\mathbf A}_1',
\end{equation}
where we assume ${\mathbf f}_t$ and ${\mathbf z}_s$ are uncorrelated for any $t$ and $s$. Therefore,
\begin{equation}\label{M1:x}
{\mathbf M}_1={\mathbf A}_1\left\{\sum_{k=1}^{k_0}\sum_{i=1}^{p_2}\sum_{j=1}^{p_2}[\boldsymbol{\Sigma}_{xp,ij}(k)\boldsymbol{\Sigma}_{xp,ij}(k)']\right\}{\mathbf A}_1'.
\end{equation}
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{opm}) 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{m1}) 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.
\subsection{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{mf}) that
\begin{equation}\label{byq}
{\mathbf B}_1'{\mathbf Y}_t{\mathbf Q}_1={\mathbf B}_1'{\mathbf A}_2{\mathbf E}_{22,t}{\mathbf P}_2'{\mathbf Q}_1,
\end{equation}
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,
\begin{equation}\label{b2q2}
{\mathbf B}_2'{\mathbf Y}_t={\mathbf B}_2'{\mathbf A}_1{\mathbf X}_t{\mathbf P}_1'+{\mathbf B}_2{\mathbf A}_1{\mathbf E}_{12,t}{\mathbf P}_2'\,\,\text{and}\,\, {\mathbf Y}_t{\mathbf Q}_2={\mathbf A}_1{\mathbf X}_t{\mathbf P}_1'{\mathbf Q}_2+{\mathbf A}_2{\mathbf E}_{21,t}{\mathbf P}_1'{\mathbf Q}_2.
\end{equation}
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{byq}). 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{mf}) and (\ref{byq}) that
\begin{equation}\label{omgy}
\mbox{Cov}({\mathbf y}_{i,t},\textnormal{vec}({\mathbf B}_1'{\mathbf Y}_t{\mathbf Q}_1))=\boldsymbol{\Omega}_{y_i}({\mathbf Q}_1\otimes{\mathbf B}_1)={\mathbf A}_2\boldsymbol{\Omega}_{e_{22,i}p}({\mathbf P}_2'{\mathbf Q}_1\otimes {\mathbf A}_2'{\mathbf B}_1).
\end{equation}
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
\begin{equation}\label{s1}
{\mathbf S}_1:=\sum_{i=1}^{p_2}[\boldsymbol{\Omega}_{y_i}({\mathbf Q}_1\otimes{\mathbf B}_1)][\boldsymbol{\Omega}_{y_i}({\mathbf Q}_1\otimes{\mathbf B}_1)]',
\end{equation}
from which we can see, via (\ref{omgy}), 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{mf}) that
\begin{equation}\label{r:ft}
{\mathbf B}_2'{\mathbf Y}_t{\mathbf Q}_2={\mathbf B}_2'{\mathbf A}_1{\mathbf X}_t{\mathbf P}_1'{\mathbf Q}_2,
\end{equation}
and consequently,
\begin{equation}\label{ft:e}
{\mathbf X}_t=({\mathbf B}_2'{\mathbf A}_1)^{-1}{\mathbf B}_2'{\mathbf Y}_t{\mathbf Q}_2({\mathbf P}_1'{\mathbf Q}_2)^{-1},
\end{equation}
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.
\subsection{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{m1}) as follows:
\begin{equation}\label{m1hat}
\widehat{\mathbf M}_1=\sum_{k=1}^{k_0}\sum_{i=1}^{p_2}\sum_{j=1}^{p_2}\widehat\boldsymbol{\Sigma}_{y,ij}(k)\widehat\boldsymbol{\Sigma}_{y,ij}(k)',
\end{equation}
where
\begin{equation}\label{yij:hat}
\widehat\boldsymbol{\Sigma}_{y,ij}(k)=\frac{1}{n}\sum_{t=k+1}^{n}({\mathbf y}_{i,t}-\bar{{\mathbf y}}_i)({\mathbf y}_{j,t-k},-\bar{{\mathbf y}}_j)',
\end{equation}
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
\begin{equation}\label{s1:hat}
\widehat{\mathbf S}_1=\sum_{i=1}^{p_2}[\widehat\boldsymbol{\Omega}_{y_i}(\widehat{\mathbf Q}_1\otimes\widehat{\mathbf B}_1)][\widehat\boldsymbol{\Omega}_{y_i}(\widehat{\mathbf Q}_1\otimes\widehat{\mathbf B}_1)]',
\end{equation}
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
\begin{equation}\label{rec:ft}
\widehat{\mathbf X}_t=(\widehat{\mathbf B}_2'\widehat{\mathbf A}_1)^{-1}\widehat{\mathbf B}_2'{\mathbf Y}_t\widehat{\mathbf Q}_2(\widehat{\mathbf P}_1'\widehat{\mathbf Q}_2)^{-1}.
\end{equation}
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 \cite{chenxiaoyang2020}.
\subsection{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 \cite{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 \cite{onatski2010}, the ratio-based method in Lam and Yao (2012), the canonical correlation analysis in \cite{gaotsay2018a}, and the white noise testing approach in \cite{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 \cite{gaotsay2018b}. The idea of our method follows from equation (\ref{byq}) 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{m1hat}) 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, \cite{hosking1980} and \cite{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 \cite{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 \cite{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 \cite{changyaozhou2017} and \cite{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 \cite{gaotsay2018b}.
\section{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.
\subsection{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.
\begin{assumption}
The process $\{\textnormal{vec}({\mathbf Y}_t),\textnormal{vec}({\mathbf F}_t)\}$ is $\alpha$-mixing with the mixing coefficient satisfying the condition $\sum_{k=1}^\infty\alpha_p(k)^{1-2/\gamma}<\infty$ for some $\gamma>2$, where
\[\alpha_p(k)=\sup_{i}\sup_{A\in\mathcal{F}_{-\infty}^i,B\in \mathcal{F}_{i+k}^\infty}|P(A\cap B)-P(A)P(B)|,
\]
and $\mathcal{F}_i^j$ is the $\sigma$-field generated by $\{(\textnormal{vec}({\mathbf Y}_t),\textnormal{vec}({\mathbf F}_t)):i\leq t\leq j\}$.
\end{assumption}
\begin{assumption}
For any $i=1,...,r_1 r_2$ and $1\leq j\leq p_1p_2-r_1 r_2$, $E|f_{i,t}|^{2\gamma}<C_1$ and $E|z_{j,t}|^{2\gamma}<C_2$, where $f_{i,t}$ and $z_{j,t}$ are the $i$-th and $j$-th element of ${\mathbf f}_t$ and ${\mathbf z}_t$, respectively, $C_1$ and $C_2>0$ are constants, and $\gamma$ is given in Assumption 1.
\end{assumption}
Assumption 1 is standard for dependent random processes. See \cite{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
\cite{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
\begin{equation}
D(\mathcal{M}({\bf H}_1),\mathcal{M}{\bf
H}_2)=\sqrt{1-\frac{1}{r}\textrm{tr}({\bf H}_1{\bf H}_1'{\bf
H}_2{\bf H}_2')}.\label{eq:D}
\end{equation}
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'$.
\begin{theorem}
Suppose Assumptions 1-2 hold and $(r_1,r_2)$ are known and fixed. Then, for fixed $p_1$ and $p_2$,
\[D(\mathcal{M}(\widehat{\mathbf A}_1),\mathcal{M}({\mathbf A}_1))=O_p(n^{-1/2}),\quad D(\mathcal{M}(\widehat{\mathbf B}_1),\mathcal{M}({\mathbf B}_1))=O_p(n^{-1/2}),\]
\[D(\mathcal{M}(\widehat{\mathbf P}_1),\mathcal{M}({\mathbf P}_1))=O_p(n^{-1/2}),\quad D(\mathcal{M}(\widehat{\mathbf Q}_1),\mathcal{M}({\mathbf Q}_1))=O_p(n^{-1/2}) \]
and
\[\quad D(\mathcal{M}(\widehat{\mathbf B}_2),\mathcal{M}({\mathbf B}_2))=O_p(n^{-1/2}),\quad D(\mathcal{M}(\widehat{\mathbf Q}_2),\mathcal{M}({\mathbf Q}_2))=O_p(n^{-1/2}), \]
as $n\rightarrow\infty$.
Furthermore,
\[\|\widehat{\mathbf A}_1\widehat{\mathbf X}_t\widehat{\mathbf P}_1'-{\mathbf A}_1{\mathbf X}_t{\mathbf P}_1'\|_2=O_p(n^{-1/2}).\]
\end{theorem}
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.
\subsection{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 \cite{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{m-factor}) and the cross dependence of ${\mathbf Y}_t$.
\begin{assumption}
(i) ${\mathbf L}_1=(\boldsymbol{\ell}_1,...,\boldsymbol{\ell}_{r_1})$ and ${\mathbf R}_1=({\mathbf r}_1,...,{\mathbf r}_{r_2})$ such that $\|\boldsymbol\ell_i\|_2^2\asymp p_1^{1-\delta_1}$, $\|{\mathbf r}_j\|_2^2\asymp p_2^{1-\delta_1}$, $i=1,...,r_1$, $j=1,...,r_2$, and $\delta_1\in[0,1)$; (ii) For each $i=1,...,r_1$, $j=1,...,r_2$ and $\delta_1$ given in (i), $\min_{\theta_k\in \mathbb{R},k\neq i}\|\boldsymbol\ell_i-\sum_{1\leq k\leq r_1,k\neq i}\theta_k\boldsymbol\ell_k\|_2^2\asymp p_1^{1-\delta_1}$ and $\min_{\theta_k\in \mathbb{R},k\neq j}\|{\mathbf r}_j-\sum_{1\leq k\leq r_2,k\neq j}\theta_k{\mathbf r}_k\|_2^2\asymp p_2^{1-\delta_1}$ .
\end{assumption}
\begin{assumption}
(i) ${\mathbf L}_2$ and ${\mathbf R}_2$ admit a singular value decomposition ${\mathbf L}_2={\mathbf A}_2{\mathbf D}_2{\mathbf U}_2'$ and ${\mathbf R}_2={\mathbf P}_2\boldsymbol{\Lambda}_2{\mathbf V}_2'$, where ${\mathbf A}_2\in \mathbb{R}^{p_1\times v_1}$ and ${\mathbf P}_2\in \mathbb{R}^{p_2\times v_2}$ are given in Equation (\ref{mf}), ${\mathbf D}_2=\mbox{diag}(d_1,...,d_{v_1})$ and ${\mathbf U}_2\in \mathbb{R}^{v_1\times v_1}$ satisfying ${\mathbf U}_2'{\mathbf U}_2={\mathbf I}_{v_1}$, $\boldsymbol{\Lambda}_2=\mbox{diag}(\gamma_1,...,\gamma_{v_2})$, ${\mathbf V}_2\in \mathbb{R}^{v_2\times v_2}$ satisfying ${\mathbf V}_2'{\mathbf V}_2={\mathbf I}_{v_2}$; (ii) There exist finite integers $1\leq k_1< v_1$ and $1\leq k_2<v_2$ such that $d_1\asymp...\asymp d_{k_1}\asymp p_1^{(1-\delta_2)/2}$ and $\gamma_1\asymp...\asymp\gamma_{k_2}\asymp p_2^{(1-\delta_2)/2}$ for some $\delta_2\in[0,1)$ and $d_{k_1+1}\asymp...\asymp d_{v_1}\asymp 1\asymp\gamma_{k_2+1}\asymp...\asymp\gamma_{v_2}$.
\end{assumption}
\begin{assumption}
(i) For any $1\leq l_1\leq v_1$, $1\leq l_2\leq v_2$, ${\mathbf h}\in \mathbb{R}^{l_1\times l_2}$, ${\mathbf U}\in\mathbb{R}^{v_1\times l_1 }$ and ${\mathbf V}\in \mathbb{R}^{v_2\times l_2}$ with $\|{\mathbf h}\|_2=c<\infty$, ${\mathbf U}'{\mathbf U}={\mathbf I}_{l_1}$ and ${\mathbf V}'{\mathbf V}={\mathbf I}_{l_2}$, we assume $E|{\mathbf h}'\textnormal{vec}({\mathbf U}'{\mathbf Z}_{22,t}{\mathbf V})|^{2\gamma}<\infty$; (ii) $\sigma_{\min}(\boldsymbol{\Xi}_1'{\mathbf B}_2^{*}{'}{\mathbf A}_1)\geq C_3$ and $\sigma_{\min}(\boldsymbol{\Xi}_2'{\mathbf Q}_2^{*}{'}{\mathbf P}_1)\geq C_4$ for some constants $C_3, C_4>0$ and some half orthogonal matrices $\boldsymbol{\Xi}_1\in \mathbb{R}^{(p_1-v_1)\times r_1}$ and $\boldsymbol{\Xi}_2\in \mathbb{R}^{(p_2-v_2)\times r_2}$ satisfying $\boldsymbol{\Xi}_1'\boldsymbol{\Xi}_1={\mathbf I}_{r_1}$ and $\boldsymbol{\Xi}_2'\boldsymbol{\Xi}_2={\mathbf I}_{r_2}$, where $\sigma_{\min}$ denotes the minimum non-zero singular value of a matrix.
\end{assumption}
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 \cite{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 \cite{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, \cite{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.
\begin{theorem}
Suppose Assumptions 1-5 hold and $r_1$ and $r_2$ are known and fixed. As $n\rightarrow\infty$, if $p_1^{\delta_1}p_2^{\delta_1}n^{-1/2}=o(1)$, then
\[\|D(\mathcal{M}({\widehat{\mathbf A}}_1),\mathcal{M}{{\mathbf A}}_1)\|_2=O_p(p_1^{\delta_1}p_2^{\delta_1}n^{-1/2})\,\,\text{and}\,\,\|D(\mathcal{M}({\widehat{\mathbf P}}_1),\mathcal{M}{{\mathbf P}}_1)\|_2=O_p(p_1^{\delta_1}p_2^{\delta_1}n^{-1/2}),\]
and the above results also hold for $\|D(\mathcal{M}({\widehat{\mathbf B}}_1),\mathcal{M}{{\mathbf B}}_1)\|_2$ and $\|D(\mathcal{M}({\widehat{\mathbf Q}}_1),\mathcal{M}{{\mathbf Q}}_1)\|_2$.
Furthermore,
\begin{equation*}
\|D(\mathcal{M}(\widehat{\mathbf B}_2^*),\mathcal{M}({\mathbf B}_2^*))\|_2
=O_p(p_1^{\delta_2}p_2^{3\delta_2/2}n^{-1/2}+p_1^{\delta_1}p_2^{\delta_1+\delta_2}n^{-1/2}),
\end{equation*}
and
\[\|D(\mathcal{M}(\widehat{\mathbf Q}_2^*),\mathcal{M}({\mathbf Q}_2^*))\|_2=O_p(p_1^{3\delta_2/2}p_2^{\delta_2}n^{-1/2}+p_1^{\delta_1+\delta_2}p_2^{\delta_1}n^{-1/2}).\]
\end{theorem}
\begin{remark}
(i) For the consistencies of $\widehat{\mathbf A}_1$ and $\widehat{\mathbf P}_1$, we require $p_1p_2=o(n^{\frac{1}{2\delta_1}})$. When the strength $\delta_1\in[0,1/2]$, the range of the total dimensions $p_1p_2$ can be greater than $\sqrt{n}$. \\
(ii) The conditions for the consistencies of $\widehat{\mathbf B}_2^*$ and $\widehat{\mathbf Q}_2^*$ are slightly stronger since they depend on the estimation error in the first step. Specifically, we require $p_1^{\delta_2}p_2^{3\delta_2/2}n^{-1/2}=o(1)$, $p_1^{\delta_1}p_2^{\delta_1+\delta_2}n^{-1/2}=o(1)$, $p_1^{3\delta_2/2}p_2^{\delta_2}n^{-1/2}=o(1)$ and $p_1^{\delta_1+\delta_2}p_2^{\delta_1}n^{-1/2}=o(1)$. To give a better illustration, we assume $p_1\asymp p_2\asymp p$, then we have $p^2=o(n^{\frac{1}{2\delta_1}})$ for the consistency of $\widehat{\mathbf A}_1$ (also $\widehat{\mathbf P}_1$), and $p^2=o(\min\{n^{\frac{2}{5\delta_2}},n^{\frac{1}{2\delta_1+\delta_2}}\} )$ for that of $\widehat{\mathbf B}_2^*$ and $\widehat{\mathbf Q}_2^*$, which is slightly stronger than the former.
\end{remark}
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.
\begin{theorem}
Under the conditions in Theorem 2, we have
\begin{align*}
(p_1p_2)^{-1/2}\|\widehat{\mathbf A}_1\widehat{\mathbf X}_t\widehat{\mathbf P}_1'-{\mathbf A}_1{\mathbf X}_t{\mathbf P}_1'\|_2=&O_p\left(p_1^{-\delta_1/2}p_2^{-\delta_1/2}(\|D(\mathcal{M}({\widehat{\mathbf A}}_1),\mathcal{M}{{\mathbf A}}_1)\|_2\right. \notag\\
&\left.+\|D(\mathcal{M}({\widehat{\mathbf P}}_1),\mathcal{M}({{\mathbf P}}_1))\|_2)+p_1^{-\delta_2/2}\|D(\mathcal{M}(\widehat{\mathbf B}_2^*),\mathcal{M}({\mathbf B}_2^*))\|_2\right.\\
&\left.+p_2^{-\delta_2/2}\|D(\mathcal{M}(\widehat{\mathbf Q}_2^*),\mathcal{M}({\mathbf Q}_2))\|_2+p_1^{-1/2}p_2^{-1/2}\right).
\end{align*}
\end{theorem}
\begin{remark}
(i) A similar result is given in Theorem 3 of \cite{LamYaoBathia_Biometrika_2011} and Theorem 5 of \cite{gaotsay2018b}, which deal with the approximate factor model and a structured factor model, respectively. When $\delta_1=\delta_2=0$, i.e. the factors and the noise terms are all strong, the convergence rate in Theorem 3 is $O_p((p_1p_2)^{-1/2}+n^{-1/2})$, which is the optimal rate specified in Theorem 3 of \cite{Bai_Econometrica_2003} when dealing with the traditional approximate factor models.\\
(ii) It is a common issue to select the number of principle components in the literature and there are many possible approaches available.
Since it is impossible to eliminate all the noise effects in recovering the factors
and we only need to guarantee that the diverging part of the noises
is removed for large $p_1$, we may select $k_1$ in a range of possible values.
In practice, let $\widehat\mu_{1,1}\geq...\geq \widehat\mu_{1,p_1}$ be the sample eigenvalues of $\widehat{\mathbf S}_1$ and define $\widehat k_{1,L}$ as
\begin{equation}\label{kl}
\widehat k_{1,L}=\arg\min_{1\leq j\leq \widehat k_{1,U}}\{\widehat\mu_{1,j+1}/\widehat\mu_{1,j}\},
\end{equation}
and $\widehat k_{1,U}$ is a pre-specified integer. We suggest $\widehat k_{1,U}=\min\{\sqrt{p_1},\sqrt{n},p_1-\widehat r_1,5\}$. Then the estimator $\widehat k_1$ for $k_1$ can assume some value between $\widehat k_{1,L}$ and $\widehat k_{1,U}$. We can select $\widehat k_2$ in a similar manner.
\end{remark}
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 \cite{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.
\begin{assumption}
$\textnormal{vec}({{\mathbf F}_t})$, $\textnormal{vec}({\mathbf Z}_{12,t})$, $\textnormal{vec}({\mathbf Z}_{21,t})$, and $\textnormal{vec}({\mathbf Z}_{22,t})$ are sub-Gaussian random vectors.
\end{assumption}
\begin{theorem}
Assume Assumptions 1-6 hold.\\
(i) If $p_1p_2=o\left\{\min\left(n^{\frac{2}{1+3\delta_1}},n^{\frac{1}{1+2\delta_1-\delta_2}}\right)\right\}$, then the test statistic $T_n$ can consistently estimate $r_1$ and $r_2$, i.e. $P(\widehat r_1=r_1,\widehat r_2=r_2)\rightarrow 1$ as $n\rightarrow\infty$.\\
(ii) If $p_1^{1+\delta_1-\delta_2/2}p_2^{1+\delta_1-\delta_2/2}n^{-1/2}\sqrt{\log(np_1p_2)}=o(1)$, then
the test statistic $T(m)$ can consistently estimate $r_1$ and $r_2$.
\end{theorem}
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{m1hat}), 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.
\section{Numerical Properties}
\subsection{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{eq:D}) 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
\begin{equation}\label{dmeasure}
\bar{D}(\mathcal{M}({\mathbf H}_1),\mathcal{M}({\mathbf H}_2))=\sqrt{1-
\frac{1}{\max{(h_1,h_2)}}\textrm{tr}({\mathbf P}_1{\mathbf P}_2)}.
\end{equation}
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{eq:D}). We only present the simulation results for $k_0=2$ in Equation (\ref{m1hat}) to save space since other choices of $k_0$ produce similar patterns.\\
{\noindent \bf Example 1.} Consider model (\ref{m-factor}) 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 \cite{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{Table1}. 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{m-factor}), yet the performance improves significantly when the sample size increases.
\begin{table}
\caption{Empirical probabilities $P(\widehat{r}_2=r_1, \widehat{r}_2=r_2)$ for Example 1 with $(r_1,r_2)=(2,3)$ and $(k_1,k_2)=(1,2)$, where $(p_1,p_2)$ and $n$ are the dimension and the sample size, respectively. $\delta_1$ and $\delta_2$ are the strength parameters of the factors and the errors, respectively. $500$ iterations are used.}
\label{Table1}
\begin{center}
\setlength{\abovecaptionskip}{0pt}
\setlength{\belowcaptionskip}{3pt}
\begin{tabular}{c|cc|ccccc}
\hline
& &&\multicolumn{5}{c}{$n$}\\
$(\delta_1, \delta_2)$ &$(p_1,p_2)$&$p_1p_2$&$300$&$500$&$1000$&$1500$&$3000$\\
\hline
(0,0.9) &$(7,7)$&49&0.956&0.982&0.984&0.980&0.976\\
&$(10,15)$&150&0.930&0.988&0.964&0.984&0.978\\
&$(20,20)$&400&0.818&0.976&0.962&0.970&0.968\\
&$(20,30)$&600&0.834&0.986&0.976&0.964&0.972\\
\hline
(0.2,0.8) &$(7,7)$&49&0.848&0.992&0.988&0.972&0.986\\
&$(10,15)$&150&0.882&0.982&0.974&0.978&0.984\\
&$(20,20)$&400&0.742&0.964&0.972&0.982&0.972\\
&$(20,30)$&600&0.816&0.994&0.976&0.968&0.966\\
\hline
(0.5,0.5) &$(7,7)$&49&0.104&0.438&0.950&0.974&0.972\\
&$(10,15)$&150&0.304&0.710&0.946&0.974&0.980\\
&$(20,20)$&400&0.028&0.074&0.334&0.696&0.980\\
&$(20,30)$&600&0.020&0.080&0.296&0.636&0.938\\
\hline
\end{tabular}
\end{center}
\end{table}
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{fig1}(a) and (b), respectively. From Figure~\ref{fig1}, 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
\begin{equation}\label{DX}
D(\widehat{\mathbf A}_1\widehat{\mathbf X}\widehat{\mathbf P}_1',{\mathbf L}_1{\mathbf F}{\mathbf R}_1')=\frac{1}{n\sqrt{p_1p_2}}\sum_{t=1}^n\|\widehat{\mathbf A}_1\widehat{\mathbf X}_t\widehat{\mathbf P}_1'-{\mathbf L}_1{\mathbf F}_t{\mathbf R}_1\|_2.
\end{equation}
The results are shown in Figure~\ref{fig2}, 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.
\begin{figure}
\begin{center}
\subfigure[]{\includegraphics[width=0.49\textwidth]{DA-1.eps}}
\subfigure[]{\includegraphics[width=0.49\textwidth]{DP-1.eps}}
\caption{(a) Boxplots of $\bar{D}(\mathcal{M}(\widehat{\mathbf A}_1),\mathcal{M}({\mathbf L}_1))$; (b) Boxplots of $\bar{D}(\mathcal{M}(\widehat{\mathbf P}_1),\mathcal{M}({\mathbf R}_1))$. We set $(r_1,r_2)=(2,3)$, $(k_1,k_2)=(1,2)$, and $(\delta_1,\delta_2)=(0,0.9)$ in Example 1. The sample sizes are $300, 500, 1000, 1500, 3000$, respectively. $500$ iterations are used.}\label{fig1}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
{\includegraphics[width=0.7\textwidth]{AXP-1.eps}}
\caption{Boxplots of $D(\widehat{\mathbf A}_1\widehat{\mathbf X}\widehat{\mathbf P}_1,{\mathbf L}_1{\mathbf F}{\mathbf R}_1')$ defined in (\ref{DX}) when $(r_1,r_2)=(2,3)$, $(k_1,k_2)=(1,2)$, and $(\delta_1,\delta_2)=(0,0.9)$ in Example 1. The sample sizes are $300, 500, 1000, 1500, 3000$, respectively. $500$ iterations are used.}\label{fig2}
\end{center}
\end{figure}
To see the advantages of the proposed method, we compare it with that of \cite{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
\begin{equation}\label{ratio}
\widehat r_i=\arg\min_{1\leq j\leq p_i/2}\{\widehat\lambda_{i,j+1}/\widehat\lambda_{i,j}\},\,\, i=1,2.
\end{equation}
Figure~\ref{fig3}(a)-(b) present the boxplots of $\widehat r_1$ and $\widehat r_2$, respectively. We see from Figure~\ref{fig3} 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 \cite{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{Table1}.
\begin{figure}
\begin{center}
\subfigure[]{\includegraphics[width=0.49\textwidth]{WLC-r1.eps}}
\subfigure[]{\includegraphics[width=0.49\textwidth]{WLC-r2.eps}}
\caption{(a) Boxplots of $\widehat r_1$ by the ratio-based method of \cite{wang2018}; (b) Boxplots of $\widehat r_2$ by the ratio-based method of \cite{wang2018}. We set $(r_1,r_2)=(2,3)$, $(k_1,k_2)=(1,2)$, and $(\delta_1,\delta_2)=(0.5,0.5)$ in Example 1. The sample sizes are $300, 500, 1000, 1500, 3000$, respectively. $500$ iterations are used.}\label{fig3}
\end{center}
\end{figure}
Finally, we compare our method with the one of \cite{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{Table2} 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.
\begin{table}
\caption{The $D(\widehat{\mathbf A}_1\widehat{\mathbf X}\widehat{\mathbf P}_1',{\mathbf L}_1{\mathbf F}{\mathbf R}_1')$ defined in (\ref{DX}) when $(r_1,r_2)=(2,3)$, $(k_1,k_2)=(1,2)$, and $(\delta_1,\delta_2)=(0.5,0.5)$ in Example 1.
The sample sizes used are $n=300, 500, 1000, 1500, 3000$. Standard errors are given in the parentheses and $500$ iterations are used. GT denotes the proposed method and WLC is the one in \cite{wang2018}.}
\label{Table2}
\footnotesize{
\begin{center}
\setlength{\abovecaptionskip}{0pt}
\setlength{\belowcaptionskip}{3pt}
\begin{tabular}{c|c|ccccc}
\hline
&&\multicolumn{5}{c}{$n$}\\
$(p_1,p_2)$&Method&$300$&$500$&$1000$&$1500$&$3000$\\
\hline
$(7,7)$>& 0.862(0.199)& 0.726(0.494)& 0.460(0.164)& 0.447(0.320)& 0.416(0.096)\\
&WLC&1.178(0.036)&1.178(0.029)&1.182(0.026)&1.183(0.025)&1.179(0.047)\\
\hline
$(10,15)$>& 0.652(0.354)& 0.403(0.185)& 0.290(0.310)& 0.254(0.208)& 0.229(0.132)\\
&WLC&0.891(0.060)&0.886(0.066)&0.862(0.084)&0.783(0.148)&0.549(0.165)\\
\hline
$(20,20)$>& 0.530(0.167)& 0.437(0.108)& 0.301(0.147)& 0.191(0.126)& 0.103(0.043)\\
&WLC&0.696(0.011)&0.695(0.010)&0.686(0.036)&0.648(0.072)&0.485(0.117)\\
\hline
$(20,30)$>& 0.485(0.181)& 0.394(0.102)& 0.278(0.122)& 0.181(0.124)& 0.098(0.069)\\
&WLC&0.662(0.010)&0.663(0.008)&0.662(0.005)&0.662(0.005)&0.651(0.044)\\
\hline
\end{tabular}
\end{center}}
\end{table}
\subsection{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 \cite{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 \cite{sharpe1964} and \cite{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{100series} with $p_1=p_2=10$ and $n=678$.
\begin{figure}
\begin{center}
{\includegraphics[width=0.75\textwidth]{ME-INV-data.eps}}
\caption{Time series plots of Fama-French 10 by 10 monthly excess return series based on Size and Investment from July 1963 to December 2019.}\label{100series}
\end{center}
\end{figure}
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 \cite{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{loading}, 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 \cite{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.
\begin{table}
\caption{Fama-French series: Size and Investment (Inv) loading matrices after being multiplied by $30$. The two-dimensional loading vectors are ordered for sizes (S1--S10) and Investment (Inv1--Inv10) from small to large and from low to high, respectively.}
\label{loading}
\begin{center}
\begin{tabular}{c|cccccccccc}
\hline
Size Factor & S1 & S2 & S3 & S4 & S5 & S6 & S7 & S8 & S9 & S10 \\ \hline
Row 1 & \cellcolor[gray]{0.5}{-18} & \cellcolor[gray]{0.6}{-12} &
\cellcolor[gray]{0.8}{-10} & \cellcolor[gray]{0.8}{-9} &
\cellcolor[gray]{0.8}{-8}
& \cellcolor[gray]{0.8}{-7} &\cellcolor[gray]{0.8}{ -7} & \cellcolor[gray]{0.8}{-6}& \cellcolor[gray]{0.8}{-5}&\cellcolor[gray]{0.8} {-3} \\
Row 2 & \cellcolor[gray]{0.5}{21} & {1} &
\cellcolor[gray]{0.8}{-6} & \cellcolor[gray]{0.8}{-6} &
\cellcolor[gray]{0.8}{-6}
& \cellcolor[gray]{0.8}{-10} &\cellcolor[gray]{0.8}{ -7} & \cellcolor[gray]{0.8}{-8}& \cellcolor[gray]{0.8}{-8}& \cellcolor[gray]{0.8}{-7} \\
\hline
\hline
Inv Factor & Inv1 & Inv2 & Inv3 & Inv4 & Inv5 & Inv6 & Inv7 & Inv8 & Inv9 & Inv10 \\ \hline
Column 1 & {\cellcolor[gray]{0.8} -11} & {\cellcolor[gray]{0.8} -10}
& {\cellcolor[gray]{0.8} -9} & {\cellcolor[gray]{0.8} -8}
& {\cellcolor[gray]{0.8} -8} & {\cellcolor[gray]{0.8} -8} & {\cellcolor[gray]{0.8}-9 }& {\cellcolor[gray]{0.8}-9} & {\cellcolor[gray]{0.8}-10} & {\cellcolor[gray]{0.8} -12 }\\
Column 2 & {\cellcolor[gray]{0.4} 27} & { -3} & {-1} & 1 & { \cellcolor[gray]{0.8}-6}
& {\cellcolor[gray]{0.8}-6} & {\cellcolor[gray]{0.8} -5}
& {\cellcolor[gray]{0.8} -5} & { 1}
& {\cellcolor[gray]{0.8} -7}\\
\hline
\end{tabular}
\end{center}
\end{table}
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{fig4}, 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{fig5}(a) as well as their corresponding spectrum in Figure~\ref{fig5}(b). From Figure~\ref{fig5}, 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 \cite{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.
\begin{figure}
\begin{center}
{\includegraphics[width=0.7\textwidth]{s-eigens.eps}}
\caption{(a) The 10 eigenvalues of $\widehat{\mathbf S}_1$; (b) The plot of ratios of consecutive eigenvalues of $\widehat{\mathbf S}_1$; (c) The 10 eigenvalues of $\widehat{\mathbf S}_2$; (d) The plot of ratios of consecutive
eigenvalues of $\widehat{\mathbf S}_2$}.\label{fig4}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\subfigure[]{\includegraphics[width=0.45\textwidth]{ext-factors.eps}}
\subfigure[]{\includegraphics[width=0.45\textwidth]{ft-spectrum.eps}}
\caption{(a) The time series plots of the extracted $2\times 2$ common factors; (b) the corresponding spectrum of the factor processes}.\label{fig5}
\end{center}
\end{figure}
Next we examine and compare the forecasting performance of the extracted factors via the proposed method (denoted by GT) and those by \cite{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 \cite{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:
\begin{equation}\label{fef}
\text{FE}_F(h)=\frac{1}{120-h+1}\sum_{\tau=558}^{678-h}\frac{1}{\sqrt{p_1p_2}}\|\widehat{\mathbf Y}_{\tau+h}-{\mathbf Y}_{\tau+h}\|_F,
\end{equation}
and
\begin{equation}\label{fe2}
\text{FE}_2(h)=\frac{1}{120-h+1}\sum_{\tau=558}^{678-h}\frac{1}{\sqrt{p_1p_2}}\|\widehat{\mathbf Y}_{\tau+h}-{\mathbf Y}_{\tau+h}\|_2,
\end{equation}
where $p_1=p_2=10$. Table~\ref{Table4} reports the 1-step to 4-step ahead forecast errors of Equations (\ref{fef}) and (\ref{fe2}) 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{Table4} 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.
\begin{table}[h]
\caption{The 1-step to 4-step ahead out-of-sample forecast errors of various methods for Example 2. GT denotes the proposed method, WLC denotes the forecasting errors based on the extracted factor by the method in \cite{wang2018}, and SAR denotes a scalar AR model to each individual return series.
Boldface numbers denote the smallest error for a given forecast horizon.}
\label{Table4}
\begin{center}
\setlength{\abovecaptionskip}{0pt}
\setlength{\belowcaptionskip}{3pt}
\begin{tabular}{c|cccccccc}
\hline
&\multicolumn{3}{c}{FE$_{F}(h)$} &&&\multicolumn{3}{c}{FE$_{2}(h)$}\\
\cline{2-4}\cline{7-9}
Step $h$>&WLC&SAR&&>&WLC&SAR\\
\cline{1-4}\cline{7-9}
1&{\bf 4.51}&4.61&4.60&&&{\bf 4.04}&4.14&4.14\\
2&{\bf 4.47}&4.51&4.60&&& {\bf 3.98}&4.04&4.14\\
3&{\bf 4.48}&4.51&4.60&&&{\bf 4.00}&4.04&4.14\\
4&{\bf 4.47}&4.49&4.57&&&{\bf 3.98}&4.01&4.11\\
\hline
\end{tabular}
\end{center}
\end{table}
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.
\section{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 \cite{TiaoTsay_1989} and \cite{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.