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.
77,570 characters
Reduced-Rank Autoregressive Models for Matrix Time Series
\maketitle
\noindent {\bf Abstract.} Matrix time series is a series of matrix data observed over time. Analytical tools for such time series is needed in many applications in finance, economics, engineering and many other fields. To avoid the use of vectorization of the matrices which loses the column and row information, and the vector autoregression framework in traditional time series analysis, \cite{chen2021autoregressive} proposed the Matrix Autoregressive (MAR) Model. The model maintains and utilizes the matrix structure, leading to a substantial dimensional reduction and admitting explicit interpretations, comparing with the vector autoregressive model on the vectorized data. However, the MAR model still encounters difficulties in dealing with large dimensional matrix time series as the coefficient matrices in MAR models are also large. In this paper we propose to achieve further dimension reduction through reduced-rank constraints of the coefficient matrices in the MAR model. Estimation and rank determination procedures are studied. Theoretical investigation and empirical examples show that the reduced-rank constraint can achieve higher statistical efficiency than the MAR model.
\bigskip
\noindent {\bf Keywords:} Forecasting; Matrix time series; Rank determination; Reduced-rank regression
\newpage
\section{Introduction}
\label{sec:intro}
\nocite{anderson:2013,horn2012matrix,chen2012sparse,chen2021autoregressive,chen2019factor}
Observations in matrix and tensor (multi-dimensional array) forms have been generated and collected more and more abundantly in many fields including biological/medical research, economics, engineering, finance, signal processing, social sciences etc. In response to the urgent need of analytical tools for analyzing such type of data in various applications, many optimization and statistical methods/procedures have been proposed and studied. Similar to the use of matrix decomposition for analysis of vector observations, tensor decomposition and estimation methods play a principal role in analyzing matrix/tensor data \citep{cichocki2015tensor,cichocki2009nonnegative,de2000best,de2000multilinear,sidiropoulos2017tensor,anandkumar2014tensor,de2008tensor}.
In many applications, the matrices are observed through time, and hence form a matrix-valued time series. Although it is possible and perhaps convenient to treat time as another mode and apply the tensor methods to such a three-way tensor, the time dimension is intrinsically different, and the temporal dependence requires careful modeling and analysis to aid practitioners on acquiring a diagnostic understanding of the dynamics and making reliable forecasts. It has been witnessed that multilinear models can reduce the dimension and improve the estimation stability for matrix/tensor data \citep{ding2018matrix,raskutti2019convex, zhao:2014,zhou:2013}.
For dependent data, there has been many recent works on factor models of matrix/tensor time series, see \cite{wang:2018,chen2019factor,han2020tensor,han2020number,gao2021two} among others. On the other hand, \cite{hoff2015multilinear} pioneered in suggesting the multilinear model for longitudinal tensor data. \cite{chen2021autoregressive} proposed the matrix autoregressive model (MAR), which retains the matrix form of the data and specifies the autoregressive relationship through a bilinear matrix product.
Besides the interpretations adherent to its matrix and bilinear form, the MAR also reduces the model complexity significantly, comparing to the approach of concatenating the matrix observation into a long vector and then using the traditional vector autoregressive (VAR) model \citep{hannan:1970,lutkepohl:2005,tsay:2014,tiao:1981}.
When the matrix observations are themselves of large dimensions, the MAR model still involves a large number of parameters. It is desirable and sometimes necessary to reduce the dimension even further. In this paper, we propose the reduced-rank matrix autoregressive model (RRMAR), which assumes the form of MAR, but requires in addition that the coefficient matrices have ranks smaller than their dimensions.
We note that another natural approach to reducing the model complexity is to impose sparsity on the MAR. This approach is closely related to recent works on sparse VAR models \citep{loh2011high,basu2015regularized,melnyk2016estimating,davis2016sparse,han:2015,kock:2015,lin2017regularized,nicholson2017varx}. In addition, \cite{basu2019low} and \cite{lin2020regularized} considered additional low rank constraints on the coefficient matrices, \cite{hall2018learning} introduced the generalized VAR model, and \cite{han2020sparse} focused on the nonlinear sparse VAR model. \cite{ghosh2019high} and \cite{ghosh2020strong} considered the high dimensional VAR from a Bayesian perspective.
In contrast to the aforementioned works based on sparsity, the thrust of the present paper hinges upon the low rank structure of the coefficient matrices in MAR. As will be elaborated in Section~\ref{sec:rrmar}, the low rank matrices in RRMAR continue to admit natural interpretations, and lead to a greater dimension reduction compared to MAR. It also relates to and provides a generative mechanism for the dynamic matrix factor models of \cite{wang:2018}, and has a close connection to the hierarchical factor models.
The proposed model and estimation procedure are related to the reduced-rank regression \citep{anderson:1951,izenman:1975,reinsel:1998}.
We consider two estimators: one based on least squares (RR.LS) and one based on maximum likelihood (RR.CC). It is worth pointing out that they correspond to two different algorithms for vector reduced-rank regression: one minimizing the trace of the sample covariance matrix, and the other the determinant, where the latter also corresponds to the canonical correlation analysis (see for example \cite{reinsel:1998} for more details). The likelihood-based RR.CC is indeed the maximum likelihood estimator if the covariance tensor of the error matrix has the form of a tensor product of two covariance matrices. Even when this assumption does not hold, the RR.CC nevertheless can be viewed as an estimator obtained together with a regularized estimation of the covariance tensor, and can still potentially lead to superior performances over the RR.LS.
Here we shall emphasize two significant differences between our model and the classical reduced-rank regression. First, the observations are in matrix form, and the model takes a bilinear form. Second, the algorithms requires running reduced-rank least squares/maximum likelihood iteratively. We develop central limit theorems for the estimators of the coefficient matrices, as well as their singular vectors. The bilinear form of the matrix model also makes the analysis substantially different from the vector case.
The estimation procedures depend on the ranks of the two coefficient matrices in the RRMAR model. We propose to use information criterion based procedures to identify the ranks of these matrices. Since two ranks are to be determined, a thorough search over all possible pairs of ranks can be very costly, so we also introduce procedures to select the two ranks separately. Asymptotic consistency of these selection procedures are established.
The rest of this article is organized as follows. The RRMAR model is introduced in Section~\ref{sec:rrmar}, together with its basic properties, interpretations and connections with other models. In Section~\ref{sec:est} we propose two estimators, RR.LS and RR.CC. Asymptotic distributions are provided for both of them in Section~\ref{sec:asymp}, as well as corresponding estimators of the leading singular vectors of the coefficient matrices. The model/rank selection procedures based on information criterion are introduced in Section~\ref{sec:rank}, with their consistency properties. We use an extensive numerical study and an example in finance to demonstrate the performances of the proposed models and estimators in Section~\ref{sec:num}. All the proofs are collected in the Appendix.
\subsection{Notations}
We gather the notations and the definitions of some special matrices in this section.
We use $\|\cdot\|_F$ to denote the Frobenius norm of a matrix, and $\rho(\cdot)$ the spectral radius. We use $\otimes$ to denote the Kronecker product, and $\circ$ the (point-wise) Hadamard product of two matrices. The notation $\boldsymbol{1}_k$ stands for a $k$-dimensional vector with all entries equal to one. For any matrix $\boldsymbol{M}$, we use $\boldsymbol{M}[i,]$ and $\boldsymbol{M}[,j]$ to denote its $i$-th row and $j$-th column respectively. The column space of $\boldsymbol{M}$ is denoted by $\mathrm{col}(\boldsymbol{M})$. The matrix vectorization, denoted by $\mathrm{vec}(\cdot)$, turns a matrix into a vector by stacking its columns.
For any positive integer $p$, let $\boldsymbol{e}_{p,j}\in\mathbb{R}^p$ be the $j$-th base vector whose $j$-the entry is 1, and others zero. For any two positive integers $p$ and $q$, let $\boldsymbol{J}_{p,q}$ be the $(pq)\times(pq)$ permutation matrix defined as
\begin{equation}
\label{eq:J}
\boldsymbol{J}_{p,q} = \left[\boldsymbol{I}_q\otimes\boldsymbol{e}_{p,1},\boldsymbol{I}_q\otimes\boldsymbol{e}_{p,2},\ldots,\boldsymbol{I}_q\otimes\boldsymbol{e}_{p,p}\right].
\end{equation}
The permutation $\boldsymbol{J}_{p,q}$ does the following: for any $p\times q$ matrix $\boldsymbol{M}$, $ \boldsymbol{J}_{p,q}\mathrm{vec}(\boldsymbol{M}')=\mathrm{vec}(\boldsymbol{M}).$
In other words, $\boldsymbol{J}_{p,q}$ connects the vectorization of a matrix and its transpose.
Let $\boldsymbol{L}_p$ be the $p^2\times p$ matrix whose $j$-th column is given by $\boldsymbol{e}_{p,j}\otimes\boldsymbol{e}_{p,j}$, i.e.
\begin{equation}
\label{eq:L}
\boldsymbol{L}_p=[\boldsymbol{e}_{p,1}\otimes\boldsymbol{e}_{p,1},\boldsymbol{e}_{p,2}\otimes\boldsymbol{e}_{p,2},\ldots,\boldsymbol{e}_{p,p}\otimes\boldsymbol{e}_{p,p}].
\end{equation}
For any $p\times p$ matrix $\boldsymbol{M}=(m_{jk})$, the following operation extracts its diagonals:
\begin{equation*}
\boldsymbol{L}_p'\mathrm{vec}(\boldsymbol{M})=(m_{11},\ldots,m_{pp})',
\end{equation*}
and furthermore,
\begin{equation*}
\mathrm{vec}^{-1}[\boldsymbol{L}_p\boldsymbol{L}_p'\mathrm{vec}(\boldsymbol{M})]=\mbox{diag}(\boldsymbol{M}),
\end{equation*}
where $\mbox{diag}(\boldsymbol{M})$ is the $p\times p$ diagonal matrix keeping $\boldsymbol{M}$'s diagonal elements.
\section{Reduced-Rank MAR Model}
\label{sec:rrmar}
The {\it reduced-rank matrix autoregressive model} (RRMAR) takes the
form
\begin{equation}
\label{eq:marrr}
\boldsymbol{X}_t=\boldsymbol{A}_1\boldsymbol{X}_{t-1}\boldsymbol{A}_2'+\boldsymbol{E}_t,
\end{equation}
where $\boldsymbol{A}_i$ are $d_i\times d_i$ autoregressive coefficient
matrices {of ranks $k_i \leq d_i$}, and $\boldsymbol{E}_t\in\mathbb{R}^{d_1\times d_2}$ is a matrix
white noise. It {is}
the same as the MAR model proposed
by \cite{chen2021autoregressive},
except for the additional low rank assumption that $\operatorname{rank}(\boldsymbol{A}_i)=k_i \leq d_i$, for
$i=1,2$. It is worth observing that the number of parameters to
determine $\boldsymbol{A}_i$ under the rank constraint is $d_i^2-(d_i-k_i)^2=(2d_i-k_i)k_i$ \citep[see for example][]{camba2003tests, reinsel:1998}, as
opposed to $d_i^2$ for the unconstrained $\boldsymbol{A}_i$, and the former can
be much smaller if $k_i\ll d_i$. We also assume that $\|\boldsymbol{A}_1\|_F=1$,
so that $\boldsymbol{A}_1$ and $\boldsymbol{A}_2$ are identified up to a sign change. To guarantee that the model \eqref{eq:marrr} is causal and stationary, we require that $\rho(\boldsymbol{A}_1)\cdot\rho(\boldsymbol{A}_2)<1$.
To better understand the implication of the low rank assumption, we
write $\boldsymbol{A}_i=\boldsymbol{A}_{il}\boldsymbol{A}_{ic}'$, where $\boldsymbol{A}_{il}$ and $\boldsymbol{A}_{ic}$ are
both $d_i\times k_i$ full rank matrices. The model \eqref{eq:marrr} is
then written as
\begin{equation}
\label{eq:marrr_fac}
\boldsymbol{X}_t=\boldsymbol{A}_{1l}\,\boxed{\boldsymbol{A}_{1c}'\boldsymbol{X}_{t-1}\boldsymbol{A}_{2c}}\,\boldsymbol{A}_{2l}'+\boldsymbol{E}_t.
\end{equation}
The boxed part $\boldsymbol{F}_{t}:=\boldsymbol{A}_{1c}'\boldsymbol{X}_{t-1}\boldsymbol{A}_{2c}$ is a
$k_1\times k_2$ matrix, which can be viewed as a composite and much
smaller version of the $d_1\times d_2$ matrix $\boldsymbol{X}_{t-1}$. The conditional
expectation of $\boldsymbol{X}_t$ given $\boldsymbol{X}_{t-1}$ is then given by loading on
$\boldsymbol{F}_{t}$ from left by $\boldsymbol{A}_{1l}$, and from right by
$\boldsymbol{A}_{2l}'$. The RRMAR model therefore provides a generating mechanism
for the matrix factor model $\boldsymbol{X}_t=\boldsymbol{A}_{l1}\boldsymbol{F}_t\boldsymbol{A}_{l2}' +\boldsymbol{E}_t$, introduced in \cite{wang:2018}, in which the factor process $\boldsymbol{F}_t$ is assumed to be latent and unobserved. In the RRMAR model in \eqref{eq:marrr_fac}, $\boldsymbol{F}_t$ depends on $\boldsymbol{X}_{t-1}$ hence is observed given the parameters. Due to this connection, we call $\boldsymbol{A}_{ic}$ the composition matrix, and $\boldsymbol{A}_{il}$ the loading matrix.
The RRMAR model is also related to the hierarchical factor models in the econometrics literature \citep{moench2013dynamic}. Specifically, let $\boldsymbol{F}^*_t = \boldsymbol{A}_{1c}'\boldsymbol{X}_{t-1}\boldsymbol{A}_{2}'$, then
\begin{equation*}
\boldsymbol{X}_t = \boldsymbol{A}_{1l}\boldsymbol{F}^*_t + \boldsymbol{E}_t,
\end{equation*}
which means that the $j$-th column of $\boldsymbol{X}_t$ follows a factor model with loading $\boldsymbol{A}_{1l}$, and factors $\boldsymbol{F}^*_t[,j]$, the $j$-th column of $\boldsymbol{F}^*_t$. In the next layer, we have
\begin{equation*}
\boldsymbol{F}_t^{*'} = \boldsymbol{A}_{2l} \boldsymbol{F}_t',
\end{equation*}
which says that the $d_2\times k_1$ factor matrix $\boldsymbol{F}_t^{*'}$ is further driven by a smaller factor matrix $\boldsymbol{F}_t'$ {(defined by the boxed part in \eqref{eq:marrr_fac})}, the $j$-th column of $\boldsymbol{F}_t^{*'}$ corresponding to the loading $\boldsymbol{A}_{2l}$, and factors $\boldsymbol{F}_t'[,j]$. Therefore, the model \eqref{eq:marrr} also gives a generating mechanism of a special instance of the hierarchical factor model.
\nocite{giannone2008nowcasting,diebold2008global}
\begin{remark}
{Again we emphasize that RRMAR is strictly an AR model, with the conditional mean of $\boldsymbol{X}_t$ given the past information solely depending on $\boldsymbol{X}_{t-1}$. As in all traditional AR models, we assume the noise process $\boldsymbol{E}_t$ is white (in time), though the elements in $\boldsymbol{E}_t$ are allowed to have (strong) correlations among them. It is possible to extend the model to have an ARMA form to allow non-white error processes, but such an ARMA model is extremely difficult to analyze and typically less useful in practice even for vector time series, due to various ambiguities. Although RRMAR can be written in a factor model form, strictly speaking, it is not a factor model, as it is {\it generative} in which $\boldsymbol{F}_t$ is based on the past information $\boldsymbol{X}_{t-1}$. On the other hand, a typical factor model is often {\it descriptive}, with the factor process $\boldsymbol{F}_t$ being latent and its estimator is typically a linear combination of the current observation $\boldsymbol{X}_t$.}
\end{remark}
There are many potential extensions of model \eqref{eq:marrr}. The first is the extension to RRMAR$(p)$ model:
\begin{equation}
\label{eq:marrrp}
\boldsymbol{X}_t=\sum_{j=1}^p\boldsymbol{A}_{j1}\boldsymbol{X}_{t-j}\boldsymbol{A}_{j2}'+\boldsymbol{E}_t,
\end{equation}
where all $\boldsymbol{A}_{ji}$ are of low ranks. The RRMAR$(p)$ model also extends the reduced-rank vector autoregressive models \citep{velu1986reduced, camba2003tests,al2019testing}. The second extension is more subtle. Although only the lag-1 observation $\boldsymbol{X}_{t-1}$ is involved on the right hand side of \eqref{eq:marrr}, there can be multiple terms in the form
\begin{equation}
\label{eq:marrrk}
\boldsymbol{X}_t=\sum_{j=1}^r\boldsymbol{A}_{1}^{(j)}\boldsymbol{X}_{t-1}\left(\boldsymbol{A}_{2}^{(j)}\right)'+\boldsymbol{E}_t.
\end{equation}
To see this extension more clearly, we take vectorization on both sides of \eqref{eq:marrrk}:
\begin{equation*}
\mathrm{vec}(\boldsymbol{X}_t) = \left(\sum_{j=1}^r \boldsymbol{A}_{2}^{(j)}\otimes\boldsymbol{A}_{1}^{(j)}\right) \mathrm{vec}(\boldsymbol{X}_{t-1}) + \boldsymbol{E}_t.
\end{equation*}
It is seen from the preceding equation that the RRMAR model \eqref{eq:marrr} amounts to restricting the coefficient matrix of the VAR(1) model to the form of a Kronecker product, and the model \eqref{eq:marrrk} is more flexible by representing the coefficient matrix as a sum of $r$ Kronecker products.
The extension to AR($p$) in \eqref{eq:marrrp} is quite straightforward, and the alternating estimation algorithms and rank determination procedures proposed below can be easily extended for this model, with heavier computational cost.
It can help to make the model more parsimonious by requiring (i) the matrices $\boldsymbol{A}_{ji}$ to have the same ranks across $1\leq j\leq p$, or (ii) the rank of $\boldsymbol{A}_{ji}$ decreases as the lag $j$ increases, or (iii) the matrices $\boldsymbol{A}_{ji}$ have the same column/row spaces across $j$.
The extension to \eqref{eq:marrrk} is more intricate,
as certain identifiability constraints are needed. We shall focus on model \eqref{eq:marrr} in the main part of the paper, and treat the extension \eqref{eq:marrrp} briefly in Appendix II, but leave the extension \eqref{eq:marrrk} for future studies. Finally, we add that the two extensions \eqref{eq:marrrp} and \eqref{eq:marrrk} can be combined to give a more comprehensive model.
\section{Estimation}
\label{sec:est}
Suppose a matrix time series $\{\boldsymbol{X}_t\}$ of length $T$ is observed. To
estimate the coefficient matrices $\boldsymbol{A}_i$, we propose to use the
alternating reduced-rank regression, updating one, while holding the
other fixed. Specifically, suppose $\boldsymbol{A}_2$ is given, we discuss how
to estimate $\boldsymbol{A}_1$. Recall that $\boldsymbol{A}[,j]$ denotes the $j$-th column
of a matrix $\boldsymbol{A}$. We also make the convention that $\boldsymbol{A}'[,j]$ denotes the
$j$-th column of $\boldsymbol{A}'$, i.e. the $j$-th row of $\boldsymbol{A}$ as a column
vector. The $j$-th column of the model equation \eqref{eq:marrr} is
\begin{equation*}
\boldsymbol{X}_t[,j]=\boldsymbol{A}_1\,\boxed{\boldsymbol{X}_{t-1}\boldsymbol{A}_2'[,j]} + \boldsymbol{E}_t[,j].
\end{equation*}
Since $\boldsymbol{A}_2$ is fixed, the preceding equation can be viewed as the
reduced-rank regression involving $(T-1)d_2$ sample units, where each
column $\boldsymbol{X}_t[,j]$ is a response vector, the boxed vector is the
covariate, and $\boldsymbol{A}_1$ is the coefficient matrix. In Section~\ref{sec:ils} we consider the estimation of $\boldsymbol{A}_1$ by least squares. On the other hand, under normality, the
classical reduced-rank regression minimizes the determinant of the
sample covariance matrix of the error vectors, under the rank
constraint, which is related and in fact equivalent to canonical
correlation analysis \citep{anderson:2013,reinsel:1998}. In Section~\ref{sec:icc} we introduce a special covariance structure of $\boldsymbol{E}_t$, under which we seek to estimate $\boldsymbol{A}_1$ by the Gaussian MLE.
{In this section we focus on the estimation of the coefficient matrices given the ranks $k_1$ and $k_2$. The determination of the ranks will be discussed in Section 5.}
\subsection{Alternating least squares}
\label{sec:ils}
The least squares estimators are solutions of
\begin{equation}
\label{eq:rrmin}
\min_{\boldsymbol{A}_1,\,\boldsymbol{A}_2: \; \operatorname{rank}(\boldsymbol{A}_1)=k_1,\,\operatorname{rank}(\boldsymbol{A}_2)=k_2}\sum_{t=2}^T\left\|\boldsymbol{X}_t-\boldsymbol{A}_1{\boldsymbol{X}_{t-1}\boldsymbol{A}_2'}\right\|_F^2
\end{equation}
We denote the least squares estimators by $\hat\boldsymbol{A}_i^{\hbox{\tiny ls}}$, and will refer to them as the RR.LS estimators.
Suppose $\boldsymbol{A}_2$ is given, the problem \eqref{eq:rrmin} becomes minimizing the trace of the sample covariance
matrix {of the residuals} under the rank constraint on $\boldsymbol{A}_1$:
\begin{align*}
& \min_{\boldsymbol{A}_1: \; \operatorname{rank}(\boldsymbol{A}_1)=k_1}\sum_{t=2}^T\left\|\boldsymbol{X}_t-\boldsymbol{A}_1{\boldsymbol{X}_{t-1}\boldsymbol{A}_2'}\right\|_F^2 \\
\Longleftrightarrow\quad &
\min_{\boldsymbol{A}_1: \; \operatorname{rank}(\boldsymbol{A}_1)=k_1}\operatorname{tr}\left[\sum_{t=2}^T\sum_{j=1}^{d_2}\left(\boldsymbol{X}_t[,j]-\boldsymbol{A}_1{\boldsymbol{X}_{t-1}\boldsymbol{A}_2'[,j]}\right)\left(\boldsymbol{X}_t[,j]-\boldsymbol{A}_1{\boldsymbol{X}_{t-1}\boldsymbol{A}_2'[,j]}\right)'\right].\nonumber
\end{align*}
Let $\boldsymbol{S}_{xx}=\sum_t \boldsymbol{X}_{t-1}\boldsymbol{A}_2'\boldsymbol{A}_2\boldsymbol{X}_{t-1}'$,
$\boldsymbol{S}_{yx}=\sum_t \boldsymbol{X}_t\boldsymbol{A}_2\boldsymbol{X}_{t-1}'$, and
$\boldsymbol{U}:=[U_1, U_2,\ldots, U_{k_1}]$, where
$U_j$ is the $j$-th leading normalized eigenvector of
$\boldsymbol{S}_{yx}\boldsymbol{S}_{xx}^{-1}\boldsymbol{S}_{xy}$. Then $\boldsymbol{A}_1$ can be updated as
\begin{equation*}
\check \boldsymbol{A}_1^{\hbox{\tiny ls}}=\boldsymbol{U}\boldsymbol{U}'\boldsymbol{S}_{yx}\boldsymbol{S}_{xx}^{-1},
\end{equation*}
see for example Equation (2.15) of \cite{reinsel:1998}. Given $\boldsymbol{A}_1$, an update of $\boldsymbol{A}_2$ can be similarly obtained. We therefore use the alternating least squares to find the minimizer of \eqref{eq:rrmin}.
\subsection{Alternating canonical correlation analysis}
\label{sec:icc}
The classical reduced-rank regression has also been situated under normality, leading to the Gaussian MLE of the coefficient matrix. To introduce the MLE for the RRMAR model, we need to assume that the covariance matrix $\Sigma_e$ of $\mathrm{vec}(\boldsymbol{E}_t)$ takes the form of a product
\begin{equation}
\label{eq:ocov}
\Sigma_e=\Sigma_2\otimes\Sigma_1,
\end{equation}
where $\Sigma_1$ and $\Sigma_2$ are $d_1\times d_1$ and $d_2\times d_2$ positive definite matrices respectively. This is equivalent to assuming $\boldsymbol{E}_t = \Sigma_1^{1/2}\boldsymbol{Z}_t\Sigma_2^{1/2}$, where $\boldsymbol{Z}$ has iid standard normal entries. This assumption allows us to separate the row and column dependence within the error matrix, with $\Sigma_1$ and $\Sigma_2$ corresponding to the row-wise and column-wise correlations among the entries of $\boldsymbol{E}_t$ respectively. This type of covariance model has been proposed and studied in the literature as the ``transposable" \citep{allen2010transposable}, ``array normal" \citep{hoff2011separable}, ``separable" \citep{tsiligkaridis2013covariance,zhou2014gemini}, and ``Kronecker product" \citep{hafner2020estimation,linton2019estimation} covariance structure. We refer the readers to \cite{hoff2011separable} and \cite{linton2019estimation} for a more detailed account on the history of the separable covariance matrix. \cite{chen:2019a} also considered the MAR model under this covariance structure.
\begin{remark}
{To introduce the Gaussian MLE, we have assumed that ${\boldsymbol{Z}}_t$ in the expression $\boldsymbol{E}_t=\Sigma_1^{1/2}{\boldsymbol{Z}}_t\Sigma_2^{1/2}$ has iid standard normal entries. Note that it is also possible to impose additional structures on ${\boldsymbol{Z}}_t$. For example, ${\boldsymbol{Z}}_t$ can have a rank one structure ${\boldsymbol{Z}}_t=\boldsymbol{z}_{1t}\boldsymbol{z}_{2t}'$ where $\boldsymbol{z}_{1t}$ and $\boldsymbol{z}_{1t}$ are random vectors of dimensions $d_1$ and $d_2$, respectively, with all elements independent and of unit-variance. In this case, $\boldsymbol{E}_t$ is a rank-one error matrix, with separated row noises and column noises. One difficulty of using this structure is that it implies that
$\boldsymbol{X}_t-\boldsymbol{A}_1\boldsymbol{X}_{t-1}\boldsymbol{A}_2'$ is of rank one for all $t$. Additional error terms may be needed. We leave such an extension to the future research.}
\end{remark}
Under the normality {and error structure \eqref{eq:ocov}}, the log likelihood of the RRMAR model is (up to some additive constants)
\begin{equation}
\label{eq:loglik}
-(T-1)(d_2\log|\Sigma_1|+d_1\log|\Sigma_2|)
-\sum_{t=2}^T\operatorname{tr}\left[ \Sigma_1^{-1}(\boldsymbol{X}_t-\boldsymbol{A}_1\boldsymbol{X}_{t-1}\boldsymbol{A}_2')\Sigma_2^{-1}(\boldsymbol{X}_t-\boldsymbol{A}_1\boldsymbol{X}_{t-1}\boldsymbol{A}_2')'\right].
\end{equation}
We will introduce an alternating algorithm to find the MLE, maximizing \eqref{eq:loglik} alternatively over one pair $(\boldsymbol{A}_i,\Sigma_i)$ while holding the other fixed. As will be seen, each iteration can be viewed as a reduced-rank regression, which is equivalent to the canonical correlation analysis \citep{reinsel:1998}.
We therefore denote the minimizer of \eqref{eq:loglik} by
$\hat\boldsymbol{A}_i^{\hbox{\tiny cc}}$ and $\hat\Sigma_i$, and refer to them as the RR.CC estimators.
We now describe how to estimate $\boldsymbol{A}_1$ and $\Sigma_1$ when $\boldsymbol{A}_2$ and $\Sigma_2$ are known. Under assumption \eqref{eq:ocov}, we can rewrite the model as
\begin{equation*}
\left(\boldsymbol{X}_t\Sigma_2^{-1/2}\right)[,j] = \boldsymbol{A}_1 \left(\boldsymbol{X}_{t-1}\boldsymbol{A}_2'\Sigma_2^{-1/2}\right)[,j] + \left(\boldsymbol{E}_t\Sigma_2^{-1/2}\right)[,j].
\end{equation*}
Note that the columns of the transformed error matrix $\boldsymbol{E}_t\Sigma_2^{-1/2}$ are iid $N({\boldsymbol{0}},\Sigma_1)$. If we let $\boldsymbol{y}_{tj}=\left(\boldsymbol{X}_t\Sigma_2^{-1/2}\right)[,j]$, \: $\boldsymbol{x}_{tj}=\left(\boldsymbol{X}_{t-1}\boldsymbol{A}_2'\Sigma_2^{-1/2}\right)[,j]$ and $\boldsymbol{\epsilon}_{tj}=\left(\boldsymbol{E}_t\Sigma_2^{-1/2}\right)[,j]$, then the preceding equation can be viewed as a reduced-rank regression with i.i.d. errors:
\begin{equation}
\label{eq:rrr2}
\boldsymbol{y}_{tj} = \boldsymbol{A}_1\boldsymbol{x}_{tj} + \boldsymbol{\epsilon}_{tj}, \quad 2\leq t\leq T,\;1\leq j\leq d_2.
\end{equation}
The MLE of $\boldsymbol{A}_1$ based on \eqref{eq:rrr2} with i.i.d. normal errors has been well studied in the classical reduced-rank regression. Here we only define necessary notations to introduce the final expression of the MLE. We refer the readers to the classical texts \cite{anderson:2013} and \cite{reinsel:1998} for more details.
Let
\begin{align*}
\tilde\boldsymbol{S}_{xx} & =\sum_t\sum_j \boldsymbol{x}_{tj}\boldsymbol{x}_{tj}'=\sum_t \boldsymbol{X}_{t-1}\boldsymbol{A}_2'\Sigma_2^{-1}\boldsymbol{A}_2\boldsymbol{X}_{t-1}',\\
\tilde\boldsymbol{S}_{yx} & =\sum_t\sum_j \boldsymbol{y}_{tj}\boldsymbol{x}_{tj}'=\sum_t \boldsymbol{X}_t\Sigma_2^{-1}\boldsymbol{A}_2\boldsymbol{X}_{t-1}'.
\end{align*}
The least squares estimator (with no rank constraint) of $\boldsymbol{A}_1$ based on \eqref{eq:rrr2} is then given by $\tilde\boldsymbol{A}_1 = \tilde\boldsymbol{S}_{yx}\tilde\boldsymbol{S}_{xx}^{-1}$. To get the MLE of $\boldsymbol{A}_1$ under the constraint $\operatorname{rank}(\boldsymbol{A}_1)=k_1$, let
\begin{equation*}
\tilde\Sigma_{\epsilon\epsilon} = \sum_t\sum_j\left(\boldsymbol{y}_{tj}-\tilde\boldsymbol{A}_1\boldsymbol{x}_{tj}\right)\left(\boldsymbol{y}_{tj}-\tilde\boldsymbol{A}_1\boldsymbol{x}_{tj}\right)' = \sum_t\left(\boldsymbol{X}_t-\tilde\boldsymbol{A}_1\boldsymbol{X}_{t-1}\boldsymbol{A}_2'\right)\Sigma_2^{-1}\left(\boldsymbol{X}_t-\tilde\boldsymbol{A}_1\boldsymbol{X}_{t-1}\boldsymbol{A}_2'\right)'.
\end{equation*}
Take $\tilde\boldsymbol{U}:=[\tilde U_1, \tilde U_2,\ldots, \tilde U_{k_1}]$, where $\tilde U_j$ is the $j$-th leading unit eigenvector of
$\tilde\Sigma_{\epsilon\epsilon}^{-1/2}\tilde\boldsymbol{S}_{yx}\tilde\boldsymbol{S}_{xx}^{-1}\tilde\boldsymbol{S}_{xy}\tilde\Sigma_{\epsilon\epsilon}^{-1/2}$. Then $\boldsymbol{A}_1$ is
updated as
\begin{equation}
\label{eq:mle-A1}
\check \boldsymbol{A}_1^{\hbox{\tiny cc}}=\tilde\Sigma_{\epsilon\epsilon}^{1/2} \tilde \boldsymbol{U} \tilde \boldsymbol{U}'\tilde\Sigma_{\epsilon\epsilon}^{-1/2}\tilde\boldsymbol{S}_{yx}\tilde\boldsymbol{S}_{xx}^{-1}.
\end{equation}
For the derivation of this update, see for example Equation (2.15) of \cite{reinsel:1998}.
Subsequently, the covariance matrix $\Sigma_1$ is updated as
\begin{equation*}
\check\Sigma_1 = \frac{1}{T-1} \sum_t\left(\boldsymbol{X}_t-\check\boldsymbol{A}_1^{\hbox{\tiny cc}}\boldsymbol{X}_{t-1}\boldsymbol{A}_2'\right)\Sigma_2^{-1}\left(\boldsymbol{X}_t-\check\boldsymbol{A}_1^{\hbox{\tiny cc}}\boldsymbol{X}_{t-1}\boldsymbol{A}_2'\right)'.
\end{equation*}
Given $\boldsymbol{A}_1$ and $\Sigma_1$, an update of $\boldsymbol{A}_2$ and $\Sigma_2$ can be similarly obtained. Therefore, we use the alternating algorithm to find the minimizer of \eqref{eq:loglik}.
\begin{remark}
Note that each step of the algorithms reduces the corresponding objective functions, hence the algorithms are likely to converge to a local/global minimum.
To ensure the estimators do not oscillating among equivalent solutions, proper constraints need to be enforced. We have assumed that $\boldsymbol{A}_1$ is always normalized so that $\|\boldsymbol{A}_1\|_F=1$, and will require the same for its estimates in the algorithm. Similarly, for the MLE we require $\|\Sigma_1\|_F=1$ and $\|\hat\Sigma_1\|_F=1$. On the other hand, since the noise covariance matrix is full rank and the ranks $k_1$ and $k_2$ are the correct ranks, there is no identifiability issue related to rank deficiency.
\end{remark}
\begin{remark}
Since the objective functions are not convex, care
needs to be taken to reach the global minimum. Using multiple random initial values is one possible approach.
Another approach is to use
the projection estimator of $\boldsymbol{A}_1$ and $\boldsymbol{A}_2$ in \cite{chen2021autoregressive} (by ignoring the rank constraints), as the initial values of both alternating algorithms. The projection estimator of $\boldsymbol{A}_1$ and $\boldsymbol{A}_2$ are derived from a global (one-step) estimation of $\boldsymbol{A}_2\otimes\boldsymbol{A}_1$, though they are not efficient.
\end{remark}
\section{Asymptotics}
\label{sec:asymp}
The asymptotic analysis is substantially different from the classical
reduced-rank regression, due to the alternating nature of the
estimation. For example, the gradient condition for the LSE $\hat\boldsymbol{A}_1^{\hbox{\tiny ls}}$ is
$\hat\boldsymbol{A}_1^{\hbox{\tiny
ls}}=\hat\boldsymbol{U}\hat\boldsymbol{U}'\hat\boldsymbol{S}_{yx}\hat\boldsymbol{S}_{xx}^{-1}$, where
$\hat \boldsymbol{U}$, $\hat\boldsymbol{S}_{xx}$ and $\hat\boldsymbol{S}_{yx}$ are defined as the
$\boldsymbol{U}$, $\boldsymbol{S}_{xx}$ and $\boldsymbol{S}_{yx}$ in Section~\ref{sec:est} with the
modification that all $\boldsymbol{A}_2$ therein need to be replaced by
$\hat\boldsymbol{A}_2^{\hbox{\tiny ls}}$. In other words, the asymptotic
behaviors of $\hat\boldsymbol{A}_1^{\hbox{\tiny ls}}$ and
$\hat\boldsymbol{A}_2^{\hbox{\tiny ls}}$ are intertwined.
The asymptotics for $\boldsymbol{A}_i^{\hbox{\tiny ls}}$ and $\boldsymbol{A}_i^{\hbox{\tiny cc}}$ involve heavy notations. First of all,
recall that we assume $\|\boldsymbol{A}_1\|_F=1$ for the parameter identifiability, so we rescale
$\hat\boldsymbol{A}_i^{\hbox{\tiny ls}}$ and $\hat\boldsymbol{A}_i^{\hbox{\tiny cc}}$ so that
$\|\hat\boldsymbol{A}_1^{\hbox{\tiny ls}}\|_F=1$ and $\|\hat\boldsymbol{A}_1^{\hbox{\tiny cc}}\|_F=1$.
In the table below, we list notations that appear in both Theorem~\ref{thm:ls} and Theorem~\ref{thm:cc}, but with different definitions in these theorems.
\begin{center}
\begin{tabular}{l|l|l}
\hline
Notations & Theorem~\ref{thm:ls} & Theorem~\ref{thm:cc} \\\hline
$\Gamma_1$ & $\mathbb{E}(\boldsymbol{X}_t\boldsymbol{A}_2'\boldsymbol{A}_2\boldsymbol{X}_t')$ & $\mathbb{E}(\boldsymbol{X}_t\boldsymbol{A}_2'\Sigma_2^{-1}\boldsymbol{A}_2\boldsymbol{X}_t')$ \\
$\Gamma_2$ & $\mathbb{E}(\boldsymbol{X}_t'\boldsymbol{A}_1'\boldsymbol{A}_1\boldsymbol{X}_t)$ & $\mathbb{E}(\boldsymbol{X}_t'\boldsymbol{A}_1'\Sigma_1^{-1}\boldsymbol{A}_1\boldsymbol{X}_t)$\\
$\mathbb{P}_i$ & \hbox{orthogonal projection to} $\mathrm{col}(\boldsymbol{A}_i)$ & \hbox{orthogonal projection to} $\mathrm{col}(\Sigma_i^{-1/2}\boldsymbol{A}_i)$\\
$\mathcal{P}_i$ & $\mathbb{P}_i$ & $\Sigma_i^{-1/2}\mathbb{P}_i\Sigma_i^{1/2}$ \\
$\boldsymbol{\alpha}_1$ & $\mathrm{vec}(\boldsymbol{A}_1)$ & Same \\
$\boldsymbol{\gamma}_1$ & $(\boldsymbol{\alpha}_1',{\boldsymbol{0}}')'$ & Same \\
$\boldsymbol{W}_t$ & $[(\boldsymbol{A}_2\boldsymbol{X}_t')\otimes\boldsymbol{I}_{d_1},\boldsymbol{I}_{d_2}\otimes(\boldsymbol{A}_1\boldsymbol{X}_t)]'$ & {Same} \\
$\boldsymbol{H}$ & $\mathbb{E}(\boldsymbol{W}_t\boldsymbol{W}_t')+\boldsymbol{\gamma}_1\boldsymbol{\gamma}_1'$ & $\mathbb{E}(\boldsymbol{W}_t\Sigma_e^{-1}\boldsymbol{W}_t')+\boldsymbol{\gamma}_1\boldsymbol{\gamma}_1'$ \\\hline
\end{tabular}
\end{center}
Define
\begin{equation*}
\boldsymbol{Q}_t:=\begin{pmatrix}
\boldsymbol{X}_t\boldsymbol{A}_2'\otimes \mathcal{P}_1+[\Gamma_1\boldsymbol{A}_1'(\boldsymbol{A}_1\Gamma_1\boldsymbol{A}_1')^+\boldsymbol{A}_1\boldsymbol{X}_t\boldsymbol{A}_2']\otimes(\boldsymbol{I}-\mathcal{P}_1) \\
\mathcal{P}_2\otimes \boldsymbol{X}_t'\boldsymbol{A}_1'+(\boldsymbol{I}-\mathcal{P}_2)\otimes [\Gamma_2\boldsymbol{A}_2'(\boldsymbol{A}_2\Gamma_2\boldsymbol{A}_2')^+\boldsymbol{A}_2\boldsymbol{X}_t'\boldsymbol{A}_1']
\end{pmatrix},
\end{equation*}
where $\boldsymbol{M}^+$ denotes the Moore-Penrose inverse of $\boldsymbol{M}$. Note that $\boldsymbol{Q}_t$ will appear in both Theorem~\ref{thm:ls} and Theorem~\ref{thm:cc}. Although it seems to have the same definition in both theorems, the two versions actually differ because the $\Gamma_i$ and $\mathcal{P}_i$ involved have different definitions.
\begin{theorem}
\label{thm:ls}
Assume that $\{\boldsymbol{E}_t\}$ are i.i.d. with mean zero and finite second
moments. Assume that
$0<\operatorname{rank}{\boldsymbol{A}_i}=k_i\leq d_i$, $\rho(\boldsymbol{A}_1)\rho(\boldsymbol{A}_2)<1$, where $\rho(\cdot)$ denotes the spectral radius of a matrix, and
$\Sigma_e$ is non-singular.
Then
\begin{equation*}
\sqrt{T}\begin{pmatrix}
\mathrm{vec}\left[\hat\boldsymbol{A}_1^{\hbox{\tiny ls}}-\boldsymbol{A}_1\right] \\
\mathrm{vec}\left[(\hat\boldsymbol{A}_2^{\hbox{\tiny ls}})'-\boldsymbol{A}_2'\right]
\end{pmatrix}
\Rightarrow N({\boldsymbol{0}}, \Xi^{\hbox{\tiny ls}}),
\end{equation*}
where
\begin{equation}
\label{eq:cov_ls}
\Xi^{\hbox{\tiny ls}} :=\boldsymbol{H}^{-1}\mathbb{E}(\boldsymbol{Q}_t\Sigma_e\boldsymbol{Q}_t')\boldsymbol{H}^{-1}.
\end{equation}
\end{theorem}
\begin{theorem}
\label{thm:cc}
Assume that $\{\boldsymbol{E}_t\}$ are i.i.d. with mean zero and finite second
moments. Assume that
$0<\operatorname{rank}{\boldsymbol{A}_i}=k_i\leq d_i$, $\rho(\boldsymbol{A}_1)\rho(\boldsymbol{A}_2)<1$, and
$\Sigma_e$ is of the form \eqref{eq:ocov}, and is non-singular.
Then
\begin{equation*}
\sqrt{T}\begin{pmatrix}
\mathrm{vec}\left[\hat\boldsymbol{A}_1^{\hbox{\tiny cc}}-\boldsymbol{A}_1\right] \\
\mathrm{vec}\left[(\hat\boldsymbol{A}_2^{\hbox{\tiny cc}})'-\boldsymbol{A}_2'\right]
\end{pmatrix}
\Rightarrow N({\boldsymbol{0}}, \Xi^{\hbox{\tiny cc}}),
\end{equation*}
where
\begin{equation}
\label{eq:cov_cc}
\Xi^{\hbox{\tiny cc}} :=\boldsymbol{H}^{-1}\mathbb{E}(\boldsymbol{Q}_t\Sigma_e^{-1}\boldsymbol{Q}_t')\boldsymbol{H}^{-1}.
\end{equation}
\end{theorem}
Besides the usual regularity conditions, we also assume in Theorem~\ref{thm:ls} and Theorem~\ref{thm:cc} that certain matrices have distinct eigenvalues. While this assumption holds for generic positive definite matrices, it is even easier to be fulfilled by the aforementioned matrices since they involve $\Gamma_i$.
As mentioned before, the RRMAR model is a MAR model with additional low rank constraints on the coefficient matrices $\boldsymbol{A}_i$. If the estimation is carried out without the rank constraints, then the procedure of \cite{chen2021autoregressive} applies and so does its asymptotic result ({e.g.} Theorem~4 {in \cite{chen2021autoregressive}}). In fact, if the estimation of $\boldsymbol{A}_i$ is given by the MLE without imposing the low rank constraints, the asymptotic covariance matrix would take the same form as $\Xi^{\hbox{\tiny cc}}$, by setting $\mathcal P_i=\boldsymbol{I}$ in the definition of $\boldsymbol{Q}_t$. Denote this covariance matrix by $\tilde \Xi$. The following theorem asserts that the MLE $\hat\boldsymbol{A}_i^{\hbox{\tiny cc}}$ under the RRMAR model are asymptotically more efficient.
\begin{theorem}
\label{thm:efficiency}
Under the assumptions of Theorem~\ref{thm:cc}, it holds that $\tilde\Xi\succeq \Xi^{\hbox{\tiny cc}}$, i.e. the difference $\tilde\Xi - \Xi^{\hbox{\tiny cc}}$ is positive semi-definite.
\end{theorem}
\begin{remark}
Since $\Sigma_e$ in Theorem~1 can be arbitrary, as long as it is non-singular, the LSE does not correspond to the MLE, and thus there is no similar result to Theorem~\ref{thm:efficiency} regarding the comparison of the LSE under the RRMAR model and the MAR model without rank constraints. On the other hand, if the entries of $\boldsymbol{E}_t$ are IID, we can show a similar result to Theorem~\ref{thm:efficiency} for the LSE, which is not presented here since it is too special.
\end{remark}
\medskip
We now consider the asymptotics
of the composition and loading matrices $\boldsymbol{A}_{il}$ and $\boldsymbol{A}_{ic}$ in \eqref{eq:marrr_fac}. The following discussion works the same for either $\hat\boldsymbol{A}_i^{\hbox{\tiny ls}}$ or $\hat\boldsymbol{A}_i^{\hbox{\tiny cc}}$. Therefore, we will use the unified notations $\hat\boldsymbol{A}_i$ and $\Xi$, dropping the superscripts $^{\hbox{\tiny ls}}$ and $^{\hbox{\tiny cc}}$. Since $\boldsymbol{A}_{ic}$ and $\boldsymbol{A}_{il}$ cannot be identified as seen from $\boldsymbol{A}_{il}\boldsymbol{A}_{ic}'=\boldsymbol{A}_{il}\boldsymbol{M}\boldsymbol{M}^{-1}\boldsymbol{A}_{ic}'$ for any invertible $k_i\times k_i$ matrix $\boldsymbol{M}$, we consider instead the singular value decomposition (SVD) of $\boldsymbol{A}_i$. Write $\boldsymbol{A}_{i}={\boldsymbol{U}}_i\boldsymbol{D}_i\boldsymbol{V}_i'$, where both $\boldsymbol{U}_i$ and $\boldsymbol{V}_i$ are $d_i\times k_i$ ortho-normal matrices. Denote the $j$-th diagonal element of $\boldsymbol{D}_i$ by $d_{ij}$, and define $\boldsymbol{d}_i=(d_{i1},\ldots,d_{i,k_i})'$. Comparing \eqref{eq:marrr_fac}, we see that $\boldsymbol{V}_i$ corresponds to the composition matrix $\boldsymbol{A}_{ic}$, $\boldsymbol{U}_i$ corresponds to the loading matrix $\boldsymbol{A}_{il}$, and $\boldsymbol{D}_i$ can be absorbed into either $\boldsymbol{A}_{il}$ or $\boldsymbol{A}_{ic}$.
Let $\hat\boldsymbol{A}_i =\hat\boldsymbol{U}_i\hat\boldsymbol{D}_i(\hat\boldsymbol{V}_i)'$ be the SVD of $\hat\boldsymbol{A}_i$. Since $\hat\boldsymbol{A}_i\hat\boldsymbol{A}_i' = \hat\boldsymbol{U}_i\hat\boldsymbol{D}_i^2(\hat\boldsymbol{U}_i)'$, the asymptotic distribution of $\hat\boldsymbol{U}_i$ can be obtained based on that of $\hat\boldsymbol{A}_i\hat\boldsymbol{A}_i'$. Similarly, the asymptotic distribution of $\hat\boldsymbol{V}_i$ can be derived from that of $\hat\boldsymbol{A}_i'\hat\boldsymbol{A}_i$.
Note that the asymptotic covariance matrix of $\mathrm{vec}(\hat\boldsymbol{A}_i)$ (for $i=1,2$) is a submatrix of $\Xi$ and can be extracted from \eqref{eq:cov_ls} or \eqref{eq:cov_cc}. Following that, we let $\Xi_{i1}$ be the asymptotic covariance matrix of $\mathrm{vec}(\hat\boldsymbol{A}_i\hat\boldsymbol{A}_i')$, which can be obtained through the expansion $$\hat\boldsymbol{A}_i\hat\boldsymbol{A}_i' = \boldsymbol{A}_i\boldsymbol{A}_i'+(\hat\boldsymbol{A}_i-\boldsymbol{A}_i)\boldsymbol{A}_i'+\boldsymbol{A}_i(\hat\boldsymbol{A}_i-\boldsymbol{A}_i)'+o_P(T^{-1/2}).$$ More specifically, when $i=1$,
\begin{equation*}
\Xi_{11}=\left[\boldsymbol{A}_1\otimes\boldsymbol{I}_{d_1}+(\boldsymbol{I}_{d_1}\otimes\boldsymbol{A}_1)\boldsymbol{J}_{d_1,d_1}\right]\left\{ \Xi[1:d_1^2,1:d_1^2] \right\}\left[\boldsymbol{A}_1\otimes\boldsymbol{I}_{d_1}+(\boldsymbol{I}_{d_1}\otimes\boldsymbol{A}_1)\boldsymbol{J}_{d_1,d_1}\right]',
\end{equation*}
where $\Xi[1:d_1^2,1:d_1^2]$ is the upper left $d_1^2\times d_1^2$ block of $\Xi^{\hbox{\tiny ls}}$ or $\Xi^{\hbox{\tiny cc}}$, and the matrix $\boldsymbol{J}_{d_1,d_1}$ is defined in \eqref{eq:J}.
The asymptotic covariance matrix of $\hat\boldsymbol{A}_1'\hat\boldsymbol{A}_1$, denoted by $\Xi_{12}$, has a similar expression. When $i=2$, the matrices $\Xi_{21}$ and $\Xi_{22}$, related to $\hat\boldsymbol{A}_2$, are also defined similarly.
Define the matrix ${\boldsymbol{R}}_{i1}$ as
\begin{equation*}
{\boldsymbol{R}}_{i1}=(\boldsymbol{I}_{k_i}\otimes\boldsymbol{U}_i,\boldsymbol{I}_{k_i}\otimes \boldsymbol{U}_i^{\perp})
\begin{pmatrix}
(\boldsymbol{D}_i^2\otimes\boldsymbol{I}_{k_i}-\boldsymbol{I}_{k_i}\otimes\boldsymbol{D}_i^2+\boldsymbol{L}_{k_i}\boldsymbol{L}_{k_i}')^{-1}(\boldsymbol{I}_{k_i^2}-\boldsymbol{L}_{k_i}\boldsymbol{L}_{k_i}')(\boldsymbol{U}_i'\otimes\boldsymbol{U}_i') \\
(\boldsymbol{D}_i^{-2}\boldsymbol{U}_i')\otimes(\boldsymbol{U}_i^\perp)'
\end{pmatrix},
\end{equation*}
and define ${\boldsymbol{R}}_{i2}$ similarly, but replacing $\boldsymbol{U}_i$ with $\boldsymbol{V}_i$.
Note that even when $\boldsymbol{A}_i$ has distinct singular values, the columns of $\boldsymbol{U}_i$ and $\boldsymbol{V}_i$ are only identified up to sign changes. We shall adopt the following convection to identify $\boldsymbol{U}_i$: the first nonzero element of each column of $\boldsymbol{U}_i$ is positive. Since $\boldsymbol{A}_1$ and $\boldsymbol{A}_2$ are also only identified up to sign changes, we make one more requirement to identify $\boldsymbol{V}_i$: the first nonzero element of the first column of $\boldsymbol{V}_2$ is positive. Subsequently, we also require the estimators $\hat\boldsymbol{U}_i$ and $\hat\boldsymbol{V}_i$ to satisfy these identifiability conditions.
Now, as a consequence of Theorem~\ref{thm:ls} and Theorem~\ref{thm:cc}, we have the following result regarding $\hat\boldsymbol{U}_i$ and $\hat\boldsymbol{V}_i$.
\begin{corollary}
\label{cor:svd}
Assume the conditions of Theorem~\ref{thm:ls} or Theorem~\ref{thm:cc} hold, and that the singular values of $\boldsymbol{A}_1$ are distinct, and so are those of $\boldsymbol{A}_2$. For each of $i=1,2$, it holds that
\begin{equation*}
\sqrt{T}\mathrm{vec}(\hat\boldsymbol{U}_i-\boldsymbol{U}_i)\Rightarrow
N\left({\boldsymbol{0}}, {\boldsymbol{R}}_{i1} \Xi_{i1} {\boldsymbol{R}}_{i1}'\right),
\end{equation*}
and
\begin{equation*}
\sqrt{T}\mathrm{vec}(\hat\boldsymbol{V}_i-\boldsymbol{V}_i)\Rightarrow
N\left({\boldsymbol{0}}, {\boldsymbol{R}}_{i2} \Xi_{i2} {\boldsymbol{R}}_{i2}'\right),
\end{equation*}
and
\begin{equation*}
\sqrt{T}(\hat\boldsymbol{d}_{i} -\boldsymbol{d}_i) \Rightarrow N\left({\boldsymbol{0}},\frac{1}{4}\boldsymbol{D}_i^{-1}\boldsymbol{L}_{k_i}'(\boldsymbol{U}_i'\otimes\boldsymbol{U}_i')\Xi_{i1}(\boldsymbol{U}_i\otimes\boldsymbol{U}_i)\boldsymbol{L}_{k_i}\boldsymbol{D}_i^{-1}\right).
\end{equation*}
\end{corollary}
\begin{remark}
Let $\boldsymbol{G}_i$ be the $k_i\times k_i$ matrix defined as
\begin{equation*}
\boldsymbol{G}_i[j,k] = \left\{\begin{array}{ll}
(d^2_{ik}-d^2_{ij})^{-1} & \hbox{when } j\neq k, \\
0 & \hbox{when } j=k.
\end{array}\right.
\end{equation*}
Furthermore, let $\boldsymbol{U}_i^\perp$ be the $d_i\times (d_i-k_i)$ orthonormal matrix such that $(\boldsymbol{U}_i,\boldsymbol{U}_i^\perp)$ is an orthogonal matrix. The asymptotic distribution of $\hat\boldsymbol{U}_i$ can also be obtained from the following equation
\begin{equation*}
\hat\boldsymbol{U}_i-\boldsymbol{U}_i = (\boldsymbol{U}_i,\boldsymbol{U}_i^{\perp})
\begin{pmatrix}
\boldsymbol{G}_i\circ \left[\boldsymbol{U}_i'(\hat\boldsymbol{A}_i\hat\boldsymbol{A}_i' - \boldsymbol{A}_i\boldsymbol{A}_i')\boldsymbol{U}_i\right] \\
(\boldsymbol{U}_i^{\perp})'(\hat\boldsymbol{A}_i\hat\boldsymbol{A}_i' - \boldsymbol{A}_i\boldsymbol{A}_i')\boldsymbol{U}_i\boldsymbol{D}_i^{-2}
\end{pmatrix}
+ o_P\left(\frac{1}{\sqrt{T}}\right),
\end{equation*}
where $\circ$ denotes the entry-wise Hadamard or Schur product of two matrices.
\end{remark}
\begin{remark}
The joint distribution of $\hat\boldsymbol{U}_i$ and $\hat\boldsymbol{V}_i$ can also be derived. However, we choose not to spell the details out here for two reasons: the notations are already very complicated, and more importantly, there does not seem to be any direct applications of such a joint distribution.
\end{remark}
\section{Identification of the Rank}
\label{sec:rank}
In most applications the ranks of $\boldsymbol{A}_i$ are unknown, and it is
important to determine them
from the data. This problem has been considered for multivariate reduced-rank regression by \cite{anderson:1951} and \cite{anderson:2013}, and for reduced-rank autoregressive model by \cite{kohn1979asymptotic}, \cite{reinsel:1998}, \cite{tiao1989model} and \cite{tsay1985use}, among others. For high dimensional reduced-rank regression based on independent samples, penalized least squares can select the ranks along with the estimation, where the penalty is based on nuclear norm \citep{yuan:2007,negahban:2011}, $\ell_0$ norm
\citep{bunea:2011}, or Schatten-$q$ quasi-norm
\citep{rohde:2011}. \cite{basu2019low} and \cite{lin2020regularized} considered low rank VAR and tensor models,
combining least squares estimation with the nuclear norm
penalty.
We propose to use an information criterion to select the ranks. For a given pair of ranks $(r_1,r_2)$, it is defined as
\begin{equation}
\label{eq:bic}
\begin{aligned}
\mathrm{EBIC}(r_1,r_2) = & \log\left[\frac{1}{Td_1d_2}\sum_{t=2}^T\|\boldsymbol{X}_t-\hat\boldsymbol{A}_1^{\hbox{\tiny ls}}\boldsymbol{X}_{t-1}(\hat\boldsymbol{A}_2^{\hbox{\tiny ls}})'\|_F^2\right] \\
& + \frac{1}{Td_1d_2}\cdot[\log(Td_2)\cdot r_1(2d_1-r_1)+\log(Td_1)\cdot r_2(2d_2-r_2)].
\end{aligned}
\end{equation}
This can be viewed as an extended version of the Bayesian Information Criterion \citep{schwarz1978estimating}, so we use the acronym EBIC. Here the likelihood is calculated for the model where the entries of $\boldsymbol{E}_t$ are iid $N(0,\sigma^2)$, so it is best viewed as a “quasi”-likelihood. It is not precisely derived according to the posterior probability under the Bayesian framework \citep{haughton1988choice}. Instead, we combine the quasi-log-likelihood and a penalty term where the numbers of parameters are multiplied by the logarithm of the sample sizes. Instead of simply counting the number of parameters, the effective number of parameters \citep{mukherjee2015degrees,yuan2016degrees} can also be used in the EBIC. However, we choose the current version for simplicity. The selected pair of ranks $(\hat k_1,\hat k_2)$ minimizes the EBIC over all the pairs $(r_1,r_2)$ such that $1\leq r_1\leq r_{1\!\max}$ and $1\leq r_2\leq r_{2\!\max}$, where $r_{1\!\max}$ and $r_{2\!\max}$ are pre-determined maximum ranks of $\boldsymbol{A}_1$ and $\boldsymbol{A}_2$, respectively. If no information is available and the dimensions $d_1$ and $d_2$ are not too large, we can simply use $r_{1\!\max}=d_1$ and $r_{2\!\max}=d_2$. The following Theorem \ref{thm:bic} confirms that the EBIC in \eqref{eq:bic} does achieve the consistency. Its empirical performances are also outstanding, as will be shown in Section~\ref{sec:simulation}.
Since there are two ranks to be determined, a direct search via \eqref{eq:bic} over all possible pairs of ranks can be very costly when both $r_{1\!\max}$ and $r_{2\!\max}$ are large. We also consider selecting these two ranks separately. Specifically, the selected ranks are
\[
\hat{r}_1=\arg\min_{r_1}\mathrm{EBIC}(r_1,r_{2\!\max})\quad\quad \hat{r}_2=\arg\min_{r_2}\mathrm{EBIC}(r_{1\!\max},r_{2}).
\]
\begin{theorem}\label{thm:bic}
Assume that $\{\boldsymbol{E}_t\}$ are i.i.d. with mean zero and finite second moments.
Also assume that $0<\operatorname{rank}{\boldsymbol{A}_i}=k_i\leq d_i$, $\rho(\boldsymbol{A}_1)\rho(\boldsymbol{A}_2)<1$, and
$\Sigma_e$ is non-singular. Then both the joint $\mathrm{EBIC}$ and the separate $\mathrm{EBIC}_i$ select the true ranks consistently, given that $k_i\leq r_{i\!\max}$.
\end{theorem}
\begin{remark}
{Theorem~\ref{thm:bic} continues to hold if the penalty {term} in EBIC \eqref{eq:bic} {(the second term)} is scaled by any positive constant $c$. {The finite sample performance of EBIC may depend on such a constant.} The resampling based tuning method introduced by \cite{hallin2007} can be {adopted}
here to determine $c$. The benefit of using a data-driven penalty can be significant when the dimension is high.}
\end{remark}
\begin{remark}
{If the second term in \eqref{eq:bic} is replaced by {$2(Td_1d_2)^{-1}[r_1(2d_1-r_1)+r_2(2d_2-r_2)]$}, the criterion is similar to AIC. Although it may not be consistent, it often works well in small sample \citep{shao1997asymptotic,brockwell:1991}. }
\end{remark}
\section{Numerical Studies}
\label{sec:num}
\subsection{Simulations}
\label{sec:simulation}
In this section, we investigate the finite sample performance of the proposed estimation and rank determination procedures for the RRMAR models under various simulation setups. The simulation study consists of three parts. The first part is designed to compare the empirical behavior of the proposed alternating least square estimator $\hat\boldsymbol{A}_i^{\hbox{\tiny ls}}$, labelled as RR.LS in the figures, and the alternating MLE $\hat\boldsymbol{A}_i^{\hbox{\tiny cc}}$, labelled as RR.CC. The true ranks are taken as known. The least squares estimator without rank constraints (labelled as LSE) in \cite{chen2021autoregressive} is also included as a benchmark for comparison. In the second part, we report the coverage probabilities of the confidence intervals constructed based on Theorem~\ref{thm:ls}, Theorem~\ref{thm:cc} and Corollary~\ref{cor:svd}. The third part examines the rank determination based on the EBIC proposed in Section~\ref{sec:rank}. We also experiment with rank selection by rolling forecasting.
For given dimensions $d_i$ and ranks $k_i$, the observed
data $\boldsymbol{X}_t$ are simulated according to model \eqref{eq:marrr}. The matrix $\boldsymbol{A}_1$ is generated according to $\boldsymbol{A}_1=\boldsymbol{Q}_1\Lambda\boldsymbol{Q}_2'$, where the entries of the $k_1\times k_2$ diagonal matrix $\Lambda$ are sampled from the uniform distribution over the interval $[0.5,1.5]$, and the $d_1\times k_1$ orthonormal matrices $\boldsymbol{Q}_1$ and $\boldsymbol{Q}_2$ are generated randomly from the Haar distribution.
The matrix $\boldsymbol{A}_2$ is generated in the same way. The two matrices $\boldsymbol{A}_1$ and $\boldsymbol{A}_2$ are then rescaled so that $\rho:=\rho(\boldsymbol{A}_1)\rho(\boldsymbol{A}_2)<1$ and $||\boldsymbol{A}_1||_F=1$.
Throughout all the simulation studies, two different settings of the covariance structure of the innovation matrix $\boldsymbol{E}_t$ are considered:
\begin{enumerate}
\item[(I)] The covariance matrix $\Sigma_e=\operatorname{Cov}(\mathrm{vec}(\boldsymbol{E}_t))$ is randomly generated according to $\Sigma_e= \boldsymbol{Q}\Lambda\boldsymbol{Q}'$, where the entries of the diagonal matrix $\Lambda$ are equally spaced over $[1,10]$, and $\boldsymbol{Q}$ is a random orthogonal matrix generated from the Haar distribution.
\item [(II)] The covariance matrix $\Sigma_e$ takes the form \eqref{eq:ocov}, where each of $\Sigma_1$ and $\Sigma_2$ is generated in the same way as the $\Sigma_e$ in Setting I, except that the diagonal entries of $\Lambda$ are equally spaced over $[1,5]$.
\end{enumerate}
For a particular simulation setting with multiple repetitions, the
matrices $\boldsymbol{A}_i$ and $\Sigma_e$ are fixed.
\begin{remark}
We randomize the population parameters in the simulation, while controlling the key parameters (e.g.
$\rho(\boldsymbol{A}_1)\rho(\boldsymbol{A}_2)$ and the entries of the diagonal matrix $\Lambda$ in the noise covariance structure). The main reason is that the theoretical results show that the estimation performance mainly depend on these controlled parameters. Due to the large number of ``free parameters", it is difficult to construct a specific
design of the parameter set, and the individual parameter is of less importance and interests. Also one won't be able to cover all possible combinations.
In a way, the randomly generated population parameter set can be viewed as a ``representative" set.
\end{remark}
In the first experiment, for each configuration of sample size $T$, dimensions $d_i$ and ranks $k_i$, we repeat the simulation 100 times, and show a box plot of the estimation error
\begin{equation*}
\log(\|\hat\boldsymbol{A}_2\otimes\hat\boldsymbol{A}_1-\boldsymbol{A}_2\otimes\boldsymbol{A}_1\|_F^2).
\end{equation*}
The spectral radius is fixed at $\rho=.75$. Figure~\ref{fig:simu1} and Figure~\ref{fig:simu2} use the {boxplots\footnote{{Each boxplot shows the distribution of the log errors from 100 repeated experiments in each setting. The box in the middle shows the range of $1$-st quartile to $3$-rd quantile, and the stems extends to the minimum of 1.5 interquartile range (IQR) and the maximum/minimum of the data. The (rare) dots show the observations outside the 1.5 IRQ range (often considered as outliers).}}} to compare LSE, RR.LS and RR.CC, under Settings (I) and (II) respectively. It is seen from both figures that the advantage of RR.LS and RR.CC over LSE gets bigger as the dimensions grow higher. On the other hand, for fixed dimensions, this advantage becomes smaller as the ranks increase.
From Figure~\ref{fig:simu1} we also find that, under Setting (I) of the error covariance matrix, RR.CC performs similarly as RR.LS does, even though the covariance matrix $\Sigma_e$
does not have the form \eqref{eq:ocov}, which is assumed for RR.CC. On the other hand, Figure~\ref{fig:simu2} clearly demonstrates the advantage of RR.CC over RR.LS under setting (II), when $\Sigma_e$ does bear the form \eqref{eq:ocov}. {In addition, it is seen from Figure 2 that the differences (in log error) between LSE and RR.LS/RR.CC remain roughly the same for difference sample sizes. This confirms the results in Theorem 3, which states that the gain in efficiency (in terms of reduction in variance of the estimators) is proportional to $T$, hence the absolute difference between the full model in Chen, Xiao and Yang (2019) and the reduced-rank model becomes smaller as $T$ increases, but their ratio or log difference remain roughly the same.}
\begin{figure}[h]
\centering
\includegraphics[width=2.1in]{figures/6_4_1_1_01}~~\includegraphics[width=2.1in]{figures/6_4_3_2_01}~~\includegraphics[width=2.1in]{figures/6_4_5_3_01}\\
\includegraphics[width=2.1in]{figures/9_6_1_1_01}~~\includegraphics[width=2.1in]{figures/9_6_3_2_01}~~\includegraphics[width=2.1in]{figures/9_6_5_3_01}\\
\includegraphics[width=2.1in]{figures/15_10_1_1_01}~~\includegraphics[width=2.1in]{figures/15_10_3_2_01}~~\includegraphics[width=2.1in]{figures/15_10_5_3_01}\\
\caption{Comparison of LSE, RR.LS and RR.CC. The three panels in each figure correspond to sample sizes 200, 400 and 1000 respectively. The errors are generated according to Setting I.}
\label{fig:simu1}
\end{figure}
\begin{figure}[h]
\centering
\includegraphics[width=2.1in]{figures/6_4_1_1_02}~~\includegraphics[width=2.1in]{figures/6_4_3_2_02}~~\includegraphics[width=2.1in]{figures/6_4_5_3_02}\\
\includegraphics[width=2.1in]{figures/9_6_1_1_02}~~\includegraphics[width=2.1in]{figures/9_6_3_2_02}~~\includegraphics[width=2.1in]{figures/9_6_5_3_02}\\
\includegraphics[width=2.1in]{figures/15_10_1_1_02}~~\includegraphics[width=2.1in]{figures/15_10_3_2_02}~~\includegraphics[width=2.1in]{figures/15_10_5_3_02}\\
\caption{Comparison of LSE, RR.LS and RR.CC. The three panels in each figure correspond to sample sizes 200, 400 and 1000 respectively. The errors are generated according to Setting II.}
\label{fig:simu2}
\end{figure}
In the second part, we consider the coverage probabilities of the confidence intervals based on Theorem~\ref{thm:ls}, Theorem~\ref{thm:cc} and Corollary~\ref{cor:svd}. In this experiment the true ranks are fixed at $k_1=3$ and $k_2=2$. We run simulations 1000 times for sample size $T=200,400,1000$, and consider the cases of
dimension $(d_1,d_2)=(6,4),(9,6),(15,10)$ and $\rho=0.25,0.5,0.75$. For RR.LS, the
error covariance matrix settings (I) and (II) are considered. For RR.CS, we consider two `correct' settings (II') and (II), where in (II') we use $\Sigma_e=\boldsymbol{I}_{d_2}\otimes \boldsymbol{I}_{d_1}$.
The confidence intervals of the entries of the matrices $\boldsymbol{A}_1,\boldsymbol{A}_2, \boldsymbol{U}_1,\boldsymbol{V}_1, \boldsymbol{U}_2,\boldsymbol{V}_2$ are constructed. Table \ref{tab:ci} shows the percentage that the true parameters fall within their corresponding marginal 95\% confidence intervals. Each percentage records the average empirical coverage over all involved matrix entries.
It can be seen from the table that the coverage is quite accurate, especially when the sample size is large ($T=1000$). The empirical coverage probabilities are closer to the nominal ones under Setting (I) for RR.LS and Setting (II') for RR.CS than those under Setting (II). This is probably due to the fact that the singulars values of $\Sigma_e$ under Setting (II) are more spread out, making $\Sigma_e$ more different from the scalar matrix.
\begin{table}[htbp]
\centering
\small
\begin{tabular}{c|cc|rrr | rrr | rrr | rrr}
\toprule
& & & \multicolumn{6}{c|}{RR.LS} & \multicolumn{6}{c}{RR.CC} \\
\hline
& & {Setting} &\multicolumn{3}{c|}{I} & \multicolumn{3}{c|}{II} &\multicolumn{3}{c|}{II$'$} & \multicolumn{3}{c}{II} \\
\hline
& \multicolumn{2}{r|}{$T$} & $200$ & $400$ & $1000$ & $200$ & $400$ & $1000$ & $200$ & $400$ & $1000$ & $200$ & $400$ & $1000$\\
&$\rho$ &$(d_1,d_2)$ & &&&&&&&&&&&\\
\hline
\multirow{9}{*}{$(\boldsymbol{A}_1,\boldsymbol{A}_2)$} &\multirow{3}{*}{$0.75$}
& $(6,4)$ & 93.5 & 93.8 & 94.3 & 92.1 & 92.3 & 92.5 & 93.9 & 94.2 & 95.0 & 91.0 & 91.3 & 91.3 \\
& & $(9,6)$ & 94.3 & 94.6 & 95.3 & 92.2 & 92.6 & 92.8 & 94.0 & 94.4 & 95.0 & 91.8 & 92.1 & 92.4 \\
& & $(15,10)$ & 94.5 & 94.7 & 95.5 & 92.5 & 92.8 & 93.3 & 93.8 & 94.1 & 94.9 & 92.2 & 92.5 & 93.0 \\
\cline{2-15}
& \multirow{3}{*}{$0.5$}
& $(6,4)$ & 93.5 & 93.8 & 94.5 & 91.4 & 91.9 & 92.1 & 93.7 & 94.0 & 94.9 & 90.5 & 91.0 & 91.3 \\
& & $(9,6)$ & 94.3 & 94.6 & 95.4 & 91.7 & 92.4 & 92.6 & 93.8 & 94.3 & 95.0 & 91.2 & 92.0 & 92.2 \\
& & $(15,10)$ & 94.4 & 94.7 & 95.4 & 92.3 & 92.8 & 93.2 & 93.7 & 94.1 & 94.9 & 92.0 & 92.5 & 92.9 \\
\cline{2-15}
& \multirow{3}{*}{$0.25$}
& $(6,4)$ & 92.8 & 93.7 & 94.6 & 88.0 & 90.1 & 91.3 & 92.6 & 93.6 & 94.7 & 87.2 & 89.5 & 90.8 \\
& & $(9,6)$ & 93.3 & 94.4 & 95.3 & 88.8 & 90.9 & 91.9 & 92.8 & 93.8 & 94.8 & 88.4 & 90.7 & 91.8 \\
& & $(15,10)$ & 93.8 & 94.5 & 95.3 & 90.6 & 92.0 & 92.9 & 93.1 & 93.8 & 94.8 & 90.4 & 91.8 & 92.6 \\
\hline
\multirow{9}{*}{$(\boldsymbol{U}_1,\boldsymbol{V}_1)$} &\multirow{3}{*}{$0.75$}
& $(6,4)$ & 94.4 & 94.2 & 94.2 & 92.4 & 92.7 & 92.5 & 94.4 & 94.3 & 94.5 & 91.3 & 91.5 & 91.1 \\
& & $(9,6)$ & 95.0 & 95.0 & 95.1 & 92.4 & 92.8 & 92.6 & 94.5 & 94.4 & 94.8 & 92.0 & 92.2 & 92.3 \\
& & $(15,10)$ & 94.9 & 95.2 & 95.3 & 92.7 & 92.8 & 93.3 & 94.2 & 94.5 & 94.8 & 92.8 & 93.0 & 93.3 \\
\cline{2-15}
& \multirow{3}{*}{$0.5$}
& $(6,4)$ & 94.6 & 94.2 & 94.1 & 92.9 & 92.5 & 92.2 & 94.5 & 94.3 & 94.2 & 91.8 & 91.5 & 91.1 \\
& & $(9,6)$ & 95.1 & 95.0 & 95.1 & 92.7 & 93.0 & 92.6 & 94.8 & 94.5 & 94.9 & 92.2 & 92.4 & 92.3 \\
& & $(15,10)$ & 94.9 & 95.0 & 95.4 & 92.8 & 92.7 & 93.0 & 94.1 & 94.4 & 94.8 & 92.7 & 92.6 & 92.8 \\
\cline{2-15}
& \multirow{3}{*}{$0.25$}
& $(6,4)$ & 95.2 & 94.8 & 94.7 & 93.0 & 92.6 & 92.4 & 95.0 & 94.6 & 94.4 & 92.3 & 91.7 & 91.6 \\
& & $(9,6)$ & 95.3 & 95.3 & 95.3 & 92.5 & 92.8 & 92.6 & 94.9 & 94.8 & 94.9 & 92.0 & 92.4 & 92.2 \\
& & $(15,10)$ & 95.0 & 95.0 & 95.1 & 92.4 & 92.4 & 92.6 & 94.3 & 94.2 & 94.5 & 92.3 & 92.3 & 92.6 \\
\hline
\multirow{9}{*}{$(\boldsymbol{U}_2,\boldsymbol{V}_2)$} &\multirow{3}{*}{$0.75$}
& $(6,4)$ & 94.1 & 94.3 & 94.2 & 91.6 & 91.7 & 91.1 & 94.5 & 94.5 & 94.3 & 91.2 & 91.3 & 90.6 \\
& & $(9,6)$ & 94.7 & 94.6 & 95.5 & 92.0 & 91.9 & 92.5 & 94.4 & 94.2 & 94.8 & 91.9 & 92.1 & 92.9 \\
& & $(15,10)$ & 95.1 & 95.6 & 95.3 & 92.6 & 93.2 & 93.1 & 94.3 & 95.0 & 94.7 & 93.2 & 93.7 & 93.6 \\
\cline{2-15}
& \multirow{3}{*}{$0.5$}
& $(6,4)$ & 94.3 & 94.4 & 94.6 & 90.9 & 91.4 & 91.0 & 94.1 & 94.7 & 94.6 & 90.8 & 91.1 & 90.7 \\
& & $(9,6)$ & 94.8 & 94.7 & 95.3 & 91.8 & 91.8 & 92.3 & 94.4 & 94.2 & 94.5 & 91.5 & 91.9 & 92.4 \\
& & $(15,10)$ & 95.0 & 95.5 & 95.2 & 92.4 & 93.1 & 92.8 & 94.2 & 94.9 & 94.6 & 92.7 & 93.4 & 93.1 \\
\cline{2-15}
& \multirow{3}{*}{$0.25$}
& $(6,4)$ & 93.7 & 93.7 & 94.8 & 88.8 & 89.8 & 90.6 & 93.3 & 93.9 & 94.7 & 88.3 & 89.9 & 91.0 \\
& & $(9,6)$ & 93.8 & 94.7 & 95.3 & 89.4 & 91.0 & 91.9 & 93.5 & 94.1 & 94.4 & 89.2 & 91.2 & 92.1 \\
& & $(15,10)$ & 94.6 & 95.6 & 95.1 & 91.0 & 92.5 & 92.4 & 93.7 & 94.8 & 94.5 & 91.5 & 92.9 & 92.7 \\
\bottomrule
\end{tabular}
\caption{Empirical coverage probabilities (in percentage) of the 95\% confidence intervals.}
\label{tab:ci}
\end{table}
The third part of the simulation considers the performance of the rank determination procedure using the joint EBIC \eqref{eq:bic} and the separate EBIC. Simulations are conducted under various configurations of the sample size, dimensions, ranks and signal strength $\rho$, and the empirical probabilities of selecting the correct ranks out of 100 repetitions
are recorded. When the true ranks are $k_1=k_2=1$, or when the signal strength $\rho$ is not too small ($\rho\geq 0.1$), both selection procedures are able to determine the ranks
perfectly. Therefore, we choose to report in Table~\ref{tab:bic} only the results for the configurations that are more challenging, with $\rho=.15$ when the true ranks are $(3,2)$, and $\rho=.25$ when the true ranks are $(5,3)$. A closer look of the simulation results (not shown in the table) reveals that, in these very low signal to noise ratio cases, both EBIC procedures tend to select ranks smaller than the true ranks.
However, larger sampling sizes significantly enhance the performance. Moreover, when the autocorrelation strength $\rho$ is larger than those reported in Table~\ref{tab:bic}, both procedures make nearly perfect choices of the ranks for all configurations and covariance settings. We also note that the performances of the joint and separate procedures are almost the same.
It is also observed that the performance under the error covariance setting (II) is worse than that under Setting (I).
This is due to the design of $\Sigma_e$ in these two settings. The eigenvalues of $\Sigma_e$ spread over $[1,10]$ in Setting (I), and over $[1,25]$ in Setting (II). Therefore, both the estimation and model selection are more challenging under Setting (II).
We also experiment with using rolling forecasting to choose the ranks. We consider the range $1\le r_i \le \min\{d_i,k_i+2\}$, $i=1,2$, as the candidate set of $\operatorname{rank}(\boldsymbol{A}_i)$. For each configuration of $(\rho,T,k_1,k_2,d_1,d_2)$, we choose $T/4$ as the rolling forecast origin, calculate the entry-wise squared forecast error (SFE) of the one-step ahead prediction of $\boldsymbol{X}_{s+1}$, $T/4\le s\le T$, then take the average over the $d_1d_2$ series and over the time,
\[
\hbox{MSFE}(r_1,r_2)=\frac{1}{d_1d_2(3T/4)}\sum_{s=T/4}^{T-1}\sum_{i=1}^{d_1}\sum_{j=1}^{d_2}
\left|\hat \boldsymbol{X}_{s+1}^{(r_1,r_2)}[i,j]-\boldsymbol{X}_{s+1}[i,j]\right|^2.
\]
The estimated ranks $(\hat k_1,\hat k_2)$ is the pair $(r_1,r_2)$ with the smallest MSFE$(r_1,r_2)$. Table~\ref{tab:forecast} shows the proportion of the correct selection out of 100 repetitions.
It is seen that, although rolling forecast criterion still performs very well in most cases, it has a much higher variability than EBIC. It performs better than EBIC in the case $\rho=.25$ and $(k_1,k_2)=(5,3)$, when the sample size is small.
\begin{table}[H]
\centering
\begin{tabular}{cc r rrr r rrr}
\toprule
\multirow{2}{*}{}&\multirow{2}{*}{} && \multicolumn{3}{c}{$(r_1,r_2)=(3,2)$, $\rho=.15$} && \multicolumn{3}{c}{$(r_1,r_2)=(5,3)$, $\rho=.25$} \\
\cmidrule(lr){4-6}\cmidrule(lr){8-10}
& $(d_1,d_2)$ && 200 & 400 & \phantom{100}1000 && 200 & 400 & \phantom{100}1000 \\
\midrule
\multirow{3}{*}{I}
& $(6,4)$ && (.01, .01) & (.67, .68) & (1, 1) && (.36, .36) & (.96, .96) & (1, 1) \\
& $(9,6)$ && (.15, .15) & (.90, .92) & (1, 1) && (.07, .07) & (.71, .71) & (1, 1) \\
& $(15,10)$ && (1, 1) & (1, 1) & (1, 1) && (.00, .00) & (.47, .48) & (1, 1) \\
\midrule
\multirow{3}{*}{II}
& $(6,4)$ && (.02, .04) & (.38, .39) & (.99, .99) && (.25, .23) & (.92, .92) & (1, 1) \\
& $(9,6)$ && (.00, .00) & (.45, .47) & (1, 1) && (.10, .09) & (.70, .69) & (1, 1) \\
& $(15,10)$ && (.98, .99) & (1, 1) & (1, 1) && (.00, .00) & (.03, .04) & (1, 1) \\
\bottomrule
\end{tabular}
\caption{Empirical probabilities of the correct rank selection by the EBIC. I, II stand for different covariance structures of $\boldsymbol{E}_t$. For each cell, two numbers correspond to the joint and separate selections respectively. The second row shows the sample sizes.}
\label{tab:bic}
\end{table}
\begin{table}[H]
\centering
\begin{tabular}{ccc r rrr r rrr r rrr}
\toprule
\multirow{2}{*}{}&\multirow{2}{*}{}&\multirow{2}{*}{} && \multicolumn{3}{c}{$(1,1)$} && \multicolumn{3}{c}{$(3,2)$} && \multicolumn{3}{c}{$(5,3)$} \\
\cmidrule(lr){5-7}\cmidrule(lr){9-11}\cmidrule(lr){13-15}
&& $(d_1,d_2)$ && 200 & 400 & 1000 && 200 & 400 & 1000 && 200 & 400 & 1000 \\
\midrule
\multirow{6}{*}{$\rho=.5$} &\multirow{3}{*}{I}
& $(6,4)$ && 0.94 & 0.98 & 0.98 && 0.86 & 0.91 & 0.92 && 0.68 & 0.74 & 0.78 \\
&& $(9,6)$ && 0.99 & 1 & 1 && 0.97 & 0.99 & 0.99 && 0.94 & 0.98 & 1 \\
&& $(15,10)$ && 1 & 1 & 1 && 0.99 & 1 & 1 && 1 & 1 & 1 \\
\cline{2-15}
&\multirow{3}{*}{II}
& $(6,4)$ && 0.95 & 0.97 & 0.99 && 0.87 & 0.88 & 0.92 && 0.74 & 0.77 & 0.79 \\
&& $(9,6)$ && 0.99 & 0.99 & 1 && 0.96 & 0.98 & 0.98 && 0.94 & 0.95 & 0.96 \\
&& $(15,10)$ && 1 & 1 & 1 && 0.99 & 1 & 1 && 1 & 1 & 1 \\
\midrule
\multirow{6}{*}{$\rho=.25$}&\multirow{3}{*}{I}
& $(6,4)$ && 0.95 & 0.98 & 0.98 && 0.81 & 0.88 & 0.90 && 0.62 & 0.68 & 0.74 \\
&& $(9,6)$ && 1 & 1 & 1 && 0.96 & 0.98 & 0.99 && 0.94 & 0.98 & 0.98 \\
&& $(15,10)$ && 1 & 1 & 1 && 1 & 1 & 1 && 0.96 & 1 & 1 \\
\cline{2-15}
&\multirow{3}{*}{II}
& $(6,4)$ && 0.96 & 0.97 & 0.98 && 0.44 & 0.79 & 0.93 && 0.64 & 0.76 & 0.82 \\
&& $(9,6)$ && 0.99 & 1 & 1 && 0.69 & 0.98 & 0.98 && 0.95 & 0.95 & 0.95 \\
&& $(15,10)$ && 1 & 1 & 1 && 0.99 & 1 & 1 && 0.94 & 1 & 1 \\
\bottomrule
\end{tabular}
\caption{Empirical probabilities of the correct rank selection by rolling forecast. I, II stands for different covariance structures of $\boldsymbol{E}_t$}
\label{tab:forecast}
\end{table}
\subsection{Example}
\label{sec:example}
We use the RRMAR model to study the eight key short term economic indicators (116 quarters, from 1991 Q1 to 2019 Q4) from ten countries. The data is downloaded from Organisation for Economic Co-operation and Development (OECD, {\tt https://www.oecd.org/}). The 8 indicators are Consumer Price Index (CPI, growth rate), GDP (growth rate), 3-month interbank Interest Rate (IR3, difference), Long Term government bond yield (IRLT, difference), International Trade total Export Value (ITEX, growth rate) and Import Value (ITIM, growth rate), Total Industrial PRoduction excludingn construction (PRTI, growth rate) and Total Manufacturing PRoduction (PRTM, growth rate). The 10 countries are Australia (AUS), Austria (AUT), Canada (CAN), France (FRA), Germany (DEU), Netherlands (NLD), Norway (NOR), Sweden (SWE), United Kingdom (GBR) and United States (USA). All the 80 series have been centered before attempting the model. We also standardize each indicator across all the countries, i.e. the 10 series corresponding to each indicator have an overall standard deviation 1.
The EBIC \eqref{eq:bic} selects the ranks as $k_1=1$ and $k_2=4$ for this data set. Using the RR.CC approach, the
estimated leading singular vectors $\hat\boldsymbol{U}_i$ and $\hat\boldsymbol{V}_i$ of ${\boldsymbol{A}_{i}}$ ($i=1,2$), and their corresponding
estimated standard errors are shown in Tables~\ref{table:A1} and \ref{table:A2}.
Entries which are not significant at 10\% level are shown in light gray color.
\begin{table}[H]
\begin{center}
\begin{tabular}{rrrrrrrrrrr}
\hline
& AUS & AUT & CAN & DEU & FRA & GBR & NLD & NOR & SWE & USA\\
$\hat\boldsymbol{U}_1'$ & 0.25 & 0.30 & 0.35 & 0.35 & 0.31 & 0.28 & 0.30 & 0.28 & 0.40 & 0.32 \\
{s.e.} & 0.02 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 & 0.01 & 0.02 & 0.01 & 0.01 \\ \hline
$\hat\boldsymbol{V}_1'$ & {\color{lightgray}0.09} & 0.49 & {\color{lightgray}0.05} & {\color{lightgray}0.08} & {\color{lightgray}-0.04} & 0.61 & -0.27 & {\color{lightgray}0.01} & 0.27 & 0.47\\
{s.e.} & 0.10 & 0.12 & 0.13 & 0.15 & 0.17 & 0.10 & 0.11 & 0.07 & 0.1 & 0.13 \\ \hline
\end{tabular}
\caption{Estimated singular vectors of the coefficient matrix $\boldsymbol{A}_1$, with their corresponding estimated standard errors.}\label{table:A1}
\end{center}
\end{table}
\begin{table}[H]
\begin{center}
\begin{tabular}{rrrrrrrrr}
\hline
& CPI & GDP & IR3 & IRLT & ITEX & ITIM & PRTI & PRTM\\
$\hat\boldsymbol{U}_2'[1,]$& 0.25 & 0.44 & 0.15 & 0.21 & 0.29 & 0.29 & 0.45 & 0.56\\
{s.e.} & 0.05 & 0.03 & 0.05 & 0.04 & 0.03 & 0.03 & 0.03 & 0.02 \\ \hline
$\hat\boldsymbol{V}_2'[1,]$ & {\color{lightgray} -0.05} & 0.84 & -0.43 & 0.21 & 0.17 & {\color{lightgray} 0.14} & {\color{lightgray} -0.03} & {\color{lightgray} 0.08}\\
{s.e.} & 0.06 & 0.04 & 0.05 & 0.07 & 0.10 & 0.10 & 0.08 & 0.09\\ \hline
\end{tabular}
\caption{Estimated leading singular vectors of the coefficient matrix $\boldsymbol{A}_2$, with their corresponding estimated standard errors.}\label{table:A2}
\end{center}
\end{table}
The implication of using $k_1=1$ and $k_2=4$ is that the observations in the previous quarter form a 4 composite indexes $\boldsymbol{f}_{t-1}:=\boldsymbol{V}_1'\boldsymbol{X}_{t-1}\boldsymbol{V}_2$, and the conditional expectation $\mathbb{E}(\boldsymbol{X}_t\mid\boldsymbol{X}_{t-1})$ is given by $d_{11}\cdot \boldsymbol{U}_1 \boldsymbol{f}_{t-1}\boldsymbol{D}_2\boldsymbol{U}_2'$, where $d_{11}$ is the largest singular value of $\boldsymbol{A}_1$, and $\boldsymbol{D}_2$ is the $4\times 4$ diagonal matrix containing the singular values of $\boldsymbol{A}_2$, see \eqref{eq:marrr_fac} as well. We report the estimated leading singular vectors of $\boldsymbol{U}_2$ and $\boldsymbol{V}_2$ in Table~\ref{table:A2}. It is very interesting to observe that when the indicators are combined to form the first element of $\boldsymbol{f}_{t-1}$ using $\hat\boldsymbol{V}_2[,1]$, GDP is most dominant, followed by IR3, IRLT and ITEX, while CPI, ITIM, PRTI, PRTM have less importance. It is also worth noting that GDP has a positive coefficient, and IR has a negative one in $\hat{\boldsymbol{V}}_2[,1]$. When the countries are combined using $\hat\boldsymbol{V}_1$ (see Table~\ref{table:A1}), 5 countries play more important roles (i.e. the 5 significant entries in $\hat\boldsymbol{V}_1$), and both USA and GBR are among them. At time $t$, all indicators from all countries significantly load on $\boldsymbol{f}_{t-1}$.
\begin{table}[H]
\begin{center}
\begin{tabular}{lrrrrrrrrrr}\hline
&iAR(1)&VAR(1) & PROJ&LSE&MLE&RR.LS & RR.CC\\\hline
MSE & 0.5258 & 3.6988 & 1.6362 & 0.5734 & 0.5187 & 0.5816 & 0.5016 \\
\# par & 80 & 6,400 & 163 & 163 & 163 & 102 & 102 \\\hline
\end{tabular}
\caption{Out-sample prediction performance comparison of various models for the matrix series of 8 indicators from 10 OECD countries.}\label{table:FFpred}
\end{center}
\end{table}
The model selected by the EBIC does not leading to the best rolling forecast performance. For this purpose, we consider the model with $k_1=5$ and $k_2=2$, which leads to almost the smallest mean sqaured rolling forecast error, but is still much more parsimonious than the full rank MAR model. The mean squared errors of the one-step rolling forecast of the last 8 years are summarized in Table~\ref{table:FFpred}, in which we compare the following seven methods.
\begin{enumerate}
\item {\bf iAR(1):\ } Fit an AR(1) model to each individual series.
\item {\bf VAR(1):\ } Fit a VAR(1) model to $\mathrm{vec}(\boldsymbol{X}_t)$.
\item {\bf PROJ, LSE, MLE:\ } Fit the MAR(1) model (without rank constraint) to $\boldsymbol{X}_t$ using projection, least squares and MLE methods. See \cite{chen2021autoregressive} for details.
\item {\bf RR.LS:\ } Reduced-rank MAR(1) model, fitted by least squares.
\item {\bf RR.CC:\ } Reduced-rank MAR(1) model, fitted by MLE under the assumption \eqref{eq:ocov}.
\end{enumerate}
From Table~\ref{table:FFpred}, it is seen that VAR(1) model involves a $80\times 80$ coefficient matrix and significantly overfits the data, with the worst out-sample prediction performance. The MAR and RRMAR models with MLE, have better performance than fitting each individual series separately (iAR(1)).
Comparing to the MAR model without rank constraint, the reduced-rank model estimated by MLE (RR.CC) has the smallest rolling forecast error with less parameters (163 vs 102).
\section{Conclusion}
\label{sec:con}
We introduce the reduced-rank matrix autoregressive model, which relies on an autoregressive term involving bilinear coefficient matrices, and assumes rank deficiency of the coefficient matrices.
Comparing with the MAR model without the low rank structure \citep{chen2021autoregressive}, the RRMAR model involves a greatly reduced number of parameters and leads to more efficient estimation. On the other hand, we use the MAR estimates as the warm-start initial values for the estimation of the RRMAR model.
Both LSE and MLE are studied, where the latter is considered under an additional assumption that the covariance tensor of the error matrix is separable. We propose to use extended BIC to select the ranks of the coefficient matrices. Our numerical analysis suggests that even if the separability assumption on the covariance tensor does not hold, MLE still has reasonable and almost equally good performance, comparing with LSE. On the other hand, MLE can perform much better when that assumption does stand. Therefore, we would recommend the use of MLE in practice.
There are a number of directions to extend the study of the reduced-rank autoregressive model. For example the conditional mean can involve multiple terms of the form $\sum_{j=1}^J \boldsymbol{A}_{j1}\boldsymbol{X}_{t-1}\boldsymbol{A}_{j2}'$, and multiple lagged terms $\boldsymbol{X}_{t-1},\ldots,\boldsymbol{X}_{t-p}$. The model can be extended for tesor time series as well. More importantly, the asymptotic analysis has been carried out for the fixed dimensional case in the current paper. It is interesting and important to study the model under the high dimensional paradigm. In particular, we would like to understand: (i) what are the convergence rates of $\hat\boldsymbol{A}_i$; and (ii) how to obtain initial estimates of $\boldsymbol{A}_i$ to start the alternating algorithm. To select the ranks, either the information criterion based procedure can be adapted to account for the high dimensionality, or the singular (eigen-)value based approach \citep{lam:2012,wang:2018} can be employed. The relationship between the reduced-rank tensor autoregressive model and the dynamic tensor factor model \citep{chen2019factor} is also worth exploring.
{In this paper the theoretical results are obtained under the fix data dimension assumption. It is an interesting and important problem to expand the theoretical results to high and diverging dimension setting. It is a challenging problem and may require additional structure of the model. Due to the stationary condition $\rho(\boldsymbol{A}_2\otimes \boldsymbol{A}_1)<1$ needed for the autoregressive model, the signal to noise ratio is constrained, different from typical regression models. This is similar to the simple AR(1) model $x_t=\phi x_{t-1}+e_t$, in which the signal to noise ratio is always $\phi^2/(1-\phi^2)$ no matter how large or small the noise variance is. In order to achieve consistency results for the diverging dimensional setting, the reduced-rank structure is not sufficient. It seems that additional sparsity structure or other type of structure is needed. We are currently investigating this problem. }
\bigskip
\noindent {\bf Acknowledgement.}\ \ We thank the Editor, Associate Editor and Referees for their constructive comments and suggestions.
\bibliographystyle{apalike}
\bibliography{mybib}
\clearpage