Factor Models of Matrix-Valued Time Series: Nonstationarity and Cointegration

Factor Models of Matrix-Valued Time Series: Nonstationarity and Cointegration

\centerline{\bf Abstract}

In this paper, we consider the nonstationary matrix-valued time series with common stochastic trends. Unlike the traditional factor analysis which flattens matrix observations into vectors, we adopt a matrix factor model in order to fully explore the intrinsic matrix structure in the data, allowing interaction between the row and column stochastic trends, and subsequently improving the estimation convergence. It also reduces the computation complexity in estimation. The main estimation methodology is built on the eigenanalysis of sample row and column covariance matrices when the nonstationary matrix factors are of full rank and the idiosyncratic components are temporally stationary, and is further extended to tackle a more flexible setting when the matrix factors are cointegrated and the idiosyncratic components may be nonstationary. Under some mild conditions which allow the existence of weak factors, we derive the convergence theory for the estimated factor loading matrices and nonstationary factor matrices. In particular, the developed methodology and theory are applicable to the general case of heterogeneous strengths over weak factors. An easy-to-implement ratio criterion is adopted to consistently estimate the size of latent factor matrix. Both simulation and empirical studies are conducted to examine the numerical performance of the developed model and methodology in finite samples.

{\em Keywords}: common stochastic trends, eigenanalysis, matrix error-correction models, matrix factor models, ratio criterion.

Introduction

\setcounter{equation}{0}

There has been increasing attention in modeling large-scale matrix-valued data over the past two decades with applications in various fields such as environmental science, economics, finance, public health and social networks. Typical examples include the international trade flow among a large number of countries CC23, the Netflix movie-rating challenge with customers and movies being matrix rows and columns, respectively HTW15, and friendships between pupils in schools CFIK20. A traditional method for tackling matrix-valued data is to convert them into multivariate data, i.e., matrices are flattened into vectors, and then adopt statistical tools which have been well studied for multivariate data (with possible high dimension). However, this naive method neglects the intrinsic matrix nature of the data and may lose important information contained in the matrix structure, leading to sup-optimal estimation convergence theory CF23, HKTY23.

There have been extensive studies on estimation and inference of independent matrix-valued data, see, for example, NW11, RT11, LT12, HTW15 and the references therein. In practice, matrix-valued data are often collected over time with significant serial correlation which may be useful to the subsequent causality and forecasting analysis. There have been notable developments in models, methodology and theory for matrix-variate time series. CXY21 study the estimation of autoregressive models for matrix time series via a bilinear form; CTC19, WLC19, CHYY23, CF23, HKTY23 and HYZC24 propose various matrix factor models based on eigenanalysis using either contemporaneous covariance or dynamic auto-covariance matrices; and HCZY24 propose simultaneous decorrelation of matrix time series to achieve an effective dimension reduction.

Most of the aforementioned literature requires a stationarity assumption on matrix time series, facilitating theoretical derivation of standard estimation and inferential theory. However, this assumption is too restrictive and often rejected when we test some empirical matrix time series collected over a long time span. A transformation mechanism may be required to convert nonstationary matrix time series to stationary ones. For example, a first- or second-order difference is taken to transform the multinational macroeconomic time series to the stationary matrix-variate process CXY21, CF23, HKTY23. In practice, the nonstationarity may be due to the existence of common stochastic trends. Conventional transformation mechanisms such as differencing can result in loss of important trending information.

The present paper introduces a general matrix factor model framework for nonstationary matrix time series, where the common stochastic trending patterns of matrix observations are captured by the latent common component through a low-rank structure. In particular, we study the following two scenarios: (i) the nonstationary matrix factors are of full rank \footnote{Here, by “full rank” we mean that the covariance matrix of the (normalized) row or column factor process is of full rank, which indicates the absence of cointegration among the factors.} and the idiosyncratic components are temporally stationary, and (ii) the matrix factors are cointegrated and the idiosyncratic components may be nonstationary. The nonstationary matrix factors in scenario (i) capture the common stochastic trends whereas the cointegrated matrix factors in scenario (ii) are generated by the matrix error-correction model which extends the classic vector error-correction model to matrix time series. The developed matrix factor model is a natural extension of the conventional approximate factor model which has been extensively studied in the literature for stationary and nonstationary vector time series CR83, BN02, B04, BN04, FLM13, BLL21. The advantage of the matrix factor model over the vector factor model (by flattening matrix observations into vector ones) is that the intrinsic matrix structure can be fully explored and the interaction between the row and column stochastic trends is allowed, which would lead to improvement of the estimation convergence for unknown elements in the matrix factor model. It also reduces the computation complexity in estimation (see the last paragraph in Section (ref) below).

For scenario (i) with full-rank I(1) (integrated of order one) matrix factors and stationary idiosyncratic components, the main estimation technique is built on the principal component analysis (PCA) using either a sample row or a sample column covariance matrix, which is called mPCA throughout the paper. This estimation methodology is also valid when the nonnegative definite matrix is constructed via aggregation of auto-covariance matrices over lags LYB11, LY12, CTC19, WLC19, in which case the idiosyncratic components are often assumed to be temporally uncorrelated. In the present paper, due to stronger signals from the stochastic trends in common components, we relax this restriction, allowing the idiosyncratic components in the matrix factor model to be serially correlated. The convergence theory in Section (ref) reveals that the developed estimators achieve the so-called “super-fast" convergence rates as in the recent work by CGGT25. The remarkable rate improvement is due to the following two reasons. First, the common components in the proposed matrix factor model contain latent stochastic trends, resulting in much stronger signals in the estimation equations than those in the stationary setting. This significantly improve the estimation convergence for the factor loading matrices compared with that for stationary matrix time series WLC19, CF23. Second, the developed estimation methodology makes full use of the intrinsic matrix structure and thus improves the convergence rates over the traditional method via the vectorized factor model B04, BLL21.

For scanario (ii) with cointegrated matrix-valued factors and mixed stationary and nonstationary idiosyncratic components, we extend BN04's PANIC (Panel Analysis of Nonstationarity in Idiosyncratic and Common components) technique from vector time series to matrix time series and call the proposed method as mPANIC. Through the eigenanalysis of the differenced matrix time series, we obtain the row and column factor loading estimation which has typical stationary convergence rates as shown in Theorem (ref). Furthermore, assuming a “sparsity" restriction which limits the number of nonstationary entries in the matrix idiosyncratic components, we derive the uniform convergence of the estimated cointegrated factor matrix.

Unlike the existing literature on matrix factor models with strong factors or loadings CF23, HKTY23, CGGT25, we not only consider weak factors in the developed matrix factor model but also allow heterogeneous factor strength. The existence of weak factors requires an adaptive estimation formula with random normalizers for estimating the common factor matrix. As in BN23 for vectorized factor model estimation theory, the weak and varying factor strength further slows down both the mean squared and the uniform convergence rates to be developed in Sections (ref) and (ref) and makes it more challenging to correctly estimate the row and column factor numbers. Our asymptotic theory reveals that the proposed eigenvalue ratio criterion can consistently estimate the factor dimension without any restriction for the minimum factor strength in scenario (i), but requires the weakest factor strength to exceed a fixed threshold as in F22 to achieve consistency in scanario (ii).

The finite-sample Monte-Carlo simulation study shows that both the mPCA and mPANIC estimates converge as the sample size increases, the mPCA (which directly tackles the nonstationary matrix time series) has faster convergence than the mPANIC, and the estimation errors increase as the factor strength decreases. The empirical application to international trade flows reveals that the estimated factor is connected to the so-called global demand index and the patterns of the estimated trending factor may be related to some major global events. In addition, the estimated row and column factor loading matrices can be used to capture the structure of bilateral trade linkages across some major economies.

The rest of the paper is organized as follows. Section (ref) introduces the matrix factor model with discussion on model identifiability, and the mPCA and mPANIC estimation methods. Section (ref) presents the convergence properties of the developed estimators with assumptions and remarks. They are extended to the case of heterogeneous factor strength in Section (ref). Section (ref) reports both the simulation and empirical studies. Section (ref) concludes the paper. The supplemental document contains the proofs of the main asymptotic results and some technical lemmas with proofs. Throughout the paper, we let ${\sf vec}(\cdot)$ be the vectorization of a matrix and ${\sf vec}^{-1}(\cdot)$ inverse of the vectorization operator. We use $\otimes$ and $\odot$ to denote the Kronecker and Hadamard products, respectively. Let $\Vert\cdot\Vert$ be the Euclidean norm of a vector or operator norm of a matrix, and $\Vert\cdot\Vert_F$ the matrix Frobenius norm. Let $\psi_k(\mathbf{A})$ denote the $k$-th largest eigenvalue of a positive semidefinite matrix $\mathbf{A}$. For notational simplicity, we write “$\max\{a,b\}$" as “$a\vee b$", “$\min\{a,b\}$" as “$a\wedge b$", and “with probability approaching one" as “{\em w.p.a.1}".

Model and methodology

\setcounter{equation}{0}

In this section, we introduce the factor model for nonstationary matrix time series and then propose two types of PCA-based estimation techniques depending on whether the nonstationary matrix factors are of full rank or cointegrated, and whether the idiosyncratic components are stationary.

Model formulation

Suppose that we observe a sequence of matrix time series ${\boldsymbol X}_t$, $t=1,\ldots,T$, where \[ {\boldsymbol X}_t=\left(

\right)_{p_1\times p_2}, \] where both the dimensions $p_1$ and $p_2$ diverge to infinity together with $T$. Consider the following matrix factor model formulation:

where ${\boldsymbol R}$ and ${\boldsymbol C}$ are row and column factor loading matrices with sizes $p_1\times r_1$ and $p_2\times r_2$, respectively, ${\boldsymbol F}_t\in{\mathscr R}^{r_1\times r_2}$ is a matrix of nonstationary factors, and ${\boldsymbol E}_t\in{\mathscr R}^{p_1\times p_2}$ is a matrix of idiosyncratic components. Throughout the paper, we assume that both $r_1$ and $r_2$ are fixed.

It is well known that the latent factors ${\boldsymbol F}_t$ and loading matrices ${\boldsymbol R}$ and ${\boldsymbol C}$ are not identifiable, i.e., for any two invertible matrices ${\boldsymbol V}_1$ and ${\boldsymbol V}_2$ with sizes $r_1\times r_1$ and $r_2\times r_2$, respectively, model ((ref)) still holds when the triplet $({\boldsymbol R}, {\boldsymbol F}_t, {\boldsymbol C})$ are replaced by $({\boldsymbol R}{\boldsymbol V}_1, {\boldsymbol V}_1^{-1}{\boldsymbol F}_t{\boldsymbol V}_2^{-1}, {\boldsymbol C}{\boldsymbol V}_2^{^\intercal})$. However, the column spaces spanned by ${\boldsymbol R}$ and ${\boldsymbol C}$, ${\cal M}({\boldsymbol R})$ and ${\cal M}({\boldsymbol C})$, are uniquely determined. In fact, assuming that ${\boldsymbol R}$ and ${\boldsymbol C}$ are of full column rank, we have the following QR decomposition: ${\boldsymbol R}={\boldsymbol Q}_R{\boldsymbol W}_R$ and ${\boldsymbol C}={\boldsymbol Q}_C{\boldsymbol W}_C$, where ${\boldsymbol Q}_R$ and ${\boldsymbol Q}_C$ are $p_1\times r_1$ and $p_2\times r_2$ matrices with orthogonal columns, and ${\boldsymbol W}_R$ and ${\boldsymbol W}_C$ are $r_1\times r_1$ and $r_2\times r_2$ upper triangular matrices. Writing ${\boldsymbol W}_t={\boldsymbol W}_R{\boldsymbol F}_t{\boldsymbol W}_C^{^\intercal}$, a transformed factor matrix of integrated processes, by ((ref)), we readily have that \[ {\boldsymbol X}_t={\boldsymbol Q}_R{\boldsymbol W}_t{\boldsymbol Q}_C^{^\intercal}+{\boldsymbol E}_t. \] It is clear that ${\cal M}({\boldsymbol R})$ is the same as ${\cal M}({\boldsymbol Q}_R)$ and ${\cal M}({\boldsymbol C})$ is the same as ${\cal M}({\boldsymbol Q}_C)$.

On the other hand, by flattening matrix observations to vector ones, we can rewrite the matrix factor model ((ref)) as the following vectorized factor model: \[ {\sf vec}({\boldsymbol X}_t)=\left({\boldsymbol C}\otimes{\boldsymbol R}\right){\sf vec}({\boldsymbol F}_t)+{\sf vec}({\boldsymbol E}_t). \] The conventional factor model has been extensively studied in the literature for vector time series CR83, BN02, FLM13. In particular, model ((ref)) can be seen as a natural extension of the classic approximation factor model for vector nonstationary time series B04,BN04,BLL21, fully exploring the intrinsic matrix structure and allowing interaction between the row and column stochastic trends in ${\boldsymbol F}_t$. It also extends factor models for stationary matrix time series WLC19, CXY21, CF23, CHYY23, HKTY23 to the nonstationary setting.

We are interested in the following two scenarios for latent factors and idiosyncratic components: (i) the nonstationary matrix factors are of full rank and the idiosyncratic components are stationary over time, and (ii) the matrix factors are cointegrated and the idiosyncratic components are either stationary or nonstationary. The nonstationarity of matrix time series in scenario (i) is fully captured by the common stochastic trends contained in matrix factors, whereas appropriate rotation of the cointegrated factors in scenario (ii) can separate out stationary and nonstationary factors. In the following two subsections, we provide detailed model settings and propose two types of PCA-based estimation techniques for these two scenarios, respectively.

mPCA estimation

Assume that ${\boldsymbol F}_t$ is generated by an integrated I(1) matrix-valued process:

where ${\boldsymbol U}_t\in{\mathscr R}^{r_1\times r_2}$ is a stationary and I(0) (integrated of order zero) matrix-valued process. In addition, the matrix idiosyncratic component ${\boldsymbol E}_t\in{\mathscr R}^{p_1\times p_2}$ is stationary and I(0). Under this model assumption, the matrix common component ${\boldsymbol Z}_t$ captures the stochastic trending patterns of matrix time series via the latent low-rank structure in ((ref)). Throughout this subsection, we assume that ${\boldsymbol F}_t$ is a matrix of full-rank I(1) random elements. Our main interest lies in estimation of the factor loading matrices ${\boldsymbol R}$ and ${\boldsymbol C}$ and construction of the matrix ${\boldsymbol F}_t$ for common stochastic trends. For this, we next introduce the mPCA approach via eigenanalysis of a nonnegative definite matrix which is either a sample row or column covariance matrix.

Assume $r_1$ and $r_2$ are known for the timebeing. Their consistent estimation can be constructed via a ratio criterion, see ((ref)) and ((ref)) below. Without loss of generality, assume that ${\boldsymbol X}_t$ has zero mean. We start with the estimation method for the row factor loading matrix ${\boldsymbol R}$ and define the sample row covariance matrix with size $p_1\times p_1$:

With the eigenanalysis of $\widehat{\boldsymbol\Omega}_{R}$, we obtain $\widehat{\boldsymbol R}=(\widehat{R}_{1},\ldots,\widehat{R}_{r_1})$ with $\widehat{R}_{k}$ being the eigenvector of $\widehat{\boldsymbol\Omega}_{R}$ corresponding to the $k$-th largest eigenvalue. The column space ${\cal M}({\boldsymbol R})$ is estimated by ${\cal M}(\widehat{\boldsymbol R})$. Equivalently $\widehat{\boldsymbol R}$ is the estimate of ${\boldsymbol R}$ subject to appropriate normalization and rotation (with the rotation matrix defined in Section (ref)). Defining the sample column covariance matrix: \[ \widehat{\boldsymbol\Omega}_{C}=\frac{1}{T}\sum_{t=1}^{T}{\boldsymbol X}_{t}^{^\intercal}{\boldsymbol X}_{t}, \] we may estimate the column factor loading matrix ${\boldsymbol C}$ in exactly the same way by conducting the eigenanalysis of $\widehat{\boldsymbol\Omega}_{C}$, and denote the resulting estimate by $\widehat{\boldsymbol C}$ which consists of the eigenvectors of $\widehat{\boldsymbol\Omega}_{C}$ corresponding to the $r_2$ largest eigenvalues. With the estimates $\widehat{\boldsymbol R}$ and $\widehat{\boldsymbol C}$, the matrix ${\boldsymbol F}_t$ of common stochastic trends can be estimated by

where $\widetilde\lambda_{R,1}$ is the maximum eigenvalue of $T^{-1}\widehat{\boldsymbol\Omega}_{R}$, $\widetilde{\boldsymbol V}_R$ is an $r_1\times r_1$ diagonal matrix with the diagonal entries being the $r_1$ largest eigenvalues of $T^{-1}\widehat{\boldsymbol\Omega}_R$ (in a descending order) whereas $\widetilde{\boldsymbol V}_C$ is an $r_2\times r_2$ diagonal matrix with the diagonal entries being the $r_2$ largest eigenvalues of $T^{-1}\widehat{\boldsymbol\Omega}_C$. The adaptive estimation formula in ((ref)) is due to possible existence of weak factors with the extra adaptive factors $\widetilde\lambda_{R,1}^{1/2}$, $\widetilde{\boldsymbol V}_R^{-1/2}$ and $\widetilde{\boldsymbol V}_C^{-1/2}$ related to the (varying) factor strength, see Assumptions (ref)(ii) and (ref)(i). When both the row and column factors are assumed to be strong, the estimation formula in ((ref)) can be simplified to $\widehat{\boldsymbol F}_t= (p_1p_2)^{-1/2} \widehat{\boldsymbol R}^{^\intercal} {\boldsymbol X}_t\widehat{\boldsymbol C}$.

The numbers of row and column factors, $r_1$ and $r_2$, are often unknown in practice. We next estimate them using the ratio criterion proposed by LY12 and AH13. Let $\{\widehat{\lambda}_{R,i}:\ i=1,\ldots,p_1\}$ and $\{\widehat{\lambda}_{C,i}:\ i=1,\ldots,p_2\}$ be sets of eigenvalues (arranged in a descending order) of $\widehat{\boldsymbol\Omega}_{R}$ and $\widehat{\boldsymbol\Omega}_{C}$, respectively. We estimate $r_1$ and $r_2$ by

and

respectively, where $K_1$ and $K_2$ are user-specified upper bounds.

In the construction of the mPCA above, we only need to conduct the eigenanalysis for matrices of sizes $p_i \times p_i$ ($i=1, 2$). This is in contrast to the PCA for flatterned ${\boldsymbol X}_t$ which involves an eigenanalysis for a matrix of size $(p_1 p_2)\times (p_1 p_2)$. For moderately large $p_1$ and/or $p_2$, the computational complexity is substantially reduced.

mPANIC estimation

We next remove the full-rank assumption on ${\boldsymbol F}_t$ and suppose that ${\boldsymbol F}_t$ is cointegrated and generated by the following error-correction model:

where $\Delta{\boldsymbol F}_t={\boldsymbol F}_t-{\boldsymbol F}_{t-1}$, ${\boldsymbol A}_1={\boldsymbol\alpha}_1{\boldsymbol\beta}_1^{^\intercal}$, ${\boldsymbol A}_2={\boldsymbol\alpha}_2{\boldsymbol\beta}_2^{^\intercal}$, ${\boldsymbol\alpha}_j$ and ${\boldsymbol\beta}_j$ are full-rank $r_j\times k_j$ matrices, $1\leq k_j\leq r_j$, $j=1,2$, and $\{{\boldsymbol V}_t\}$ is a sequence of $r_1\times r_2$ white noise matrices. Model ((ref)) is a natural extension of the vector error-correction (VEC) model to matrix time series. LX24 and HRW25 consider estimating the matrix error-correction model with extra lag and constant terms for observed matrix time series. In contrast, the factor process $\{{\boldsymbol F}_t\}$ is latent in the present paper. Flattening the factor matrix to a vector, we may rewrite ((ref)) as the following classic VEC model formulation: \[ {\sf vec}(\Delta {\boldsymbol F}_t)=({\boldsymbol A}_2\otimes {\boldsymbol A}_1){\sf vec}({\boldsymbol F}_{t-1})+{\sf vec}({\boldsymbol V}_t) \] with \[ ({\boldsymbol A}_2\otimes {\boldsymbol A}_1)=\left({\boldsymbol\alpha}_2{\boldsymbol\beta}_2^{^\intercal}\right)\otimes \left({\boldsymbol\alpha}_1{\boldsymbol\beta}_1^{^\intercal}\right)=({\boldsymbol\alpha}_2\otimes {\boldsymbol\alpha}_1)({\boldsymbol\beta}_2\otimes {\boldsymbol\beta}_1)^{^\intercal}=:{\boldsymbol\alpha}{\boldsymbol\beta}^{^\intercal}. \] Letting ${\boldsymbol\beta}_\bot$ be the orthogonal complement of ${\boldsymbol\beta}$ and defining \[ {\boldsymbol P}=({\boldsymbol P}_1,\ {\boldsymbol P}_2),\quad {\boldsymbol P}_1={\boldsymbol\beta}_\bot\left({\boldsymbol\beta}_\bot^{^\intercal}{\boldsymbol\beta}_\bot\right)^{-1/2},\quad {\boldsymbol P}_2={\boldsymbol\beta}\left({\boldsymbol\beta}^{^\intercal}{\boldsymbol\beta}\right)^{-1/2}, \] we may rotate ${\sf vec}({\boldsymbol F}_{t})$ via \[ {\boldsymbol P}^{^\intercal}{\sf vec}({\boldsymbol F}_{t})=\left[

\right], \] to separate out the full-rank I(1) random vector ${\boldsymbol P}_1^{^\intercal}{\sf vec}({\boldsymbol F}_{t})$ and the stationary I(0) random vector ${\boldsymbol P}_2^{^\intercal}{\sf vec}({\boldsymbol F}_{t})$, see CP09. Hence, the proposed structure is comparable to the matrix factor model setting in CGGT25 which includes both stationary and nonstationary matrix-valued factors and propose a projection-based three-stage estimation method.

Let $e_{t,(i,j)}$ be the $(i,j)$ entry of ${\boldsymbol E}_{t}$ and suppose that

where $|\rho_{ij}|< 1$ or $\rho_{ij}=1$, $L$ denotes the lag operator, and $\{\varepsilon_{t,(i,j)}\}$ is a sequence of stationary I(0) random variables for each $(i,j)$. Note that $e_{t,(i,j)}$ is nonstationary if $\rho_{ij}=1$ and is stationary if $|\rho_{ij}|<1$. Similar assumptions can be found in BN04 and BLL21 for elements of idiosyncratic error vectors.

Neither the mPCA method proposed in Section (ref) nor the projection-based estimation in CGGT25 is applicable when ${\boldsymbol F}_t$ and ${\boldsymbol E}_t$ satisfy ((ref)) and ((ref)), respectively. We next adopt a matrix version of BN04's PANIC estimation technique (or mPANIC). Taking the first-order difference on both sides of ((ref)), we obtain

which becomes a factor model for stationary matrix-valued time series WLC19, CF23. As in Section (ref), we define the sample row covariance matrix using $\Delta{\boldsymbol X}_t$:

and then obtain $\overline{\boldsymbol R}=(\overline{R}_{1},\ldots,\overline{R}_{r_1})$ with $\overline{R}_{k}$ being the eigenvector of $\overline{\boldsymbol\Omega}_{R}$ corresponding to the $k$-th largest eigenvalue. The estimate of column factor loading matrix ${\boldsymbol C}$ is constructed in a similar way, denoted by $\overline{\boldsymbol C}=(\overline{C}_{1},\ldots,\overline{C}_{r_2})$. The cointegrated factor matrix ${\boldsymbol F}_t$ is estimated by

where $\overline\lambda_{R,1}$ is the maximum eigenvalue of $\overline{\boldsymbol\Omega}_{R}$, $\overline{\boldsymbol V}_R$ is an $r_1\times r_1$ diagonal matrix with the diagonal entries being the $r_1$ largest eigenvalues of $\overline{\boldsymbol\Omega}_R$ (in a descending order) whereas $\overline{\boldsymbol V}_C$ is an $r_2\times r_2$ diagonal matrix with the diagonal entries being the $r_2$ largest eigenvalues of $\overline{\boldsymbol\Omega}_C$. Finally, the numbers of row and column factors, $r_1$ and $r_2$, are determined via the ratio criterion:

and

where $\{\overline{\lambda}_{R,k}:\ k=1,\ldots,p_1\}$ and $\{\overline{\lambda}_{C,k}:\ k=1,\ldots,p_2\}$ are sets of eigenvalues (arranged in a descending order) of $\overline{\boldsymbol\Omega}_{R}$ and $\overline{\boldsymbol\Omega}_{C}$, respectively.

Large-sample theory

\setcounter{equation}{0}

In this section, we give some regularity conditions and establish the convergence properties for the mPCA and mPANIC estimators developed in Section (ref). The factor strength is assumed to be homogeneous over (weak) factors in this section and the extension to the case of heterogeneous factor strength will be explored in Section (ref).

Convergence theory for mPCA estimation

We start with the introduction of some notation. Let \[ {\mathbf f}_t={\sf vec}({\boldsymbol F}_t),\ \ {\mathbf u}_t={\sf vec}({\boldsymbol U}_t),\ \ {\mathbf e}_t={\sf vec}({\boldsymbol E}_t). \] Denote the $i$-th row vectors of ${\boldsymbol R}$ and ${\boldsymbol C}$ by ${\boldsymbol R}_{i\bullet}$ and ${\boldsymbol C}_{i\bullet}$, respectively. Let $e_{t,(i,j)}$ be the $(i,j)$ entry of ${\boldsymbol E}_{t}$. We require the following conditions to derive the convergence properties of mPCA estimation.

\setcounter{assumption}{0}

\setcounter{remark}{0}

Due to the model identifiability issue discussed in Section (ref), we can only consistently estimate the row/column factor loading matrices up to appropriate rotation. Define two rotation matrices:

where $\widehat{\boldsymbol V}_R$ is an $r_1\times r_1$ diagonal matrix with the diagonal entries being the $r_1$ largest eigenvalues of $(p_1^{\alpha_R}p_2^{\alpha_C}T)^{-1}\widehat{\boldsymbol\Omega}_R$ (in a descending order) and $\widehat{\boldsymbol V}_C$ is an $r_2\times r_2$ diagonal matrix with the diagonal entries being the $r_2$ largest eigenvalues of $(p_1^{\alpha_R}p_2^{\alpha_C}T)^{-1}\widehat{\boldsymbol\Omega}_C$. Lemma (ref) in Appendix B of the supplement shows that both $\widehat{\boldsymbol H}_R$ and $\widehat{\boldsymbol H}_C$ are invertible {\em w.p.a.1}, and derives their respective limits. The following theorem gives the mean squared convergence rates for the (normalized) mPCA estimates.

\setcounter{theorem}{0}

We next present the uniform convergence rate of the estimated common trends via mPCA.

We require a further condition to establish the consistency of eigenvalue ratio estimation.

The following theorem establishes the consistency property of the ratio criterion.

Convergence theory for mPANIC estimation

We next derive the convergence properties for the mPANIC estimator developed in Section (ref).

As in Section (ref), we define two rotation matrices:

where $\check{\boldsymbol V}_R$ is an $r_1\times r_1$ diagonal matrix with the diagonal entries being the $r_1$ largest eigenvalues of $(p_1^{\alpha_R}p_2^{\alpha_C})^{-1}\overline{\boldsymbol\Omega}_R$ (in a descending order), and $\check{\boldsymbol V}_C$ is defined similarly with $r_1$ and $\overline{\boldsymbol\Omega}_R$ replaced by $r_2$ and $\overline{\boldsymbol\Omega}_C$, respectively. The following theorem gives the mean squared convergence rates for (normalized) mPANIC estimates of row and column factor loadings.

Write

and \[ {\boldsymbol\Gamma}^\circ=\left(\Gamma_{ij}^\circ\right)_{p_1\times p_2}\quad \text{with}\quad \Gamma_{ij}^\circ={\sf I}(|\rho_{ij}|<1)\quad \text{and}\quad{\boldsymbol\Gamma}^\dag=\left(\Gamma_{ij}^\dag\right)_{p_1\times p_2}\quad \text{with}\quad \Gamma_{ij}^\dag={\sf I}(\rho_{ij}=1), \] where ${\sf I}(\cdot)$ denotes the indicator function. It follows from ((ref)) that ${\boldsymbol E}_t$ is decomposed into a sum of I(0) matrix idiosyncratic component ${\boldsymbol E}_t^\circ$ and I(1) component ${\boldsymbol E}_t^\dag$. In addition, we impose the following structural restriction on ${\boldsymbol\Gamma}^\dag$, i.e.,

where ${\boldsymbol\Gamma}_1^\dag$ is an $s_R\times s_C$ block matrix with all the entries being ones, and ${\boldsymbol O}_{k\times l}$ is a null matrix with size $k\times l$. The condition ((ref)) limits the number of nonstationary entries in the idiosyncratic matrix, which is similar to the commonly-used sparsity restriction in high-dimensional econometrics.

We next state the uniform convergence property of the cointegrated matrix factor estimation.

Extension to heterogeneous factor strength

\setcounter{equation}{0}

The convergence theory developed in Section (ref) assumes that the row/column (weak) factor strengths are homogeneous with the constant values $\alpha_R$ and $\alpha_C$, see Assumption (ref)(ii). This assumption, however, may be too restrictive for real-world applications. In this section, we further extend the methodology and theory, allowing for heterogeneous factor strength.

Convergence of mPCA with heterogeneous factor strength

Due to heterogeneity in the weak factor strength, we replace Assumptions (ref)(ii)(iii), (ref)(iii) and (ref) by the following conditions.

\setcounter{assumption}{0}

\setcounter{remark}{0}

Due to the heterogeneous factor strength, a new normalization matrix is required to properly scale the sample moments of row/column factor loadings, see Assumption (ref)(i). Accordingly, the rotation matrices also need to be modified, i.e.,

where $\widehat{\boldsymbol V}_R^\ast = p_2^{-\alpha_{C,1}}\mathbf{B}_R^{-2}\widetilde{\boldsymbol V}_R$, $\widehat{\boldsymbol V}_C^\ast = p_1^{-\alpha_{R,1}}\mathbf{B}_C^{-2}\widetilde{\boldsymbol V}_C$, $\widetilde{\boldsymbol V}_R$ and $\widetilde{\boldsymbol V}_C$ are defined as in (ref). Lemma (ref) in Appendix B of the supplement shows that both $\widehat{\boldsymbol H}_R$ and $\widehat{\boldsymbol H}_C$ are invertible {\em w.p.a.1}. The following theorem extends Theorems (ref)--(ref) to the case of varying factor strength.

\setcounter{theorem}{0}

Convergence of mPANIC with heterogeneous factor strength

Similarly to Section (ref), we replace Assumptions (ref)(iii), (ref), and (ref) by the following assumptions.

Define

where $\check{\boldsymbol V}_R^\ast = p_2^{-\alpha_{C,1}}\mathbf{B}_R^{-2}\overline{\boldsymbol V}_R$, $\check{\boldsymbol V}_C^\ast = p_1^{-\alpha_{R,1}}\mathbf{B}_C^{-2}\overline{\boldsymbol V}_C$, $\overline{\boldsymbol V}_R$ and $\overline{\boldsymbol V}_C$ are defined as in (ref). The following theorem extends Theorems (ref)--(ref) to the case of heterogeneous factor strength.

Numerical studies

\setcounter{equation}{0}

In this section we provide both Monte-Carlo simulation and empirical studies to examine the finite-sample performance of the proposed model and methods.

Simulation studies

To generate common factors, we adopt model (ref) for integrated full-rank $\boldsymbol{F}_t$, and (ref) for cointegrated $\boldsymbol{F}_t$, where we set $r_1=r_2=2$ and the initial value $\boldsymbol{F}_0 = 0$. In (ref), entries of $\boldsymbol{U}_t$ are mutually independent AR(1) processes defined by $$ u_{t,(i,j)} = 0.3 u_{t-1,(i,j)} + \varepsilon_{t,(i,j)}^u \quad \text{with}\quad u_{0,(i,j)}=0,\quad \varepsilon_{t,(i,j)}^u\stackrel{i.i.d.}\sim {\sf N}(0,1). $$ In ((ref)), entries of $\boldsymbol{V}_t$ are mutually independent sequences of i.i.d. standard normal variables, ${\boldsymbol\alpha}_1 = [-0.1,0.1]^{^\intercal}$, ${\boldsymbol\beta}_1 = [1,-1]^{^\intercal}$, ${\boldsymbol\alpha}_2 = [0.1,-0.1]^{^\intercal}$, ${\boldsymbol\beta}_2 = [1,-1]^{^\intercal}$. To generate factor loadings, we set $\mathbf{R} = \mathbf{U}_R\mathbf{B}_R$ and $\mathbf{C} = \mathbf{V}_C\mathbf{B}_C$, where $\mathbf{U}_R$ and $\mathbf{V}_C$ are random orthonormal matrices of dimensions $p_1 \times r_1$ and $p_2 \times r_2$, respectively, which are obtained via the QR decomposition of random matrices with i.i.d. standard normal entries. We consider three settings for the factor loadings with different levels of factor strength:

Weak factors (or factor loadings) exist in cases (ii) and (iii). The idiosyncratic error matrix $\boldsymbol{E}_t$ is generated by $$ \boldsymbol{E}_t = 0.3 \boldsymbol{E}_{t-1} + {\boldsymbol\Gamma}_R^{1/2} {\boldsymbol \Xi}_t {\boldsymbol\Gamma}_C^{1/2},$$ where $\boldsymbol{E}_0 = {\bf O}_{p_1\times p_2}$, $\{{\boldsymbol \Xi}_t\}$ is a sequence of matrices containing i.i.d. standard normal elements, ${\boldsymbol\Gamma}_R$ and ${\boldsymbol\Gamma}_C$ are $p_1\times p_1$ and $p_2\times p_2$ matrices with the $(i,j)$-entry being $0.5^{|i-j|}$.

For each setting, we repeat the simulation $1000$ times. Since both common factors and factor loadings are estimated only up to some rotation matrices, we report the root mean squared error between the estimated and true projection matrices (onto the corresponding column spaces). For instance, we compute the root mean squared error of row factor loading matrix estimation by $$ \text{RMSE}(\widehat{\mathbf{R}}) = \frac{1}{1000}\sum_{i=1}^{1000} \left\|\widehat{\mathbf{R}}^{(i)}\widehat{\mathbf{R}}^{(i)^\intercal} - \mathbf{U}_R\mathbf{U}_R^{^\intercal} \right\|_F, $$ where $\widehat{\mathbf{R}}^{(i)}$ denotes the estimate of $\mathbf{R}$ in the $i$-th replication, and $\mathbf{U}_R$ the matrix of left singular vectors of $\mathbf{R}$. $\text{RMSE}(\widehat{\mathbf{C}})$ and $\text{RMSE}(\widehat{\mathbf{F}})$ are defined similarly. When implementing the eigenvalue ratio estimators (ref), (ref) (ref) and (ref), we set $K_1=K_2=10$.

Table (ref) presents the simulation results for both mPCA and mPANIC under various settings, where $\text{Mean}(\widehat{r}_1)$ and $\text{Mean}(\widehat{r}_2)$ represent the average estimated number of row and column factors over 1000 simulation replications, $\text{CP}(\widehat{r}_1)$ and $\text{CP}(\widehat{r}_2)$ denote the corresponding proportions of correct estimations. Several key insights can be drawn from the simulation results. First, we observe that estimation errors for factor loadings and common factors decrease as the sample size increases, confirming the consistency of both the mPCA and mPANIC methods. Second, mPCA consistently yields smaller estimation errors than mPANIC, as it directly targets the nonstationary time series, resulting in faster convergence. Third, as the weakest factor strength (i.e., $\alpha_{R,2}$ and $\alpha_{C,2}$) decreases, estimation errors increase and the accuracy in estimating the number of factors deteriorates. When factor strengths are $1$ and $0.8$ in case (ii), both the mPCA and mPANIC methods perform well in recovering the true number of factors, with $\text{CP}(\widehat{r}_1)$ and $\text{CP}(\widehat{r}_2)$ close to 1. However, when the weakest factor strength is $0.6$ in case (iii), $\text{CP}(\widehat{r}_1)$ and $\text{CP}(\widehat{r}_2)$ under mPANIC are close to 0, whereas mPCA achieves a higher CP of around 0.9, which is consistent with the implication of Theorems (ref)(iii) and (ref)(iii). In fact, to ensure consistent selection of the factor number in mPANIC, the current simulation design requires both $\alpha_{R,2}$ and $\alpha_{C,2}$ to be larger than $0.75$.

An empirical application

In this section, we apply the proposed model and methodology to analyze international trade flows. By modeling possibly nonstationary trade volumes directly, we aim to examine the evolution pattern of global trade and the international trade network, and further assess the impact of major events such as the global financial crisis and the COVID-19 pandemic. Our empirical analysis is based on monthly bilateral export volumes of commodity goods among 23 countries and regions, obtained from the International Monetary Fund’s Direction of Trade Statistics. The data cover free-on-board export values denominated in US dollars over a 273-month period from January 2000 to September 2022. The raw trade volume data are log-transformed, whereby the first differences correspond to the growth rates of international trade flows. The countries and regions include Australia (AUS), Canada (CAN), Mainland China (CHN), Denmark (DNK), Finland (FIN), France (FRA), Germany (DEU), Hong Kong China (HKG), Indonesia (IDN), Ireland (IRL), Italy (ITA), Japan (JPN), Korea (KOR), Malaysia (MYS), Mexico (MEX), Netherlands (NLD), New Zealand (NZL), Singapore (SGP), Spain (ESP), Sweden (SWE), Thailand (THA), United Kingdom (GBR), and United States (USA).

We first estimate the size of matrix factors. The proposed eigenvalue ratio criterion yields $\widehat{r}_1 = \widehat{r}_2 = 1$ under model (ref), and $\overline{r}_1 = 2,\ \overline{r}_2 = 1$ under model (ref). Figure (ref) displays the sorted eigenvalues of the sample row (left) and column (right) covariance matrices for model (ref) (top panels) and model (ref) (bottom panels). As shown in Figure (ref), the eigenvalue gaps are more pronounced under model (ref), facilitating the identification of the number of factors. In contrast, the eigenvalues under model (ref), where the data are differenced, are more tightly clustered, making it more challenging to determine the true factor dimensionality. Furthermore, for the two estimated factors under model (ref), we apply the model-free method proposed in ZRY19 to test for cointegration. The result, however, indicates that these two factors are not cointegrated. Therefore, we proceed with model (ref) in the subsequent analysis.

Figure (ref) plots the time series of the estimated common factor. From a modeling perspective, the strong trending pattern in the factor reflects its nonstationary feature. Model (ref) explicitly accommodates integrated components, thereby capturing persistent global movements in international trade. In this context, the estimated factor can be interpreted as a global demand index: bilateral trade volumes fall worldwide when it drops; and international trade expands globally when it rises. Thus, the model effectively identifies the dominant source of co-movement in international trade flows, highlighting its relevance for capturing global shocks in a parsimonious and interpretable manner. Furthermore, the dynamic behavior of the estimated factor aligns closely with several major global events. The plot shows a persistent upward trend beginning in 2002, which may be attributed to China's accession to the World Trade Organization and the subsequent acceleration in its trade activities with major economies (e.g., HO16). Notably, it exhibits a sharp drop around 2009, corresponding to the so-called “Great Trade Collapse” due to the 2008 Global Financial Crisis. The estimated factor also drops sharply in early 2020 due to the outbreak of COVID-19 pandemic, followed by a quick rebound, reflecting the V-shaped pattern of global trade disruption and recovery. In addition, a smaller but visible decline is observed around 2015--2016, coinciding with the documented slowdown in global trade during that period (e.g., CMR20).

Figure (ref) displays the heatmap of the normalized matrix $\widehat{\boldsymbol R}\widehat{\boldsymbol C}^{^\intercal}$, where each element has been rescaled to lie between 0 and 1 for ease of comparison. This matrix captures the structure of bilateral trade linkages across 23 economies, conditional on the global common factor estimated from monthly international trade data. Each entry $(i,j)$ in $\widehat{\boldsymbol R}\widehat{\boldsymbol C}^{^\intercal}$ reflects the relative intensity of trade between exporter $i$ and importer $j$ in response to global shocks, such as aggregate demand or systemic disruptions. As $\widehat{\boldsymbol R}$ and $\widehat{\boldsymbol C}$ represent country-specific sensitivities to the global factor, this matrix effectively summarizes how global economic conditions propagate through the trade network. The heatmap reveals that the strongest bilateral linkage arises between China (CHN) and the United States (USA), highlighting their pivotal roles as global trade hubs. This is consistent with their status as the two largest trading nations over the past decades.

Conclusion

\setcounter{equation}{0}

We have introduced a general matrix factor model to tackle large-scale trending matrix time series, making full case of the intrinsic matrix data structure and aiming to improve the estimation convergence. This is substantially different from the conventional vectorized factor model framework which transforms the matrix time series observations into vectors and often results in loss of sample information. We propose two types of PCA-based estimation techniques: mPCA and mPANIC, according to different nonstationary features in matrix common and idiosyncratic components. Under various model assumptions, we establish the convergence theory for the estimated factor loading matrices and nonstationary factor matrices and the consistency of the factor number estimation via the eigenvalue ratio criterion. In particular, we allow the existence of weak factors in the proposed factor model structure, where the factor strength can be heterogeneous. The obtained mean squared convergence rates (for factor loading matrix estimation) and uniform convergence rates (for factor matrix estimation) depend on the matrix size $(p_1,p_2)$, time series length $T$ and the weak factor strength. The Monte-Carlo simulation study justifies the estimation convergence in finite samples. The empirical application to international trade flows shows that the estimated factor can be viewed as the global trade index and the estimated row and column loading matrices are useful to describe bilateral trade linkages between major economies (such as China and the United States).