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.
135,416 characters
Identification and Estimation for Matrix Time Series CP-factor Models
\if11
{
\spacingset{1.25}
\title{\bf \Large Identification and Estimation for Matrix Time Series CP-factor Models}
\author[1,2]{Jinyuan Chang}
\author[1]{Yue Du}
\author[1]{Guanglin Huang}
\author[3]{Qiwei Yao}
\affil[1]{\it \small Joint Laboratory of Data Science and Business
Intelligence, Southwestern University of Finance and Economics, Chengdu, China}
\affil[2]{\it \small State Key Laboratory of Mathematical Sciences, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing, China}
\affil[3]{\it \small Department of Statistics, The London School of Economics and Political Science, London, U.K.}
\setcounter{Maxaffil}{0}
\date{\vspace{-5ex}}
\maketitle
} \fi
\if01
{
\bigskip
\bigskip
\bigskip
\begin{center}
{
\Large \bf Identification and Estimation for Matrix Time Series CP-factor Models
}
\end{center}
\medskip
} \fi
\spacingset{1.5}
\begin{abstract}
We propose a new method for identifying and estimating the CP-factor models for matrix time series. Unlike the generalized eigenanalysis-based method of \cite{chang2023modelling} for which the convergence rates of the associated estimators may suffer from small eigengaps as the asymptotic theory is based on some matrix perturbation analysis, the proposed new method enjoys faster convergence rates which are free from any eigengaps. It achieves this by turning the problem into a joint diagonalization of several matrices whose elements are determined by a basis of a linear system, and by choosing the basis carefully to avoid near co-linearity (see Proposition \ref{pro:theta-unique} and Section \ref{sec:theta-hat-est}). Furthermore, unlike \cite{chang2023modelling} which
requires the two factor loading matrices to be full-ranked, the proposed new method can handle rank-deficient factor loading matrices.
Illustration with both simulated and real matrix time series data shows the advantages of the proposed new method.
\end{abstract}
\noindent {\sl Keywords}: CP-decomposition; dimension-reduction; matrix time series; non-orthogonal joint diagonalization.
\spacingset{1.69}
\setlength{\abovedisplayskip}{0.2\baselineskip}
\setlength{\belowdisplayskip}{0.2\baselineskip}
\setlength{\abovedisplayshortskip}{0.2\baselineskip}
\setlength{\belowdisplayshortskip}{0.2\baselineskip}
\section{Introduction}
The modern capacity for data collection has resulted in an abundance of time series data, with those in high-dimensional matrix format increasingly prevalent across diverse fields such as economics, finance, engineering, environmental sciences, medical research, network traffic monitoring, image processing and others. The demand of modeling and forecasting high-dimensional matrix time series brings the opportunities with challenges. Let ${\mathbf Y}_t = (y_{i,j,t})$ be a $p \times q$ matrix recorded at time $t$, where $y_{i,j,t}$ represents the value of, for example, the $j$-th variable on the $i$-th individual at time $t$. A popular approach to model ${\mathbf Y}_t$ in the existing literature is via the so-called Tucker decomposition, namely the matrix Tucker-factor model. See, for example, \cite{wang2019factor}, \cite{chen2019constrained}, \cite{chen_chen_2022}, and \cite{han2024tensor}. It represents a high-dimensional matrix time series as a linear combination of a lower-dimensional matrix process. The Tucker decomposition can be viewed as a natural extension of the factor model for vector time series considered in \cite{lam2012factor} and \cite{Chang2015}. Similarly we can only identify the factor loading spaces (the linear spaces spanned by the columns of the factor loading matrices) in the matrix Tucker-factor model while the factor loading matrices themselves are not uniquely defined. Parallel to the approaches based on Tucker decomposition,
\cite{chang2023modelling} and \cite{han2024cp}
consider to model ${\mathbf Y}_t$ via the so-called canonical polyadic (CP) decomposition, namely the matrix CP-factor model. It provides a more comprehensive dimensionality reduction as the dynamic structure of a matrix time series is driven by a vector process rather than a matrix process. Furthermore the factor loading matrices in the matrix CP-factor model can be identified uniquely up to the column reflection and permutation indeterminacy under some regularity conditions.
The CP-factor model for matrix time series ${\mathbf Y}_t$
admits the form
\begin{equation}\label{eq:abm}
{\mathbf Y}_t = {\mathbf A} {\mathbf X}_{t} {\mathbf B}^{{{\mathrm{\scriptscriptstyle \top} }}} + \boldsymbol{\varepsilon}_t\,, ~~~~ t\geq1\,,
\end{equation}
where ${\mathbf X}_t={\rm diag}(\mathbf{x}_t)$ with $\mathbf{x}_t =(x_{t,1},\ldots, x_{t,d})^{\mathrm{\scriptscriptstyle \top} }$ being a $d \times 1$ time series, $\boldsymbol{\varepsilon}_t$ is a $p\times q$ matrix white noise, and ${\mathbf A}=({\mathbf a}_1,\ldots,{\mathbf a}_d)$ and ${\mathbf B}=({\mathbf b}_1,\ldots,{\mathbf b}_d)$ are, respectively, $p\times d$ and $q \times d$ constant matrices which are called
factor loading matrices. See, for example, \cite{chang2023modelling}.
Without loss of generality, we assume $|{\mathbf a}_\ell|_2=1= |{\mathbf b}_\ell|_2$ for each $\ell = 1,\ldots,d$.
For matrix CP-factor model \eqref{eq:abm}, we cannot observe $({\mathbf A},{\mathbf B},{\mathbf X}_t,\boldsymbol{\varepsilon}_t)$ and only assume $1 \le d < \min(p,q)$ is an unknown fixed integer.
Based on the assumption ${\rm rank}({\mathbf A}) =d ={\rm rank}({\mathbf B}) $, \cite{chang2023modelling} proposes a one-pass estimation procedure for $(d,{\mathbf A},{\mathbf B})$ which
identifies $({\mathbf A}, {\mathbf B})$ uniquely up to the column reflection and permutation indeterminacy. In contrast to the standard alternating least squares method and its variations \citep{han2022tensor, han2024cp},
the estimation procedure proposed in \cite{chang2023modelling} is based on solving some generalized eigenequations and requires no iterations.
Note that the incoherence conditions imposed in \cite{han2022tensor} and \cite{han2024cp} also require both ${\mathbf A}$
and ${\mathbf B}$ to be full-ranked. In fact those conditions imply that both $\{ {\mathbf a}_\ell\}_{\ell=1}^d$ and $\{ {\mathbf b}_\ell\}_{\ell=1}^d$ are two sets of near-orthogonal vectors. We do not require such an incoherence condition in this paper.
In this paper, we investigate the identification issue of the CP-factor model \eqref{eq:abm} for matrix time series without imposing the condition ${\rm rank}({\mathbf A})=d={\rm rank}({\mathbf B})$. Let
\begin{equation*}
{\rm rank}({\mathbf A})=d_1~~\textrm{and}~~{\rm rank}({\mathbf B})=d_2\,.
\end{equation*}
Then $1\le d_1, d_2 \le d$.
As the CP-decomposition for 3-way tensors often exhibits rank-deficient factor loading matrices
\citep{kolda2009tensor}, i.e., in model \eqref{eq:abm} it may hold that
$\max(d_1, d_2)<d$.
We identify the condition under which ${\mathbf A}$ and ${\mathbf B}$ are uniquely identifiable up to the column reflection and permutation indeterminacy. Our setting allows all scenarios in terms of the relationships among $d_1,\, d_2 $ and $ d$.
The proposed new estimation procedure consists of several steps (see Section \ref{sec:estimation}). The key idea is to transform the $p\times q$ matrix CP-factor model (\ref{eq:abm}) to a $(d_1d_2)$-vector factor model, and then to identify the columns of ${\mathbf A}$ and ${\mathbf B}$ by a joint diagonalization of several symmetric matrices whose elements are determined by a basis, and in fact any basis, of a linear system (see Proposition \ref{pro:theta-unique}). Therefore, we can choose an appropriate basis to avoid near co-linearity such that our estimator enjoys faster convergence rate than those eigenanalysis-based estimators (see Section \ref{sec:theta-hat-est}).
Note that the convergence rates of the eigenanalysis-based estimators are derived based on some matrix perturbation analysis, and may suffer from the adverse
impact of
eigen-gap (i.e., the minimum pairwise gap among a set of eigenvalues). Our newly proposed estimator is free from this adversity.
For example, the convergence rate of the estimator of
\cite{chang2023modelling} can be formulated as the product of the rate of our new estimator and the inverse of an eigen-gap (See Remark \ref{rek:chang}). Note that the eigen-gap typically diminishes to 0 when $p$ or/and $q$ diverge to infinity.
The rest of the paper is organized as follows. Section \ref{sec:pre} gives preliminaries of the matrix CP-factor model \eqref{eq:abm}. A general identification strategy for the matrix CP-factor model is presented in Section \ref{sec:identify-ab}. Section \ref{sec:estimation} provides a one-pass estimation procedure for $(d_1,d_2,d,{\mathbf A},{\mathbf B})$. Section \ref{sec:prediction} gives a unified prediction approach for the matrix CP-factor model. We investigate the associated theoretical properties of the proposed method in Section \ref{sec:asymptotics}. Numerical results with simulation studies and real data analysis are given in Section \ref{section:simulation}. The \textsf{R}-function \texttt{CP\_MTS} for implementing our newly proposed method is available publicly in the \texttt{HDTSA} package \citep{chang2024hdtsa}. All technical proofs and some additional simulation studies are relegated in the supplementary material.
\textit{Notation}. For a positive integer $m$, write $[m] = \{1, \ldots , m\}$, and denote by ${\mathbf I}_m$ the $m \times m$ identity matrix. Denote by $I(\cdot)$ the indicator function. For an $m_1 \times m_2$ matrix ${\mathbf H} = (h_{i,j} )_{m_1 \times m_2}$, let $\mathcal{R}({\mathbf H})=\max\{k:{\text{any}\ k\ \text{columns of the matrix ${\mathbf H}$ are linearly independent}}\}$, and denote by $\mathcal{M}({\mathbf H})$ the linear space spanned by the columns of ${\mathbf H}$.
Let $\|{\mathbf H}\|_2$, $\|{\mathbf H}\|_\text{F}$, ${\rm rank}({\mathbf H})$, $\lambda_{i}({\mathbf H})$, and $\sigma_{i}({\mathbf H})$ be, respectively, the spectral norm, Frobenius norm, rank, $i$-th largest eigenvalue, and $i$-th largest singular value of matrix ${\mathbf H}$.
Specifically, if $m_2 = 1$, we use $|{\mathbf H}|_1 = \sum_{i=1}^{m_1}|h_{i,1} |$ and $|{\mathbf H}|_2=(\sum_{i=1}^{m_1}h_{i,1}^2)^{1/2}$
to denote, respectively, the $L_1$-norm and $L_2$-norm of the $m_1$-dimensional vector ${\mathbf H}$.
Also, denote by ${\mathbf H}^{{{\mathrm{\scriptscriptstyle \top} }}}$ and ${\mathbf H}^{+}$, respectively, the transpose and the Moore-Penrose inverse of ${\mathbf H}$.
The operator ${\rm diag}(\cdot)$ stacks a vector into a square diagonal matrix. Let $\otimes$ denote the Kronecker product, and
$\odot$ denote the Khatri-Rao product such that $\check{{\mathbf H}} \odot \tilde{{\mathbf H}} = (\check{{\mathbf h}}_1\otimes\tilde{{\mathbf h}}_1,\ldots,\check{{\mathbf h}}_m\otimes\tilde{{\mathbf h}}_m)$ for any matrices $\check{{\mathbf H}} = (\check{{\mathbf h}}_1,\ldots,\check{{\mathbf h}}_m)$ and $\tilde{{\mathbf H}} = (\tilde{{\mathbf h}}_1,\ldots,\tilde{{\mathbf h}}_m)$.
Moreover, for any two sequences of positive numbers $\{\tau_{k}\}$ and $\{\tilde{\tau}_{k}\}$, we write $\tau_{k} \asymp \tilde{\tau}_{k}$ if $\tau_{k}/\tilde{\tau}_{k}=O(1)$ and $\tilde{\tau}_{k}/\tau_{k}=O(1)$ as $k\rightarrow\infty$, and write $\tau_k \ll \tilde{\tau}_k$ or $ \tilde{\tau}_k \gg \tau_k$ if $\lim\sup_{k\to \infty} \tau_k/\tilde{\tau}_k=0$. To simplify our presentation, for a matrix ${\mathbf H} =(h_{i,j})_{m_1\times m_2}$, we write $\vec{\mathbf H}$ or ${\rm vec}({\mathbf H})$ as an $(m_1m_2)$-dimensional vector with the $\{(j-1)m_1 +i\}$-th element being $h_{i,j}$, and for a tensor $\mathcal{H}=(h_{i,j,k,l})_{m_1\times m_2\times m_3 \times m_4}$, we write $\vec{\mathcal{H}}$ as an $(m_1m_2m_3m_4)$-dimensional vector with the $\{(i-1)m_2m_3m_4 + (j-1)m_3m_4 + (k-1)m_4 + l\}$-th element being $h_{i,j,k,l}$.
\section{Preliminary}\label{sec:pre}
Recall that, in the matrix CP-factor model \eqref{eq:abm}, ${\mathbf A}$ and ${\mathbf B}$ are, respectively, $p\times d$ and $q\times d$ matrices with ${\rm rank}({\mathbf A})=d_1$ and ${\rm rank}({\mathbf B})=d_2$, and $d_1,d_2 \in [d]$.
Model \eqref{eq:abm} can be equivalently represented as
\begin{equation*}
\vec{\mathbf Y}_t=({\mathbf B} \odot {\mathbf A})\mathbf{x}_t + \vec\boldsymbol{\varepsilon}_t\,,~~~~t\geq1\,,
\end{equation*}
where ${\mathbf B} \odot {\mathbf A} =({\mathbf b}_1\otimes {\mathbf a}_1,\ldots, {\mathbf b}_{d}\otimes {\mathbf a}_{d})$ and $\mathbf{x}_t= (x_{t,1},\ldots, x_{t,d})^{{\mathrm{\scriptscriptstyle \top} }}$.
Condition \ref{cd:ra}(i) below holds naturally.
If ${\rm rank}({\mathbf B} \odot {\mathbf A}) = \tilde{d} < d$, the matrix ${\mathbf B} \odot {\mathbf A}$ has $\tilde{d}$ linearly independent columns that span its column space. Therefore, we can find $\{{\mathbf b}_{\ell_1}\otimes {\mathbf a}_{\ell_1},\ldots, {\mathbf b}_{\ell_{\tilde{d}}}\otimes {\mathbf a}_{\ell_{\tilde{d}}}\}$ with some distinct $\ell_1,\ldots,\ell_{\tilde{d}} \in [d]$ such that they provide a basis for $\mathcal{M}({\mathbf B} \odot {\mathbf A})$. The remaining columns of ${\mathbf B} \odot {\mathbf A}$ can be expressed as linear combinations of this set of basis vectors. Then ${\mathbf B} \odot {\mathbf A} = ({\mathbf b}_{\ell_1}\otimes {\mathbf a}_{\ell_1},\ldots, {\mathbf b}_{\ell_{\tilde{d}}}\otimes {\mathbf a}_{\ell_{\tilde{d}}}) \tilde{{\mathbf C}}$ for some $\tilde{d} \times d$ matrix $\tilde{{\mathbf C}}$. Since $({\mathbf A},{\mathbf B},{\mathbf X}_t)$ are unobserved, we can reformulate $\vec{\mathbf Y}_t$ in a new form $\vec{\mathbf Y}_t=(\tilde{{\mathbf B}} \odot \tilde{{\mathbf A}})\tilde{\mathbf{x}}_t + \vec\boldsymbol{\varepsilon}_t$ with $\tilde{{\mathbf A}}=({\mathbf a}_{\ell_1},\ldots,{\mathbf a}_{\ell_{\tilde{d}}})$, $\tilde{{\mathbf B}}= (\tilde{{\mathbf b}}_{\ell_1},\ldots,\tilde{{\mathbf b}}_{\ell_{\tilde{d}}})$ and $\tilde{\mathbf{x}}_t=\tilde{{\mathbf C}}\mathbf{x}_t$. In this new form, the newly defined factor loading matrices $\tilde{{\mathbf A}}\in \mathbb{R}^{p \times \tilde{d}}$ and $\tilde{{\mathbf B}} \in \mathbb{R}^{q \times \tilde{d}}$ satisfy $\textup{rank}(\tilde{{\mathbf B}}\odot\tilde{{\mathbf A}}) = \tilde{d}$.
On the other hand, since $\boldsymbol{\varepsilon}_t$ is a matrix white noise, Condition \ref{cd:ra}(ii) holds automatically.
\begin{cd}\label{cd:ra}
{\rm(i)} ${\rm rank}({\mathbf B} \odot {\mathbf A})=d$. {\rm(ii)} $\mathbb{E}(\boldsymbol{\varepsilon}_t)=\bf{0}$ for any $t\geq1$, $\mathbb{E}(\boldsymbol{\varepsilon}_t \otimes \boldsymbol{\varepsilon}_s)=\bf{0}$ for all $t\ne s$, and $\mathbb{E}(x_{t,\ell}\boldsymbol{\varepsilon}_s)=\bf{0}$ for any $\ell \in[d]$ and $t\le s$.
\end{cd}
When $d_1=d_2=d$, \cite{chang2023modelling} provides a one-pass estimator for $({\mathbf A},{\mathbf B})$ by solving some generalized eigenequations defined by the matrices
\begin{equation}\label{eq:sigma_yxi}
\boldsymbol{\Sigma}_{{\mathbf Y},\xi}(k) = \frac{1}{n-k}\sum_{t = k+1}^{n}\mathbb{E}[ \{ {\mathbf Y}_t -\mathbb{E}(\bar{{\mathbf Y}})\}\{ \xi_{t-k} -\mathbb{E}(\bar{\xi})\} ]\,,~~~~k\geq1\,,
\end{equation}
where $\bar{{\mathbf Y}} = n^{-1}\sum_{t = 1}^{n}{\mathbf Y}_t$, $\xi_t$ is a scalar defined as a linear combination of the elements of ${\mathbf Y}_t$, and $\bar{\xi} = n^{-1}\sum_{t = 1}^{n}\xi_t$. For example, we can select $\xi_{t}$ as the first principal component of
$\vec {\mathbf Y}_t$.
Recall ${\mathbf A}^+$ and $ {\mathbf B}^+$ are, respectively,
the Moore-Penrose inverse of ${\mathbf A}$ and ${\mathbf B}$. The key requirement underlying the results of \cite{chang2023modelling} is ${\mathbf A}^+{\mathbf A} = {\mathbf B}^+{\mathbf B}={\mathbf I}_d$, which only holds when $d_1=d_2=d$.
Hence, the estimation method of \cite{chang2023modelling} is not applicable when $\min(d_1, d_2) <d$. Note that the CP-decomposition for 3-way tensors can often exhibit rank-deficient factor loading matrices \citep{kolda2009tensor}, i.e., in model (\ref{eq:abm}) it may hold that $\min(d_1, d_2) < d$ or even $\max(d_1, d_2) < d$.
In this paper, we consider a new approach which
identifies $(d,{\mathbf A},{\mathbf B})$ without the condition $d_1=d_2=d$. Furthermore we propose a
unified and more efficient one-pass estimation
for $({\mathbf A},{\mathbf B})$ regardless they are rank-deficient or not.
\section{Identification of
$({\mathbf A},{\mathbf B})$}\label{sec:identify-ab}
We need to identify in model \eqref{eq:abm} the order $d$ and the factor loading pairs $({\mathbf a}_1, {\mathbf b}_1),\ldots,({\mathbf a}_d,{\mathbf b}_d)$.
To carry out this task, we first introduce a reduced model for a $d_1 \times d_2$ matrix time series, and then identify $d$ and the CP-factor loadings for the reduced model via (i) a factor model for a vector time series, and (ii) a non-orthogonal joint diagonalization of $d$ symmetric matrices.
\subsection{A reduced model}
For a prescribed integer $K>1$ and
$\boldsymbol{\Sigma}_{{\mathbf Y}, \xi}(k)$ specified in \eqref{eq:sigma_yxi}, define
\begin{equation}\label{eq:M1M2}
{\mathbf M}_1=\sum_{k=1}^{K}\boldsymbol{\Sigma}_{{\mathbf Y}, \xi}(k)\boldsymbol{\Sigma}_{{\mathbf Y}, \xi}(k)^{{{\mathrm{\scriptscriptstyle \top} }}}~~\textrm{and}~~ {\mathbf M}_2=\sum_{k=1}^{K}\boldsymbol{\Sigma}_{{\mathbf Y}, \xi}(k)^{{{\mathrm{\scriptscriptstyle \top} }}}\boldsymbol{\Sigma}_{{\mathbf Y}, \xi}(k)\,.
\end{equation}
Furthermore, due to ${\mathbf X}_t={\rm diag}(\mathbf{x}_t)$ with $\mathbf{x}_t=(x_{t,1},\ldots,x_{t,d})^{{{\mathrm{\scriptscriptstyle \top} }}}$, we let
$$ {\mathbf G}_k={\rm diag}({\mathbf g}_k)=\frac{1}{n-k}\sum_{t=k+1}^{n}\mathbb{E}[\{{\mathbf X}_t-\mathbb{E}(\bar{{\mathbf X}})\}\{\xi_{t-k}-\mathbb{E}(\bar{\xi})\}] \,,~~~~k\in[K]\,,$$
where $\bar{{\mathbf X}} = n^{-1}\sum_{t = 1}^{n}{\mathbf X}_t$. It follows from \eqref{eq:abm} and Condition \ref{cd:ra} that
\begin{equation*}
{\mathbf M}_1 = {\mathbf A} \bigg(\sum_{k = 1}^K{\mathbf G}_k {\mathbf B}^{{\mathrm{\scriptscriptstyle \top} }} {\mathbf B} {\mathbf G}_k\bigg) {\mathbf A}^{{{\mathrm{\scriptscriptstyle \top} }}} ~~ \text{and} ~~ {\mathbf M}_2 = {\mathbf B} \bigg(\sum_{k = 1}^K{\mathbf G}_k {\mathbf A}^{{\mathrm{\scriptscriptstyle \top} }} {\mathbf A} {\mathbf G}_k\bigg) {\mathbf B}^{{{\mathrm{\scriptscriptstyle \top} }}}\,.
\end{equation*}
Let ${\mathbf G}=\sum_{k=1}^K {\mathbf g}_{k}{\mathbf g}_{k}^{{\mathrm{\scriptscriptstyle \top} }}$. Proposition \ref{pro:m1-rank-con} shows that $d_1$ and $d_2$ can be identified, repsectively, by ${\rm rank}({\mathbf M}_1)$ and ${\rm rank}({\mathbf M}_2)$.
\begin{proposition}\label{pro:m1-rank-con}
Let Condition \ref{cd:ra} hold and all the main diagonal elements of ${{\mathbf G}}$
are non-zero. The following two assertions hold.
\begin{enumerate}[(i)]
\item If $\max\{ \mathcal{R}({{\mathbf G}}) + d_2,\mathcal{R}({\mathbf B}^{\mathrm{\scriptscriptstyle \top} } {\mathbf B}) + {\rm rank}({{\mathbf G}}) \}>d$, then ${\rm rank}({\mathbf M}_1)=d_1$.
\item If $\max\{ \mathcal{R}({{\mathbf G}}) + d_1,\mathcal{R}({\mathbf A}^{\mathrm{\scriptscriptstyle \top} }{\mathbf A}) + {\rm rank}({{\mathbf G}}) \}>d$, then ${\rm rank}({\mathbf M}_2)=d_2$.
\end{enumerate}
\end{proposition}
The conditions required in Proposition \ref{pro:m1-rank-con} are mild. Notice that all the main diagonal elements of ${\mathbf G}$ are positive if all the components of some ${\mathbf g}_k$ are non-zero with $k\in [K]$, which implies $\mathcal{R}({\mathbf G})\geq1$. For the scenario $d_1=d_2=d$, Proposition \ref{pro:m1-rank-con} holds automatically. When $\min(d_1, d_2) < d$, we suppose that $d_2 \le d_1$ without loss of generality. For the scenario $d_2 <d_1 = d$, we only need to identify $d_2$. Proposition \ref{pro:m1-rank-con}(ii) holds automatically in this scenario, which implies $d_2$ could be identified trivially. For the scenario $\max(d_1, d_2) < d$, by Condition \ref{cd:ra}(i), we know ${\mathbf a}_\ell\neq {\mathbf 0}$ and ${\mathbf b}_\ell\neq {\mathbf 0}$ for each $\ell\in[d]$, which implies $\mathcal{R}({\mathbf A}^{{\mathrm{\scriptscriptstyle \top} }}{\mathbf A})\geq1$ and $\mathcal{R}({\mathbf B}^{{\mathrm{\scriptscriptstyle \top} }}{\mathbf B})\geq1$. Proposition \ref{pro:rankwith-xt} proposes some sufficient conditions such that $ {\rm rank}({{\mathbf G}})=d$, which make Proposition \ref{pro:m1-rank-con} hold automatically. Define $$\boldsymbol{\Sigma}_{\mathbf{x}}(k) = \frac{1}{n-k}\sum_{t=k+1}^{n} \mathbb{E}[ \{\mathbf{x}_{t}-\mathbb{E}(\bar{\mathbf{x}})\} \{\mathbf{x}_{t-k}-\mathbb{E}(\bar{\mathbf{x}})\}^{{\mathrm{\scriptscriptstyle \top} }} ]\,,~~~~k\in[K]\,,$$ where $\bar{\mathbf{x}} =n^{-1}\sum_{t=1}^{n}\mathbf{x}_t$. Write $\xi_{t} = \boldsymbol{\omega}^{{\mathrm{\scriptscriptstyle \top} }}\vec{{\mathbf Y}}_{t} $ and $\boldsymbol{\Sigma}_{\mathbf{x},K}=\{\vec{\boldsymbol{\Sigma}}_{\mathbf{x}} (1) , \ldots, \vec{\boldsymbol{\Sigma}}_{\mathbf{x}} (K) \} \in \mathbb{R}^{d^2 \times K}$.
\begin{proposition}\label{pro:rankwith-xt}
Assume that $\mathbb{E}(\mathbf{x}_t \otimes \vec{\boldsymbol{\varepsilon}}_{t-k} ) ={\mathbf 0}$ for any $k\in[K]$ with $K\ge d^2$. If $\boldsymbol{\omega}^{{\mathrm{\scriptscriptstyle \top} }} ({\mathbf B} \odot {\mathbf A}) \ne {\mathbf 0}$ and ${\rm rank} (\boldsymbol{\Sigma}_{\mathbf{x},K}) =d^2 $, then ${\rm rank}({\mathbf G}) =d$.
\end{proposition}
Due to ${\mathbf a}_\ell\neq {\mathbf 0}$ and ${\mathbf b}_\ell\neq {\mathbf 0}$ for each $\ell\in[d]$,
the requirement $\boldsymbol{\omega}^{{\mathrm{\scriptscriptstyle \top} }} ({\mathbf B} \odot {\mathbf A}) \ne {\mathbf 0}$ is generally mild and can be satisfied by appropriately choosing a non-zero vector $\boldsymbol{\omega}$. If $\mathbf{x}_t$ satisfies ${\rm rank} (\boldsymbol{\Sigma}_{\mathbf{x},K}) =d^2 $, and $\mathbf{x}_t$ and $ \vec{\boldsymbol{\varepsilon}}_{t-k}$ are uncorrelated for $k\in[K]$,
Proposition \ref{pro:rankwith-xt} shows that ${\rm rank}({\mathbf G}) =d$. Combining with Proposition \ref{pro:m1-rank-con},
it is reasonable to assume Condition \ref{cd:M1M2eigenvalues}, which ensures ${\cal M}({\mathbf M}_1) = {\cal M}({\mathbf A})$ and
${\cal M}({\mathbf M}_2) = {\cal M}({\mathbf B}),$ i.e., the information on the loadings $\{{\mathbf a}_\ell\}_{\ell=1}^d$ and $\{{\mathbf b}_\ell\}_{\ell=1}^d$ is, respectively, kept in ${\mathbf M}_1$ and ${\mathbf M}_2$.
\begin{cd}\label{cd:M1M2eigenvalues}
${\rm rank}({\mathbf M}_1)=d_1$ and ${\rm rank}({\mathbf M}_2)=d_2$.
\end{cd}
Now perform the spectral decomposition for ${\mathbf M}_1$ and ${\mathbf M}_2$:
\begin{align}\label{eq:definition of PQ}
{\mathbf M}_1 = {\mathbf P}{\mathbf D}_1{\mathbf P}^{{{\mathrm{\scriptscriptstyle \top} }}}~~\textrm{and}~~
{\mathbf M}_2 = {\mathbf Q}{\mathbf D}_2{\mathbf Q}^{{{\mathrm{\scriptscriptstyle \top} }}}\,,
\end{align}
where ${\mathbf D}_1$ and ${\mathbf D}_2$ are, respectively, $d_1\times d_1$ and $d_2 \times d_2$ full-ranked diagonal matrices, ${\mathbf P}^{\mathrm{\scriptscriptstyle \top} } {\mathbf P}={\mathbf I}_{d_1}$ and ${\mathbf Q}^{\mathrm{\scriptscriptstyle \top} } {\mathbf Q} = {\mathbf I}_{d_2}$.
As ${\cal M}({\mathbf P})={\cal M}({\mathbf M}_1)= {\cal M}({\mathbf A})$ and
${\cal M}({\mathbf Q})={\cal M}({\mathbf M}_2)={\cal M}({\mathbf B})$, then
\begin{align}\label{eq:FA}
{\mathbf A} =
{\mathbf P} {\mathbf U} ~~{\rm and}~~
{\mathbf B}= {\mathbf Q} {\mathbf V}\,,
\end{align}
where ${\mathbf U}$ and ${\mathbf V}$ are, respectively,
${d_1\times d}$ and $d_2 \times d$ matrices with unit column vectors. Since ${\mathbf P}$ and ${\mathbf Q}$ are determined by
the spectral decomposition (\ref{eq:definition of PQ}), we only need to identify $({\mathbf U}, {\mathbf V})$ in order to identify
$({\mathbf A}, {\mathbf B})$. When $d=1$, we may take ${\mathbf a}_1 = {\mathbf A} = {\mathbf P}$ and ${\mathbf b}_1 = {\mathbf B} = {\mathbf Q}$. Therefore only the non-trivial case with $d\ge 2$ will be considered in the sequel.
Define a $d_1 \times d_2$ process ${\mathbf Z}_t = {\mathbf P}^{\mathrm{\scriptscriptstyle \top} } {\mathbf Y}_t {\mathbf Q}$. It follows from \eqref{eq:abm} and \eqref{eq:FA} that
\begin{align}\label{eq:zt}
{\mathbf Z}_t={\mathbf U} {\mathbf X}_t {\mathbf V}^{{{\mathrm{\scriptscriptstyle \top} }}} +\boldsymbol{\Delta}_t\,,~~~~t\geq1\,,
\end{align}
where $\boldsymbol{\Delta}_t={\mathbf P}^{{{\mathrm{\scriptscriptstyle \top} }}} \boldsymbol{\varepsilon}_t{\mathbf Q}$ is a matrix white noise. This is a reduced form of the CP-factor model \eqref{eq:abm} for the matrix time series ${\mathbf Y}_t$.
We will identify $({\mathbf U},{\mathbf V})$ based on this reduced model.
\subsection{
A vector factor model}
\label{sec:identifyUV}
Recall ${\mathbf X}_t={\rm diag}(\mathbf{x}_t)$ with $\mathbf{x}_t =(x_{t,1},\ldots, x_{t,d})^{\mathrm{\scriptscriptstyle \top} }$. It follows from \eqref{eq:zt} that
\begin{equation}
\label{eq:facm}
\vec{\mathbf Z}_t = ({\mathbf V}\odot {\mathbf U})\mathbf{x}_t + \vec\boldsymbol{\Delta}_t\,,~~~~t\geq1\,.
\end{equation}
This is the standard factor model for vector time series
considered by \cite{lam2012factor} and \cite{Chang2015}. Note that ${\mathbf B} \odot {\mathbf A} = ({\mathbf Q} \otimes {\mathbf P})({\mathbf V} \odot {\mathbf U})$, where ${\mathbf P}\in\mathbb{R}^{p\times d_1}$ and ${\mathbf Q}\in\mathbb{R}^{q\times d_2}$ with ${\rm rank}({\mathbf P})=d_1$ and ${\rm rank}({\mathbf Q})=d_2$. By Condition \ref{cd:ra}(i), we know the dimension of the factor loading space ${\cal M}({\mathbf V}\odot{\mathbf U})$ in \eqref{eq:facm} is $d$, as ${\rm rank}({\mathbf V}\odot {\mathbf U}) = {\rm rank}({\mathbf B}\odot{\mathbf A}) =d$.
Using the techniques developed in \cite{Chang2015}, we can identify $d$ and
${\cal M}({\mathbf V}\odot{\mathbf U})$ uniquely based on an eigenanalysis. More precisely, we can find a $(d_1 d_2)\times d$ matrix ${\mathbf W}$, with
${\mathbf W}^{\mathrm{\scriptscriptstyle \top} } {\mathbf W} ={\mathbf I}_d$, such that
\begin{equation}\label{eq:iduv}
{\mathbf V}\odot {\mathbf U} \equiv ({\mathbf v}_1 \otimes {\mathbf u}_1, \ldots, {\mathbf v}_d\otimes {\mathbf u}_d) = {\mathbf W}\boldsymbol{\Theta}\,,
\end{equation}
where $\boldsymbol{\Theta}$ is an unknown $d\times d$ invertible matrix with unit column vectors. Since the $d$ columns of ${\mathbf W}$ are the orthogonal basis of ${\cal M}({\mathbf V}\odot{\mathbf U})$, we can select ${\mathbf W}$ in \eqref{eq:iduv} as an arbitrary $(d_1d_2)\times d$ matrix such that ${\mathbf W}^{{\mathrm{\scriptscriptstyle \top} }}{\mathbf W}={\mathbf I}_d$ and $\mathcal{M}({\mathbf W})=\mathcal{M}({\mathbf U}\odot{\mathbf V})$. In \eqref{eq:iduv}, different selections of ${\mathbf W}$ will lead to different $\boldsymbol{\Theta}$. As we will show in Section \ref{sec:theta-iden}, for any given ${\mathbf W}$, the associated rotation matrix $\boldsymbol{\Theta}$ can be uniquely identified up to the column reflection and permutation indeterminacy. Write
\begin{equation}
\label{eq:c=wtheta}
{\mathbf C}\equiv (\vec{\mathbf C}_1, \ldots, \vec{\mathbf C}_d) = {\mathbf W}\boldsymbol{\Theta} \equiv (\vec {\mathbf W}_1, \ldots, \vec{\mathbf W}_d)\boldsymbol{\Theta}\,,
\end{equation}
where ${\mathbf C}_\ell$ and ${\mathbf W}_\ell$ are $d_1\times d_2$ matrices.
We put the columns of both ${\mathbf C}$ and ${\mathbf W}$ in the form of vectorized $d_1\times d_2$ matrices for some technical convenience which will be obvious soon. It follows from
(\ref{eq:iduv}) and \eqref{eq:c=wtheta} that $\vec{\mathbf C}_\ell = {\mathbf v}_\ell\otimes {\mathbf u}_\ell $, which implies ${\mathbf u}_\ell {\mathbf v}_\ell^{\mathrm{\scriptscriptstyle \top} } = {\mathbf C}_\ell$.
Given ${\mathbf W}$ and its associated rotation matrix $\boldsymbol{\Theta}$, the $(d_1d_2)\times d$ matrix ${\mathbf C}$ specified in \eqref{eq:c=wtheta} is uniquely identified, which can be used to identify $({\mathbf U},{\mathbf V})$. See Proposition \ref{pro:rank-1} for details.
\begin{proposition}
\label{pro:rank-1}
Let Conditions \ref{cd:ra} and
\ref{cd:M1M2eigenvalues} hold. Then
matrices
${\mathbf C}_1, \ldots, {\mathbf C}_d$ specified in \eqref{eq:c=wtheta} are all of rank 1 with the nonzero singular value equal to 1, and $({\mathbf u}_\ell, {\mathbf v}_\ell)$ are the unit singular vectors of ${\mathbf C}_\ell$ for each $\ell \in [d]$.
\end{proposition}
\subsection{A non-orthogonal joint diagonalization}\label{sec:theta-iden}
For given ${\mathbf W}$ in \eqref{eq:iduv}, Proposition \ref{pro:rank-1} implies that the task of identifying $({\mathbf U},{\mathbf V})$ boils down to identifying $\boldsymbol{\Theta}$ specified in \eqref{eq:iduv} such that ${\mathbf C}_1,\ldots,{\mathbf C}_d$ defined in \eqref{eq:c=wtheta} satisfying ${\rm rank}({\mathbf C}_\ell)=1$ for each $\ell \in [d]$.
By \eqref{eq:c=wtheta}, it holds that
\begin{equation}
\label{eq:bi-representation}
{\mathbf C}_\ell = \sum_{i=1}^d \theta_{i,\ell} {\mathbf W}_i ~~
{\rm and} ~~
{\mathbf W}_\ell = \sum_{i=1}^d \theta^{i,\ell} {\mathbf C}_i\,, ~~~~ \ell\in[d]\,,
\end{equation}
where $\theta_{i,j}$ and $\theta^{i,j}$ denote, respectively, the $(i,j)$-th elements of $\boldsymbol{\Theta}$ and
$\boldsymbol{\Theta}^{-1}$.
For any two matrices ${\mathbf D}=(d_{i,j})$ and ${\mathbf F}=(f_{i,j})$ of the same
size, define $\boldsymbol{\Psi}({\mathbf D}, {\mathbf F})$ to be a 4-way tensor with the
$(i,j,k,\ell)$-th element
$
d_{i,k}f_{j,\ell} + d_{j,\ell}f_{i,k} - d_{i,\ell} f_{j,k} - d_{j,k}f_{i,\ell}.
$
By Theorem 2.1 of \cite{de2006link},
for any matrix ${\mathbf D}\ne \bf0$, ${\rm rank}({\mathbf D})=1$ if and only if $\boldsymbol{\Psi}({\mathbf D}, {\mathbf D})= \bf0$.
Hence, by \eqref{eq:bi-representation}, for given ${\mathbf W}$ in \eqref{eq:iduv}, we know $\boldsymbol{\Theta} = (\theta_{i,j})$ is the solution of
\begin{equation}\label{eq: quadratic solution}
{\bf0} = \boldsymbol{\Psi}({\mathbf C}_\ell, {\mathbf C}_\ell) = \sum_{i,j=1}^d \theta_{i,\ell} \theta_{j, \ell}
\boldsymbol{\Psi}({\mathbf W}_i, {\mathbf W}_j)\,, ~~~~ \ell\in[d]\,.
\end{equation}
This is a set of quadratic equations.
Consider a $(d_1^2d_2^2)\times d(d+1)/2$
matrix
\begin{align}\label{eq:omega-D}
\boldsymbol{\Omega}=\big(\vec \boldsymbol{\Psi}({\mathbf W}_1, {\mathbf W}_1), \ldots, \vec \boldsymbol{\Psi}({\mathbf W}_1, {\mathbf W}_d),
\vec\boldsymbol{\Psi}({\mathbf W}_2, {\mathbf W}_2), \ldots, \vec\boldsymbol{\Psi}({\mathbf W}_d, {\mathbf W}_d)\big)\,.
\end{align}
Proposition \ref{pro:rankomega} is instrumental in solving those quadratic equations.
\begin{proposition}\label{pro:rankomega}
Let $d\ge 2$ and Condition \ref{cd:M1M2eigenvalues} hold. The following three assertions hold.
\begin{enumerate}[(i)]
\item ${\rm rank}(\boldsymbol{\Omega}) \le d(d-1)/2$.
\item ${\rm rank}(\boldsymbol{\Omega})=d(d-1)/2$ if and only if
the ${d(d-1)/2}$ vectors $\vec \boldsymbol{\Psi}({\mathbf C}_1, {\mathbf C}_2), \ldots, \vec \boldsymbol{\Psi}({\mathbf C}_1, {\mathbf C}_d),$
$\vec\boldsymbol{\Psi}({\mathbf C}_2, {\mathbf C}_3), \ldots, \vec\boldsymbol{\Psi}({\mathbf C}_{d-1}, {\mathbf C}_d)$ are linearly independent.
\item Let ${\rm ker}(\boldsymbol{\Omega}) = \{{\mathbf h} \in \mathbb{R}^{d(d+1)/2} : \boldsymbol{\Omega} {\mathbf h} = \mathbf{0}\}$. Then ${\rm dim}\{{\rm ker}(\boldsymbol{\Omega})\} = d$ if and only if ${\rm rank}(\boldsymbol{\Omega})=d(d-1)/2$.
\end{enumerate}
\end{proposition}
Now assume ${\rm rank}(\boldsymbol{\Omega}) ={d(d-1)/2}$. Let
${\mathbf h}_m =
(h_{1,1}^m, \ldots, h_{1,d}^m, h_{2,2}^m, \ldots, h_{d,d}^m)^{\mathrm{\scriptscriptstyle \top} }$, $m\in [d]$,
be a set of basis vectors of ${\rm ker}(\boldsymbol{\Omega})$.
Recall $\boldsymbol{\Psi}({\mathbf C}_{\ell}, {\mathbf C}_{\ell})=\bf0$ for any $\ell\in[d]$. By (\ref{eq:bi-representation}), it holds that
\begin{align*}
\bf0 &= \sum_{1\le i\le j\le d} h^m_{i,j} \boldsymbol{\Psi}({\mathbf W}_i, {\mathbf W}_j)
= \sum_{1\le i\le j\le d} h^m_{i,j} \sum_{k,\ell=1}^d \theta^{k,i} \theta^{\ell,j}\boldsymbol{\Psi}({\mathbf C}_k, {\mathbf C}_\ell)\\
&= \sum_{1\le k<\ell \le d} \boldsymbol{\Psi}({\mathbf C}_k, {\mathbf C}_\ell) \sum_{1\le i\le j\le d}
(\theta^{k,i} \theta^{\ell,j} + \theta^{k,j} \theta^{\ell,i})h^m_{i,j}\,.
\end{align*}
By Proposition \ref{pro:rankomega}(ii), we have
\begin{equation}
\label{eq:off-diagonal}
\sum_{1\le i\le j\le d}
(\theta^{k,i} \theta^{\ell,j} + \theta^{k,j} \theta^{\ell,i}) h^m_{i,j} =0 ~~\mbox{for all}~~ 1\le k < \ell \le d\,.
\end{equation}
Let ${\mathbf H}_{m}$ be a $d\times d$ matrix with the $(i,i)$-th element being $h_{i,i}^m$ for any $i$, and the $(i,j)$-th and $(j,i)$-th elements being $h_{i,j}^m/2$ for any $i<j$. Based on \eqref{eq:off-diagonal}, we know $\boldsymbol{\Gamma}_m \equiv \boldsymbol{\Theta}^{-1} {\mathbf H}_m (\boldsymbol{\Theta}^{-1})^{\mathrm{\scriptscriptstyle \top} }$ is a diagonal matrix, i.e.,
we can find $\boldsymbol{\Theta}^{-1}$ which diagonalizes jointly ${\mathbf H}_m = \boldsymbol{\Theta} \boldsymbol{\Gamma}_m \boldsymbol{\Theta}^{\mathrm{\scriptscriptstyle \top} }$ for each $m\in[d]$.
It is also clear from \eqref{eq:off-diagonal} that
the diagonal property is independent of the norms
of row vectors of $\boldsymbol{\Theta}^{-1}$. Hence all the columns of $\boldsymbol{\Theta}$ can be set as unit vectors. Proposition \ref{pro:theta-unique} shows that $\boldsymbol{\Theta}$ is invariant with respect to the choice of the basis vectors
for ${\rm ker}(\boldsymbol{\Omega})$.
The available algorithms for this joint diagonalization
include the joint approximate diagonalization of \cite{pham2001blind}, the fast Frobenius diagonalization of \cite{ziehe2004fast}, and the quadratic diagonalization of \cite{vollgraf2006quadratic}.
\begin{proposition}\label{pro:theta-unique}
Let Condition \ref{cd:M1M2eigenvalues} hold. For a given $(d_1d_2) \times d$ matrix ${\mathbf W}$ in \eqref{eq:iduv} such that ${\mathbf W}^{{\mathrm{\scriptscriptstyle \top} }}{\mathbf W}={\mathbf I}_d$ and $\mathcal{M}({\mathbf W})=\mathcal{M}({\mathbf U}\odot{\mathbf V})$, if $\boldsymbol{\Omega}$ defined in \eqref{eq:omega-D} satisfies ${\rm rank}(\boldsymbol{\Omega})=d(d-1)/2$, then $\boldsymbol{\Theta}$ in \eqref{eq:iduv} can be uniquely identified by the non-orthogonal joint diagonalization ${\mathbf H}_m =\boldsymbol{\Theta} \boldsymbol{\Gamma}_m \boldsymbol{\Theta}^{\mathrm{\scriptscriptstyle \top} }$, $m\in [d]$,
up to the column reflection and permutation indeterminacy, and $\boldsymbol{\Theta}$ is invariant with respect to the choice of the basis vectors of
${\rm ker}(\mathbf{\Omega})$.
\end{proposition}
Proposition
\ref{pro:Psitouniqueness} provides a sufficient condition under which ${\rm rank}(\boldsymbol{\Omega})=d(d-1)/2$. Such sufficient condition holds automatically when $d_1=d_2=d$, as then $\mathcal{R}({\mathbf A}) = \mathcal{R}({\mathbf B})=d$. When $d_1\neq d_2$, we assume $d_2<d_1$ without loss of generality. For the scenario $d_2<d_1=d$, since $\mathcal{R}({\mathbf A})=d$, Proposition \ref{pro:Psitouniqueness} indicates that ${\rm rank}(\boldsymbol{\Omega})= d(d-1)/2$ if $\mathcal{R}({\mathbf B}) \ge 2$. Actually, the requirement $\mathcal{R}({\mathbf B}) \ge 2$ is necessary for the identification of $({\mathbf A}, {\mathbf B})$ when $d_2<d_1=d$. Recall $\vec {\mathbf Y}_{t}=({\mathbf B}\odot {\mathbf A})\mathbf{x}_{t} + \vec \boldsymbol{\varepsilon}_{t}$. If $\mathcal{R}({\mathbf B})=1$, since $|{\mathbf b}_{\ell}|_2=1$ for each $\ell\in[d]$, we can assume ${\mathbf b}_{2}={\mathbf b}_{1}$ without loss of generality.
Let $\tilde{{\mathbf B}}={\mathbf B}$ and $\tilde{{\mathbf A}}=(\tilde{{\mathbf a}}_1,\ldots, \tilde{{\mathbf a}}_d)$, where $\tilde{{\mathbf a}}_{\ell}={\mathbf a}_{\ell}$ for any $\ell \ge 2$, and $\tilde{{\mathbf a}}_1=c_1{\mathbf a}_1+c_2{\mathbf a}_2$ for some nonzero constants $c_1, c_2$ such that $|\tilde{{\mathbf a}}_1|_2=1$.
Select $\boldsymbol{\Xi}=(\xi_{i,j})$ with $\xi_{1,1}=c_1$, $\xi_{2,1}=c_2$, $\xi_{i,i}=1$ for any $ 2 \le i \le d$, and $\xi_{i,j}=0$ otherwise. Then ${\mathbf Y}_t$ can be also formulated by another matrix CP-factor model $ \vec{\mathbf Y}_{t} = (\tilde{{\mathbf B}} \odot \tilde{{\mathbf A}}) \boldsymbol{\Xi}^{-1}\mathbf{x}_{t} + \vec \boldsymbol{\varepsilon}_{t}$.
\begin{proposition}\label{pro:Psitouniqueness}
Let $d\ge 2$, and Conditions \ref{cd:ra} and
\ref{cd:M1M2eigenvalues} hold.
Then ${\rm rank}(\boldsymbol{\Omega})= d(d-1)/2$ provided that
$\mathcal{R}({\mathbf A}) + d_2 \ge d + 2$ and $\mathcal{R}({\mathbf B}) + d_1 \ge d + 2$. \end{proposition}
By Propositions \ref{pro:rank-1} and \ref{pro:theta-unique}, if ${\rm rank}(\boldsymbol{\Omega})=d(d-1)/2$, then ${\mathbf U}$ and ${\mathbf V}$ specified in \eqref{eq:FA} can be uniquely defined up to the column reflection and permutation indeterminacy, which implies ${\mathbf A}$ and ${\mathbf B}$ can be uniquely defined up to the column reflection and permutation indeterminacy. For $d\geq2$, Proposition
\ref{pro:Omega-impossible} shows that the requirement ${\rm rank}(\boldsymbol{\Omega})=d(d-1)/2$ is necessary for identifying $({\mathbf A},{\mathbf B})$, and it is impossible to obtain the consistent estimators for $({\mathbf A}, {\mathbf B})$ without such requirement.
\begin{proposition}
\label{pro:Omega-impossible}
Let $d\ge 2$. Consider the following parameter space for the matrix CP-factor model \eqref{eq:abm}:
\begin{align*}
&\mathcal{U} = \big\{({\mathbf A},{\mathbf B}): {\mathbf A}=({\mathbf a}_1,\ldots, {\mathbf a}_d)~\textrm{and}~{\mathbf B}=({\mathbf b}_1,\ldots, {\mathbf b}_d)~\textrm{with}~|{\mathbf a}_\ell|_2 =1= |{\mathbf b}_\ell|_2 \\
&~~~~~~~~~~~~~~~~~~~~~~~~\, \textrm{for each}~\ell \in [d], ~\textrm{and}~ {\rm rank}(\boldsymbol{\Omega}) < d(d-1)/2 ~\textrm{with}~\boldsymbol{\Omega}~\textrm{defined as}~\eqref{eq:omega-D} \big\}\,.
\end{align*}
Write $\mathcal{G}=\{(\breve{{\mathbf A}}, \breve{{\mathbf B}}):\breve{{\mathbf A}}=(\breve{{\mathbf a}}_1, \ldots, \breve{{\mathbf a}}_d) \in \mathbb{R}^{p\times d}, \breve{{\mathbf B}}=(\breve{{\mathbf b}}_1, \ldots,\breve{{\mathbf b}}_d) \in \mathbb{R}^{q\times d} \}$ for the class of all measurable estimators of $({\mathbf A},{\mathbf B})$ based on the data $\{{\mathbf Y}_t\}_{t = 1}^n$. Under Conditions \ref{cd:ra} and
\ref{cd:M1M2eigenvalues}, it holds that
\begin{equation*}
\inf_{(\breve{{\mathbf A}},\breve{{\mathbf B}}) \in \mathcal{G}} \sup_{({\mathbf A},{\mathbf B}) \in
\mathcal{U}} \mathbb{P} \bigg[ \max \{\mathscr{D}(\breve{{\mathbf A}},{\mathbf A}), \mathscr{D}(\breve{{\mathbf B}},{\mathbf B})\} \ge \frac{1}{8} \bigg] \ge \frac{1}{2} \,,
\end{equation*}
where $\mathscr{D}(\breve{{\mathbf A}},{\mathbf A}) = \max_{\ell \in[d]} | \breve{{\mathbf a}}_{\ell} - {\mathbf a}_{\ell}|_2$ and $\mathscr{D}(\breve{{\mathbf B}},{\mathbf B}) = \max_{\ell \in[d]} | \breve{{\mathbf b}}_{\ell} - {\mathbf b}_{\ell}|_2$.
\end{proposition}
\section{Estimation}\label{sec:estimation}
Based on Section \ref{sec:identify-ab}, we can estimate $d$ and
$({\mathbf a}_{\ell}, {\mathbf b}_{\ell})$ for $\ell \in[d]$ via the following five steps:
\begin{enumerate}[{\it Step 1.}]
\item Based on \eqref{eq:definition of PQ}, we can obtain the estimates for $d_1, \, d_2$, ${\mathbf P}$ and ${\mathbf Q}$, denoted by $\hat{d}_1$, $\hat{d_2}$, $\hat{{\mathbf P}}$ and ${\hat{{\mathbf Q}}}$, respectively.
\item Based on \eqref{eq:facm}, we can obtain the estimates for $d$ and ${\mathbf W}$ (the orthogonal basis of $\mathcal{M}({\mathbf V}\odot{\mathbf U})$) with replacing ${\mathbf Z}_t$ by $\hat{{\mathbf Z}}_t = \hat{{\mathbf P}}^{{\mathrm{\scriptscriptstyle \top} }}{\mathbf Y}_{t}\hat{{\mathbf Q}} $. Denote by $\hat{d}$ and $\hat{{\mathbf W}} $ the associated estimators.
\item With replacing ${\mathbf W}$ involved in \eqref{eq:iduv} by $\hat{{\mathbf W}}$, we can use the joint diagonalization algorithm mentioned in Section \ref{sec:theta-iden} to obtain $\hat{\boldsymbol{\Theta}}$, the estimate of $\boldsymbol{\Theta}$ involved in \eqref{eq:iduv}.
\item Let $\hat{{\mathbf C}} =({\rm vec}(\hat{{\mathbf C}}_1), \ldots, {\rm vec}(\hat{{\mathbf C}}_{\hat{d}})) =\hat{{\mathbf W}} \hat{\boldsymbol{\Theta}}$. For each $\ell\in[\hat{d}]$, we select $\hat{{\mathbf u}}_{\ell}$ and $\hat{{\mathbf v}}_{\ell}$, respectively, as the unit eigenvectors corresponding to the largest eigenvalues of $\hat{{\mathbf C}}_\ell\hat{{\mathbf C}}_\ell^{{\mathrm{\scriptscriptstyle \top} }}$ and $\hat{{\mathbf C}}_\ell^{{\mathrm{\scriptscriptstyle \top} }}\hat{{\mathbf C}}_\ell$. Based on Proposition \ref{pro:rank-1}, we can estimate $({\mathbf u}_\ell,{\mathbf v}_\ell)$ by $(\hat{{\mathbf u}}_\ell,\hat{{\mathbf v}}_\ell)$ for each $\ell\in[\hat{d}]$.
\item Based on
\eqref{eq:FA}, we can estimate ${\mathbf A} $ and ${\mathbf B}$, respectively, by $\hat{{\mathbf A}}=\hat{{\mathbf P}}(\hat{{\mathbf u}}_1, \ldots, \hat{{\mathbf u}}_{\hat{d}})$ and $\hat{{\mathbf B}}=\hat{{\mathbf Q}}(\hat{{\mathbf v}}_1, \ldots, \hat{{\mathbf v}}_{\hat{d}})$.
\end{enumerate}
Steps 4 and 5 are straightforward. More details of Steps 1--3 are given, respectively, in Sections \ref{sec:phqh-est}--\ref{sec:theta-hat-est}. Especially Step 3 involves a further rotation to improve the convergence rate of the estimation. All the estimation is based on observations $\{{\mathbf Y}_t\}_{t=1}^{n}$.
\subsection{Estimating $d_1, \, d_2, \, {\mathbf P}$ and ${\mathbf Q}$}\label{sec:phqh-est}
Let $\xi_t$ be a prescribed linear combination of ${\mathbf Y}_t$ (e.g. the first principal component of
$\vec {\mathbf Y}_t$), and $K> 1$ be a prescribed integer. Based on \eqref{eq:M1M2}, we put
\begin{align}
\label{eq:m1h-m2h}
\hat{{\mathbf M}}_1=&~\sum_{k=1}^{K}T_{\delta_1}\{\hat{\boldsymbol{\Sigma}}_{{\mathbf Y},\xi}(k)\}T_{\delta_1}
\{\hat{\boldsymbol{\Sigma}}_{{\mathbf Y},\xi}(k)^{{{\mathrm{\scriptscriptstyle \top} }}}\} \,, \notag\\ \hat{{\mathbf M}}_2=&~\sum_{k=1}^{K}T_{\delta_1}\{\hat{\boldsymbol{\Sigma}}_{{\mathbf Y},\xi}(k)^{{{\mathrm{\scriptscriptstyle \top} }}}\}T_{\delta_1}\{\hat{\boldsymbol{\Sigma}}_{{\mathbf Y},\xi}(k)\}\,,
\end{align}
where $T_{\delta_1}(\cdot)$ is a truncation operator with the threshold level $\delta_1\geq0$, i.e., $T_{\delta_1}({\mathbf S})=(s_{i,j}I(|s_{i,j}|\ge \delta_1))$ for any matrix ${\mathbf S} =(s_{i,j})$, and
\begin{equation}
\label{eq:esthatsig}
\hat{\boldsymbol{\Sigma}}_{{\mathbf Y},\xi}(k) =\frac{1}{n-k}\sum_{t=k+1}^{n}({\mathbf Y}_t
-\bar{{\mathbf Y}})(\xi_{t-k} -\bar{\xi})\,,~~~~ k\in[K]\,.
\end{equation}
We set $\delta_1>0$ in \eqref{eq:m1h-m2h} when $pq\ge n$. Note that $\hat {\mathbf M}_1$ is a $p\times p$ matrix, and $\hat {\mathbf M}_2$ is a $q\times q$ matrix.
By Condition \ref{cd:M1M2eigenvalues}, we can estimate $d_1$ and $d_2$ by the eigenvalue-ratio based method \citep{Chang2015} as follows:
\begin{equation}
\label{eq:d1h}
\hat{d}_1=\arg\min_{j\in[
p]}\frac{\lambda_{j+1}(\hat{{\mathbf M}}_1)+c_{1,n}}
{\lambda_{j}(\hat{{\mathbf M}}_1)+c_{1,n}}~~ {\rm and}
~~
\hat{d}_2=\arg\min_{j\in[q]}\frac{\lambda_{j+1}(\hat{{\mathbf M}}_2)+c_{2,n}}{\lambda_{j}(\hat{{\mathbf M}}_2)+c_{2,n}}
\end{equation}
for some $c_{1,n}, c_{2,n} \to 0^+$ as $n\rightarrow\infty$. The proposed eigenvalue-ratio based method here is an extension of that in \cite{lam2012factor}. Adding $c_{1,n}$ and $c_{2,n}$ is to avoid the technical difficulties associated with handling potential ``0/0'' cases and can lead to consistent estimates for $d_1$ and $d_2$. See Theorem \ref{thm:rank} in Section \ref{sec:asymptotics} for details. In contrast, the eigenvalue-ratio based method proposed in \cite{lam2012factor} without adding $c_{1,n}$ and $c_{2,n}$ only ensures that the numbers of factors are not underestimated, without providing consistency.
Perform the spectral decomposition for the non-negative definite matrices $\hat {\mathbf M}_1$ and $\hat {\mathbf M}_2$. Let $\hat {\mathbf P}$ be the $p\times \hat d_1$ matrix of which the columns are the $\hat d_1$ orthonormal eigenvectors of $\hat {\mathbf M}_1$ corresponding to its $\hat d_1$ largest eigenvalues,
and $\hat {\mathbf Q}$ be the $q\times \hat d_2$ matrix of which the columns are the $\hat d_2$ orthonormal eigenvectors of $\hat {\mathbf M}_2$ corresponding to its $\hat d_2$ largest eigenvalues. Now we are ready to reduce the original $p\times q$ process ${\mathbf Y}_t$ to the $\hat d_1 \times \hat d_2$ process
\begin{equation*}
\hat {\mathbf Z}_t = \hat {\mathbf P}^{\mathrm{\scriptscriptstyle \top} } {\mathbf Y}_t \hat {\mathbf Q} \,,~~~~t\geq1\,.
\end{equation*}
\subsection{Estimating $d$ and ${\mathbf W}_1, \ldots, {\mathbf W}_d$}\label{sec:wh-est}
Based on \eqref{eq:facm} and \eqref{eq:iduv}, we can reformulate \eqref{eq:facm} as $\vec {\mathbf Z}_{t} = {\mathbf W} \mathbf{x}_{t}^{*} +\vec \boldsymbol{\Delta}_{t}$ with $\mathbf{x}_{t}^{*} = \boldsymbol{\Theta}\mathbf{x}_{t}$.
Hence, we can estimate a factor loading matrix ${\mathbf W} = (\vec {\mathbf W}_1, \ldots, \vec{\mathbf W}_d)$ based on the method proposed in \cite{lam2011estimation}, \cite{lam2012factor} and \cite{Chang2015}. To do this, we put
\begin{align}\label{eq:mh}
\hat{{\mathbf M}} = \sum_{k=1}^{\tilde{K}} \hat{\boldsymbol{\Sigma}}_{\vec{{\mathbf Z}}}(k)\hat{\boldsymbol{\Sigma}}_{\vec{{\mathbf Z}}}(k)^{{{\mathrm{\scriptscriptstyle \top} }}}
\end{align}
with a prescribed integer $\tilde{K}\ge 1$ and
\begin{align}\label{eq:zh}
\hat{\boldsymbol{\Sigma}}_{\vec{{\mathbf Z}}}(k)= (\hat{{\mathbf Q}}^{{{\mathrm{\scriptscriptstyle \top} }}} \otimes \hat{{\mathbf P}}^{{{\mathrm{\scriptscriptstyle \top} }}}) T_{\delta_2}\{\hat{\boldsymbol{\Sigma}}_{\vec{{\mathbf Y}}}(k)\} (\hat{{\mathbf Q}} \otimes \hat{{\mathbf P}})\,,~~~~k\in[\tilde{K}]\,,
\end{align}
where $T_{\delta_2}(\cdot) $ is a truncation operator with the threshold level $\delta_{2} \ge 0$, and
\begin{align*}
\hat{\boldsymbol{\Sigma}}_{\vec{{\mathbf Y}}}(k)=\frac{1}{n-k}\sum_{t=k+1}^{n}(\vec{{\mathbf Y}}_{t}-\bar{\vec{{\mathbf Y}}})(\vec{{\mathbf Y}}_{t-k}-\bar{\vec{{\mathbf Y}}})^{{\mathrm{\scriptscriptstyle \top} }} \,, ~~~~k\in[\tilde{K}]\,,
\end{align*}
with $\bar{\vec{{\mathbf Y}}}= n^{-1}\sum_{t=1}^{n}\vec{{\mathbf Y}}_{t}$.
Analogous to \eqref{eq:d1h}, we can estimate $d$ as
\begin{equation}
\label{eq:dh}
\hat{d}=\bigg\{\arg\min_{j\in [\hat d_1 \hat d_2]}\frac{\lambda_{j+1}(\hat{{\mathbf M}})+c_{3,n}}{\lambda_{j}(\hat{{\mathbf M}})+c_{3,n}}\bigg\}I(\hat{d}_1\hat{d}_2 \ge 2) + I(\hat{d}_1\hat{d}_2 =1)
\end{equation}
for some $c_{3,n} \to 0^+$ as $n\to \infty$.
Furthermore we let
$\hat {\mathbf W} \equiv ( {\rm vec}(\hat {\mathbf W}_1), \ldots, {\rm vec}(\hat
{\mathbf W}_{\hat d}))$ be the $(\hat d_1 \hat d_2)\times \hat d$ matrix of which the columns are
the $\hat d$ orthonormal eigenvectors of $\hat {\mathbf M}$ corresponding to its largest $\hat d$ eigenvalues.
\begin{remark}\label{rek:two-stage}
We can also consider an alternative two-stage procedure to estimate $d$ and ${\mathbf W}$ in Step 2. Notice that $\vec{\mathbf Y}_t=({\mathbf B} \odot {\mathbf A})\mathbf{x}_t + \vec\boldsymbol{\varepsilon}_t$ for $t \ge 1$. We can firstly obtain the estimates of $d$ and $\mathbf{T}$ (the orthogonal basis of $\mathcal{M}({\mathbf B} \odot {\mathbf A})$), denoted by $\hat{d}$ and $\hat{\mathbf{T}}$, based on the method proposed in \cite{lam2011estimation}, \cite{lam2012factor} and \cite{Chang2015}. Recall ${\mathbf V} \odot {\mathbf U} = ({\mathbf Q} \otimes {\mathbf P})^{{\mathrm{\scriptscriptstyle \top} }}({\mathbf B} \odot {\mathbf A})$ and ${\mathbf W}$ is an orthogonal basis of $\mathcal{M}({\mathbf V} \odot {\mathbf U})$. Based on $(\hat{{\mathbf P}},\hat{{\mathbf Q}})$, the estimates of ${\mathbf P}$ and ${\mathbf Q}$ obtained in Step 1, we can then estimate ${\mathbf W}$ by $ (\hat{{\mathbf Q}} \otimes \hat{{\mathbf P}})^{{\mathrm{\scriptscriptstyle \top} }}\hat{\mathbf{T}}$. Figure \ref{fig: compare-W} in the supplementary material shows that, although this alternative two-stage approach yields estimation errors nearly identical to those of our proposed method, it is considerably more computationally expensive when $p$ and $q$ are large. This is because the first stage of this alternative approach requires an eigen-decomposition of a $(pq) \times (pq)$ matrix defined based on the sample auto-covariance matrices of $\{\vec{\mathbf{Y}}_t\}_{t=1}^{n}$, whereas Step 2 of our proposed method only involves an eigen-decomposition of a $(\hat{d}_1\hat{d}_2) \times (\hat{d}_1\hat{d}_2)$ matrix. This indicates that our proposed method can significantly reduce the computational complexity, especially in high-dimensional settings.
\end{remark}
\subsection{Estimating $\mathbf{\Theta}$ via joint diagonalization}\label{sec:theta-hat-est}
For 4-way tensor $\boldsymbol{\Psi}(\cdot,\cdot)$ defined in Section \ref{sec:theta-iden}, we define a $(\hat d_1^2\hat d_2^2)\times \hat{d}(\hat d+1)/2$ matrix $\hat \boldsymbol{\Omega}$ as follows:
\begin{align}\label{eq:omega-h}
\hat{\boldsymbol{\Omega}} = \big(\vec \boldsymbol{\Psi}(\hat{\mathbf W}_1, \hat{\mathbf W}_1), \ldots, \vec \boldsymbol{\Psi}(\hat{\mathbf W}_1, \hat{\mathbf W}_{ \hat{d} }),
\vec\boldsymbol{\Psi}(\hat{\mathbf W}_2, \hat{\mathbf W}_2), \ldots, \vec\boldsymbol{\Psi}(\hat{\mathbf W}_{ \hat{d} }, \hat{\mathbf W}_{ \hat{d} }) \big)\,,
\end{align}
which is an estimate of $\boldsymbol{\Omega}$ defined as in \eqref{eq:omega-D}. Let
\begin{equation}
\label{eq:init-basis}
\tilde {\mathbf h}_m=
(\tilde h_{1,1}^m, \ldots, \tilde h_{1, \hat d}^m, \tilde h_{2,2}^m, \ldots, \tilde h_{\hat d,\hat d}^m)^{\mathrm{\scriptscriptstyle \top} }\,, ~~~~ m\in [\hat{d}]\,,
\end{equation}
be the right-singular vectors of $\hat \boldsymbol{\Omega}$ corresponding to the $\hat d$ smallest singular values. Such selected $\{\tilde{{\mathbf h}}_m\}_{m=1}^{\hat{d}}$ provides the estimate for a basis of ${\rm ker}(\boldsymbol{\Omega})$. By Proposition \ref{pro:theta-unique}, an estimator for $\boldsymbol{\Theta}$ can be obtained by the joint diagonalization of $\tilde {\mathbf H}_1, \ldots, \tilde{{\mathbf H}}_{\hat{d}}$, which are constructed in the same manner as ${\mathbf H}_m$ with ${\mathbf h}_m$ replaced by $\tilde {\mathbf h}_m$. See the statement below \eqref{eq:off-diagonal}.
Though $\boldsymbol{\Theta}$ can be uniquely identified by any set of basis $\{ {\mathbf h}_m\}_{m=1}^{d}$ of ${\rm ker}(\boldsymbol{\Omega})$ (see Proposition \ref{pro:theta-unique}), the accuracy of its estimator depends on the choice of $\{ {\mathbf h}_m\}_{m=1}^{d}$ sensitively. Motivated by Proposition \ref{pro:rotation} at the end of this section, a good choice is to rotate
the basis vectors $\{\tilde{{\mathbf h}}_m\}_{m=1}^{\hat{d}}$ in \eqref{eq:init-basis} first. More specifically, let $(\hat{{\mathbf h}}_1, \ldots, \hat{{\mathbf h}}_{\hat{d}}) = (\tilde{{\mathbf h}}_1,\ldots, \tilde{{\mathbf h}}_{\hat{d}}) \hat \boldsymbol{\Pi}$ with
\begin{align*}
\hat \boldsymbol{\Pi} = \{2(\hat {\mathbf \Upsilon}_0^{\mathrm{\scriptscriptstyle \top} } \hat{\mathbf \Upsilon}_2)( \hat{\mathbf \Upsilon}_1^{\mathrm{\scriptscriptstyle \top} } \hat{\mathbf \Upsilon}_2 + \hat{\mathbf \Upsilon}_2^{\mathrm{\scriptscriptstyle \top} } \hat{\mathbf \Upsilon}_1)^{-1}( \hat{\mathbf \Upsilon}_2^{\mathrm{\scriptscriptstyle \top} } \hat{\mathbf \Upsilon}_0) \}^{-1/2}\,,
\end{align*}
where
\begin{gather*}
\hat{\mathbf \Upsilon}_0 =
\big(\textup{vec}(\tilde{\mathbf H}_1),\ldots,\textup{vec}(\tilde{\mathbf H}_{\hat{d}})\big) \, , ~~
\hat{\mathbf \Upsilon}_1 =
\big(\textup{vec}(\tilde{\mathbf H}^{-1}\tilde{\mathbf H}_1),\ldots,\textup{vec}(\tilde{\mathbf H}^{-1}\tilde{\mathbf H}_{\hat d})\big)\,,\\
\hat{\mathbf \Upsilon}_2 = \big(\textup{vec}(\tilde{\mathbf H}_1\tilde{\mathbf H}^{-1}),\ldots,\textup{vec}(\tilde{\mathbf H}_{\hat{d}}\tilde{\mathbf H}^{-1})\big)\,, ~~
\tilde {\mathbf H} = \sum_{m=1}^{\hat d} \phi_m \tilde {\mathbf H}_m
\end{gather*}
for some $\hat d$-dimensional vector $(\phi_1, \ldots, \phi_{\hat d})^{\mathrm{\scriptscriptstyle \top} } \ne {\mathbf 0}$ such that $\tilde{{\mathbf H}}$ is invertible. Section \ref{sec:simulation setting up} specifies how to select $(\phi_1, \ldots, \phi_{\hat d})^{\mathrm{\scriptscriptstyle \top} }$ in practice. Define $\hat{{\mathbf H}}_{1},\ldots, \hat{{\mathbf H}}_{\hat{d}}$ in the same manner as ${\mathbf H}_{m}$ but with replacing ${\mathbf h}_m$ by $\hat{{\mathbf h}}_m$. Utilizing the fast Frobenius diagonalization algorithm introduced by \cite{ziehe2004fast}, we can obtain the non-orthogonal joint diagonalizer $\boldsymbol{\Phi}$ for $\hat{{\mathbf H}}_1, \ldots, \hat{{\mathbf H}}_{\hat{d}}$ such that all the columns of $\boldsymbol{\Phi}^{-1}$ are unit vectors. Then, $\boldsymbol{\Theta}$ involved in \eqref{eq:iduv} can be estimated by $\hat{\boldsymbol{\Theta}} =\boldsymbol{\Phi}^{-1}$.
Now we give some illustrations on how the set of basis $\{{\mathbf h}_m\}_{m=1}^d$ of ${\rm ker}(\boldsymbol{\Omega})$ used to identify $\boldsymbol{\Theta}$ affects the convergence rate of the associated estimator of $\boldsymbol{\Theta}$ based on the non-orthogonal joint diagonalization. Recall
\[
{\mathbf H}_m= \boldsymbol{\Theta} \,{\rm diag}(\gamma_{1,m} ,
\ldots, \gamma_{d,m}) \,\boldsymbol{\Theta}^{\mathrm{\scriptscriptstyle \top} }, ~~~~ m\in[d]\,.
\]
By Theorem 3 of \cite{afsari2008sensitivity}, the convergence
rate of the estimator for $\boldsymbol{\Theta}$ based on the fast Frobenius diagonalization algorithm is bounded by
\[
\frac{\eta({{\mathbf h}_1}, \ldots, {\mathbf h}_d)}{1-\rho^2({{\mathbf h}}_1, \ldots, {\mathbf h}_d) } \times K_n(d,\boldsymbol{\Theta}) \,,
\]
where $K_{n}(d, \boldsymbol{\Theta})$ is a universal quantity only depending on $(n,d, \boldsymbol{\Theta})$, and
\begin{align}\label{eq:alpha and rho}
\rho({{\mathbf h}}_1, \ldots, {\mathbf h}_d) &= \max_{k,\ell\in[d]:\,k\ne \ell}\frac{|\sum_{m =
1}^d\gamma_{k,m}\gamma_{\ell,m}|}{(\sum_{m =
1}^d\gamma_{k,m}^2)^{1/2} (\sum_{m = 1}^d\gamma_{\ell,m}^2)^{1/2}}
\, , \notag\\
\eta({{\mathbf h}}_1, \ldots, {\mathbf h}_d) &= \max_{k,\ell\in[d]:\,k\ne\ell}\bigg(\frac{1}{\sum_{m = 1}^d\gamma_{\ell,m}^2} + \frac{1}{\sum_{m = 1}^d\gamma_{k,m}^2}\bigg)\,.
\end{align}
Ideally we should choose $\{ {\mathbf h}_m\}_{m=1}^{d}$ such that
$\rho({\mathbf h}_1, \ldots, {\mathbf h}_d)=0$ and $\eta({\mathbf h}_1, \ldots, {\mathbf h}_d)$ as
small as possible. For any given set of basis $\{{\mathbf h}_m\}_{m=1}^d$ of ${\rm ker}(\boldsymbol{\Omega})$, Proposition \ref{pro:rotation} indicates that we should replace $\{{\mathbf h}_m\}_{m=1}^d$ by its rotation $\{{\mathbf h}_m^*\}_{m=1}^d$ such that
$({\mathbf h}_1^*, \ldots, {\mathbf h}_d^{*}) = ({\mathbf h}_1, \ldots, {\mathbf h}_d) \boldsymbol{\Pi} $ with
\begin{align*}
\boldsymbol{\Pi} = \{2({\mathbf \Upsilon}_0^{\mathrm{\scriptscriptstyle \top} }{\mathbf \Upsilon}_2)({\mathbf \Upsilon}_1^{\mathrm{\scriptscriptstyle \top} }{\mathbf \Upsilon}_2 + {\mathbf \Upsilon}_2^{\mathrm{\scriptscriptstyle \top} }{\mathbf \Upsilon}_1)^{-1}({\mathbf \Upsilon}_2^{\mathrm{\scriptscriptstyle \top} }{\mathbf \Upsilon}_0) \}^{-1/2} \,,
\end{align*}
where
\begin{gather}\label{eq:bUpsilon_k}
{\mathbf \Upsilon}_0 =
\big(\textup{vec}({\mathbf H}_1),\ldots,\textup{vec}({\mathbf H}_d)\big) \, , ~~
{\mathbf \Upsilon}_1 =
\big(\textup{vec}({\mathbf H}^{-1}{\mathbf H}_1),\ldots,\textup{vec}({\mathbf H}^{-1}{\mathbf H}_d)\big)\,, \notag\\
{\mathbf \Upsilon}_2 = \big(\textup{vec}({\mathbf H}_1{\mathbf H}^{-1}),\ldots,\textup{vec}({\mathbf H}_d{\mathbf H}^{-1})\big)\,, ~~ {\mathbf H} = \sum_{m=1}^d \phi_m {\mathbf H}_m
\end{gather}
for some $d$-dimensional vector $(\phi_1, \ldots, \phi_d)^{\mathrm{\scriptscriptstyle \top} } \ne {\mathbf 0}$ such that ${\mathbf H}$ is invertible.
\begin{proposition}
\label{pro:rotation}
$\rho({\mathbf h}^*_1, \ldots, {\mathbf h}_d^*) =0$ and $\eta({\mathbf h}^*_1, \ldots, {\mathbf h}_d^*) =2$.
\end{proposition}
\section{Prediction}\label{sec:prediction}
Given observations $\{{\mathbf Y}_t\}_{t=1}^n$, we can also use the matrix CP-factor model \eqref{eq:abm} to forecast the future values ${\mathbf Y}_{n+h}$ for $h\ge 1$. More specifically, we can predict ${\mathbf Y}_{n+h}$ by recovering the latent process $\{{\mathbf X}_t\}_{t=1}^n$. Let
$\hat{{\mathbf L}} = \hat{{\mathbf B}} \odot \hat{{\mathbf A}} $ with $\hat{{\mathbf A}}\in\mathbb{R}^{p\times\hat{d}}$ and $\hat{{\mathbf B}}\in\mathbb{R}^{q\times\hat{d}}$ being, respectively, the estimates of the factor loading matrices ${\mathbf A}$ and ${\mathbf B}$ in the matrix CP-factor model \eqref{eq:abm}. If ${\rm rank}(\hat{{\mathbf L}})=\hat{d}$, we can recover ${\mathbf X}_t$ by $\hat{{\mathbf X}}_t = \text{diag}(\hat{\mathbf{x}}_t)$ with $\hat{\mathbf{x}}_t=\hat{{\mathbf L}}^{+}{\mathbf Y}_t = (\hat{x}_{t,1},\ldots,\hat{x}_{t,\hat{d}})^{\mathrm{\scriptscriptstyle \top} }$. In order to predict ${\mathbf Y}_{n+h}$, we only need to fit a $\hat{d}$-dimensional multivariate time series model for $\{\hat{\mathbf{x}}_t\}^n_{t=1}$. Then we can predict ${\mathbf Y}_{n+h}$ by $\hat{{\mathbf Y}}_{n+h} = \hat{{\mathbf A}}\tilde{\hat{{\mathbf X}}}_{n+h}\hat{{\mathbf B}}^{\mathrm{\scriptscriptstyle \top} }$ with $\tilde{\hat{{\mathbf X}}}_{n+h} = \textup{diag} (\tilde{\hat{\mathbf{x}}}_{n+h})$, where $\tilde{\hat{\mathbf{x}}}_{n+h}$ is the $h$-step ahead forecast of $\hat{\mathbf{x}}_{n+h}$ based on the fitted model for $\{\hat{\mathbf{x}}_t\}^n_{t=1}$. \cite{chang2023modelling} uses this idea to predict ${\mathbf Y}_{n+h}$ under the assumption $d_1 = d_2 = d$ based on the CP-refined estimate considered there for $({\mathbf A},{\mathbf B})$. Since the CP-refined estimate \citep{chang2023modelling} does not work if the assumption $d_1=d_2=d$ is not satisfied, we cannot select $(\hat{{\mathbf A}},\hat{{\mathbf B}})$ as the CP-refined estimate to recover ${\mathbf X}_t$ in these cases. When $d_1=d_2=d$ is not satisfied, if the factor loading matrices ${\mathbf A}$ and ${\mathbf B}$ can be uniquely identified, we can select $(\hat{{\mathbf A}},\hat{{\mathbf B}})$ as our newly proposed estimate specified in Section \ref{sec:estimation}, and use the same idea to predict ${\mathbf Y}_{n+h}$.
As shown in Proposition \ref{pro:Omega-impossible}, if ${\rm rank}(\boldsymbol{\Omega}) < d(d-1)/2$, the factor loading matrices ${\mathbf A}$ and ${\mathbf B}$ cannot be uniquely identified, which implies that we cannot recover ${\mathbf X}_t$ successfully. Hence, above mentioned strategy for predicting ${\mathbf Y}_{n+h}$ does not always work. A natural question is that whether we can propose a unified prediction procedure for ${\mathbf Y}_{n+h}$ based on the matrix CP-factor model \eqref{eq:abm} without any assumption on the relationship among $d_1$, $d_2$ and $d$. By \eqref{eq:FA} and \eqref{eq:zt}, we have
\begin{align}\label{eq:newmodel}
{\mathbf Y}_t={\mathbf P}{\mathbf Z}_t{\mathbf Q}^{{\mathrm{\scriptscriptstyle \top} }}+\underbrace{\boldsymbol{\varepsilon}_t - {\mathbf P}{\mathbf P}^{{\mathrm{\scriptscriptstyle \top} }}\boldsymbol{\varepsilon}_t{\mathbf Q}{\mathbf Q}^{{\mathrm{\scriptscriptstyle \top} }}}_{{\rm white~noise}}
\end{align}
with ${\mathbf Z}_t = {\mathbf P}^{{\mathrm{\scriptscriptstyle \top} }}{\mathbf Y}_t{\mathbf Q}$. In order to predict ${\mathbf Y}_{n+h}$, we only need to predict ${\mathbf Z}_{n+h}$. For $(\hat{{\mathbf P}},\hat{{\mathbf Q}},\hat{{\mathbf W}})$ specified in Sections \ref{sec:phqh-est} and \ref{sec:wh-est}, we define
\begin{align*}
\hat{\mathbf{x}}^*_t = \hat{{\mathbf W}} ^{{\mathrm{\scriptscriptstyle \top} }}{\rm vec}(\hat{{\mathbf P}}^{{\mathrm{\scriptscriptstyle \top} }}{\mathbf Y}_{t}\hat{{\mathbf Q}})\,,~~~~ t\in[n]\,.
\end{align*}
Proposition \ref{pro:PQW} in Section \ref{sec:asymptotics} indicates that such defined $\hat{d}$-dimensional vector $\hat{\mathbf{x}}_t^*$ provides a recovery of $\mathbf{E}_3\mathbf{x}_t^*$, where $ \mathbf{x}_t^* = \boldsymbol{\Theta}\mathbf{x}_{t} $, and $\mathbf{E}_3$ is an orthogonal matrix specified in Proposition \ref{pro:PQW}. Hence, we can fit a $\hat{d}$-dimensional vector time series model for $\{\hat{\mathbf{x}}_t^*\}^n_{t = 1}$.
Let $\tilde{\hat{\mathbf{x}}}_{n+h}^*$ be the $h$-step ahead forecast of $\hat{\mathbf{x}}_{n+h}^*$. By \eqref{eq:facm} and Proposition \ref{pro:PQW}, we know $\hat{{\mathbf W}}\tilde{\hat{\mathbf{x}}}_{n+h}^*$ provides a prediction of $({\mathbf E}_2 \otimes {\mathbf E}_1)\vec{{\mathbf Z}}_{n+h}$, where ${\mathbf E}_1$ and ${\mathbf E}_2$ are two orthogonal matrices specified in Proposition \ref{pro:PQW}. Let $\hat{{\mathbf Z}}_{n+h}$ satisfy ${\rm vec}(\hat{{\mathbf Z}}_{n+h})= \hat{{\mathbf W}} \tilde{\hat{\mathbf{x}}}_{n+h}^*$. Applying Proposition \ref{pro:PQW} again, by \eqref{eq:newmodel}, we know
$\hat{{\mathbf Y}}_{n+h} = \hat{{\mathbf P}}\hat{{\mathbf Z}}_{n+h}\hat{{\mathbf Q}}^{{\mathrm{\scriptscriptstyle \top} }}$ provides a prediction of ${\mathbf Y}_{n+h}$. This new prediction idea only depends on the calculation of three matrices $\hat{{\mathbf P}}$, $\hat{{\mathbf Q}}$ and $\hat{{\mathbf W}}$. As we have discussed in Sections \ref{sec:phqh-est} and \ref{sec:wh-est}, determining $\hat{{\mathbf P}}$, $\hat{{\mathbf Q}}$ and $\hat{{\mathbf W}}$ only involves the spectral decomposition of $\hat{{\mathbf M}}_1$, $\hat{{\mathbf M}}_2$ and $\hat{{\mathbf M}}$, respectively, which does not require any additional assumption on the relationship among $d_1$, $d_2$ and $d$. Hence, our newly proposed prediction strategy provides a unified prediction procedure for ${\mathbf Y}_{n+h}$ based on the matrix CP-factor model \eqref{eq:abm} regardless of the relationship among $d_1$, $d_2$ and $d$. When the linear dynamic structure is concerned for the latent process ${\mathbf X}_t$, our numerical studies in Section \ref{sec:sim-est} indicate that if the factor loading matrices ${\mathbf A}$ and ${\mathbf B}$ can be uniquely identified, the finite-sample performance of our newly proposed prediction method is almost identical to the prediction idea considered in \cite{chang2023modelling} with selecting $(\hat{{\mathbf A}},\hat{{\mathbf B}})$ as our proposed estimate of $({\mathbf A},{\mathbf B})$ specified in Section \ref{sec:estimation}. However, if the factor loading matrices ${\mathbf A}$ and ${\mathbf B}$ cannot be uniquely identified, our newly proposed prediction method outperforms that of \cite{chang2023modelling}.
\section{Asymptotic properties}\label{sec:asymptotics}
As we do not impose the stationarity on $\{{\mathbf Y}_t\}$, we use the concept of ``$\alpha$-mixing'' to characterize
the serial dependence of $\{{\mathbf Y}_t\}$ with the $\alpha$-mixing coefficients defined as
\begin{align}\label{eq:mixinga}
\alpha(k)=\sup_{r}\sup_{A\in \mathcal{F}_{-\infty}^{r},B\in \mathcal{F}_{r+k}^{\infty} } |\mathbb{P}(AB)-\mathbb{P}(A)\mathbb{P}(B)|\,, ~~~~ k\ge 1\,,
\end{align}
where $\mathcal{F}_{r}^{s}$ is the $\sigma$-filed generated by $\{{\mathbf Y}_t: r\le t\le s\}$. Write
\begin{align*}
\boldsymbol{\Sigma}_{\vec{{\mathbf Y}}}(k) = \frac{1}{n-k}\sum_{t=k+1}^{n} \mathbb{E}[\{\vec{{\mathbf Y}}_{t}-\mathbb{E}(\bar{\vec{{\mathbf Y}}})\}\{\vec{{\mathbf Y}}_{t-k}-\mathbb{E}(\bar{\vec{{\mathbf Y}}})\} ^{{\mathrm{\scriptscriptstyle \top} }}]\,,~~~~ k\ge 1\,,
\end{align*}
where $\bar{\vec{{\mathbf Y}}}= n^{-1}\sum_{t=1}^{n}\vec{{\mathbf Y}}_{t}$. Define ${\mathbf M} = \sum_{k=1}^{\tilde{K}} \boldsymbol{\Sigma}_{\vec{{\mathbf Z}}}(k) \boldsymbol{\Sigma}_{\vec{{\mathbf Z}}}(k)^{{\mathrm{\scriptscriptstyle \top} }} $ with $\tilde{K}$ given in \eqref{eq:mh} and
\begin{align*}
\boldsymbol{\Sigma}_{\vec{{\mathbf Z}}}(k) =\frac{1}{n-k}\sum_{t=k+1}^{n} \mathbb{E} [\{ \vec{{\mathbf Z}}_t - \mathbb{E}(\bar{\vec{ {\mathbf Z}}} ) \} \{ \vec{{\mathbf Z}}_{t-k}- \mathbb{E}(\bar{\vec{{\mathbf Z}}} )\}^{{\mathrm{\scriptscriptstyle \top} }} ]\,,~~~~k\ge1\,,
\end{align*}
where $\bar{\vec{{\mathbf Z}}}=n^{-1}\sum_{t=1}^{n}\vec {\mathbf Z}_{t}$. Following the arguments in \cite{Chang2015}, we can identify $d$ as $d={\rm rank}({\mathbf M})$, and select $\vec{\mathbf W}_1,\ldots, \vec{\mathbf W}_d$ involved in \eqref{eq:c=wtheta} as the $d$ orthonormal eigenvectors of ${\mathbf M}$ corresponding to the $d$ non-zero eigenvalues $\lambda_1({\mathbf M}) \ge \cdots \ge \lambda_d({\mathbf M}) > 0 $, i.e., $\vec{\mathbf W}_{\ell}$ is the eigenvector associated with the eigenvalue $\lambda_{\ell}({\mathbf M})$ for $\ell\in[d]$. We need the following regularity conditions in our theoretical analysis.
\begin{cd}\label{cd:bounded-value}
{\rm(i)} The nonzero singular values of ${\mathbf B} \odot {\mathbf A}$ are uniformly bounded away from zero. {\rm(ii)} The nonzero eigenvalues of ${\mathbf M}_1$, ${\mathbf M}_2$ and ${\mathbf M}$ are uniformly bounded away from zero.
\end{cd}
\begin{cd}\label{cd:tail}
{\rm (i)} There exist some universal constants $K_1>0$, $K_2>0$ and $r_1\in(0,2]$ such that
$
\mathbb{P}(|y_{i,j,t}|>x)\le K_1 \exp(-K_2x^{r_1})$ and $\mathbb{P}(|\xi_t|>x)\le K_1 \exp(-K_2x^{r_1})$
for any $x>0$, $i\in[p] $, $j\in[q]$ and $t\in[n]$. {\rm (ii)} There exist some universal constants $K_3>0$, $K_4>0$ and $r_2\in (0,1]$ such that the $\alpha$-mixing coefficients $\alpha(k)$ defined as in \eqref{eq:mixinga} satisfy
$
\alpha(k)\le K_3\exp(-K_4k^{r_2})$
for $k\ge 1$.
\end{cd}
\begin{cd}\label{cd:sgm_yxi}
{\rm (i)} There exists a universal constant $K_5>0$ such that
$ \|\boldsymbol{\Sigma}_{{\mathbf Y},\xi}(k)\|_2\le K_5 $ for any $k\in[K]$, and $ \|\boldsymbol{\Sigma}_{\vec{{\mathbf Y}}}(k)\|_2 \le K_5$ for any $k \in [\tilde{K}]$.
{\rm (ii)} Write $\boldsymbol{\Sigma}_{{\mathbf Y},\xi}(k)= (\sigma_{y,\xi,i,j}^{(k)})_{p\times q}$ and $\boldsymbol{\Sigma}_{\vec{{\mathbf Y}}} (k) = (\sigma_{i,j}^{(k)})_{pq\times pq}$. There exists a universal constant $\iota \in[0,1)$ such that
$
\sum_{j_1=1}^{q}|\sigma_{y,\xi,i_1,j_1}^{(k)}|^{\iota} \le s_1$, $\sum_{i_1=1}^{p}|\sigma_{y,\xi,i_1,j_1}^{(k)}|^{\iota} \le s_2$, $\sum_{j_2=1}^{pq}|\sigma_{i_2,j_2}^{(k)}|^{\iota} \le s_3$ and $\sum_{i_2=1}^{pq}|\sigma_{i_2,j_2}^{(k)}|^{\iota} \le s_4$
for any $i_1 \in [p]$, $j_1 \in [q]$ and $i_2,j_2 \in [pq]$, where $s_1$, $s_2$, $s_3$ and $s_4$ may, respectively, diverge together with $p$ and $q$.
\end{cd}
Condition \ref{cd:bounded-value} is used to
simplify the presentation for the results. Our technical proofs indeed allow the nonzero singular values of ${\mathbf B} \odot {\mathbf A}$, and the nonzero eigenvalues of ${\mathbf M}_1$, ${\mathbf M}_2$ and ${\mathbf M}$ decay to zero as $p$ and/or $q$ grow to infinity.
Condition \ref{cd:tail} is also used in \cite{chang2023modelling}, which is a common assumption in the literature on ultrahigh-dimensional data analysis. See \cite{chang2023modelling} for the discussion of their validity. We impose Condition \ref{cd:sgm_yxi}(i) just for simplifying the presentation. Our technical proofs indeed allow $\max_{k \in [K]} \|\boldsymbol{\Sigma}_{{\mathbf Y},\xi}(k)\|_2$ and $ \max_{k \in [K]}\|\boldsymbol{\Sigma}_{\vec{{\mathbf Y}}}(k)\|_2$ to diverge as $p$ and/or $q$ grow to infinity. Condition \ref{cd:sgm_yxi}(ii) imposes some sparsity requirement on $\boldsymbol{\Sigma}_{{\mathbf Y},\xi}(k)$ and $\boldsymbol{\Sigma}_{\vec{{\mathbf Y}}}(k)$.
Under some sparsity condition on ${\mathbf A}$ and ${\mathbf B}$, applying the technique used to derive Lemma 5 of
\cite{Chang2018}, we can show that Condition \ref{cd:sgm_yxi}(ii) holds for certain $(s_1,s_2,s_3,s_4)$. Let
\begin{equation}\label{eq:Pi1 and Pi2}
\Pi_{1,n}=(s_1s_2)^{1/2}\bigg\{\frac{\log(pq)}{n}\bigg\}^{(1-\iota)/2} ~~ \text{and} ~~ \Pi_{2,n}= (s_3s_4)^{1/2} \bigg\{\frac{\log(pq)}{n}\bigg\}^{(1-\iota)/2} \,.
\end{equation}
Theorem \ref{thm:rank} shows that the eigenvalue-ratio based estimators $\hat{d}_1$, $\hat{d}_2$ and $\hat{d}$ provide consistent estimates for $d_1$, $d_2$ and $d$, respectively.
\begin{theorem}\label{thm:rank}
Let Conditions \ref{cd:ra}--\ref{cd:sgm_yxi} hold. Select the threshold levels in \eqref{eq:m1h-m2h} and \eqref{eq:zh} as
\[
\delta_1=\breve{C}\sqrt{\frac{\log(pq)}{n}}~~\textrm{and}~~ \delta_2=\tilde{C}\sqrt{\frac{\log(pq)}{n}}
\]
for some sufficiently large constants $\breve{C}, \tilde{C}>0$. For any $(c_{1,n},c_{2,n}, c_{3,n})$ given in \eqref{eq:d1h} and \eqref{eq:dh} satisfying $\Pi_{1,n} \ll c_{1,n},c_{2,n} \ll 1$ and $\max(\Pi_{1,n}, \Pi_{2,n}) \ll c_{3,n} \ll 1$, it holds that
\[
\mathbb{P}(\hat{d_1}=d_1) \to 1\,,~\mathbb{P}(\hat{d_2}=d_2) \to 1~~\textrm{and}~~\mathbb{P}(\hat{d}=d) \to 1\]
as $n \to \infty$,
provided that $\Pi_{1,n}+\Pi_{2,n}\ll1$ and $\log(pq)\ll n^c$ for some constant $c \in (0,1)$ depending only on $r_1$ and $r_2$.
\end{theorem}
Proposition \ref{pro:PQW} states the asymptotic performance of $\hat{{\mathbf P}}$, $\hat{{\mathbf Q}}$ and $\hat{{\mathbf W}}$.
\begin{proposition}\label{pro:PQW}
Let Conditions \ref{cd:ra}--\ref{cd:sgm_yxi} hold. Select the threshold levels in \eqref{eq:m1h-m2h} and \eqref{eq:zh} as
\[
\delta_1=\breve{C}\sqrt{\frac{\log(pq)}{n}}~~\textrm{and}~~ \delta_2=\tilde{C}\sqrt{\frac{\log(pq)}{n}}
\]
for some sufficiently large constants $\breve{C}, \tilde{C}>0$. Assume that $\Pi_{1,n} + \Pi_{2,n}\ll1$ and $\log(pq)\ll n^c$ for some constant $c \in (0,1)$ depending only on $r_1$ and $r_2$. If $(\hat{d}_1,\hat{d}_2)=(d_1,d_2)$, there exist some orthogonal matrices ${\mathbf E}_1 \in \mathbb{R}^{d_1 \times d_1}$ and ${\mathbf E}_2 \in \mathbb{R}^{d_2 \times d_2}$ such that
\[
\|\hat{{\mathbf P}}{\mathbf E}_1 - {\mathbf P} \|_2=O_{\rm p}(\Pi_{1,n})=\|\hat{{\mathbf Q}}{\mathbf E}_2 - {\mathbf Q} \|_2\,.
\]
Furthermore, if $(\hat{d}_1,\hat{d}_2,\hat{d})=(d_1,d_2,d)$, there exists an orthogonal matrix ${\mathbf E}_3 \in \mathbb{R}^{d \times d}$ such that
\[
\|({\mathbf E}_2 \otimes {\mathbf E}_1)^{{{\mathrm{\scriptscriptstyle \top} }}} \hat{{\mathbf W}}{\mathbf E}_3 -{\mathbf W} \|_2=O_{\rm p}(\Pi_{1,n} + \Pi_{2,n})\,.
\]
\end{proposition}
If the nonzero eigenvalues of ${\mathbf M}$ are distinct, ${\mathbf E}_3$ will be a diagonal matrix with its diagonal elements being $1$ or $-1$. For the trivial case $d=1$, if $(\hat{d}_1,\hat{d}_2,\hat{d})=(d_1,d_2,d)$, we have ${\mathbf E}_1 =\pm 1$ and $ {\mathbf E}_2 = \pm 1$. Following the discussion below \eqref{eq:FA}, it holds in this trivial case that $|\hat{{\mathbf A}}{\mathbf E}_1 - {\mathbf A} |_2=O_{\rm p}(\Pi_{1,n})=|\hat{{\mathbf B}}{\mathbf E}_2 - {\mathbf B} |_2 $ provided that $\Pi_{1,n} \ll 1$ and $\log(pq)\ll n^c$ for some constant $c \in (0,1)$ depending only on $r_1$ and $r_2$. Note that $\mathbb{P}\{(\hat{d}_1,\hat{d}_2,\hat{d})=(d_1,d_2,d)\} \to 1$ as $n\to \infty$. Hence, in the trivial case $d=1$, $({\mathbf A},{\mathbf B})$ can be consistently estimated up to the reflection indeterminacy.
For the non-trivial case $d\ge 2$, the convergence rates of the estimation errors for $ {\mathbf A} $ and $ {\mathbf B} $ will be shown in Theorem \ref{thm:hat_a_hat_b}.
To present Theorem \ref{thm:hat_a_hat_b}, we need to introduce some notation first. For $({\mathbf E}_1,{\mathbf E}_2,{\mathbf E}_3)$ specified in Proposition \ref{pro:PQW}, let $\breve{{\mathbf W}} = ({\mathbf E}_2 \otimes {\mathbf E}_1) {\mathbf W} {\mathbf E}_3^{{\mathrm{\scriptscriptstyle \top} }}$. Define $\breve{\boldsymbol{\Omega}}$ in the same manner as $\boldsymbol{\Omega}$ given in \eqref{eq:omega-D} but with replacing ${\mathbf W}$ by $\breve{{\mathbf W}}$.
Following the discussions of Propositions \ref{pro:theta-unique}--\ref{pro:Omega-impossible}, the requirement ${\rm rank}(\boldsymbol{\Omega}) =d(d-1)/2$ is crucial for the identification of $({\mathbf A}, {\mathbf B})$ when $d\ge 2$.
As shown in Section \ref{sec:sub-theta} in the supplementary material for the proof of Lemma \ref{lemma:theta}, we know ${\rm rank}(\breve{\boldsymbol{\Omega}}) ={\rm rank}(\boldsymbol{\Omega})$. By Proposition \ref{pro:rankomega}(i), it holds that ${\rm rank}(\boldsymbol{\Omega}) = d(d-1)/2$ if and only if $\lambda_{d(d-1)/2}(\breve{\boldsymbol{\Omega}}^{{\mathrm{\scriptscriptstyle \top} }}\breve{\boldsymbol{\Omega}}) >0$.
Note that $\breve{\boldsymbol{\Omega}}^{{\mathrm{\scriptscriptstyle \top} }}\breve{\boldsymbol{\Omega}}$ is a $\{d(d+1)/2\} \times \{d(d+1)/2\}$ matrix. We require the following mild condition in our theoretical analysis.
\begin{cd}\label{cd:eigen-bomega}
$\lambda_{d(d-1)/2}(\breve{\boldsymbol{\Omega}}^{{\mathrm{\scriptscriptstyle \top} }}\breve{\boldsymbol{\Omega}})$ is uniformly bounded away from zero.
\end{cd}
Write $\hat{{\mathbf A}} = (\hat{{\mathbf a}}_1,\ldots,\hat{{\mathbf a}}_{\hat{d}})$ and $\hat{{\mathbf B}} = (\hat{{\mathbf b}}_1,\ldots,\hat{{\mathbf b}}_{\hat{d}})$, where $(\hat{{\mathbf A}},\hat{{\mathbf B}})$ are specified in Section \ref{sec:estimation}. Recall ${\mathbf A} = ({\mathbf a}_1,\ldots,{\mathbf a}_{d})$ and ${\mathbf B} = ({\mathbf b}_1,\ldots,{\mathbf b}_{d})$. Theorem \ref{thm:hat_a_hat_b} indicates that the columns of $\hat{{\mathbf A}} $ and $\hat{{\mathbf B}}$ are, respectively, consistent to those of ${\mathbf A}$ and ${\mathbf B}$ up to the reflection and permutation indeterminacy.
\begin{theorem}\label{thm:hat_a_hat_b}
Let $d\ge2$ and Conditions \ref{cd:ra}--\ref{cd:eigen-bomega} hold. Select the threshold levels in \eqref{eq:m1h-m2h} and \eqref{eq:zh} as
\[
\delta_1=\breve{C}\sqrt{\frac{\log(pq)}{n}}~~\textrm{and}~~ \delta_2=\tilde{C}\sqrt{\frac{\log(pq)}{n}}
\]
for some sufficiently large constants $\breve{C}, \tilde{C}>0$. If $(\hat{d}_1, \hat{d}_2, \hat{d}) = (d_1,d_2,d)$, there exists a permutation of $(1,\ldots,d)$, denoted by $(j_1,\ldots, j_d)$, such that
\begin{align*}
\max_{\ell \in [d]}| {\kappa}_{1,\ell}\hat{{\mathbf a}}_{j_\ell }-{\mathbf a}_\ell |_2= O_{\rm p} (\Pi_{1,n} + \Pi_{2,n}) = \max_{\ell \in[d]}| {\kappa}_{2,\ell}\hat{{\mathbf b}}_{j_\ell }-{\mathbf b}_\ell|_2
\end{align*}
with some $ {\kappa}_{1,\ell}, {\kappa}_{2,\ell} \in \{1,-1\}$, provided that $ \Pi_{1,n}+\Pi_{2,n} \ll1$ and $\log(pq)\ll n^c$ for some constant $c \in (0,1)$ depending only on $r_1$ and $r_2$.
\end{theorem}
\begin{remark}\label{rek:chang}
The convergence rates of the estimates for ${\mathbf a}_{\ell}$ and ${\mathbf b}_{\ell}$ suggested in \cite{chang2023modelling} are, respectively, $(1+\vartheta_{\ell
}^{-1})\cdot O_{\rm p}(\tilde{\Pi}_{1,n} + \tilde{\Pi}_{2,n})$ and $\{1+(\vartheta_{\ell
}^{*})^{-1}\}\cdot O_{\rm p}(\tilde{\Pi}_{1,n} + \tilde{\Pi}_{2,n})$, where $\vartheta_\ell$ and $\vartheta^*_\ell$ are the eigen-gaps defined as in Equation (38) of \cite{chang2023modelling}, $\tilde{\Pi}_{1,n} = \Pi_{1,n}$, and $\tilde{\Pi}_{2,n}=(\tilde{s}_3\tilde{s}_4)^{1/2} \{n^{-1}\log(pq)\}^{(1-\iota)/2} $. Here, $\tilde{s}_3$ and $\tilde{s}_4$ control the sparsity of the matrix
\begin{align*}
\boldsymbol{\Sigma}_{\mathring{{\mathbf Y}}}(k) = \frac{1}{n-k} \sum_{t=k+1}^{n} \mathbb{E}[\{{\mathbf Y}_{t}-\mathbb{E}(\bar{{\mathbf Y}})\} \otimes {\rm vec}\{{\mathbf Y}_{t-k}-\mathbb{E}(\bar{{\mathbf Y}})\} ] =:\big(\sigma_{\mathring{y}, r,s}^{(k)}\big)_{(p^2q)\times q}
\end{align*}
in the sense that $\sum_{s=1}^{q}|\sigma_{\mathring{y},r,s}^{(k)}|^{\iota} \le \tilde{s}_3$ and $\sum_{r=1}^{p^2q}|\sigma_{\mathring{y},r,s}^{(k)}|^{\iota} \le \tilde{s}_4$ for any $r\in[p^2q]$ and $s\in[q]$. Recall $\Pi_{2,n}=(s_3s_4)^{1/2} \{n^{-1}\log(pq)\}^{(1-\iota)/2} $ with $(s_3, s_4)$ specified in Condition \ref{cd:sgm_yxi}(ii). By direct calculation, we have $s_3 \le p \tilde{s}_3$ and $\tilde{s}_4 \le p s_4$. Under some mild conditions, it holds that $s_3s_4 \asymp \tilde{s}_3\tilde{s}_4$, which implies $\tilde{\Pi}_{2,n} \asymp \Pi_{2,n}$. Hence, if $\vartheta_\ell$ and $\vartheta^*_\ell$ are uniformly bounded away from zero, Theorem \ref{thm:hat_a_hat_b} indicates that our new estimators share the same convergence rates of those proposed in \cite{chang2023modelling}. If $\vartheta_{\ell} \to 0$ or $\vartheta_{\ell}^{*} \to 0$, our new estimators will have faster convergence rates than the estimators considered in \cite{chang2023modelling}.
\end{remark}
\begin{remark}
The model considered in \cite{han2024cp} for order 2 tensor is in the same form as our CP-factor model \eqref{eq:abm}. Therefore, the two estimation procedures proposed in \cite{han2024cp}, the composite PCA (denoted by cPCA) and the High-Order Projection Estimators (denoted by HOPE), can also be used to estimate the loading matrices ${\mathbf A}$ and ${\mathbf B}$ in our CP-factor model \eqref{eq:abm}, where cPCA is a one-pass estimation and HOPE is an iterative refinement initialized at the cPCA solution. \cite{han2024cp} assumes each latent factor $x_{t,\ell}=w_\ell f_{t,\ell}$ where $\{f_{t,\ell}\}_{t\ge 1}$ is stationary with $\mathbb{E}(f_{t,\ell}^2)=1$, and $w_\ell$ represents the signal strength. Under the model setting of \cite{han2024cp}, the latent factor process $\{x_{t,\ell}\}_{t\ge 1}$ is stationary for each $\ell\in[d]$. Moreover, \cite{han2024cp} also assumes $\mathbb{E}(f_{t-h,\ell_1}f_{t,\ell_2}) = 0$ for all $\ell_1 \neq \ell_2$ and $h \ge 1$, which implies $\mathbb{E}(x_{t-h,\ell_1}x_{t,\ell_2}) = 0$ for all $\ell_1 \neq \ell_2$ and $h \ge 1$. However, these assumptions imposed on the latent factors are not necessary in our proposed method. Write $\delta = \| ({\mathbf B} \odot {\mathbf A})^{{\mathrm{\scriptscriptstyle \top} }}({\mathbf B} \odot {\mathbf A}) - \mathbf{I}_d \|_2$, $\psi_{\ell}=w_{\ell}^2\mathbb{E}(f_{t-h,\ell} f_{t,\ell})$ with some fixed lag $h\geq1$, and $\psi_{*} =\min_{\ell \in[d+1]}(\psi_{\ell-1} -\psi_{\ell})$ with $\psi_{0}=\infty$ and $\psi_{d+1}=0$. To simplify the comparison between the theoretical results of \cite{han2024cp} and our proposed method, we ignore the permutation indeterminacy among the estimators. Theorem 1 of \cite{han2024cp} shows that the cPCA estimators $\hat{{\mathbf a}}^{\textup{cpca}}_1,\ldots,\hat{{\mathbf a}}^{\textup{cpca}}_d$, $\hat{{\mathbf b}}^{\textup{cpca}}_1,\ldots,\hat{{\mathbf b}}^{\textup{cpca}}_d$ satisfy
\begin{align*} \label{eq: cpca error bound}
& \max_{\ell\in[d]}\{1 - ({\mathbf a}_\ell^{{\mathrm{\scriptscriptstyle \top} }}\hat{{\mathbf a}}^{\textup{cpca}}_\ell)^2\}^{1/2} + \max_{\ell\in[d]}\{1 - ({\mathbf b}_\ell^{{\mathrm{\scriptscriptstyle \top} }}\hat{{\mathbf b}}^{\textup{cpca}}_\ell)^2\}^{1/2}\\
&~~~~~~~~\lesssim \bigg(1+\frac{2\psi_1}{\psi_{*}}\bigg) \delta + \psi_{*}^{-1}\bigg\{\max_{\ell\in[d]}w_{\ell}^2\sqrt{\frac{\log n}{n}} + \bigg(1+\max_{\ell\in[d]}w_{\ell} \bigg)\sqrt{\frac{pq}{n}}\bigg\}
\end{align*}
with probability at least $1 - (nd)^{-C_1} - e^{-pq}$, where $C_1$ is a positive constant. Theorem 2 of \cite{han2024cp} shows that, after a sufficient number of iterations, the HOPE estimators $\hat{{\mathbf a}}^{\textup{iso}}_1,\ldots,\hat{{\mathbf a}}^{\textup{iso}}_d$, $\hat{{\mathbf b}}^{\textup{iso}}_1,\ldots,\hat{{\mathbf b}}^{\textup{iso}}_d$ satisfy
\begin{align*}
\max_{\ell\in[d]}\{1 - ({\mathbf a}_\ell^{{\mathrm{\scriptscriptstyle \top} }}\hat{{\mathbf a}}^{\textup{iso}}_\ell)^2\}^{1/2}+ \max_{\ell\in[d]}\{1 - ({\mathbf b}_\ell^{{\mathrm{\scriptscriptstyle \top} }}\hat{{\mathbf b}}^{\textup{iso}}_\ell)^2\}^{1/2} &\lesssim (\psi_{d}^{-1} +\psi_{d}^{-1/2} ) \sqrt{\frac{\max(p,q)}{n}}
\end{align*}
with probability at least $1 - (nd)^{-C_2} - e^{-p} - e^{-q}$, provided that the cPCA estimators satisfy certain convergence rates, where $C_2$ is a positive constant. For our proposed estimators $\hat{{\mathbf a}}_1,\ldots,\hat{{\mathbf a}}_d,\hat{{\mathbf b}}_1,\ldots,\hat{{\mathbf b}}_d$, due to $1 - (\hat{{\mathbf a}}_\ell^{{\mathrm{\scriptscriptstyle \top} }}{\mathbf a}_\ell)^2 \le | \kappa_{1,\ell} \hat{{\mathbf a}}_{\ell} - {\mathbf a}_\ell |_2^2 $ and $ 1 - (\hat{{\mathbf b}}_\ell^{{\mathrm{\scriptscriptstyle \top} }}{\mathbf b}_\ell)^2 \le | \kappa_{2,\ell} \hat{{\mathbf b}}_{\ell} - {\mathbf b}_\ell |_2^2 $ for $\kappa_{1,\ell},\kappa_{2,\ell} \in \{1,-1\}$, then
\begin{align*}
\max_{\ell\in[d]}\{1 - ({\mathbf a}_\ell^{{\mathrm{\scriptscriptstyle \top} }}\hat{{\mathbf a}}_\ell)^2\}^{1/2} + \max_{\ell\in[d]}\{1 - ({\mathbf b}_\ell^{{\mathrm{\scriptscriptstyle \top} }}\hat{{\mathbf b}}_\ell)^2\}^{1/2} &\lesssim \Pi_{1,n} + \Pi_{2,n}
\end{align*}
with probability approaching one, where $\Pi_{1,n}$ and $\Pi_{2,n}$ are specified in \eqref{eq:Pi1 and Pi2}. Hence, the two estimation procedures proposed in \cite{han2024cp} can only work for $pq \ll n$, while our proposed method allows $p,q\gg n$. More importantly, in order to obtain the consistency of the cPCA estimators, we need to require ${\mathbf B}\odot {\mathbf A}$ to be very close to an orthonormal matrix ($\delta\rightarrow0$ as $n\rightarrow\infty$). However, such requirement may be too restrictive in practice. The larger $\delta$ is, or the smaller $\psi_*$ is, the worse convergence rate of the cPCA estimators will be. Since the HOPE estimators are obtained through an iterative refinement method initialized with the cPCA estimators, the HOPE estimators will perform poorly if the cPCA estimators have large estimation errors. However, the convergence rate of our proposed method does not depend on these quantities.
\end{remark}
Theorem \ref{thm:hat_a_hat_b} requires $(\hat{d}_1, \hat{d}_2, \hat{d}) = (d_1,d_2,d)$. By Theorem \ref{thm:rank}, we have $\mathbb{P}\{(\hat{d}_1,\hat{d}_2,\hat{d})=(d_1,d_2,d)\}\rightarrow1$ as $n\rightarrow\infty$. Hence, such requirement is reasonable in our theoretical analysis. More generally, without assuming $(\hat{d}_1, \hat{d}_2, \hat{d}) =(d_1, d_2,d)$, we can consider to measure the difference between ${\mathbf A} = ({\mathbf a}_1,\ldots,{\mathbf a}_d)$ and $\hat{{\mathbf A}} = (\hat{{\mathbf a}}_1,\ldots,\hat{{\mathbf a}}_{\hat{d}})$ by
\begin{equation}\label{eq:rho_A}
\varpi^2({\mathbf A},\hat {\mathbf A}) = \max_{\ell \in [d]}\min_{j \in [\hat{d}]}(1 - |\hat{{\mathbf a}}_j^{{\mathrm{\scriptscriptstyle \top} }} {\mathbf a}^{ }_\ell|^2)\,.
\end{equation}
Also, we can measure the difference between ${\mathbf B} = ({\mathbf b}_1,\ldots,{\mathbf b}_d)$ and $\hat{{\mathbf B}} = (\hat{{\mathbf b}}_1,\ldots,\hat{{\mathbf b}}_{\hat{d}})$ by
\begin{equation}\label{eq:rho_B}
\varpi^2({\mathbf B},\hat {\mathbf B}) = \max_{\ell \in [d]}\min_{j \in [\hat{d}]}(1 - |\hat{{\mathbf b}}_j^{{\mathrm{\scriptscriptstyle \top} }} {\mathbf b}^{ }_\ell|^2)\,.
\end{equation}
Consider the event $\mathcal{G} =\{(\hat{d}_1, \hat{d}_2, \hat{d}) =(d_1, d_2,d)\}$. Due to $|\hat{{\mathbf a}}_{j}|_{2} = 1=|{\mathbf a}_{\ell}|_{2}$ and $| \kappa_{1,\ell}\hat{{\mathbf a}}_{j_\ell} -{\mathbf a}_\ell |^2_2 \ge 2 - 2 | \hat{{\mathbf a}}_{j_\ell}^{{\mathrm{\scriptscriptstyle \top} }} {\mathbf a}_\ell |$ for any $ \kappa_{1,\ell} \in \{1, -1\}$, restricted on $\mathcal{G} $, Theorem \ref{thm:hat_a_hat_b} indicates that
$1 - | \hat{{\mathbf a}}_{j_\ell}^{{\mathrm{\scriptscriptstyle \top} }}{\mathbf a}_\ell |^2 \le 2( 1 - |\hat{{\mathbf a}}_{j_\ell}^{{\mathrm{\scriptscriptstyle \top} }}{\mathbf a}_\ell |) = O_{\rm p}(\Pi_{1,n}^2 + \Pi_{2,n}^2)$
provided that $ \Pi_{1,n} + \Pi_{2,n} \ll 1$ and $\log(pq)\ll n^c$ for some constant $c \in (0,1)$ depending only on $r_1$ and $r_2$. Hence, restricted on $\mathcal{G} $, for any $\epsilon >0$, there exists some constant $C_{\epsilon}>0$ such that $\mathbb{P}\{\varpi^2({\mathbf A},\hat {\mathbf A}) > C_{\epsilon} (\Pi_{1,n}^2 + \Pi_{2,n}^2)\,|\, \mathcal{G} \} \le \epsilon$. Together with Theorem \ref{thm:rank}, we have
\begin{align*}
&\mathbb{P}\{\varpi^2({\mathbf A},\hat {\mathbf A}) > C_{\epsilon} (\Pi_{1,n}^2 + \Pi_{2,n}^2) \} \\
&~~~~~~~\le \mathbb{P}\{\varpi^2({\mathbf A},\hat {\mathbf A}) > C_{\epsilon} (\Pi_{1,n}^2 + \Pi_{2,n}^2) \,|\, \mathcal{G} \} \, \mathbb{P}(\mathcal{G} )+ \mathbb{P}(\mathcal{G}^{\rm c})\\
&~~~~~~~\le \mathbb{P}\{\varpi^2({\mathbf A},\hat {\mathbf A}) > C_{\epsilon} (\Pi_{1,n}^2 + \Pi_{2,n}^2) \,|\, \mathcal{G} \} + \mathbb{P}(\hat{d_1}\ne d_1) + \mathbb{P}(\hat{d_2}\ne d_2)+ \mathbb{P}(\hat{d}\ne d)\\
&~~~~~~~\le \epsilon + o(1)\to \epsilon
\end{align*}
as $n\rightarrow\infty$, which implies $\varpi^2({\mathbf A},\hat {\mathbf A}) =O_{\rm p}(\Pi_{1,n}^2 + \Pi_{2,n}^2)$. Also, we can show $\varpi^2({\mathbf B},\hat {\mathbf B})=O_{\rm p}(\Pi_{1,n}^2 + \Pi_{2,n}^2)$.
\section{Numerical studies}\label{section:simulation}
In this section, we will evaluate the finite-sample performance of our proposed method by simulation and real data analysis. The simulation setup is given in Section \ref{sec:simulation setting up}, and the analysis of the simulation results is presented in Section \ref{sec:sim-est}. The real data analysis is given in Section \ref{sec:application}.
\subsection{Setting up}\label{sec:simulation setting up}
Let ${\mathbf A}^\dag \equiv (a^\dag_{i,j})_{p \times d}$ and ${\mathbf B}^\dag \equiv (b^\dag_{i,j})_{q \times d}$ with the elements drawn from the uniform distribution on $[-3,3]$ independently satisfying ${\rm rank}({\mathbf A}^\dag) =d = {\rm rank}({\mathbf B}^\dag)$. Define ${\mathbf P} \in \mathbb{R}^{p\times d_1}$ and ${\mathbf Q}\in\mathbb{R}^{q\times d_2}$ such that the columns of ${\mathbf P}$ and ${\mathbf Q}$ are, respectively, the $d_1$ and $d_2$ left-singular vectors corresponding to the $d_1$ and $d_2$ largest singular values of ${\mathbf A}^{\dag}$ and ${\mathbf B}^{\dagger}$.
Let ${\mathbf U}^* = {\mathbf P}^{\mathrm{\scriptscriptstyle \top} } {\mathbf A}^\dagger =({\mathbf u}_1^*,\ldots,{\mathbf u}_d^* )$ and ${\mathbf V}^* = {\mathbf Q}^{\mathrm{\scriptscriptstyle \top} } {\mathbf B}^\dagger =({\mathbf v}_1^*,\ldots,{\mathbf v}_d^* )$. Derive ${\mathbf U} = ({\mathbf u}_1,\ldots,{\mathbf u}_d )$ and ${\mathbf V} = ({\mathbf v}_1,\ldots,{\mathbf v}_d )$ with ${\mathbf u}_j = {\mathbf u}_j^*/ |{\mathbf u}_j^*|_2$ and ${\mathbf v}_j = {\mathbf v}_j^*/ | {\mathbf v}_j^* |_2$ for any $j \in [d]$. Write $\mathbf{x}^*_j = (x^*_{1,j},\ldots,x^*_{n,j})^{{\mathrm{\scriptscriptstyle \top} }}$ and let $\mathbf{x}^*_1,\ldots,\mathbf{x}^*_d$ be $d$ independent AR(1) processes with independent $\mathcal{N}(0,1)$ innovations, and the autoregressive coefficients drawn from the uniform distribution on $[-0.95,-0.6]\cup [0.6,0.95]$. Let ${\mathbf X}_t = \text{diag}( x_{t,1},\ldots, x_{t,d} )$ with $ x_{t,j} = x^*_{t,j}| {\mathbf v}_j^* |_2 | {\mathbf u}_j^* |_2 $ for each $t\in [n]$. The elements of the error term $\boldsymbol{\varepsilon}_t$ are drawn from $\mathcal{N}(0,1)$ independently.
Finally, we generate ${\mathbf Y}_t = {\mathbf A} {\mathbf X}_t {\mathbf B}^{{\mathrm{\scriptscriptstyle \top} }} + \boldsymbol{\varepsilon}_t$ for any $t \in [n]$ with ${\mathbf A}={\mathbf P}{\mathbf U}$ and ${\mathbf B}={\mathbf Q}{\mathbf V}$. We set $n \in \{300, 600, 900\}$, $d \in \{3,5,7\}$ and $p,q$ taking values between 10 and 160. We consider three different scenarios for $(d,d_1,d_2)$:
\begin{enumerate}
\item[(R1)] Let $d_1=d_2=d$. In this scenario, ${\mathbf A}$ and ${\mathbf B}$ are full rank.
\item[(R2)] Let $d_1 = d - 1$ and $d_2 = d$. In this scenario, only ${\mathbf B}$ is full rank.
\item[(R3)] Let $d_1 = d_2 = d - 1$. In this scenario, both ${\mathbf A}$ and ${\mathbf B}$ are not full rank.
\end{enumerate}
We follow \cite{chang2023modelling} to specify $\xi_t$ involved in \eqref{eq:esthatsig}. Let ${\mathbf Y}=(\vec{\mathbf Y}_1,\ldots, \vec {\mathbf Y}_{n})^{{\mathrm{\scriptscriptstyle \top} }}$. Perform the principal component analysis for ${\mathbf Y}$ and select $\xi_t$ as the average of the first $m$ principal components corresponding to the eigenvalues which count for at least 99\% of the total variations. Let $\hat\sigma_0^2=(npq)^{-1}\| {\mathbf Y}\|_{\rm F}^2$. We set $\delta_1 = \delta_2 = \hat\sigma_0 \{n^{-1}\log(pq)\}^{1/2}$ in \eqref{eq:m1h-m2h} and \eqref{eq:zh}, and set $c_{1,n} = c_{2,n} = c_{3,n} = \hat\sigma_0 n^{-1}$ in \eqref{eq:d1h} and \eqref{eq:dh}. We also choose $K = 20$ and $\tilde{K} = 10$ with $K$ and $\tilde{K}$ given in \eqref{eq:m1h-m2h} and \eqref{eq:mh}, respectively. Here, using a relatively large value for $K$ is to ensure that ${\mathbf M}_1$ and ${\mathbf M}_2$ defined in \eqref{eq:M1M2} satisfy $\textup{rank}({\mathbf M}_1) = d_1$ and $\textup{rank}({\mathbf M}_2) = d_2$. These two requirements are essential for our proposed method. See Propositions \ref{pro:m1-rank-con} and \ref{pro:rankwith-xt}. As shown in \eqref{eq:mh}, $\tilde{K}$ is the number of lags used in the methods of \cite{lam2011estimation}, \cite{lam2012factor} and \cite{Chang2015} to estimate the linear space spanned by the columns of the factor loading matrix in the standard factor model. In practice, a small $\tilde{K}$ (i.e., $1\leq \tilde{K}\leq 10$) is enough and the estimation results are generally robust to the specific choice of $\tilde{K}$. See our sensitivity analysis with respect to the tuning parameters $K$ and $\tilde{K}$ in Figures \ref{fig:K-estimation}--\ref{fig:Ktilde-rank} of the supplementary material for more details.
As mentioned in Section \ref{sec:theta-hat-est}, we need to select an appropriate constant vector
$\boldsymbol{\phi} = (\phi_1, \ldots, \phi_{\hat{d}})^{{\mathrm{\scriptscriptstyle \top} }}$ to ensure that $\tilde{{\mathbf H}} = \sum_{i=1}^{\hat{d}} \phi_i \tilde{{\mathbf H}}_i$ is an invertible matrix with $\tilde{{\mathbf H}}_i$ defined below \eqref{eq:init-basis}. Let $\phi_{i} = I\{\sigma_{\hat{d}}(\tilde{{\mathbf H}}_{i}) = \max_{j\in[\hat{d}]}\sigma_{\hat{d}}(\tilde{{\mathbf H}}_{j}) > 0\}$ for any $i\in[\hat{d}]$. If $|\boldsymbol{\phi}|_1 = 0$, we randomly generate a unit vector $\boldsymbol{\phi}$ such that $\sigma_{\hat{d}}(\tilde{{\mathbf H}}) > 0$. If $|\boldsymbol{\phi}|_1 \ge 1$, we arbitrarily keep one non-zero element in $\boldsymbol{\phi}$ and set all other elements to zero. The simulation results show that our proposed procedure based on such selected $\boldsymbol{\phi}$ exhibits good finite-sample performance. We also compare our proposed method with the refined method (denoted by CP-refined) introduced by \cite{chang2023modelling}, and the cPCA and the HOPE methods proposed by \cite{han2024cp} with the recommended tuning parameter $h = 1$ therein.
All simulations are implemented in \textsf{R}. Our proposed method is available in \textsf{R}-package \texttt{HDTSA}, which is implemented by calling the \textsf{R}-function \texttt{CP\_MTS} with setting \texttt{method = `CP.Unified'}. The CP-refined method of \cite{chang2023modelling} can also be implemented by calling the \textsf{R}-function \texttt{CP\_MTS} with setting \texttt{method = `CP.Refined'}. All simulation results are based on 2000 replications.
\subsection{Simulation results }\label{sec:sim-est}
We first consider the finite-sample performance of the estimation $(\hat{d}_1, \hat{d}_2,\hat{d})$ given in \eqref{eq:d1h} and \eqref{eq:dh}.
Note that the CP-refined method of \cite{chang2023modelling} is developed under the assumption $d_1=d_2=d$. To fairly compare our proposed method and the CP-refined method, we compare the relative frequency estimate of ${\mathbb{P}}_{d}: ={\mathbb{P}}(\hat{d} = d)$ with $\hat{d}$ specified in \eqref{eq:dh}, and the relative frequency estimate of ${\mathbb{P}}_c :={\mathbb{P}}( \hat{d} = d )$ with $\hat{d} $ estimated by the CP-refined method.
Note that ${\mathbb{P}}_{1,2,d}: ={\mathbb{P}}\{(\hat{d}_1 , \hat{d}_2
, \hat{d} ) = (d_1,d_2,d)\}\leq \mathbb{P}_d$ with $(\hat{d}_1, \hat{d}_2,\hat{d})$ estimated by our proposed method. Table \ref{table:rf-all} indicates that (i) our proposed method outperforms the CP-refined method across Scenarios R1--R3, and (ii) $(d_1,d_2,d)$ can be consistently estimated by our proposed method. To conserve space, we omit the results for $p<q$ in Scenarios R1 and R3, as the symmetry in the data-generating process leads to results that are nearly identical to those obtained when $p> q$.
Figure \ref{fig:est-case1} reports the averages of the estimation errors $\varpi^2({\mathbf A},\hat{{\mathbf A}})$ and $\varpi^2({\mathbf B},\hat{{\mathbf B}})$ defined in \eqref{eq:rho_A} and \eqref{eq:rho_B} based on 2000 repetitions across different scenarios.
Our proposed method consistently outperforms all competing methods, except in Scenario R1 with $p > q$, where it performs comparably to the HOPE method in estimating ${\mathbf B}$. In contrast, the estimation errors of the CP-refined method are very large in Scenarios R2 and R3, which indicates that the CP-refined method does not work for the matrix CP-factor model \eqref{eq:abm} with rank-deficient factor loading matrices ${\mathbf A}$ and ${\mathbf B}$. Also, in Scenarios R2 and R3, the HOPE method offers no notable improvement over the cPCA method and even underperforms the cPCA method in some settings, suggesting that the iterative method HOPE is ineffective when the factor loading matrices ${\mathbf A}$ and ${\mathbf B}$ are rank-deficient. Additionally, all methods lose efficiency when $(d,d_1,d_2) = (3,2,2)$ since $({\mathbf A},{\mathbf B})$ cannot be identified uniquely, but our proposed method still yields the smallest estimation errors.
The averages and standard deviations of the estimation errors $\varpi^2({\mathbf A},\hat{{\mathbf A}})$ and $\varpi^2({\mathbf B},\hat{{\mathbf B}})$ based on 2000 repetitions are summarized in Tables \ref{table:varpi-case1}--\ref{table:varpi-case3} in the supplementary material.
Next, we evaluate the finite-sample performance of our proposed prediction method introduced in Section \ref{sec:prediction}. We generate a sequence $\{{\mathbf Y}_t\}_{t=1}^{n+m+1}$ defined in Section \ref{sec:simulation setting up} with $m=20$. For any $s \in[m]$, we apply our proposed prediction method to the data $\{{\mathbf Y}_t\}^{n+s-1}_{t=s}$
and then, respectively, obtain the one-step forecast of ${\mathbf Y}_{n+s}$ (denoted by $\hat{{\mathbf Y}}^{(1)}_{n+s}$) and the two-step forecast of ${\mathbf Y}_{n+s+1}$ (denoted by $\hat{{\mathbf Y}}^{(2)}_{n+s+1}$).
We also consider the prediction method introduced in \cite{chang2023modelling} to obtain the one-step ahead forecast of ${\mathbf Y}_{n+s}$ and the two-step ahead forecast of ${\mathbf Y}_{n+s+1}$ using $(\hat{{\mathbf A}}, \hat{{\mathbf B}})$ estimated from the data $\{{\mathbf Y}_t\}^{n+s-1}_{t=s}$ for each $s \in [m]$, where $(\hat{{\mathbf A}}, \hat{{\mathbf B}})$ can be selected as either (i) our proposed estimate of $({\mathbf A}, {\mathbf B})$ specified in Section \ref{sec:estimation}, or (ii) the CP-refined estimate of $({\mathbf A}, {\mathbf B})$ given in \cite{chang2023modelling}.
Here, for the obtained univariate time series, we fit it by an autoregressive (AR) model with the order determined by the Akaike information criterion (AIC). For the obtained multivariate time series, we fit it by a vector autoregressive (VAR) model with the order determined by the AIC. Based on 2000 repetitions,
Figure \ref{fig:fore-case1} plots the averages of the one-step ahead $${\textup{RMSE} } := \frac{1}{m\sqrt{pq}} \sum_{s=1}^{m} \| \hat{{\mathbf Y}}^{(1)}_{n+s} - {\mathbf Y}_{n+s} \|_\text{F}\,.$$
It can be observed that (i) in all cases, the finite-sample performance of our newly proposed prediction method is better than the prediction method introduced in \cite{chang2023modelling} with selecting $(\hat{{\mathbf A}},\hat{{\mathbf B}})$ as the CP-refined estimate, (ii) in the cases expect $(d,d_1,d_2) = (3,2,2)$, the averages of the one-step ahead $\textup{RMSE}$ of our newly proposed prediction method are almost identical to those of the prediction method introduced in \cite{chang2023modelling} with selecting $(\hat{{\mathbf A}},\hat{{\mathbf B}})$ as our proposed estimate, and (iii) in the case $(d,d_1,d_2) = (3,2,2)$, our newly proposed prediction method outperforms the prediction method introduced in \cite{chang2023modelling} with selecting $(\hat{{\mathbf A}},\hat{{\mathbf B}})$ as our proposed estimate. Note that the factor loading matrices ${\mathbf A}$ and ${\mathbf B}$ cannot be uniquely identified in the case $(d,d_1,d_2) = (3,2,2)$. Hence, we can conclude that (i) when $({\mathbf A},{\mathbf B})$ can be uniquely identified, the prediction method introduced in \cite{chang2023modelling} with selecting $(\hat{{\mathbf A}},\hat{{\mathbf B}})$ as our proposed estimate works quite well, which has almost identical performance as our newly proposed prediction method; and (ii) our newly proposed prediction method works very well regardless of whether $({\mathbf A},{\mathbf B})$ can be uniquely identified or not. The results of two-step ahead forecasting are similar to that of one-step ahead forecasting. See Figure \ref{fig:fore2-case1} in the supplementary material for details.
\subsection{Real data analysis}\label{sec:application}
In this section, we illustrate the proposed method for the matrix CP-factor model \eqref{eq:abm} by using the Fama-French $10 \times 10$ return series. We collect the monthly returns from January 1964 to December 2021, which contains 69600 observations for total 696 months. The data are downloaded from \url{http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html}. The portfolios are formed by the intersections of 10 levels of size, denoted by (${\rm S}_{1},\ldots,{\rm S}_{10}$), and 10 levels of the book equity to market equity ratio (BE), denoted by $({\rm BE}_{1}, \ldots, {\rm BE}_{10}) $. The data contain a small number of missing values in the early years and we transform them to zeros. Since all the 100 series are clearly related to the overall market condition, following \cite{wang2019factor}, we decide to remove the influence of market effects before empirical analysis. Two filtering approaches are considered: (i) (CAPM filtering) fitting a standard CAPM model \citep{fama1973risk} to each of the series to remove the market effect, (ii) (Demean filtering) subtracting the corresponding monthly excess market return from each of the series. The market return data are obtained from the same website above. Based on each filtering approach, we finally obtain 100 market-adjusted return series.
The 100 market-adjusted return series can be represented as a $10\times10$ matrix time series ${\mathbf Y}_t = (y_{i,j,t})$ for $t\in[696]$ (i.e., $p=q=10$, $n=696$), where $y_{i,j,t}$ is the market-adjusted return at the $i$-th level of size ${\rm S}_{i}$ and the $j$-th level of the BE-ratio ${\rm BE}_{j}$ at time $t$. Figure \ref{fig:app-timeseries} shows the time series plots of the market-adjusted return series $\{y_{i,j,t}\}_{t=1}^n$ based on the CAPM filtering for $i,j \in[10]$. The rows in Figure \ref{fig:app-timeseries} correspond to the ten levels of size and the columns correspond to the ten levels of the BE-ratio. All series are stationary because they reject the null hypothesis of Augmented Dickey-Fuller test at 5\% significance level.
We evaluate the post-sample forecasting performance of our proposed method introduced in Section \ref{sec:prediction} by performing the one-step and two-step ahead rolling forecasts for the 240 monthly readings in the last twenty years (2002--2021). To do this, we first use the data $\{{\mathbf Y}_t\}_{t=1}^{456}$ to determine the rank parameters $(d,d_1,d_2)$. With the tuning parameters selected as those in Section \ref{sec:simulation setting up}, our proposed method obtains $(\hat{d}, \hat{d}_1, \hat{d}_2) = (2, 2, 1)$, which aligns with the conventional scree plots of $\hat{{\mathbf M}}_1$ and $\hat{{\mathbf M}}_2$ given in Figures \ref{fig: app-eigenall}(a) and \ref{fig: app-eigenall}(b), respectively. We adopt $(\hat{d},\hat{d}_1,\hat{d}_2) = (2,2,1)$ in the rolling forecasts. For each $s\in [240]$, we apply our proposed prediction method to the data $\{{\mathbf Y}_t\}_{t=s}^{455+s}$ and then obtain the one-step forecast of ${\mathbf Y}_{456+s}$, denoted by $\hat{{\mathbf Y}}^{(1)}_{456+s} = (\hat{y}^{(1)}_{i,j,456+s}) $. For the two-step ahead forecast, we apply our proposed prediction method to the data $\{{\mathbf Y}_{t}\}_{t=s}^{454+s}$, and the two-step ahead forecast $\hat{{\mathbf Y}}^{(2)}_{456+s} = (\hat{y}^{(2)}_{i,j,456+s}) $ can be obtained by plug-in the one-step forecast into the fitted model. More specifically, for each $s\in[240]$, we fit the obtained 2-dimensional time series by a VAR model with the order determined by the AIC. For comparison, we can also fit $\{{\mathbf Y}_{t}\}_{t = s}^{455+s}$ and $\{{\mathbf Y}_{t}\}_{t = s}^{454+s}$ by the following methods and obtain the associated one-step and two-step ahead forecasts:
\begin{itemize}
\item (CP-refined) The CP-refined method of \cite{chang2023modelling} with the pre-determined parameter $K = 10$ therein. The associated rank in this method is estimated as $\hat{d} = 1$ based on $\{{\mathbf Y}_t\}_{t=1}^{456}$ and then fixed in the rolling forecasts. Motivated by the scree plot in Figure \ref{fig: app-eigenall}(c), we also consider $\hat{d} = 2$ as an alternative. For $\hat{d} = 1$, we fit the obtained univariate time series by an AR model with the order determined by the AIC. For $\hat{d} = 2$, we fit the obtained 2-dimensional time series by a VAR model with the order determined by the AIC. The methods with $\hat{d} = 1$ and $\hat{d} = 2$ are referred to as CP-refined(1) and CP-refined(2), respectively.
\item (cPCA, HOPE) The composite PCA and High-Order Projection Estimators in \cite{han2024cp} with the recommended tuning parameter $h = 1$ therein. Following the same rank specification strategy as in the CP-refined method, we consider both $\hat{r} = 1$ and $\hat{r} = 2$ for the associated rank in these two methods. For $\hat{r} = 1$, we fit the obtained univariate time series by an AR model with the order determined by the AIC. For $\hat{r} = 2$, we fit the obtained 2-dimensional time series by a VAR model with the order determined by the AIC. The methods with $\hat{r} = 1$ are referred to as cPCA(1) and HOPE(1), while those with $\hat{r} = 2$ are denoted as cPCA(2) and HOPE(2).
\item (FAC) The matrix Tucker-factor model with the FAC method proposed by \cite{wang2019factor} with the pre-determined parameter $h_0 = 1$ as suggested therein. The associated ranks in this model are estimated as $(\hat{k}_1,\hat{k}_2) = (1,1)$ by the ratio estimators suggested therein based on $\{{\mathbf Y}_{t}\}_{t = 1}^{456}$, and are fixed in the rolling forecasts. Motivated by the scree plots in Figures \ref{fig: app-eigenall}(d) and \ref{fig: app-eigenall}(e), we also consider an alternative setting with $(\hat{k}_1,\hat{k}_2) = (2,1)$. For $(\hat{k}_1,\hat{k}_2) = (1,1)$, we fit the obtained univariate time series by an AR model with the order determined by the AIC. For $(\hat{k}_1,\hat{k}_2) = (2,1)$, we fit the obtained 2-dimensional time series by a VAR model with the order determined by the AIC. The methods with $(\hat{k}_1,\hat{k}_2) = (1,1)$ and $(\hat{k}_1,\hat{k}_2) = (2,1)$ are referred to as FAC(1,1) and FAC(2,1), respectively.
\item (TOPUP, TIPUP) The Time series Outer-Product Unfolding Procedure and the Time series Inner-Product Unfolding Procedure proposed by \cite{han2024tensor} for the matrix Tucker-factor model. The associated ranks in this model are estimated as $(\hat{k}_1,\hat{k}_2) = (2,2)$ by the information criterion considered in \cite{han2022rank} based on $\{{\mathbf Y}_{t}\}_{t = 1}^{456}$, and are fixed in the rolling forecasts. We fit the obtained 4-dimensional time series by a VAR model with the order determined by the AIC. The methods are implemented using the R package \texttt{tensorTS}.
\item (MAR) The matrix-AR(1) model of \cite{chen2021autoregressive}.
\item (TS-PCA) Apply the principal component analysis for time series proposed by \cite{Chang2018} to the 100-dimensional time series $\{\vec{{\mathbf Y}}_{t}\}_{t = s}^{455+s}$ and $\{\vec{{\mathbf Y}}_{t}\}_{t = s}^{454+s}$, respectively, to obtain the associated one-step and two-step ahead forecasts. The method is implemented using the R package \texttt{HDTSA}. For the obtained univariate time series, we fit it by an AR model with the order determined by the AIC. For the obtained multivariate time series, we fit it by a VAR model with the order determined by the AIC.
\item (UniAR) Fit each of 100 component time series by an AR model with the order determined by the AIC.
\end{itemize}
For each $s \in [240]$, the one-step ahead forecasting performance is evaluated by the $\textup{rRMSE}(s)$ and $\textup{rMAE}(s)$ defined as
\begin{gather*}
\text{rRMSE}(s) = \bigg\{ \frac{1}{100} \sum_{i = 1}^{10}\sum_{j = 1}^{10}|\hat{y}^{(1)}_{i,j,456+s} - y_{i,j,456+s}|^2 \bigg\}^{1/2}\,,\\
\text{rMAE}(s) = \frac{1}{100} \sum_{i = 1}^{10}\sum_{j = 1}^{10} |\hat{y}^{(1)}_{i,j,456+s} - y_{i,j,456+s}| \,.
\end{gather*}
For the two-step ahead forecast, we can evaluate it by the associated $\textup{rRMSE}(s)$ and $\textup{rMAE}(s)$ analogously. Table \ref{table:app-forecast} reports the averages of $\{\textup{rRMSE}(s)\}_{s=1}^{240}$ and $\{\textup{rMAE}(s)\}_{s=1}^{240}$, denoted by $\textup{rRMSE}$ and $\textup{rMAE}$, respectively. The standard deviations of $\{\textup{rRMSE}(s)\}_{s=1}^{240}$ and $\{\textup{rMAE}(s)\}_{s=1}^{240}$ are reported in parentheses.
As shown in Table \ref{table:app-forecast}, under CAPM filtering (Panel A), our proposed method achieves the lowest rRMSE and rMAE for both one- and two-step ahead forecasts, outperforming all competing methods. Under Demean filtering (Panel B), although the HOPE(2) achieves the lowest rRMSE and rMAE, our proposed method performs comparably and yields smaller standard deviations than the HOPE(2).
Overall, the results show that our proposed method delivers robust and accurate forecasts across different market-adjustment schemes, often outperforming alternatives in both accuracy and stability.
\begin{acks}
The authors thank Yuefeng Han for sharing code for implementing the methods proposed in \cite{han2024cp}.
\end{acks}
\section*{Funding}
J. Chang, Y. Du and G. Huang were supported in part by the National Natural Science
Foundation of China (Grant nos. 72125008 and 72495122). Q. Yao was supported in part by the U.K. Engineering and Physical Sciences Research
Council (Grant nos. EP/V007556/1 and EP/X002195/1).
\section*{Supplement Material}
{\bf Supplement to ``Identification and Estimation for Matrix Time Series CP-factor Models''.}
This supplement contains additional simulation studies and all technical proofs.
\bibliographystyle{jasa}
\spacingset{0.95}\selectfont
\bibliography{mybibfile}
\begin{landscape}
$ $\\
$ $\\
\begin{table}[htbp]
\scriptsize
\caption{
Relative frequency estimates of ${\mathbb{P}}_{1,2,d} ={\mathbb{P}}\{(\hat{d}_1 , \hat{d}_2
, \hat{d} ) = (d_1,d_2,d)\}$ and ${\mathbb{P}}_{d} = {\mathbb{P}}(\hat{d} = d)$ with $(\hat{d}_1, \hat{d}_2,\hat{d})$ estimated by our proposed method, and the relative frequency estimate of ${\mathbb{P}}_c={\mathbb{P}}( \hat{d} = d )$ with $\hat{d} $ estimated by the CP-refined method of \cite{chang2023modelling} in Scenarios R1--R3. All numbers reported below are multiplied by 100.
}
\label{table:rf-all}
\resizebox{22.9cm}{!}{
\begin{tabular}{c|c|cccccccc|ccccccccc|cccccc}
\hline\hline
\multirow{3}{*}{$d$} & \multirow{3}{*}{$n$} & \multicolumn{8}{c|}{R1} & \multicolumn{9}{c|}{R2} & \multicolumn{6}{c}{R3} \\ \cline{3-25}
& & \multicolumn{4}{c|}{$p = q$} & \multicolumn{4}{c|}{$p > q$} & \multicolumn{3}{c|}{$p = q$} & \multicolumn{3}{c|}{$p > q$} & \multicolumn{3}{c|}{$p < q$} & \multicolumn{3}{c|}{$p = q$} & \multicolumn{3}{c}{$p > q$} \\
& & $(p,q)$ & $\mathbb{P}_{1,2,d}$ & $\mathbb{P}_{d}$ & \multicolumn{1}{c|}{$\mathbb{P}_{c}$} & $(p,q)$ & $\mathbb{P}_{1,2,d}$ & $\mathbb{P}_{d}$ & $\mathbb{P}_{c}$ & $(p,q)$ & $\mathbb{P}_{1,2,d}$ & \multicolumn{1}{c|}{$\mathbb{P}_{c}$} & $(p,q)$ & $\mathbb{P}_{1,2,d}$ & \multicolumn{1}{c|}{$\mathbb{P}_{c}$} & $(p,q)$ & $\mathbb{P}_{1,2,d}$ & $\mathbb{P}_{c}$ & $(p,q)$ & $\mathbb{P}_{1,2,d}$ & \multicolumn{1}{c|}{$\mathbb{P}_{c}$} & $(p,q)$ & $\mathbb{P}_{1,2,d}$ & $\mathbb{P}_{c}$ \\ \hline
\multirow{9}{*}{3} & 300 & \multirow{3}{*}{$(20,20)$} & 96.59 & 97.04 & \multicolumn{1}{c|}{94.58} & \multirow{3}{*}{$(40,10)$} & 95.11 & 96.74 & 95.52 & \multirow{3}{*}{$(20,20)$} & 85.38 & \multicolumn{1}{c|}{77.02} & \multirow{3}{*}{$(40,10)$} & 77.51 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(10,40)$} & 86.57 & 79.05 & \multirow{3}{*}{$(20,20)$} & 88.86 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(40,10)$} & 87.44 & 0.00 \\
& 600 & & 97.35 & 97.70 & \multicolumn{1}{c|}{96.15} & & 95.16 & 97.02 & 96.01 & & 86.36 & \multicolumn{1}{c|}{79.23} & & 76.70 & \multicolumn{1}{c|}{0.00} & & 87.91 & 81.21 & & 89.86 & \multicolumn{1}{c|}{0.00} & & 89.09 & 0.00 \\
& 900 & & 97.65 & 98.00 & \multicolumn{1}{c|}{97.25} & & 96.41 & 97.88 & 97.07 & & 84.23 & \multicolumn{1}{c|}{77.54} & & 78.55 & \multicolumn{1}{c|}{0.00} & & 89.29 & 82.90 & & 90.97 & \multicolumn{1}{c|}{0.00} & & 89.32 & 0.00 \\ \cline{2-25}
& 300 & \multirow{3}{*}{$(40,40)$} & 98.75 & 98.80 & \multicolumn{1}{c|}{97.90} & \multirow{3}{*}{$(80,10)$} & 94.99 & 97.04 & 96.07 & \multirow{3}{*}{$(40,40)$} & 90.88 & \multicolumn{1}{c|}{83.56} & \multirow{3}{*}{$(80,10)$} & 78.11 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(10,80)$} & 91.00 & 83.13 & \multirow{3}{*}{$(40,40)$} & 92.83 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(80,10)$} & 90.02 & 0.00 \\
& 600 & & 99.50 & 99.55 & \multicolumn{1}{c|}{99.05} & & 96.56 & 97.67 & 97.06 & & 90.74 & \multicolumn{1}{c|}{85.34} & & 78.31 & \multicolumn{1}{c|}{0.00} & & 90.85 & 84.11 & & 93.68 & \multicolumn{1}{c|}{0.00} & & 90.51 & 0.00 \\
& 900 & & 99.80 & 99.80 & \multicolumn{1}{c|}{99.40} & & 96.97 & 98.69 & 97.88 & & 91.17 & \multicolumn{1}{c|}{85.51} & & 76.49 & \multicolumn{1}{c|}{0.00} & & 90.19 & 84.41 & & 93.52 & \multicolumn{1}{c|}{0.00} & & 91.54 & 0.00 \\ \cline{2-25}
& 300 & \multirow{3}{*}{$(80,80)$} & 99.75 & 99.80 & \multicolumn{1}{c|}{99.35} & \multirow{3}{*}{$(160,10)$} & 96.55 & 98.12 & 97.11 & \multirow{3}{*}{$(80,80)$} & 94.31 & \multicolumn{1}{c|}{88.58} & \multirow{3}{*}{$(160,10)$} & 80.50 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(10,160)$} & 92.07 & 84.92 & \multirow{3}{*}{$(80,80)$} & 95.74 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(160,10)$} & 90.92 & 0.00 \\
& 600 & & 100.00 & 100.00 & \multicolumn{1}{c|}{99.50} & & 97.78 & 99.04 & 98.38 & & 94.34 & \multicolumn{1}{c|}{89.89} & & 79.97 & \multicolumn{1}{c|}{0.00} & & 91.85 & 86.37 & & 95.08 & \multicolumn{1}{c|}{0.00} & & 91.35 & 0.00 \\
& 900 & & 100.00 & 100.00 & \multicolumn{1}{c|}{99.55} & & 97.84 & 99.14 & 98.49 & & 94.87 & \multicolumn{1}{c|}{90.66} & & 77.69 & \multicolumn{1}{c|}{0.00} & & 91.17 & 84.95 & & 95.77 & \multicolumn{1}{c|}{0.00} & & 92.01 & 0.00 \\ \hline
\multirow{9}{*}{5} & 300 & \multirow{3}{*}{$(20,20)$} & 97.49 & 98.04 & \multicolumn{1}{c|}{94.88} & \multirow{3}{*}{$(40,10)$} & 94.14 & 98.41 & 95.42 & \multirow{3}{*}{$(20,20)$} & 97.25 & \multicolumn{1}{c|}{91.30} & \multirow{3}{*}{$(40,10)$} & 89.70 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(10,40)$} & 97.12 & 93.58 & \multirow{3}{*}{$(20,20)$} & 97.99 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(40,10)$} & 96.91 & 0.00 \\
& 600 & & 98.20 & 98.55 & \multicolumn{1}{c|}{96.29} & & 94.69 & 98.62 & 97.24 & & 97.19 & \multicolumn{1}{c|}{93.32} & & 89.86 & \multicolumn{1}{c|}{0.00} & & 97.36 & 94.37 & & 98.64 & \multicolumn{1}{c|}{0.00} & & 97.68 & 0.00 \\
& 900 & & 98.10 & 98.55 & \multicolumn{1}{c|}{96.14} & & 95.63 & 98.78 & 97.51 & & 96.40 & \multicolumn{1}{c|}{91.79} & & 90.33 & \multicolumn{1}{c|}{0.00} & & 97.26 & 95.39 & & 98.49 & \multicolumn{1}{c|}{0.00} & & 97.53 & 0.00 \\ \cline{2-25}
& 300 & \multirow{3}{*}{$(40,40)$} & 99.60 & 99.60 & \multicolumn{1}{c|}{99.05} & \multirow{3}{*}{$(80,10)$} & 95.24 & 99.02 & 97.15 & \multirow{3}{*}{$(40,40)$} & 98.99 & \multicolumn{1}{c|}{96.88} & \multirow{3}{*}{$(80,10)$} & 91.13 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(10,80)$} & 97.93 & 94.61 & \multirow{3}{*}{$(40,40)$} & 99.40 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(80,10)$} & 97.66 & 0.00 \\
& 600 & & 99.45 & 99.55 & \multicolumn{1}{c|}{99.10} & & 96.12 & 99.23 & 98.16 & & 99.35 & \multicolumn{1}{c|}{97.39} & & 91.48 & \multicolumn{1}{c|}{0.00} & & 98.14 & 95.98 & & 99.55 & \multicolumn{1}{c|}{0.00} & & 97.93 & 0.00 \\
& 900 & & 99.70 & 99.70 & \multicolumn{1}{c|}{99.35} & & 96.89 & 99.49 & 98.73 & & 99.05 & \multicolumn{1}{c|}{97.39} & & 91.32 & \multicolumn{1}{c|}{0.00} & & 98.54 & 96.48 & & 99.55 & \multicolumn{1}{c|}{0.00} & & 98.08 & 0.00 \\ \cline{2-25}
& 300 & \multirow{3}{*}{$(80,80)$} & 99.90 & 99.90 & \multicolumn{1}{c|}{99.70} & \multirow{3}{*}{$(160,10)$} & 96.26 & 99.18 & 98.05 & \multirow{3}{*}{$(80,80)$} & 99.75 & \multicolumn{1}{c|}{98.50} & \multirow{3}{*}{$(160,10)$} & 92.28 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(10,160)$} & 98.29 & 96.28 & \multirow{3}{*}{$(80,80)$} & 99.80 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(160,10)$} & 98.49 & 0.00 \\
& 600 & & 99.95 & 99.95 & \multicolumn{1}{c|}{99.65} & & 96.84 & 99.23 & 98.52 & & 99.70 & \multicolumn{1}{c|}{98.85} & & 91.21 & \multicolumn{1}{c|}{0.00} & & 98.59 & 96.83 & & 100.00 & \multicolumn{1}{c|}{0.00} & & 98.04 & 0.00 \\
& 900 & & 99.95 & 99.95 & \multicolumn{1}{c|}{99.70} & & 97.01 & 99.70 & 99.24 & & 99.85 & \multicolumn{1}{c|}{99.00} & & 90.56 & \multicolumn{1}{c|}{0.00} & & 99.30 & 97.79 & & 99.95 & \multicolumn{1}{c|}{0.00} & & 97.63 & 0.00 \\ \hline
\multirow{9}{*}{7} & 300 & \multirow{3}{*}{$(20,20)$} & 97.94 & 98.60 & \multicolumn{1}{c|}{95.84} & \multirow{3}{*}{$(40,10)$} & 84.63 & 98.69 & 96.64 & \multirow{3}{*}{$(20,20)$} & 98.53 & \multicolumn{1}{c|}{94.99} & \multirow{3}{*}{$(40,10)$} & 83.97 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(10,40)$} & 97.68 & 96.36 & \multirow{3}{*}{$(20,20)$} & 99.00 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(40,10)$} & 98.23 & 0.00 \\
& 600 & & 98.59 & 99.39 & \multicolumn{1}{c|}{97.33} & & 87.07 & 99.00 & 97.84 & & 99.19 & \multicolumn{1}{c|}{95.76} & & 84.77 & \multicolumn{1}{c|}{0.00} & & 98.24 & 96.78 & & 99.70 & \multicolumn{1}{c|}{0.00} & & 98.47 & 0.00 \\
& 900 & & 98.85 & 99.35 & \multicolumn{1}{c|}{97.25} & & 89.27 & 98.95 & 98.06 & & 98.74 & \multicolumn{1}{c|}{96.38} & & 85.96 & \multicolumn{1}{c|}{0.00} & & 98.23 & 96.87 & & 99.40 & \multicolumn{1}{c|}{0.00} & & 97.93 & 0.00 \\ \cline{2-25}
& 300 & \multirow{3}{*}{$(40,40)$} & 99.85 & 99.85 & \multicolumn{1}{c|}{99.30} & \multirow{3}{*}{$(80,10)$} & 85.49 & 98.85 & 97.39 & \multirow{3}{*}{$(40,40)$} & 99.90 & \multicolumn{1}{c|}{98.60} & \multirow{3}{*}{$(80,10)$} & 83.71 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(10,80)$} & 98.13 & 97.73 & \multirow{3}{*}{$(40,40)$} & 100.00 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(80,10)$} & 98.53 & 0.00 \\
& 600 & & 99.90 & 99.90 & \multicolumn{1}{c|}{99.50} & & 89.37 & 99.32 & 98.27 & & 99.80 & \multicolumn{1}{c|}{99.05} & & 84.98 & \multicolumn{1}{c|}{0.00} & & 98.39 & 97.99 & & 99.90 & \multicolumn{1}{c|}{0.00} & & 98.13 & 0.00 \\
& 900 & & 99.85 & 99.90 & \multicolumn{1}{c|}{99.30} & & 89.15 & 99.53 & 98.81 & & 99.75 & \multicolumn{1}{c|}{98.65} & & 85.63 & \multicolumn{1}{c|}{0.00} & & 99.15 & 98.49 & & 99.95 & \multicolumn{1}{c|}{0.00} & & 98.79 & 0.00 \\ \cline{2-25}
& 300 & \multirow{3}{*}{$(80,80)$} & 100.00 & 100.00 & \multicolumn{1}{c|}{99.70} & \multirow{3}{*}{$(160,10)$} & 87.47 & 99.53 & 99.06 & \multirow{3}{*}{$(80,80)$} & 99.95 & \multicolumn{1}{c|}{99.70} & \multirow{3}{*}{$(160,10)$} & 84.41 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(10,160)$} & 99.15 & 98.24 & \multirow{3}{*}{$(80,80)$} & 100.00 & \multicolumn{1}{c|}{0.00} & \multirow{3}{*}{$(160,10)$} & 98.83 & 0.00 \\
& 600 & & 99.95 & 99.95 & \multicolumn{1}{c|}{99.85} & & 89.06 & 99.74 & 99.11 & & 99.90 & \multicolumn{1}{c|}{99.65} & & 86.38 & \multicolumn{1}{c|}{0.00} & & 99.25 & 98.79 & & 100.00 & \multicolumn{1}{c|}{0.00} & & 98.48 & 0.00 \\
& 900 & & 100.00 & 100.00 & \multicolumn{1}{c|}{99.80} & & 92.05 & 100.00 & 99.38 & & 100.00 & \multicolumn{1}{c|}{99.70} & & 87.72 & \multicolumn{1}{c|}{0.00} & & 99.35 & 99.04 & & 100.00 & \multicolumn{1}{c|}{0.00} & & 98.79 & 0.00 \\ \hline\hline
\end{tabular}
}
\end{table}
\end{landscape}
\begin{landscape}
$ $\\
\begin{figure}[htbp]
\centering
\centerline{\includegraphics[width= 23cm]{plot/est-20250716.png}}
\caption{The lineplots for the averages of estimation errors $\varpi^2({\mathbf A},\hat{{\mathbf A}})$ and $\varpi^2({\mathbf B},\hat{{\mathbf B}})$ based on 2000 repetitions in Scenarios R1--R3. The legend is defined as follows: \textup{(i)} our proposed method ($\color{black}{-\bullet-}$), \textup{(ii)} the CP-refined method of \cite{chang2023modelling} ($\color{red}{-\blacktriangle-}$), \textup{(iii)} the cPCA of \cite{han2024cp} ($\color{blue}{-\blacksquare-}$), and \textup{(iv)} the HOPE of \cite{han2024cp} ($\color{green}{-+-}$).}
\label{fig:est-case1}
\end{figure}
\end{landscape}
\begin{landscape}
$ $\\
$ $\\
\begin{figure}[htbp]
\centering
\centerline{\includegraphics[width= 23cm]{plot/forecast_20250715_1step.png}}
\caption{The lineplots for the averages of one-step ahead forecast RMSE based on 2000 repetitions in Scenarios R1--R3. The legend is defined as follows: \textup{(i)} our proposed prediction method ($\color{black}{-\bullet-}$), \textup{(ii)} the prediction method introduced in \cite{chang2023modelling} with $(\hat{{\mathbf A}},\hat{{\mathbf B}})$ selected as our proposed estimate ($\color{blue}{-\blacktriangle-}$), \textup{(iii)} the prediction method introduced in \cite{chang2023modelling} with $(\hat{{\mathbf A}},\hat{{\mathbf B}})$ selected as the CP-refined estimate ($\color{red}{-\blacksquare-}$).}
\label{fig:fore-case1}
\end{figure}
\end{landscape}
\begin{figure}[htbp]
\centerline{\includegraphics[width= 14cm]{plot/app-timeseries.png}}
\caption{The time series plots of 100 market-adjusted returns formed on different levels of size (by rows) and book equity to market equity ratio (by columns). The horizontal axis represents time and the vertical axis represents the monthly returns.}
\label{fig:app-timeseries}
\end{figure}
\begin{figure}[htbp]
\centering
\subfigure[Proposed method ($\hat{{\mathbf M}}_1$)]{\includegraphics[width=0.32\textwidth]{plot/a1.png}}
\subfigure[Proposed method ($\hat{{\mathbf M}}_2$)]{\includegraphics[width=0.32\textwidth]{plot/a2.png}}
\subfigure[CP-refined method]{\includegraphics[width=0.32\textwidth]{plot/a3.png}}
\vspace{1em}
\subfigure[FAC (column matrix)]{\includegraphics[width=0.33\textwidth]{plot/a4.png}} \hspace{2em}
\subfigure[FAC (row matrix)]{\includegraphics[width=0.33\textwidth]{plot/a5.png}}
\caption{Scree plots for our proposed method, the CP-refined method and the FAC based on $\{{\mathbf Y}_t\}_{t=1}^{456}$. The black solid line represents the eigenvalues, while the blue dashed line indicates the ratios of adjacent eigenvalues. }
\label{fig: app-eigenall}
\end{figure}
\begin{landscape}
$ $\\
$ $\\
\begin{table}[htbp]
\footnotesize
\centering
\caption{The averages and standard deviations (in parentheses) of one-step and two-step ahead forecasting errors of the market-adjusted returns from January, 2002 to December, 2021. Panel A and Panel B represent the market-adjusted returns obtained by demean filtering and CAPM filtering, respectively.
}
\label{table:app-forecast}
\begin{tabular}{ccccccccccccccc}
\hline\hline
& Proposed & CP-refined(1) & CP-refined(2) & cPCA(1) & cPCA(2) & HOPE(1) & HOPE(2) & FAC(1,1) & FAC(2,1) & TOPUP & TIPUP & MAR & TS-PCA & UniAR \\ \hline
\multicolumn{15}{l}{Panel A: CAPM filtering} \\ \hline
\multicolumn{15}{c}{one-step ahead forecast} \\
{rRMSE} & \textbf{3.4302} & 3.4485 & 3.4408 & 3.4402 & 3.4361 & 3.4423 & 3.4373 & 3.4482 & 3.4610 & 3.4597 & 3.4669 & 3.4669 & 3.4710 & 3.4895 \\
& (\textbf{1.5062}) & (1.5223) & (1.4992) & (1.5254) & (1.5154) & (1.5299) & (1.5163) & (1.5226) & (1.5271) & (1.5277) & (1.5303) & (1.5364) & (1.5028) & (1.5099) \\
{rMAE} & \textbf{2.6218} & 2.6424 & 2.6325 & 2.6340 & 2.6294 & 2.6340 & 2.6307 & 2.6433 & 2.6541 & 2.6534 & 2.6566 & 2.6600 & 2.6579 & 2.6711 \\
& (\textbf{1.0522}) & (1.0746) & (1.0468) & (1.0712) & (1.0561) & (1.0755) & (1.0569) & (1.0751) & (1.0791) & (1.0763) & (1.0764) & (1.0849) & (1.0578) & (1.0601) \\ \hline
\multicolumn{15}{c}{two-step ahead forecast} \\
rRMSE & \textbf{3.4297} & 3.4488 & 3.4378 & 3.4436 & 3.4362 & 3.4447 & 3.4370 & 3.4505 & 3.4627 & 3.4602 & 3.4641 & 3.4401 & 3.4626 & 3.4873 \\
& (\textbf{1.5003}) & (1.5238) & (1.4972) & (1.5432) & (1.5178) & (1.5477) & (1.5161) & (1.5254) & (1.5291) & (1.5303) & (1.5331) & (1.5162) & (1.4947) & (1.5111) \\
rMAE & \textbf{2.6241} & 2.6444 & 2.6317 & 2.6378 & 2.6327 & 2.6374 & 2.6328 & 2.6456 & 2.6558 & 2.6538 & 2.6552 & 2.6353 & 2.6545 & 2.6705 \\
& (\textbf{1.0526}) & (1.0767) & (1.0494) & (1.0896) & (1.0684) & (1.0942) & (1.0668) & (1.0774) & (1.0805) & (1.0789) & (1.0823) & (1.0694) & (1.0575) & (1.0640) \\ \hline
\multicolumn{15}{l}{Panel B: Demean filtering} \\ \hline
\multicolumn{15}{c}{one-step ahead forecast} \\
rRMSE & 3.4905 & 3.5146 & 3.4947 & 3.5293 & 3.4987 & 3.5255 & \textbf{3.4862} & 3.5057 & 3.5165 & 3.5268 & 3.5303 & 3.5154 & 3.5244 & 3.5470 \\
& (1.5698) & (1.5829) & (1.5725) & (1.5883) & (1.5878) & (1.5876) & (\textbf{1.5824}) & (1.5952) & (1.5884) & (1.5826) & (1.5891) & (1.6093) & (1.6005) & (1.5789) \\
rMAE & 2.6676 & 2.6910 & 2.6718 & 2.7047 & 2.6741 & 2.7008 & \textbf{2.6622} & 2.6834 & 2.6923 & 2.7022 & 2.7036 & 2.6923 & 2.6999 & 2.7143 \\
& (1.1133) & (1.1288) & (1.1250) & (1.1333) & (1.1273) & (1.1327) & (\textbf{1.1182}) & (1.1402) & (1.1365) & (1.1283) & (1.1319) & (1.1597) & (1.1586) & (1.1250) \\ \hline
\multicolumn{15}{c}{two-step ahead forecast} \\
rRMSE & 3.4869 & 3.5142 & 3.4912 & 3.5239 & 3.4920 & 3.5208 & \textbf{3.4700} & 3.5018 & 3.5105 & 3.5269 & 3.5294 & 3.4959 & 3.5124 & 3.5413 \\
& (1.5685) & (1.5863) & (1.5721) & (1.5926) & (1.5939) & (1.5943) & (\textbf{1.5954}) & (1.5954) & (1.5894) & (1.5899) & (1.5961) & (1.6057) & (1.5881) & (1.5817) \\
rMAE & 2.6674 & 2.6922 & 2.6703 & 2.7013 & 2.6721 & 2.6982 & \textbf{2.6511} & 2.6805 & 2.6885 & 2.7036 & 2.7038 & 2.6765 & 2.6919 & 2.7130 \\
& (1.1170) & (1.1358) & (1.1237) & (1.1408) & (1.1391) & (1.1422) & (\textbf{1.1376}) & (1.1403) & (1.1368) & (1.1367) & (1.1406) & (1.1539) & (1.1483) & (1.1323) \\ \hline\hline
\end{tabular}
\end{table}
\end{landscape}
\newpage