Dynamic Matrix Factor Models for High Dimensional Time Series

Dynamic Matrix Factor Models for High Dimensional Time Series

Introduction

Tensor data are widely collected in various fields and have become increasingly popular. When the data is collected over time, {one of the modes often represents time, leading to a unique treatment of this temporal dimension.} This type of data is referred to as matrix time series data when the data collected at each time point forms a matrix. Matrix time series data can be found in many applications, such as economics where economic indicators are reported from different countries on a monthly or quarterly basis. These data can be modeled as a matrix time series, where the countries are represented by rows and the economic indicators by columns. In this case, the variables in the same column will be strongly correlated due to reporting from the same country, and the variables in the same row will also be correlated due to the influence of different countries. Thus, it is more reasonable to consider the matrix time series as a whole, preserving the row-wise and column-wise relationships. Another example is a sequence of images, which can be modeled as matrices observed over time. Using matrix models is essential as the spatial structure would be lost if modeled in vector or univariate form. Recent studies on matrix-valued time series include work by wang2019 and chen2022factor, which analyze factor models for matrix and tensor time series, and studies by Hoff2015, chen2021autoregressive, and han2023rr, which explore autoregressive models in bilinear matrix form.

One of the major challenges in analyzing matrix time series is its high dimensionality. The curse of dimensionality is a notorious problem for high-dimensional data analysis, as it makes it difficult to analyze the data and extract meaningful insights. To overcome this challenge, a common approach is to adopt a factor approach. The factor model describes the co-movement of data with lower-dimensional latent factors. It is based on the idea that high-dimensional data can be reduced to lower-dimensional representations by extracting a small number of underlying factors that capture the most important information in the data.

Factor models are widely used in a variety of fields, including finance, economics, and social sciences. {Numerous studies have focused on factor models for time series, with seminal works like bai2002, bai2003, stock2005 exploring vector time series factor models.} {Extensions to matrix and tensor factor models have been introduced in wang2019, chen2022factor, han2020, han2023cp, chen2023statistical, chen2024semi, among others. In certain studies, the latent factor process is presumed to follow specific dynamic structures, typically autoregressive} Molenaar1985ADF,tiao1989,engle1995,forni2000,Jasiak2001,hallin2016,fan2013. {A critical aspect of building a factor model is to determine the number of factors. This topic has been thoroughly explored in vector factor models, such as bai2002,lam2012,hallin2007,amengual2007.} For matrix and tensor factor models, han2022rank proposed information criteria and eigenvalue ratio based methods that use specified penalty functions, {while} lam2021rank,chen2024rank proposed a correlation thresholding method without using penalty functions and is free of parameter tuning.

While there have been numerous studies on dynamic factor models for panel data or vector time series, there is still a lack of literature on matrix-valued dynamic factor models. wang2019 introduces a factor model for matrix-valued time series, but it lacks an assumed dynamic structure for the factor process, making it difficult to predict the factor model without additional specifications. Our goal is to extend the matrix factor model in wang2019,chen2022factor by including a dynamic structure for the latent matrix factor process. This approach is similar to dynamic factor models for vector time series in stock2016. Specifically, we assume that the matrix factor process follows a matrix autoregressive (MAR) structure as in Hoff2015, chen2021autoregressive, represented by ${{\boldsymbol{F} }}_t = {{\boldsymbol{A} }}_1{{\boldsymbol{F} }}_{t-1}{{\boldsymbol{A} }}_2^{\top} + {{\boldsymbol{E} }}_t$. The use of this dynamic structure brings two main benefits: (i) It allows for straightforward prediction of the factor model through a two-step computational efficient procedure. (ii) The autoregressive dynamic structure provides a clear and insightful understanding of the hidden dynamics of the factor process. It should be noted that our work differs from that of yuan2023two, which is based on maximum likelihood estimation and explores a different matrix factor model.

The rest of this paper is organized as follows. In section 2 we specify the proposed dynamic matrix factor model, with detailed discussion of the motivations using such a model. In Section 3, we propose a two step estimation procedure for dynamic matrix factor model. The prediction procedure is also discussed. The theoretical properties of the estimators are presented in Section 4. Empirical studies are presented in Section 5 to demonstrate the finite sample properties of the estimators. In Section 6, an analysis of a transport network example is shown to demonstrate the usefulness of the model. Section 7 concludes.

Dynamic Matrix Factor Models

The proposed Dynamic Matrix Factor Model (DMFM) is a combination of two components, the matrix factor model of wang2019 and the matrix autoregressive model of chen2021autoregressive. The core innovation lies in its two-component structure, which uniquely combines the strengths of the two models for dimensionality reduction and time series modeling.

Model setting

The Dynamic Matrix Factor Model (DMFM) takes the following two-component form:

where in the first part (ref), ${{\boldsymbol{X} }}_t \in \mathbb{R}^{d_1\times d_2}$ is the observed matrix time series, ${{\boldsymbol{F} }}_t\in \mathbb{R}^{r_1\times r_2}$ is the unobserved latent factor process, with $r_1 \ll d_1$ and $r_2 \ll d_2$. ${{\boldsymbol{U} }}_1\in \mathbb{R}^{r_1\times d_1}$ and ${{\boldsymbol{U} }}_2\in \mathbb{R}^{r_2\times d_2}$ are orthonormal loading matrices, and $\lambda$ singles out the signal strength. {In the typical strong factor model setting, $\lambda\asymp \sqrt{d_1 d_2}$.} We assume the ${{\boldsymbol{E} }}_t\in \mathbb{R}^{d_1\times d_2}$ is a white noise series so that ${\rm Cov}(\hbox{\rm vec}({{\boldsymbol{E} }}_{t_1}),\hbox{\rm vec}({{\boldsymbol{E} }}_{t_2}))=0$ for $t_1\neq t_2$. Contemporary correlations among the elements of ${{\boldsymbol{E} }}_t$ are allowed, with ${\rm Cov}(\hbox{\rm vec}({{\boldsymbol{E} }}_t)) = (d_1d_2)^{-1} \boldsymbol{\Sigma}_E$.

The second part ((ref)) of DMFM aims to describe the dynamics of factors after dimension reduction by the matrix factor model ((ref)). We assume that the latent factor process ${{\boldsymbol{F} }}_t$ follows an autoregressive model, with the coefficient matrices ${{\boldsymbol{A} }}_1\in \mathbb{R}^{r_1\times r_1}$ and ${{\boldsymbol{A} }}_2\in \mathbb{R}^{r_2\times r_2}$ being the left and right autoregressive coefficients. ${{\boldsymbol{\xi} }}_t\in \mathbb{R}^{r_1\times r_2}$ are innovation series which are assumed to be white noise but also with possible contemporary correlations among its elements, with ${\rm Cov}(\hbox{\rm vec}({{\boldsymbol{\xi} }}_t)) = (r_1r_2)^{-1}{{\boldsymbol{\Sigma} }}_{{\boldsymbol{\xi} }}$. In addition, we assume that ${{\boldsymbol{E} }}_t$ and ${{\boldsymbol{\xi} }}_{t+h}$ are uncorrelated for all $h\in \mathbb{Z}$, making ${{\boldsymbol{E} }}_t$ and ${{\boldsymbol{F} }}_{t+h}$ uncorrelated for all $h \in \mathbb{Z}$ as well.

The first component (ref) is a Matrix Factor Model wang2019. It projects the temporal dynamics of the high-dimensional matrix time series ${{\boldsymbol{X} }}_t$ into the temporal dynamics of the low-dimensional latent factor ${{\boldsymbol{F} }}_t$. This transformation is crucial for computational efficiency and for the extraction of meaningful latent features from the raw data. The second component (ref) is a Matrix Autoregressive Model chen2021autoregressive that provides a parametric model for the dynamics of the latent factors. It makes the DMFM a generative model with a relatively simple prediction capability which the matrix factor model (ref) alone does not have. wang2024high considered a low rank tensor AR model and uses a nuclear norm penalty to enforce the low rank structure and optimization algorithms for estimation. It is quite different from our approach.

Relation to Reduced Rank MAR Model

A reduced rank MAR Model was developed in han2023rr, which has a close connection with dynamic matrix factor model. The reduced rank MAR Model takes the form

where ${{\boldsymbol{A} }}_{il}$ and ${{\boldsymbol{A} }}_{ic}$ are both $d_i\times r_i$ full rank matrices. It is essentially the MAR model ${{\boldsymbol{X} }}_t={{\boldsymbol{A} }}_1 {{\boldsymbol{X} }}_{t-1} {{\boldsymbol{A} }}_2^{\top} +{{\boldsymbol{E} }}_t$, with ${{\boldsymbol{A} }}_i = {{\boldsymbol{A} }}_{il}{{\boldsymbol{A} }}_{ic}^{\top}$ ($i=1,2$) being a $d_i\times d_i$ matrix of rank $r_i$.

Let ${{\boldsymbol{F} }}^*_t = {{\boldsymbol{A} }}_{1c}^{\top} {{\boldsymbol{X} }}_{t-1}{{\boldsymbol{A} }}_{2c}$ and notice that \[ {{\boldsymbol{F} }}^*_t = {{\boldsymbol{A} }}_{1c}^{\top} {{\boldsymbol{X} }}_{t-1}{{\boldsymbol{A} }}_{2c} ={{\boldsymbol{A} }}_{1c}^{\top}[{{\boldsymbol{A} }}_{1l}{{\boldsymbol{A} }}_{1c}^{\top} {{\boldsymbol{X} }}_{t-2} {{\boldsymbol{A} }}_{2c}{{\boldsymbol{A} }}_{2l}^{\top}+{{\boldsymbol{E} }}_{t-1}]{{\boldsymbol{A} }}_{2c} ={{\boldsymbol{A} }}^*_{1}{{\boldsymbol{F} }}^*_{t-1}{{\boldsymbol{A} }}_{2}^{*\top}+ {{\boldsymbol{A} }}_{1c}^{\top}{{\boldsymbol{E} }}_{t-1} {{\boldsymbol{A} }}_{2c}, \] where ${{\boldsymbol{A} }}_i^*={{\boldsymbol{A} }}_{ic}^\top{{\boldsymbol{A} }}_{il}$. Hence model (ref) becomes \[ {{\boldsymbol{X} }}_t= {{\boldsymbol{A} }}_{1l}{{\boldsymbol{F} }}_t^*{{\boldsymbol{A} }}_{2l}^{\top} +{{\boldsymbol{E} }}_t, \mbox{\ \ } {{\boldsymbol{F} }}^*_t ={{\boldsymbol{A} }}^*_1 {{\boldsymbol{F} }}^*_{t-1} {{\boldsymbol{A} }}_2^{*\top} + {{\boldsymbol{A} }}_{1c}^{\top}{{\boldsymbol{E} }}_{t-1} {{\boldsymbol{A} }}_{2c}. \] This is very similar to DMFM model in (ref) and (ref), with constraints in the coefficient matrices and shared noise process in ${{\boldsymbol{E} }}_t$ and ${{\boldsymbol{A} }}_{1c}^{\top}{{\boldsymbol{E} }}_{t-1} {{\boldsymbol{A} }}_{2c}$ in the two components.

Similarly, from (ref) and (ref) we have for the DMFM model

This is in a MAR model form, with reduced rank coefficient matrices ${{\boldsymbol{U} }}_i {{\boldsymbol{A} }}_i {{\boldsymbol{U} }}_i^{\top}$, and error process $\lambda {{\boldsymbol{U} }}_1 {{\boldsymbol{\xi} }}_t {{\boldsymbol{U} }}_2^{\top}- {{\boldsymbol{U} }}_1 {{\boldsymbol{A} }}_1 {{\boldsymbol{U} }}_1^{\top} {{\boldsymbol{E} }}_{t-1} {{\boldsymbol{U} }}_2 {{\boldsymbol{A} }}_2^{\top} {{\boldsymbol{U} }}_2^{\top}+{{\boldsymbol{E} }}_t$. The noise process consists of a moving average process involving both ${{\boldsymbol{E} }}_t$ and ${{\boldsymbol{E} }}_{t-1}$, and another independent process involving ${{\boldsymbol{\xi} }}_t$. Even without the noise term ${{\boldsymbol{E} }}_t$, the model is not equivalent to the reduced rank MAR as the noise $\lambda {{\boldsymbol{U} }}_1 {{\boldsymbol{\xi} }}_t {{\boldsymbol{U} }}_2^{\top}$ is now restricted to the space determined by the coefficient matrices ${{\boldsymbol{U} }}_i$.

Estimation and Prediction

Model Formulation

Before introducing the estimation and prediction procedures, we revisit the model formulation (ref) and (ref) and make an equivalent modification to facilitate the discussion. The latent equation (ref) should be understood as the generating mechanism of a MAR(1) model with coefficient matrices ${{\boldsymbol{A} }}_i\ (i=1,2)$ and innovation covariance matrix ${{\boldsymbol{\Sigma} }}_{{\boldsymbol{\xi} }}$, which determines the scale of the factor process ${{\boldsymbol{F} }}_t$. The strength of the factor component is then determined by $\lambda$ in (ref). However, the estimation of ${{\boldsymbol{A} }}_i$ and the subsequent predictions are not affected by the scale of the factor process $\hat{{\boldsymbol{F} }}_t$ since an autoregressive model is entertained in (ref). In other words, if we absorb the strength $\lambda$ into the factors ${{\boldsymbol{F} }}_t$, both the estimation of ${{\boldsymbol{A} }}_i$ and the prediction of ${{\boldsymbol{X} }}_t$ remain unchanged. Therefore, we consider the equivalent formulation of the DMFM model:

where $\lambda$ is absorbed into ${{\boldsymbol{F} }}_t$. Note that the scale of ${{\boldsymbol{\xi} }}_t$ is also multiplied by $\lambda$, comparing to (ref).

Our subsequent discussion on the estimation and prediction will be based on (ref) and (ref). As a result, we will suppress the estimation of $\lambda$ and focus on ${{\boldsymbol{U} }}_i$ and ${{\boldsymbol{A} }}_i$. More specifically, after obtaining the two orthonomal loading matrices $\hat{{\boldsymbol{U} }}_1$ and $\hat{{\boldsymbol{U} }}_2$, we will treat $\hat{{\boldsymbol{F} }}_t=\hat{{\boldsymbol{U} }}_1^\top{{\boldsymbol{X} }}_t\hat{{\boldsymbol{U} }}_2$ as the estimated factors and proceed with the estimation of ${{\boldsymbol{A} }}_i$ using $\hat{{\boldsymbol{F} }}_t$.

The estimation error in $\hat{{\boldsymbol{F} }}_t$ involves not only the error in $\hat{{\boldsymbol{U} }}_i$, but more importantly, the noise ${{\boldsymbol{E} }}_t$ in (ref). One of the major contributions of this paper is to propose a lag-2 moment based estimator of ${{\boldsymbol{A} }}_i$ to alleviate the impact of ${{\boldsymbol{E} }}_t$, and to introduce a corresponding prediction procedure. To facilitate the discussion, we note that if the ${{\boldsymbol{U} }}_i$ are known, the “estimated" factors $\tilde{{\boldsymbol{F} }}_t:={{\boldsymbol{U} }}_1^\top{{\boldsymbol{X} }}_t{{\boldsymbol{U} }}_2$ is $\tilde{{\boldsymbol{F} }}_t={{\boldsymbol{F} }}_t+{{\boldsymbol{U} }}_1^\top{{\boldsymbol{E} }}_t{{\boldsymbol{U} }}_2$. We therefore utilize the following MAR model with the measurement error to introduce the idea and method of estimating ${{\boldsymbol{A} }}_i$.

It is immediately seen that the preceding model is in the state-space form, where ${{\boldsymbol{F} }}_t$ are latent state variable, and $\tilde{{\boldsymbol{F} }}_t$ are observed. The link between the DMFM model (ref) and (ref) and the preceding state space model is that ${{\boldsymbol{\zeta} }}_t={{\boldsymbol{U} }}_1^\top{{\boldsymbol{E} }}_t{{\boldsymbol{U} }}_2$. We use ${{\boldsymbol{\Sigma} }}_\zeta$ to denote the covariance matrix of $\hbox{\rm vec}({{\boldsymbol{\zeta} }}_t)$.

We propose a two-stage estimation procedure for the dynamic matrix factor model. In the first stage, we employ matrix factor model estimation to obtain $\hat{{\boldsymbol{U} }}_i$ and $\hat{\boldsymbol{F}}_t=\hat{{\boldsymbol{U} }}_1^\top{{\boldsymbol{X} }}_t\hat{{\boldsymbol{U} }}_2$. Subsequently, we treat the estimated factor series as observed with measurement error, and use it to estimate the parameters in the matrix autoregressive model. Such a two-stage estimation procedure is commonly used in estimating dynamic factor models stock2016,Jasiak2001,hallin2016,otto2022approximate.

For estimating the matrix factor model component (ref), we use the iTOPUP and iTIPUP methods, as introduced in chen2022factor and han2020. These methods are extended version of principal component analysis (PCA) and are designed for matrix and tensor time series. The resulting estimators are denoted as $\hat{\lambda}$, $\hat{{{\boldsymbol{U} }}}_1$ and $\hat{{{\boldsymbol{U} }}}_2$, respectively.

For estimating the MAR coefficient matrices ${{\boldsymbol{A} }}_i$ in (ref), the PROJ, LSE and MLE estimators in chen2021autoregressive can be applied directly, treating $\hat{{\boldsymbol{F} }}_t$ as the observation. These estimators are obtained by directly fitting a MAR model on $\hat{{{\boldsymbol{F} }}}_t$ from the first stage, thus ignoring the estimation error in $\hat{{\boldsymbol{F} }}_t$ (also see (ref)). However, the estimation error in $\hat{{\boldsymbol{F} }}_t$ can be substantial when the signal to noise ratio is at a low level. This is similar to the case of modeling time series observed with measurement errors. Hence, we propose a lagged estimator that reduces the influence of the errors from the factor estimation. To enhance the estimation of the factors given all the estimated parameters in both parts of DMFM, we also use Kalman Filter to obtain more accurate estimate of the factor process under the state-space formulation (ref) and (ref).

The estimation accuracy is driven by the signal-to-noise ratio (SNR) of the model. There are in fact two relevant SNRs, one related to the estimation of the loadings ${{\boldsymbol{U} }}_i$, and one related to the estimation of ${{\boldsymbol{A} }}_i$. Since we are adopting the existing method of han2020 for the estimation of ${{\boldsymbol{U} }}_i$, we will only introduce the SNR relevant to the estimation of ${{\boldsymbol{A} }}_i$. In view of (ref), we define

Factor estimation and measurement error

The first step is to determine the rank of the factor matrix. To achieve this, we utilize the method proposed in han2022rank, which suggests either the Bayesian Information Criterion (BIC) or the eigenvalue ratio (ER) method for estimating the number of factors. Once the ranks $r_1$ and $r_2$ are determined, we can estimate the loading matrices $\hat{\boldsymbol{U}}_1$ and $\hat{\boldsymbol{U}}_2$, as well as the estimated latent factor process $\hat{\boldsymbol{F}}_t$, using the iTIPUP or iTOPUP matrix factor estimation methods. It is important to note that the estimated $\hat{\boldsymbol{U}}_1$, $\hat{\boldsymbol{U}}_2$, and $\hat{\boldsymbol{F}}_t$ are not unique; however, any choice of these estimates will not impact the prediction results. In the following sections, we will delve deeper into the identification problem. The estimation error in $\hat{{{\boldsymbol{F} }}}_t$ may have significant impact on the second stage estimation. Without considering the dynamics in the factor process, han2020 established that the error rate of $\hat{{{\boldsymbol{F} }}}_t$ is of the order of $O_{\mathbb{P}} \left(\frac{\sigma\sqrt{d_{\max}}}{\lambda\sqrt{T}}+\frac{\sigma}{\lambda} \right)$; {see also (ref)}. When the signal to noise ratio is high, the error is small enough to obtain accurate estimation in the second stage. When the signal to noise ratio is not sufficiently high, we introduce a (higher) lagged estimator to alleviate the influence of estimation error from the estimation of the factors.

Estimation of the coefficient matrices in (ref)

Here we propose two moment estimators of ${{\boldsymbol{A} }}_i$ in (ref). Although the second estimator is designed to deal with un-ignorable measurement errors in $\hat{{{\boldsymbol{F} }}}_t$, in this section we simply assume the first stage factor estimation produces the factor process without error, and develop the estimator of the coefficient matrices in (ref) with the observed true ${{\boldsymbol{F} }}_t$ process. The consequences of the measurement error will become apparent later in our empirical and theoretic investigation of the estimators.

Lag-1 Moment Estimator (LSE)

We consider lag-1 Yule-Walker estimator in the form of

where $\hat{{\boldsymbol{\Gamma} }}_0$ and $\hat{{\boldsymbol{\Gamma} }}_1$ are the lag-0 and lag-1 sample autocovariance matrices of the vectorized factor series $\hbox{\rm vec}({{\boldsymbol{F} }}_t)$, i.e. $\hat{{\boldsymbol{\Gamma} }}_k=T^{-1}\sum_{t=k+1}^T \hbox{\rm vec}({{\boldsymbol{F} }}_t)\hbox{\rm vec}({{\boldsymbol{F} }}_{t-k})^\top$, for $k\geq 0$. Note that $\hat{{{\boldsymbol{\Gamma} }}}_1\hat{{{\boldsymbol{\Gamma} }}}_0^{-1}$ is the Yule-Walker estimator of the coefficient matrix ${{\boldsymbol{\Phi} }}$ in the VAR(1) model $\hbox{\rm vec}({{\boldsymbol{F} }}_t)={{\boldsymbol{\Phi} }}\hbox{\rm vec}({{\boldsymbol{F} }}_{t-1})+\hbox{\rm vec}({{\boldsymbol{\xi} }}_t)$. To estimate ${{\boldsymbol{A} }}_i$ in the MAR model (ref), a straightforward approach is to project $\hat{{{\boldsymbol{\Gamma} }}}_1\hat{{{\boldsymbol{\Gamma} }}}_0^{-1}$ to the cone of Kronecker products through

The estimator in (ref) can be viewed as a projection weighted by $\hat{{\boldsymbol{\Gamma} }}_0^{1/2}$. This weighting scheme is due to the proposition below.

The proof is given in Appendix (ref). Proposition (ref) shows that asymptotic covariance matrix of each row of $\hat{{{\boldsymbol{\Gamma} }}}_1\hat{{{\boldsymbol{\Gamma} }}}_0^{-1}$ is proportional to ${{\boldsymbol{\Gamma} }}_0$, justifying the projection weighted by $\hat{{\boldsymbol{\Gamma} }}_0^{1/2}$ in (ref).

A direct calculation shows that (ref) is equivalent to the least squares problem

which is considered in chen2021autoregressive. Specifically, the solution of (ref) can be found by iteratively updating

until convergence. To initialize, one can use the projection estimator in (ref).

For simplicity, we will refer to the solution of (ref) (equivalently (ref)) as {\bf LSE}. chen2021autoregressive considers both the LSE and the projection estimator (ref), and has shown empirically that the LSE has a superior performance.

Lag-2 Moment Estimator (L2E)

Similar to the lag-1 Yule-Walker estimator in (ref), we can estimate ${{\boldsymbol{A} }}_1$ and ${{\boldsymbol{A} }}_2$ by solving the optimization problem

The $\hat{{{\boldsymbol{\Gamma} }}}_2 \hat{{{\boldsymbol{\Gamma} }}}_1^{-1}$ can be viewed as the lag-2 Yule-Walker estimator of the VAR model for $\hbox{\rm vec}({{\boldsymbol{F} }}_t)$. Similar to the LSE, the estimator in (ref) is obtained by projecting $\hat{{{\boldsymbol{\Gamma} }}}_2 \hat{{{\boldsymbol{\Gamma} }}}_1^{-1}$ to the cone of Kronecker products using the weighting matrix $(\hat{{{\boldsymbol{\Gamma} }}}_1 \hat{{{\boldsymbol{\Gamma} }}}_0^{-1} \hat{{{\boldsymbol{\Gamma} }}}_1^{\top})^{1/2}$. The choice of the weighting matrix is motivated by the following proposition.

The proof is given in Appendix (ref). For simplicity, we will refer to the estimator in (ref) as the {\bf Lag-2 Estimator (L2E)}.

The L2E is designed to handle the measurement error in the observed ${{\boldsymbol{F} }}_t$, which naturally arises for the DMFM, since the factors are estimated as $\hat{{\boldsymbol{F} }}_t=\hat{{\boldsymbol{U} }}_1^\top{{\boldsymbol{X} }}_t\hat{{\boldsymbol{U} }}_2=\hat{{\boldsymbol{U} }}_1^\top{{\boldsymbol{U} }}_1{{\boldsymbol{F} }}_t{{\boldsymbol{U} }}_2^\top\hat{{\boldsymbol{U} }}_2+\hat{{\boldsymbol{U} }}_1^\top{{\boldsymbol{E} }}_t\hat{{\boldsymbol{U} }}_2$, where $\hat{{\boldsymbol{U} }}_1^\top{{\boldsymbol{E} }}_t\hat{{\boldsymbol{U} }}_2$ becomes the measurement error. To fix the idea, we consider the MAR model with measurement error (ref) and (ref), where $\tilde{{\boldsymbol{F} }}_t$ are assumed to be observed. If the LSE (ref) is to be applied to the $\tilde{{\boldsymbol{F} }}_t$ directly, one needs to calculate $\hat{{\boldsymbol{\Gamma} }}_1 \hat{{\boldsymbol{\Gamma} }}_0^{-1}$ first, where the sample autocovariance matrices $\hat{{\boldsymbol{\Gamma} }}_k$ are calculated for the observed $\tilde{{\boldsymbol{F} }}_t$. Recall that ${{\boldsymbol{\Sigma} }}_{{\boldsymbol{\zeta} }}$ is the covariance matrice of $\hbox{\rm vec}({{\boldsymbol{\zeta} }}_t)$. It holds that $\hat{{\boldsymbol{\Gamma} }}_0 \asymp {{\boldsymbol{\Gamma} }}_0 + {{\boldsymbol{\Sigma} }}_{{\boldsymbol{\zeta} }} + O_p(T^{-1/2})$, where ${{\boldsymbol{\Gamma} }}_0:={\rm Var}(\hbox{\rm vec}({{\boldsymbol{F} }}_t))$. The bias ${{\boldsymbol{\Sigma} }}_{{\boldsymbol{\zeta} }}$ will be translated into the LSE of ${{\boldsymbol{A} }}_i$ and will only be negligible if ${{\boldsymbol{\Sigma} }}_{{\boldsymbol{\zeta} }}$ is of an order smaller than $T^{-1/2}$ (see Theorem (ref) and (ref) for more detailed and precise theoretical analysis).

We propose the L2E to alleviate the impact of the measurement error. It is also inspired by the estimator in hannan1963 that leverages the lagged moments beyond lags 0 and 1. It effectively minimizes the impact of potential measurement errors. More specifically, for model (ref) and (ref), due to the whiteness of ${{\boldsymbol{\zeta} }}_t$, $\hat{{\boldsymbol{\Gamma} }}_k \asymp {{\boldsymbol{\Gamma} }}_k + O_p(T^{-1/2})$ for both $k=1,2$, and the L2E will be asymptotically unbiased as long as ${{\boldsymbol{\Sigma} }}_{{\boldsymbol{\zeta} }}=o(1)$ (see Theorem (ref)).

The optimization problem in (ref) can be solved as follows. Let ${{\boldsymbol{Y} }} = \hat{{\boldsymbol{\Gamma} }}_2 \hat{{\boldsymbol{\Gamma} }}_1^{-1} (\hat{{\boldsymbol{\Gamma} }}_1 \hat{{\boldsymbol{\Gamma} }}_0^{-1} \hat{{\boldsymbol{\Gamma} }}_1^\top)^{1/2}$ and ${{\boldsymbol{V} }} = (\hat{{\boldsymbol{\Gamma} }}_1 \hat{{\boldsymbol{\Gamma} }}_0^{-1} \hat{{\boldsymbol{\Gamma} }}_1^\top)^{1/2}$. And let ${{\boldsymbol{y} }}_i$ and ${{\boldsymbol{v} }}_i$ be the $i$-th column of matrix ${{\boldsymbol{Y} }}$ and ${{\boldsymbol{V} }}$, then (ref) becomes equivalent to

Now the optimization problem in (ref) can be solved by the iterated least squares method introduced in Section (ref) similarly.

Improved factor estimation using Kalman Filter

The first step of the estimation produces estimated factor process $\hat{{{\boldsymbol{F} }}}_t$, using only the factor model component (ref). Since the second component (ref) also provides useful information about ${{\boldsymbol{F} }}_t$, it would be desirable to utilize the full model in the estimation. Given the estimated ${{\boldsymbol{U} }}_1$, ${{\boldsymbol{U} }}_2$, ${{\boldsymbol{A} }}_1$ and ${{\boldsymbol{A} }}_2$, DMFM becomes a linear state space model with the latent factor process ${{\boldsymbol{F} }}_t$ as the state and ${{\boldsymbol{X} }}_t$ as the observation. Hence Kalman filter and smoother kalman1960 can be used to provide an improved estimator of ${{\boldsymbol{F} }}_t$. However, using the full DMFM as the state-space model for filtering would involve very high dimensional estimated loading matrices and the estimated covariance matrices of $\hbox{\rm vec}({{\boldsymbol{E} }}_t)$. The estimation errors make the system unstable. Here we use a simplified approximation by assuming the estimated factor process from the first step is the true factor process with independent measurement error. Such an approximation allows us to work with a much simpler state-space model of low dimensions.

Specifically, we assume the following vector version of the state space model (ref) and (ref):

where ${{\boldsymbol{F} }}_{t}$ is the underlying true factor series, and ${{\boldsymbol{\Phi} }} = {{\boldsymbol{A} }}_2 \otimes {{\boldsymbol{A} }}_1$ is the coefficient matrix of the VAR model (ref). Recall that the innovations ${{\boldsymbol{\xi} }}_t$ and ${{\boldsymbol{\zeta} }}_t$ have covariance structures ${\rm Cov}(\hbox{\rm vec}({{\boldsymbol{\xi} }}_t)) = {{\boldsymbol{\Sigma} }}_{{\boldsymbol{\xi} }}$ and ${\rm Cov}(\hbox{\rm vec}({{\boldsymbol{\zeta} }}_t)) = \boldsymbol{\Sigma}_{{{\boldsymbol{\zeta} }}}$ respectively. Equation ((ref)) is the state equation with $\hbox{\rm vec}({{\boldsymbol{F} }}_t)$ as the state, and ((ref)) is the observation equation, with $\hbox{\rm vec}(\tilde{{{\boldsymbol{F} }}_t})$ as the observations.

To estimate the covariance matrices ${{\boldsymbol{\Sigma} }}_{{\boldsymbol{\xi} }}$ and ${{\boldsymbol{\Sigma} }}_{{\boldsymbol{\zeta} }}$ needed for filtering, we obtain the sample “residual” ${{\boldsymbol{W} }}_t =\hbox{\rm vec}(\tilde{{{\boldsymbol{F} }}}_t) - {{\boldsymbol{\Phi} }}\hbox{\rm vec}(\tilde{{{\boldsymbol{F} }}}_{t-1}) =\hbox{\rm vec}({{\boldsymbol{\zeta} }}_t) - {{\boldsymbol{\Phi} }}\hbox{\rm vec}({{\boldsymbol{\zeta} }}_{t-1}) +\hbox{\rm vec}({{\boldsymbol{\xi} }}_{t})$. Since \[ {{\boldsymbol{G} }}_0 := {\rm Cov}({{\boldsymbol{W} }}_t) = {{\boldsymbol{\Sigma} }}_{{\boldsymbol{\xi} }} + {{\boldsymbol{\Phi} }} \boldsymbol{\Sigma}_{{\boldsymbol{\zeta} }}{{\boldsymbol{\Phi} }}^{\top} + \boldsymbol{\Sigma}_\zeta, \mbox{\ and \ } {{\boldsymbol{G} }}_1 := {\rm Cov}({{\boldsymbol{W} }}_t,{{\boldsymbol{W} }}_{t-1})= -{{\boldsymbol{\Phi} }} \boldsymbol{\Sigma}_{{\boldsymbol{\zeta} }}, \] we have \[ \boldsymbol{\Sigma}_{{\boldsymbol{\zeta} }} = -{{\boldsymbol{\Phi} }}^{-1} {{\boldsymbol{G} }}_1, \mbox{\ and \ } {{\boldsymbol{\Sigma} }}_{{\boldsymbol{\xi} }} = {{\boldsymbol{G} }}_0 - \boldsymbol{\Sigma}_{{\boldsymbol{\zeta} }}- {{\boldsymbol{\Phi} }} \boldsymbol{\Sigma}_{{\boldsymbol{\zeta} }} {{\boldsymbol{\Phi} }}^{\top}. \] In practice, we let the estimated factors $\hat{{\boldsymbol{F} }}_t$ to play the roles of $\tilde{{\boldsymbol{F} }}_t$, and use the L2E $\hat{{\boldsymbol{\Phi} }}=\hat{{\boldsymbol{A} }}_2\otimes\hat{{\boldsymbol{A} }}_1$ to calculate the $\hat{{\boldsymbol{W} }}_t =\hbox{\rm vec}(\hat{{{\boldsymbol{F} }}}_t) - \hat{{\boldsymbol{\Phi} }}\hbox{\rm vec}(\hat{{{\boldsymbol{F} }}}_{t-1})$, from which the sample autocovariance matrices $\hat{{\boldsymbol{G} }}_0$ and $\hat{{\boldsymbol{G} }}_1$ are obtained. A direct moment estimator of ${{\boldsymbol{\Sigma} }}_{{\boldsymbol{\zeta} }}$ would be $-\hat{{\boldsymbol{\Phi} }}^{-1}\hat{{\boldsymbol{G} }}_1$ but it may not be symmetric and positive semi-definite as a covariance matrix needs to be. Hence we project $-\left(\hat{{\boldsymbol{\Phi} }}^{-1}\hat{{\boldsymbol{G} }}_1+\hat{{\boldsymbol{G} }}_1^\top\hat{{\boldsymbol{\Phi} }}^{-\top}\right)/2$ to the cone of positive semi-definite matrices. Specifically, let $-\left(\hat{{\boldsymbol{\Phi} }}^{-1}\hat{{\boldsymbol{G} }}_1+\hat{{\boldsymbol{G} }}_1^\top\hat{{\boldsymbol{\Phi} }}^{-\top}\right)/2=\hat{{\boldsymbol{Q} }}\hat{{\boldsymbol{D} }}\hat{{\boldsymbol{Q} }}'$ be the eigenvalue decomposition, we take

where $\hat{{\boldsymbol{D} }}_+$ is the diagonal matrix obtained by thresholding the diagonal entries of $\hat{{\boldsymbol{D} }}$ from below by zero. The direct moment estimator of ${{\boldsymbol{\Sigma} }}_{{\boldsymbol{\xi} }}$ is then given by $\hat{{\boldsymbol{\Sigma} }}_{{\boldsymbol{\xi} }} = \hat{{\boldsymbol{G} }}_0 - \hat{\boldsymbol{\Sigma}}_{{\boldsymbol{\zeta} }}- \hat{{\boldsymbol{\Phi} }} \hat{\boldsymbol{\Sigma}}_{{\boldsymbol{\zeta} }} \hat{{\boldsymbol{\Phi} }}^{\top}$, which again is projected to the cone of positive semi-definite matrices to generate the final estimator $\hat{{\boldsymbol{\Sigma} }}_{{\boldsymbol{\xi} }}$.

Once $\hat{{\boldsymbol{\Phi} }}$, $\hat{{\boldsymbol{\Sigma} }}_{{\boldsymbol{\xi} }}$ and $\hat{{\boldsymbol{\Sigma} }}_{{\boldsymbol{\zeta} }}$ are obtained, the implementation of the Kalman filter based on (ref) and (ref) is standard, so we omit the details here. We remark that the projection step may produce singular $\hat{{\boldsymbol{\Sigma} }}_{{\boldsymbol{\xi} }}$ and $\hat{{\boldsymbol{\Sigma} }}_{{\boldsymbol{\zeta} }}$, although they are both positive semi-definite. In the implementation, we always use the Moore-Penrose inverse if a singular covariance matrix is encountered at any step of the algorithm.

Prediction

Denote by $\mathcal{F}_t$ the set of observations $\{{{\boldsymbol{X} }}_1,\cdots, {{\boldsymbol{X} }}_t\}$ up to time $t$, and by ${{\boldsymbol{X} }}_t(1)$ and ${{\boldsymbol{F} }}_t(1)$ the best linear predictor of ${{\boldsymbol{X} }}_{t+1}$ and ${{\boldsymbol{F} }}_{t+1}$ based on ${\cal F}_t$. Under the DMFM model (ref) and (ref),

where ${{\boldsymbol{F} }}_t(0)$ is the best linear prediction of ${{\boldsymbol{F} }}_t$ based on ${\cal F}_t$. To generate the prediction based on the estimated model, a straightforward approach is to plug in the estimates $\hat{{\boldsymbol{U} }}_i$, $\hat{{\boldsymbol{A} }}_i$ and $\hat{{\boldsymbol{F} }}_t$ into (ref). However, the behavior of such a plug-in prediction is subtle, as we will elaborate below. First of all, as to be shown in both Section (ref) and Section (ref), when the SNR is relatively small, the L2E has better performance than LSE. However, even when $\hat{{\boldsymbol{A} }}_i^{\rm (L2E)}$ are closer to the true ${{\boldsymbol{A} }}_i$ than $\hat{{\boldsymbol{A} }}_i^{\rm (LSE)}$, a direct plug-in of $\hat{{\boldsymbol{A} }}_i^{\rm (L2E)}$ into (ref) will not lead to a better prediction, because the prediction also involves $\hat{{\boldsymbol{F} }}_t$ as the prediction of the latent ${{\boldsymbol{F} }}_t(0)$ by ${\cal F}_t$, but $\hat{{\boldsymbol{F} }}_t$ is obtained using only (ref) without using the dynamics in (ref). It can be shown easily that ${{\boldsymbol{A} }}_1^{(\rm LSE)}$ and ${{\boldsymbol{A} }}_2^{(\rm LSE)}$ attempt to minimize the square linear prediction error of ${{\boldsymbol{F} }}_{t+1}$ by $\hat{{\boldsymbol{F} }}_{t}$, and hence that of ${{\boldsymbol{X} }}_{t+1}$ by $\hat{{\boldsymbol{F} }}_t$, if ${{\boldsymbol{U} }}_1$ and ${{\boldsymbol{U} }}_2$ are known. Therefore, the plug-in prediction (ref) using LSE is preferred to using L2E, even when the SNR is small, though L2E is the preferred estimator for ${{\boldsymbol{A} }}_1$ and ${{\boldsymbol{A} }}_2$ in such cases. Therefore, if a plug-in type prediction is entertained, one should always use the LSE:

However, the proposed L2E can be effectively used to enhance the prediction performance using the underlying state-space structure of the process, instead of the simple plug-in prediction. Specifically, we propose to use the Kalman filter in Section (ref) to obtain the best linear estimation of ${{\boldsymbol{F} }}_t$ based on $\{\hat{{\boldsymbol{F} }}_1,\ldots,\hat{{\boldsymbol{F} }}_t\}$ as an estimate of ${{\boldsymbol{F} }}_t(0)$, and then obtain ${{\boldsymbol{X} }}_t(1)$ using (ref). In practice, the Kalman filter is applied to $\{\hat{{\boldsymbol{F} }}_t\}$ with the estimated $\hat{{\boldsymbol{A} }}_i^{\rm (L2E)}$, $\hat{{\boldsymbol{\Sigma} }}_{{\boldsymbol{\xi} }}$ and $\hat{{\boldsymbol{\Sigma} }}_{{\boldsymbol{\zeta} }}$ to get $\hat{{\boldsymbol{F} }}_t(0)$. The prediction of ${{\boldsymbol{X} }}_{t+1}$ is then given by

It is critical in practice to decide whether the Kalman filter based prediction (ref) or the plug-in prediction (ref) should be adopted, which in turn depends on the level of the SNR. There is some space for a more thorough and detailed analysis on the choice between LSE and L2E, using potentially certain information criterion or bootstrap procedure. Due to the scope of the current paper, we will leave it to the future work. Here we propose a simple rule of thumb based on heuristics. Note that the estimate $\hat{{\boldsymbol{\Sigma} }}_{{\boldsymbol{\zeta} }}$ in (ref) is the estimated covariance matrix of the measurement error in $\hat{{\boldsymbol{F} }}_t$. If the measurement error is significantly large, there is a large chance that $\hat{{\boldsymbol{\Sigma} }}_{{\boldsymbol{\zeta} }}$ is strictly positive definite. Therefore, we would use the Kalman filter based prediction (ref) if the smallest eigenvalue of $\hat{{\boldsymbol{\Sigma} }}_{{\boldsymbol{\zeta} }}$ is positive and use the plug-in prediction (ref) otherwise.

Theoretical properties

In this section we present some theoretical properties of the proposed two-stage estimation procedures. To single out the impact of the SNR for the estimation of both ${{\boldsymbol{U} }}_i$ and ${{\boldsymbol{A} }}_i$, we return in this section the original model formulation (ref) and (ref), where $\lambda$ is made explicit. We introduce some notations first. Let $d_{\max}=\max\{d_1,d_2\}$, $\otimes$ be the Kronecker product. The matrix Frobenius norm is defined as $\|{{\boldsymbol{A} }}\|_{\rm F} = (\sum_{ij} a_{ij}^2)^{1/2}$. Define the spectral norm as $$ \|{{\boldsymbol{A} }}\|_{\rm S} = \max_{\|{{\boldsymbol{x} }}\|_2=1,\|{{\boldsymbol{y} }}\|_2= 1} \|{{\boldsymbol{x} }}^{\top} {{\boldsymbol{A} }} {{\boldsymbol{y} }}\|_2.$$

To establish the consistency of the proposed procedures, we impose the following assumptions.

Assumption (ref) is adopted from chen2022factor and han2020 on tensor factor models, which also corresponds to the white noise assumption of lam2011,lam2012. Different from the assumptions in chen2023statistical, it allows substantial contemporaneous correlation among the entries of ${{\boldsymbol{E} }}_t$. Note that the normality assumption, which ensures fast convergence rates in our analysis, is imposed for technical convenience. In fact it can be extended to the sub-Gaussian condition, and still maintains the same convergence rates of the factor loading matrices as those presented in chen2022factor,han2020. Assumption (ref) is a standard condition for matrix autoregressive models.

han2020 proposed iterative procedures, iTOPUP and iTIPUP, to estimate the factor loading matrices of matrix/tensor factor models. When the ranks $r_k$ and time lag $h_0$ are fixed, both methods delivered convergence rate $\| \widehat{{\boldsymbol{U} }}_k\widehat{{\boldsymbol{U} }}_k^{\top} -{{\boldsymbol{U} }}_k {{\boldsymbol{U} }}_k^{\top} \|_{\rm S} =O_{\mathbb{P}}(\sigma d_{\max}^{1/2}\lambda^{-1}T^{-1/2}),k=1, 2$, which also matches the statistical lower bound. Using their iterative procedures, we can further establish the theoretical properties of the estimated latent factor process $\widehat {{\boldsymbol{F} }}_t$ in (ref).

The proposition specifies the convergence rate for the estimated factors $\widehat{{\boldsymbol{F} }}_t$ (up to a transformation). It shows that, in order to estimate the factors consistently, the signal to noise ratio of the factor model ($\lambda/\sigma$) must go to infinity, in order to have sufficient information on the signal part at each time point $t$. Concerning the impact of $\lambda/\sigma$, it is not surprising that the stronger the factor strength is, the more useful information the observed process carries and the faster the estimated factors converge. Moreover, (ref) shows that the error rates for the sample auto-covariance of the estimated factors to the true sample auto-cross-moment is $o_{\mathbb{P}}(T^{-1/2})$ when $\lambda/\sigma \gg \sqrt{d_{\max}}$. When $\lambda/\sigma\gg \sqrt{d_{\max}}+\sqrt{T}$, (ref) is much smaller than the parametric rate $T^{-1/2}$. Hence, intuitively, this implies that it is a valid option to use the estimated factor processes as the true factor processes to model the dynamics of the factor under such situations. The results are expected to be the same as using the true factor process, without loss of efficiency. The statistical rates in (ref) and (ref) lay a foundation for further modeling of the estimated factor processes with vast repository of linear and nonlinear options.

The following theorem shows the rate of convergence for alternating least square estimators of the second stage MAR modelling of the latent factor process, directly using $\hat{{{\boldsymbol{F} }}}_t$ as observed, outlined in Section (ref). Due to the identifiability issue, we make the convention that $\|{{\boldsymbol{A} }}_1\|_{\rm F}=1, \|\widehat {{\boldsymbol{A} }}_1\|_{\rm F}=1.$

The proof of the theorem is given in Appendix (ref).

The first two terms of (ref) come from the estimation error of the latent factor process and the covariance matrices. Compared with (ref), the errors are more precisely controlled in the second stage estimation of the MAR component.

In addition to the consistency of the MAR estimators, we also have the following asymptotic normality under a stronger signal strength condition.

The proof of the theorem is given in Appendix (ref).

Furthermore, we have the following convergence rates and central limit theorems for the Lag-2 estimators, $\widehat {{\boldsymbol{A} }}_1^{(\rm L2E)}$ and $\widehat {{\boldsymbol{A} }}_2^{(\rm L2E)}$ in Section (ref).

Numerical studies

Simulations settings

To evaluate the performance of the proposed estimation and prediction procedures, we use the original model formulation (ref) and (ref) as the data generating mechanism. For each $(d_i,r_i)$, the loading matrices ${{\boldsymbol{U} }}_i$ is randomly generated as the first $r_i$ left singular vectors of a standard $d_i\times d_i$ Gaussian ensemble (with i.i.d. $N(0,1)$ entries). The matrix ${{\boldsymbol{A} }}_i$ is generated as ${{\boldsymbol{A} }}_i={{\boldsymbol{L} }}_i{{\boldsymbol{D} }}_i{{\boldsymbol{R} }}_i^\top$, where ${{\boldsymbol{L} }}_i$ and ${{\boldsymbol{R} }}_i$ are $r_i\times r_i$ orthogonal matrices generated in the same way as ${{\boldsymbol{U} }}_i$, and ${{\boldsymbol{D} }}_i$ is a diagonal matrix whose diagonal entries are randomly generated from $\rm{Unif}[0.5,1.5]$. The matrices ${{\boldsymbol{A} }}_1$ and ${{\boldsymbol{A} }}_2$ are then rescaled so that $\|{{\boldsymbol{A} }}_1\|_{\rm F} = 1$ and the spectral radius of ${{\boldsymbol{A} }}_2\otimes{{\boldsymbol{A} }}_1$ is $\rho$. Throughout the simulations, ${{\boldsymbol{E} }}_t$ is a $d_1 \times d_2$ standard Gaussian ensemble. The innovations ${{\boldsymbol{\xi} }}_t$ in (ref) are generated as i.i.d. Gaussian with covariance matrix ${{\boldsymbol{\Sigma} }}_{{\boldsymbol{\xi} }}={\rm Cov}(\hbox{\rm vec}({{\boldsymbol{\xi} }}_t))$. The ${{\boldsymbol{\Sigma} }}_{{\boldsymbol{\xi} }}$ is generated through its eigenvalue decomposition, similar to the generation of ${{\boldsymbol{A} }}_i$, but its $r_1r_2$ eigenvalues are equally spaced over $[1,10]$. For the original DMFM formulation (ref) and (ref), the SNR is defined as

We will adjust the value of $\lambda$ to reach a specified level of SNR. In the simulation studies, we fix the factor dimensions to be $r_1=r_2=3$ and consider various configurations of $\{d_1, d_2, T, \rho, \lambda\}$. For each configuration, all the parameters related to the data generating mechanism, including ${{\boldsymbol{U} }}_i$, ${{\boldsymbol{A} }}_i$, $\lambda$, $\rho$ and ${{\boldsymbol{\Sigma} }}_\xi$ are held fixed once they are generated. We will then simulate the DMFM model with 100 repetitions to inspect the empirical performance of the estimation in Section (ref) and of the prediction in Section (ref).

Estimation performance

The estimation of ${{\boldsymbol{U} }}_i$ is adopted from chen2022factor and han2020, and is not the focus of this paper, so we only examine the the estimation of ${{\boldsymbol{A} }}_i$ in this section.

To introduce the metric for the estimation performance, we note that the estimation of the loading ${{\boldsymbol{U} }}_i$ is actually estimating the column space ${{\boldsymbol{U} }}_i$. In other words, the $\hat{{\boldsymbol{U} }}_i$ is trying to estimate ${{\boldsymbol{U} }}_i$ after a orthogonal rotation. Although further identification conditions can be introduced to fully identify ${{\boldsymbol{U} }}_1$ and ${{\boldsymbol{U} }}_2$ up to a sign change, such an identification will only produce an equivalent estimation of ${{\boldsymbol{A} }}_i$ and does not affect the prediction of ${{\boldsymbol{X} }}_t$. Therefore, we shall skip the unnecessary additional identification. On the other hand, when ${{\boldsymbol{F} }}_t$ are estimated as $\hat{{\boldsymbol{U} }}_1^\top{{\boldsymbol{X} }}_t\hat{{\boldsymbol{U} }}_2$, the MAR model on ${{\boldsymbol{F} }}_t$ can be written as

Therefore, if an MAR model is fitted to $\hat{{\boldsymbol{F} }}_t$, we are essentially estimating

In the simulation we will use the squared error (on log scale) to evaluate the performance of different estimators:

where $\hat{{\boldsymbol{A} }}_i$ is either the LSE or the L2E.

In Figure (ref), we show the box plots of the errors over 100 repetitions (ref) for the case $d_1=d_2=8, r_1=r_2=3, \rho=.9$ and $T=1000$. The box plot is given against the SNR ranging from 0.6 to 128, as labelled at the horizontal axis. For each SNR, two box plots are shown, the one on the left for LSE, and the one on the right for L2E. The two horizontal lines show the median errors (corresponding to the two box plots at the very right side) for the oracle case where the true ${{\boldsymbol{F} }}_t$ are used to estimate ${{\boldsymbol{A} }}_i$. This is equivalent to an infinity SNR so there is no estimation error in $\hat{{\boldsymbol{F} }}_t$. We label it as $\infty$ in the plot.

There are a number of observations that can be made from Figure (ref). First, we see the general pattern that when the SNR is low, L2E performs better than the LSE, and when the SNR reaches 8, LSE starts to outperform L2E. This is consistent with the discussion made in Section (ref) and the Proposition (ref). Second, L2E can be significantly better than LSE, especially when the SNR is between 0.8 and 5. We note that the errors are plotted in the log scale, so a unit difference in Figure (ref) means the squared error is actually reduced by a factor of $1/e$. Last but not least, we see that as the SNR increases, eventually both the LSE and the L2E are getting closer to the oracle estimators when there is no estimation error. It is clear that L2E is less efficient in the oracle case, as Proposition (ref) indicates. We would like to further point out that the L2E starts to behave almost same as its oracle version when the SNR reaches 8, while the LSE requires a SNR as large as 16 to behave like the oracle LSE. This again confirms the fact the L2E is less susceptible to the estimation errors in $\hat{{\boldsymbol{F} }}_t$.

Similar plots for different simulation settings are give in the appendix.

Prediction Performance

In this section, we compare the performance of five prediction methods based on the DMFM. Throughout this section we consider one-step ahead prediction with origin $T$.

We compare the methods listed above through the prediction mean squared error

Note that in the simulation for different SNR, we fix the generating mechanism of ${{\boldsymbol{E} }}_t$, and change the value of $\lambda$ in (ref) in order to reach different levels of SNR. Consequently, there is a factor $(\hbox{SNR})^{-2}$ in the definition of PSE in (ref) so that the PSE at different SNR levels are comparable.

In Figure (ref) we compare the first 4 methods using the box plots of $\log(PSE)$ over 100 repetitions. For each SNR, 4 box plots are for LSE, L2E, V.LSE, V.L2E from left to right. We see that among LSE, V.LSE and V.L2E, the LSE prediction is consistently better than or comparable to others. Comparing LSE and L2E, we see L2E enjoys a better performance when the SNR is between 0.8 and 4. When the SNR is larger, the LSE leads to the best predictions. We also compare the LSE and L2E+ in Figure (ref), by plotting the ratios of the L2E+ PSE over the LSE PSE. These ratios are in general more likely to be below one, confirming the advantage of the recommended L2E+ prediction method, which adaptively combines the strength of LSE and L2E. When the SNR is as large as 5 or 6, L2E+ starts to coincide with LSE.

Real example: New York City Yellow Cab Taxi Trips

In this section, we apply the dynamic matrix factor model to analyze the New York City yellow cab taxi traffic data. This data is maintained by the Taxi & Limousine Commission of New York City and can be downloaded from https://www1.nyc.gov/site/tlc/about/tlc-trip-record-data.page. The portion of the dataset we use comprises more than 1 billion trip records spanning from January 1, 2009, to December 31, 2019. Most of the trips occur within Manhattan Island. Each trip record includes the time and location of pickup and drop-off, as well as payment information.

Similar to that in chen2022factor, we obtain a $19\times 19\times 750$ matrix time series, with 19 of 69 predefined zones in Manhattan (see the blue area in the map in Figure (ref) for more details) and $T=750$ business days in years 2015-2017. The zones selected are mainly the midtown areas of Manhattan, which include Penn Station and numerous busy business districts. Each entry $X_{ijt}$ is the aggregated number of ride from zone $i$ to zone $j$ in the hours between 7am to 10am (morning rush hours) on day $t$.

Figure (ref) plots the daily sum of trips across all 19 selected regions in Manhattan during 7am to 10am, revealing a decreasing trend starting from 2014. This decline can be attributed to the emergence of ride-hailing services. Additionally, seasonal trends within each year are visible, likely due to temperature fluctuations. These two factors contribute to the non-stationarity of the time series, necessitating the removal of trends using exponential smoothing and the subsequent building of models on the residual data. Specifically, using smoothing parameter $\alpha=0.1$, we obtain the exponential smoothing trend $\hat{{{\boldsymbol{B} }}}_{t} = \alpha {{\boldsymbol{Y} }}_t + (1-\alpha) \hat{{{\boldsymbol{B} }}}_{t-1}$ where ${{\boldsymbol{Y} }}_t$ is observed matrix time series. Let ${{\boldsymbol{X} }}_t={{\boldsymbol{Y} }}_t-{{\boldsymbol{B} }}_t$ as the detrended series. Figure (ref) shows the detrended series among 11 zones. The $i$-th rows shows the pick-up volume from $i$-th zone to one of the other zones, and the $j$-th column show the drop-off volume in $j$-th zone coming from one of the other zones. Each individual time series is mostly stationary after removing the trends. We can see that some zones, such as Zone 1 (the first zone from the 19 selected) have large pickup volumes, but very small drop-off volumes. These zones are mostly residential areas, with people going to work during morning rush hours. Zones such as 4 and 8, however, have large drop-off volumes but small pickup volumes, are commercial districts.

Dynamic matrix factor model and several other models are estimated and their out-sample rolling forecasting performance are comparied. The rolling period is the entire year 2017 period. For each day in year 2017, we use all the past observations to fit the model, and make a prediction of the current day and compare with the true trip counts on that day. We report \[ \mbox{RMSE} = \sqrt{(250\times d_1\times d_2)^{-1}\sum\limits_{t=499}^{749}\|\hat{{{\boldsymbol{X} }}}_{t}(1) - {{\boldsymbol{X} }}_{t+1}\|_{\rm F}^2}. \]

Here is a list of models we compared.

(a) Matrix factor model methods: The first step is to estimate a matrix factor model with factor rank $(5,5)$ and obtain the estimated factor process $\hat{{{\boldsymbol{F} }}}_t$. In the second step, various models are used to fit $\hat{{{\boldsymbol{F} }}}_t$ for prediction, including a MAR model with LSE (MFMLSE), a VAR model (MFV), individual independent AR model for each factor element (MFi), a random walk model for the factors (MFrw), and a constant mean model for the factors using the estimated mean (MFmean).

(b) Vector factor model methods: The first step is to estimate a vector factor model with factor rank $25$ using the stacked observations $\hbox{\rm vec}({{\boldsymbol{X} }}_t)$ and obtain the estimated factor process $\hat{\boldsymbol{f}}_t$. In the second step, various models are used to fit $\hat{\boldsymbol{f}}_t$ for prediction, including a VAR model (VFV), individual independent AR model for each factor element (VFi), a random walk model for the factors (VFrw), and a constant model for the factors using the estimated mean (VFmean).

(c) Autoregressive model on the observed ${{\boldsymbol{X} }}_t$: These models do not apply a factor model, but instead directly apply MAR(MAR) of chen2021autoregressive, reduced rank MAR (MARR) of han2023rr, VAR using the stacked $\hbox{\rm vec}({{\boldsymbol{X} }}_t)$ (VAR), individual AR (iAR), random walk (Xrw), and constant time series using the mean estimator (Xmean).

In Table (ref), we can see that the dynamic matrix factor model has the best prediction performance among all models. The individual AR models (MFi, VFi and iAR) are surprisingly better than some of the more complicated models. VAR performs the worst, as it tends to overfit with a $19^2\times 19^2$ coefficient matrix. This data set appears to have a large SNR, so the plug-in prediction (ref) based on the DMF performs the best. We also remark that the proposed rule at the end of Section (ref) always choose to use the plug-in prediction over the Kalman-filter-based prediction (ref).

Conclusion

In this paper, we present a dynamic matrix factor model as an extension of the matrix factor model. The model incorporates the dynamics of latent factors by using a matrix autoregressive model, making it more suitable for analyzing the dynamics of low-dimensional latent factors and making predictions. The dynamic matrix factor model is similar in form to the reduced rank matrix autoregressive model, with the main difference being the lower rank error. The model is estimated in two steps, starting with estimating the factor model and then incorporating the matrix autoregressive dynamics. To improve the estimation accuracy, we propose a lag-2 Yule-Walker estimator and use Kalman Filter to reduce the influence of errors from matrix factor estimation when the signal-to-noise ratio is low. This allows us to obtain more accurate estimates of the underlying factors and their dynamics, even in the presence of noise and other sources of measurement error. The prediction is done through the prediction of the latent factors. Our numerical analysis demonstrates that the dynamic matrix factor model performs well in terms of estimation and prediction performance by capturing the underlying dynamics of the latent factors. We also apply the dynamic matrix factor model to NYC yellow cab taxi data to demonstrate its advantages. Our results suggest that the dynamic matrix factor model is a powerful tool for analyzing the dynamics of matrix-valued time series data and making accurate predictions over time.