Identification and estimation for matrix time series CP-factor models

Identification and Estimation for Matrix Time Series CP-factor Models

\if11 { \spacingset{1.25}

\affil[1]{\it Joint Laboratory of Data Science and Business Intelligence, Southwestern University of Finance and Economics, Chengdu, China} \affil[2]{\it State Key Laboratory of Mathematical Sciences, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing, China} \affil[3]{\it Department of Statistics, The London School of Economics and Political Science, London, U.K.}

\setcounter{Maxaffil}{0}

} \fi \if01 {

} \fi

\spacingset{1.5}

{\sl Keywords}: CP-decomposition; dimension-reduction; matrix time series; non-orthogonal joint diagonalization.

\spacingset{1.69} {0.2\baselineskip} {0.2\baselineskip} {0.2\baselineskip} {0.2\baselineskip}

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, wang2019factor, chen2019constrained, chen_chen_2022, and 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 lam2012factor and 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, chang2023modelling and 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

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, 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 (ref), 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}) $, 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 han2022tensor, han2024cp, the estimation procedure proposed in chang2023modelling is based on solving some generalized eigenequations and requires no iterations. Note that the incoherence conditions imposed in han2022tensor and 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 (ref) for matrix time series without imposing the condition ${\rm rank}({\mathbf A})=d={\rm rank}({\mathbf B})$. Let

Then $1\le d_1, d_2 \le d$. As the CP-decomposition for 3-way tensors often exhibits rank-deficient factor loading matrices kolda2009tensor, i.e., in model (ref) 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)). The key idea is to transform the $p\times q$ matrix CP-factor model ((ref)) 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)). 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)). 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 chang2023modelling can be formulated as the product of the rate of our new estimator and the inverse of an eigen-gap (See Remark (ref)). 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) gives preliminaries of the matrix CP-factor model (ref). A general identification strategy for the matrix CP-factor model is presented in Section (ref). Section (ref) provides a one-pass estimation procedure for $(d_1,d_2,d,{\mathbf A},{\mathbf B})$. Section (ref) gives a unified prediction approach for the matrix CP-factor model. We investigate the associated theoretical properties of the proposed method in Section (ref). Numerical results with simulation studies and real data analysis are given in Section (ref). The R-function CP_MTS for implementing our newly proposed method is available publicly in the HDTSA package chang2024hdtsa. All technical proofs and some additional simulation studies are relegated in the supplementary material.

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}$.

Preliminary

Recall that, in the matrix CP-factor model (ref), ${\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 (ref) can be equivalently represented as

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)(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)(ii) holds automatically.

When $d_1=d_2=d$, chang2023modelling provides a one-pass estimator for $({\mathbf A},{\mathbf B})$ by solving some generalized eigenequations defined by the matrices

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 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 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 kolda2009tensor, i.e., in model ((ref)) 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.

Identification of $({\mathbf A},{\mathbf B})$

We need to identify in model (ref) 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.

A reduced model

For a prescribed integer $K>1$ and $\boldsymbol{\Sigma}_{{\mathbf Y}, \xi}(k)$ specified in (ref), define

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 (ref) and Condition (ref) that

Let ${\mathbf G}=\sum_{k=1}^K {\mathbf g}_{k}{\mathbf g}_{k}^{{\mathrm{\scriptscriptstyle \top} }}$. Proposition (ref) shows that $d_1$ and $d_2$ can be identified, repsectively, by ${\rm rank}({\mathbf M}_1)$ and ${\rm rank}({\mathbf M}_2)$.

The conditions required in Proposition (ref) 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) 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)(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)(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) proposes some sufficient conditions such that $ {\rm rank}({{\mathbf G}})=d$, which make Proposition (ref) 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}$.

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) shows that ${\rm rank}({\mathbf G}) =d$. Combining with Proposition (ref), it is reasonable to assume Condition (ref), 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$.

Now perform the spectral decomposition for ${\mathbf M}_1$ and ${\mathbf M}_2$:

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

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)), 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 (ref) and (ref) that

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 (ref) for the matrix time series ${\mathbf Y}_t$. We will identify $({\mathbf U},{\mathbf V})$ based on this reduced model.

A vector factor model

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 (ref) that

This is the standard factor model for vector time series considered by lam2012factor and 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)(i), we know the dimension of the factor loading space ${\cal M}({\mathbf V}\odot{\mathbf U})$ in (ref) is $d$, as ${\rm rank}({\mathbf V}\odot {\mathbf U}) = {\rm rank}({\mathbf B}\odot{\mathbf A}) =d$. Using the techniques developed in 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

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 (ref) 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 (ref), different selections of ${\mathbf W}$ will lead to different $\boldsymbol{\Theta}$. As we will show in Section (ref), for any given ${\mathbf W}$, the associated rotation matrix $\boldsymbol{\Theta}$ can be uniquely identified up to the column reflection and permutation indeterminacy. Write

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)) and (ref) 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 (ref) is uniquely identified, which can be used to identify $({\mathbf U},{\mathbf V})$. See Proposition (ref) for details.

A non-orthogonal joint diagonalization

For given ${\mathbf W}$ in (ref), Proposition (ref) implies that the task of identifying $({\mathbf U},{\mathbf V})$ boils down to identifying $\boldsymbol{\Theta}$ specified in (ref) such that ${\mathbf C}_1,\ldots,{\mathbf C}_d$ defined in (ref) satisfying ${\rm rank}({\mathbf C}_\ell)=1$ for each $\ell \in [d]$. By (ref), it holds that

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 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 (ref), for given ${\mathbf W}$ in (ref), we know $\boldsymbol{\Theta} = (\theta_{i,j})$ is the solution of

This is a set of quadratic equations. Consider a $(d_1^2d_2^2)\times d(d+1)/2$ matrix

Proposition (ref) is instrumental in solving those quadratic equations.

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)), it holds that

By Proposition (ref)(ii), we have

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 (ref), 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 (ref) 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) 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 pham2001blind, the fast Frobenius diagonalization of ziehe2004fast, and the quadratic diagonalization of vollgraf2006quadratic.

Proposition (ref) 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) 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}$.

By Propositions (ref) and (ref), if ${\rm rank}(\boldsymbol{\Omega})=d(d-1)/2$, then ${\mathbf U}$ and ${\mathbf V}$ specified in (ref) 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) 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.

Estimation

Based on Section (ref), we can estimate $d$ and $({\mathbf a}_{\ell}, {\mathbf b}_{\ell})$ for $\ell \in[d]$ via the following five steps:

Steps 4 and 5 are straightforward. More details of Steps 1--3 are given, respectively, in Sections (ref)--(ref). 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}$.

Estimating $d_1, \, d_2, \, {\mathbf P}$ and ${\mathbf Q}$

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 (ref), we put

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

We set $\delta_1>0$ in (ref) 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), we can estimate $d_1$ and $d_2$ by the eigenvalue-ratio based method Chang2015 as follows:

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 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) in Section (ref) for details. In contrast, the eigenvalue-ratio based method proposed in 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

Estimating $d$ and ${\mathbf W}_1, \ldots, {\mathbf W}_d$

Based on (ref) and (ref), we can reformulate (ref) 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 lam2011estimation, lam2012factor and Chang2015. To do this, we put

with a prescribed integer $\tilde{K}\ge 1$ and

where $T_{\delta_2}(\cdot) $ is a truncation operator with the threshold level $\delta_{2} \ge 0$, and

with $\bar{\vec{{\mathbf Y}}}= n^{-1}\sum_{t=1}^{n}\vec{{\mathbf Y}}_{t}$. Analogous to (ref), we can estimate $d$ as

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.

Estimating $\mathbf{\Theta}$ via joint diagonalization

For 4-way tensor $\boldsymbol{\Psi}(\cdot,\cdot)$ defined in Section (ref), we define a $(\hat d_1^2\hat d_2^2)\times \hat{d}(\hat d+1)/2$ matrix $\hat \boldsymbol{\Omega}$ as follows:

which is an estimate of $\boldsymbol{\Omega}$ defined as in (ref). Let

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), 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 (ref).

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)), the accuracy of its estimator depends on the choice of $\{ {\mathbf h}_m\}_{m=1}^{d}$ sensitively. Motivated by Proposition (ref) at the end of this section, a good choice is to rotate the basis vectors $\{\tilde{{\mathbf h}}_m\}_{m=1}^{\hat{d}}$ in (ref) 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

where

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) 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 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 (ref) 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 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

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) 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

where

for some $d$-dimensional vector $(\phi_1, \ldots, \phi_d)^{\mathrm{\scriptscriptstyle \top} } \ne {\mathbf 0}$ such that ${\mathbf H}$ is invertible.

Prediction

Given observations $\{{\mathbf Y}_t\}_{t=1}^n$, we can also use the matrix CP-factor model (ref) 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 (ref). 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}$. 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 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), and use the same idea to predict ${\mathbf Y}_{n+h}$.

As shown in Proposition (ref), 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 (ref) without any assumption on the relationship among $d_1$, $d_2$ and $d$. By (ref) and (ref), we have

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) and (ref), we define

Proposition (ref) in Section (ref) 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). 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 (ref) and Proposition (ref), 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). 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) again, by (ref), 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) and (ref), 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 (ref) 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) 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 chang2023modelling with selecting $(\hat{{\mathbf A}},\hat{{\mathbf B}})$ as our proposed estimate of $({\mathbf A},{\mathbf B})$ specified in Section (ref). However, if the factor loading matrices ${\mathbf A}$ and ${\mathbf B}$ cannot be uniquely identified, our newly proposed prediction method outperforms that of chang2023modelling.

Asymptotic properties

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

where $\mathcal{F}_{r}^{s}$ is the $\sigma$-filed generated by $\{{\mathbf Y}_t: r\le t\le s\}$. Write

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 (ref) and

where $\bar{\vec{{\mathbf Z}}}=n^{-1}\sum_{t=1}^{n}\vec {\mathbf Z}_{t}$. Following the arguments in Chang2015, we can identify $d$ as $d={\rm rank}({\mathbf M})$, and select $\vec{\mathbf W}_1,\ldots, \vec{\mathbf W}_d$ involved in (ref) 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.

Condition (ref) 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) is also used in chang2023modelling, which is a common assumption in the literature on ultrahigh-dimensional data analysis. See chang2023modelling for the discussion of their validity. We impose Condition (ref)(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)(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 Chang2018, we can show that Condition (ref)(ii) holds for certain $(s_1,s_2,s_3,s_4)$. Let

Theorem (ref) 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.

Proposition (ref) states the asymptotic performance of $\hat{{\mathbf P}}$, $\hat{{\mathbf Q}}$ and $\hat{{\mathbf W}}$.

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 (ref), 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).

To present Theorem (ref), we need to introduce some notation first. For $({\mathbf E}_1,{\mathbf E}_2,{\mathbf E}_3)$ specified in Proposition (ref), 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 (ref) but with replacing ${\mathbf W}$ by $\breve{{\mathbf W}}$. Following the discussions of Propositions (ref)--(ref), 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) in the supplementary material for the proof of Lemma (ref), we know ${\rm rank}(\breve{\boldsymbol{\Omega}}) ={\rm rank}(\boldsymbol{\Omega})$. By Proposition (ref)(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.

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). Recall ${\mathbf A} = ({\mathbf a}_1,\ldots,{\mathbf a}_{d})$ and ${\mathbf B} = ({\mathbf b}_1,\ldots,{\mathbf b}_{d})$. Theorem (ref) 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.

Theorem (ref) requires $(\hat{d}_1, \hat{d}_2, \hat{d}) = (d_1,d_2,d)$. By Theorem (ref), 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

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

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) 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), we have

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)$.

Numerical studies

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), and the analysis of the simulation results is presented in Section (ref). The real data analysis is given in Section (ref).

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)$:

We follow chang2023modelling to specify $\xi_t$ involved in (ref). 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 (ref) and (ref), and set $c_{1,n} = c_{2,n} = c_{3,n} = \hat\sigma_0 n^{-1}$ in (ref) and (ref). We also choose $K = 20$ and $\tilde{K} = 10$ with $K$ and $\tilde{K}$ given in (ref) and (ref), respectively. Here, using a relatively large value for $K$ is to ensure that ${\mathbf M}_1$ and ${\mathbf M}_2$ defined in (ref) 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) and (ref). As shown in (ref), $\tilde{K}$ is the number of lags used in the methods of lam2011estimation, lam2012factor and 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)--(ref) of the supplementary material for more details. As mentioned in Section (ref), 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 (ref). 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 chang2023modelling, and the cPCA and the HOPE methods proposed by han2024cp with the recommended tuning parameter $h = 1$ therein. All simulations are implemented in R. Our proposed method is available in R-package HDTSA, which is implemented by calling the R-function CP_MTS with setting method = `CP.Unified'. The CP-refined method of 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.

Simulation results

We first consider the finite-sample performance of the estimation $(\hat{d}_1, \hat{d}_2,\hat{d})$ given in (ref) and (ref). Note that the CP-refined method of 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 (ref), 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) 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) reports the averages of the estimation errors $\varpi^2({\mathbf A},\hat{{\mathbf A}})$ and $\varpi^2({\mathbf B},\hat{{\mathbf B}})$ defined in (ref) and (ref) 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 (ref) 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)--(ref) in the supplementary material.

Next, we evaluate the finite-sample performance of our proposed prediction method introduced in Section (ref). We generate a sequence $\{{\mathbf Y}_t\}_{t=1}^{n+m+1}$ defined in Section (ref) 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 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), or (ii) the CP-refined estimate of $({\mathbf A}, {\mathbf B})$ given in 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) 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 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 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 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 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) in the supplementary material for details.

Real data analysis

In this section, we illustrate the proposed method for the matrix CP-factor model (ref) 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 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 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) 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) 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) 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), 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)(a) and (ref)(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:

For each $s \in [240]$, the one-step ahead forecasting performance is evaluated by the $\textup{rRMSE}(s)$ and $\textup{rMAE}(s)$ defined as

For the two-step ahead forecast, we can evaluate it by the associated $\textup{rRMSE}(s)$ and $\textup{rMAE}(s)$ analogously. Table (ref) 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), 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.

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).

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.

\spacingset{0.95}\selectfont