EconBase
← Back to paper

A Hybrid Framework Combining Autoregression and Common Factors for Matrix Time Series Modeling

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

83,668 characters · 22 sections · 36 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

A Hybrid Framework Combining Autoregression and Common Factors for Matrix Time Series

{16pt}

abstractMatrix-valued time series are ubiquitous in modern economics and finance, yet modeling them requires navigating a trade-off between flexibility and parsimony. We propose the Matrix Autoregressive model with Common Factors (MARCF), a unified framework for high-dimensional matrix time series that bridges the structural gap between the Matrix Autoregression (MAR) and Matrix Factor Model (MFM). While MAR typically assumes distinct predictor and response subspaces and MFM enforces identical ones, MARCF explicitly characterizes the intersection of these subspaces. By decomposing the coefficient matrices into common, predictor-specific, and response-specific components, the framework accommodates distinct input and output structures while exploiting their overlap for dimension reduction. We develop a regularized gradient descent estimator that is scalable for high-dimensional data and can efficiently handle the non-convex parameter space. Theoretical analysis establishes local linear convergence of the algorithm and statistical consistency of the estimator under high-dimensional scaling. The estimation efficiency and interpretability of the proposed methods are demonstrated through simulations and an application to global macroeconomic forecasting.

Keywords: Matrix-valued time series, factor model, autoregression, dimension reduction.

\setlength\abovedisplayskip{4pt} \setlength\belowdisplayskip{4pt}

Introduction

Matrix-valued time series, where each observation is a matrix rather than a vector, have become ubiquitous in modern economics and finance. Prominent examples include multinational macroeconomic indicators (e.g., GDP, inflation, and interest rates observed across multiple countries) and financial panel data characterized by inherent two-way dependencies. Modeling such data presents a fundamental statistical challenge: one must capture complex temporal dynamics and cross-sectional correlations while respecting the intrinsic matrix structure to maintain estimation efficiency and interpretability.

Currently, the literature is dominated by two distinct modeling paradigms, each dealing with the high-dimensionality of matrix data through different structural assumptions. The first paradigm is the Matrix Autoregressive (MAR) framework, which extends vector autoregression to the matrix domain. For instance, RongChen2021MatrixAR introduced a bilinear MAR model, while XiaoRRMAR2024 proposed the Reduced-Rank MAR (RRMAR) model. A series of extensions have been developed to capture more complex data patterns, such as spatio-temporal MAR model Hsu2Spatio-TempMAR, envelope-based MAR model Samadi25EMAR, sparse MAR model jiang24regu_highdim_MAR, and additive MAR models Zhang24additiveMAR, among many others. The defining feature of the RRMAR framework is its flexibility: it allows the subspace containing predictive information (the predictor subspace) to be distinct from the subspace in which the future response evolves (the response subspace) huang2025supervised. While this framework allows for complex dynamics, it ignores potential structural overlaps, leading to parameter inefficiency when strong co-movements exist.

The second paradigm is the Matrix Factor Model (MFM) wang2019factor. These models assume that the high-dimensional observation is driven by a low-dimensional latent factor process. In the literature, two primary classes of assumptions on factors are commonly adopted. The first class assumes that factors are pervasive and influence most observed series, while allowing weak serial dynamic dependence in the error processes stock2002forecasting,bai2002determining,bai2008large,chen2023statistical. The second assumes that the latent factors capture all dynamic dependencies of the observed process, rendering the error process devoid of serial dependence LamAndYao2012,wang2019factor,gao2022modeling,gao2023twotrans. The literature has expanded to include constrained factors Chen24ConstaintMFM and tensor-based decompositions Chang23CPMFM. In contrast to MAR models, MFM frameworks typically enforce a strict symmetry: the loadings that map factors to the observed matrix are invariant over the lag structure. This implies that the predictor and response subspaces are identical. While highly parsimonious, this assumption is restrictive; it limits the ability of the model to capture asymmetric transmission mechanisms where the driving forces of the dynamics transform structurally over time (e.g., a shock in a specific financial sector propagating to the broad real economy).

Consequently, practitioners face a dilemma between the flexibility of MAR and the parsimony of MFM. Ideally, one would seek a hybrid approach that integrates the strengths of both paradigms. In the vector time series literature, significant efforts have been made to combine autoregressive dynamics with factor structures. A seminal example is the Factor-Augmented VAR (FAVAR) bernanke2005measuring,bai2016estimation, which augments standard VAR systems with factors extracted from large datasets to capture broad economic movements. More recently, miao2023high investigated high-dimensional VARs with common factors, employing low-rank constraints on the transition matrix to consistently estimate dynamics in large systems. While these approaches successfully blend factors and autoregressions, they typically treat factors as auxiliary regressors or impose low-rank structures on the vectorized process. Additionally, they do not explicitly model the relationship between the response and predictor spaces. A step in this direction was taken by wang2023commonFactor, who proposed a common basis VAR model that identifies the intersection between predictor and response subspaces. However, all these methods are fundamentally designed for vector-valued data. Direct application to matrix-valued series requires vectorization, which demolishes the inherent two-way structural information, such as the distinct correlations among countries (rows) and economic indicators (columns), and fails to address the specific identifiability challenges imposed by the bilinear matrix structure.

In this paper, we propose the Matrix Autoregressive model with Common Factors (MARCF), a unified framework that generalizes the RRMAR specification and bridges the gap between MAR and MFM. Instead of treating the row and column spaces of the coefficient matrices as arbitrary low-rank subspaces, we explicitly model the intersection between the predictor and response subspaces. We decompose the relevant factor spaces into common components shared between the past and present, predictor-specific components unique to the past, and response-specific components unique to the present; see Section (ref) for a formal introduction and interpretations of these components. Crucially, this decomposition applies simultaneously to both the row and column dimensions, allowing the model to capture not only symmetric dependencies but also complex interaction dynamics. For instance, the framework can model scenarios where dynamics are driven by a common factor in the cross-section (rows) interacting with a specific leading indicator in the variable dimension (columns). This capability allows the model to reveal asymmetric transmission channels, such as global country-level trends being predicted by specific financial variables.

This decomposition offers three key advantages. First, it creates a continuum between the existing paradigms: the model reduces to RRMAR when the intersection is empty ($d_i=0$) and approximates a structured MFM when the subspaces fully align ($d_i=r_i$). Second, by identifying the shared structure, the MARCF model achieves a dimension reduction of approximately $p_1d_1+p_2d_2$ relative to the baseline RRMAR, significantly improving forecasting performance in high-dimensional regimes. Third, it provides a transparent interpretation of the sources of variation, allowing practioners to distinguish between persistent trends, driving forces, and reactive components in a complex time series dataset.

To estimate the model, we develop a computationally efficient regularized gradient descent algorithm. Unlike existing methods that rely on alternating minimization RongChen2021MatrixAR, XiaoRRMAR2024, our approach avoids computationally expensive matrix inversions and scales efficiently with the dimension of the data. We tackle the complex identification issues inherent in bilinear models by introducing a balanced regularization strategy, for which we establish theoretical guarantees regarding statistical consistency, algorithmic convergence, and structural dimension selection consistency.

The remainder of the paper is organized as follows. Section (ref) introduces the MARCF framework and its relationship to RRMAR and MFM. Section (ref) details the estimation strategy, including the regularized loss function, gradient descent algorithm, initialization, and rank selection. Section (ref) presents the theoretical properties of the estimator, including computational convergence, statistical error rates, and rank selection consistency. Section (ref) reports simulation results, and Section (ref) illustrates the model performance and interpretability using a multinational macroeconomic dataset. Technical proofs and supplementary details are provided in the supplementary materials.

Model

Matrix Factor Model and Reduced-Rank MAR

Matrix-valued time series models extend classical multivariate time series methods to data with inherent two-way structure. A fundamental approach in this domain is the Matrix Factor Model (MFM). The MFM assumes that the high-dimensional observed matrix $\mathbf{Y}_t$ is driven by a low-dimensional latent process, given by

equation[equation omitted — 132 chars of source]

where $\boldsymbol{\varepsilon}_t\in\mathbb{R}^{p_1\times p_2}$ is a matrix-valued error term. As discussed in Section (ref), the literature varies in its assumptions regarding the factors and error terms. In this paper, we assume that the latent factors $\mathbf{F}_t$ drive all temporal dynamics and $\boldsymbol{\varepsilon}_t$ is a white noise process.

Another fundamental approach is the Matrix Autoregressive (MAR) model, which captures bilinear dependencies in an observed time series $\mathbf{Y}_t\in\mathbb{R}^{p_1\times p_2}$ through the recursion

equation[equation omitted — 131 chars of source]

where each $\mathbf{A}_i \in \mathbb{R}^{p_i \times p_i}$ is a coefficient matrix and $\mathbf{E}_t$ is a matrix-valued white noise process. To address the high-dimensional challenge of the parameter space, the Reduced-Rank MAR (RRMAR) model XiaoRRMAR2024 imposes low-rank constraints on the coefficient matrices, requiring $\text{rank}(\mathbf{A}_i)=r_i\ll p_i$ for $i=1,2$. We parameterize this constraint using the singular value decomposition (SVD) $\mathbf{A}_i=\mathbf{U}_i\mathbf{S}_i\mathbf{V}_i^\top$, where $\mathbf{U}_i\in\mathbb{R}^{p_i\times r_i}$ and $\mathbf{V}_i\in\mathbb{R}^{p_i\times r_i}$ have orthonormal columns, and $\mathbf{S}_i$ is a diagonal matrix of singular values. Although the SVD components are not unique, the column spaces of $\mathbf{U}_i$ and $\mathbf{V}_i$, as well as the corresponding subspace projection matrices $\mathbf{U}_i\mathbf{U}_i^\top$ and $\mathbf{V}_i\mathbf{V}_i^\top$, are uniquely defined.

The RRMAR and MFM frameworks are closely linked through their subspace representations. First, the low-rank structure of the RRMAR allows for a supervised factor analysis interpretation wang2022high,wang2023commonFactor,huang2025supervised. Specifically, by pre-multiplying (ref) by $\mathbf{U}_1^\top$ and post-multiplying by $\mathbf{U}_2$, we obtain the following low-dimensional representation of the dynamics:

equation[equation omitted — 209 chars of source]

Following huang2025supervised, we view $\mathbf{U}_1^\top\mathbf{Y}_t\mathbf{U}_2$ and $\mathbf{V}_1^\top\mathbf{Y}_t\mathbf{V}_2$ as the response factor and $\textit{predictor factor}$. The term “response factor” is justified because the RRMAR model implies

equation[equation omitted — 301 chars of source]

This takes the form of an MFM as in (ref) where the low-dimensional transformation $\mathbf{U}_1^\top\mathbf{Y}_t\mathbf{U}_2$ summarizes all predictable information in $\mathbf{Y}_t$. On the other hand, the relationship in (ref) implies that the predictor factor $\mathbf{V}_1^\top\mathbf{Y}_{t-1}\mathbf{V}_2$ contains all driving forces in the historical data necessary to predict the future response.

Consider now the case where the latent factor $\mathbf{F}_t$ in MFM follows a MAR process, as proposed by yu2024dynamic:

equation[equation omitted — 122 chars of source]

where $\mathbf{B}_i\in\mathbb{R}^{r_i\times r_i}$ are coefficient matrices and $\boldsymbol{\xi}_t$ is low-dimensional white noise. Substituting (ref) into (ref), the structured dynamics become

equation[equation omitted — 198 chars of source]

where $\mathbf{N}_t$ is a composite noise term involving $\boldsymbol{\xi}_t$ and a moving average of $\mathbf{E}_t$. Equation (ref) reveals the implicit MAR structure of this dynamic MFM, with the effective coefficient matrices $\mathbf{\Lambda}_i\mathbf{B}_i\mathbf{\Lambda}_i^\top$. Comparing this to the RRMAR decomposition in (ref), we observe a fundamental restriction in the MFM framework: the predictor subspace and the response subspace are identical (both spanned by $\mathbf{\Lambda}_i$). While this symmetry is parsimonious, it prevents the model from capturing rotational dynamics where the relevant subspaces shift between the predictor and the response. In contrast, the RRMAR model allows $\mathbf{V}_i$ and $\mathbf{U}_i$ to be distinct, offering greater flexibility but potentially ignoring structural commonalities shared by the response and predictor.

Matrix Autoregression with Common Factors

To reconcile the distinct subspace assumption of the RRMAR model with the shared subspace assumption of the MFM, we propose a hybrid framework. This framework explicitly parameterizes the intersection of the predictor and response subspaces. By controlling the degree of this overlap, we can capture shared dynamics while accommodating distinct structural features in the response and predictor spaces.

Let $\mathcal{M}(\mathbf{U}_i)$ and $\mathcal{M}(\mathbf{V}_i)$ denote the response and predictor subspaces, respectively. We assume their intersection has dimension $d_i$, where $0\leq d_i\leq r_i$. The parameters $d_1$ and $d_2$ quantify the extent of structural overlap. Given these dimensions, there exist orthogonal matrices $\mathbf{O}_{i1}$ and $\mathbf{O}_{i2}$ allowing the decomposition:

equation[equation omitted — 157 chars of source]

Here, $\mathbf{C}_i\in\mathbb{R}^{p_i\times d_i}$ represents the common subspaces shared by both responses and predictors. The matrix $\mathbf{R}_i\in\mathbb{R}^{p_i\times (r_i-d_i)}$ spans the response-specific subspace, while $\mathbf{P}_i\in\mathbb{R}^{p_i\times (r_i-d_i)}$ spans the predictor-specific subspace.

Substituting this decomposition into the general bilinear framework yields the Matrix Autoregressive model with Common Factors (MARCF):

equation[equation omitted — 232 chars of source]

In this formulation, $\mathbf{D}_i=\mathbf{O}_{i1}^\top\mathbf{S}_i\mathbf{O}_{i2}$ serves as a dense mixing matrix within the low-dimensional latent space.

This specification provides a unified continuum between the two existing paradigms. When $d_1=d_2=0$, the common component vanishes, and the model reduces to the RRMAR framework with disjoint subspaces. Conversely, when $d_i=r_i$, the specific components vanish, and the predictor and response subspaces align perfectly, approximating the structure of a dynamic MFM. By estimating $d_1$ and $d_2$ from the data, the MARCF model adaptively balances the efficiency of shared factors with the flexibility of rotational dynamics.

For interpretation, we let $\mathbf{C}_i^\top\mathbf{C}_i=\mathbf{I}_{d_i}$, $\mathbf{R}_i^\top\mathbf{R}_i=\mathbf{P}_i^\top\mathbf{P}_i = \mathbf{I}_{r_i - d_i}$, and $\mathbf{C}_i^\top\mathbf{R}_i = \mathbf{C}_i^\top\mathbf{P}_i = \mathbf{0}_{d_i\times(r_i-d_i)}$; while the detailed identification conditions are formally presented in Section (ref). By pre-multiplying (ref) by $[\mathbf{C}_1~\mathbf{R}_1]^\top$ and post-multiplying by $[\mathbf{C}_2~\mathbf{R}_2]$, we have

equation[equation omitted — 500 chars of source]

This representation extends the factor interpretation in (ref) by explicitly characterizing the common factors shared by responses and predictors, alongside their specific counterparts.

Furthermore, when $d_1$ and $d_2$ are nonzero, the proposed framework achieves significant dimension reduction. The total number of free parameters in the MARCF model is given by

equation[equation omitted — 120 chars of source]

This contrasts with $\mathrm{df}_{\mathrm{RRMAR}} = 2p_1 r_1 + 2p_2 r_2 + r_1^2 + r_2^2$ for the standard reduced-rank model. In high-dimensional settings, the MARCF model reduces the effective dimension of parameter space by approximately $p_1 d_1 + p_2 d_2$, enhancing both estimation efficiency and forecasting accuracy.

Model Interpretation

The MARCF framework offers a transparent decomposition of the sources of variation, enabling an innovative analysis of dynamic mechanisms. To illustrate, consider the multinational macroeconomic application discussed in Section (ref), where rows represent countries and columns represent economic indicators.

\paragraph{(1) Row Dynamics (Country-Level Heterogeneity):} The row decomposition identifies groups of countries with distinct roles in the global system.

itemize• Common Countries ($\mathbf{C}_1$): This identifies systemic economies, such as the G7 nations. These countries exhibit persistent influence: their economic state at time $t-1$ strongly predicts the global state at time $t$, and they remain the primary drivers of variance in the response at time $t$. • Predictor-Specific Countries ($\mathbf{P}_1$): This subspace captures “leading” economies. For instance, it may isolate emerging markets or commodity exporters. While these nations might not dominate total global variance in response, shocks to their economies, such as a sudden shift in manufacturing output, often serve as early warning signals, transmitting information that predicts future movements in the core economies. • Response-Specific Countries ($\mathbf{R}_1$): This captures “reactive” economies, possibly small open economies. These countries absorb global shocks and exhibit volatility in the response period, but their idiosyncratic fluctuations have little predictive power for the future global state.

\paragraph{(2) Column Dynamics (Variable-Level Lead-Lag Relationships):} The column decomposition separates indicators based on their temporal properties.

itemize• Common Variables ($\mathbf{C}_2$): These represent fundamental macroeconomic aggregates like industrial production or GDP. These variables are auto-correlated and central to describing the economic cycle; they predict themselves and are the primary targets of prediction. • Predictor-Specific Variables ($\mathbf{P}_2$): These correspond to financial variables or sentiment indices (e.g., stock returns, PMI). These are classic leading indicators; they contain rich information about the future direction of the economy but may be noisy or less relevant for defining the real economic state at the current moment. • Response-Specific Variables ($\mathbf{R}_2$): These capture lagging indicators, such as unemployment rates. These variables react to past economic conditions but are slow to adjust and do not necessarily drive future cycles.

\paragraph{(3) Interaction Dynamics (Cross-Sectional Spillovers):} Beyond the diagonal blocks, the MARCF model captures crucial cross-terms (e.g., corresponding to $\mathbf{C}_1^\top \mathbf{Y}_{t-1} \mathbf{P}_2$). These interactions allow for asymmetric dependencies where the source of predictive information differs structurally from the response it generates. For instance, the model can describe how predictor-specific variables (e.g., leading indicators, $\mathbf{P}_2$) in common countries ($\mathbf{C}_1$) drive future common real output ($\mathbf{C}_2$). This capability allows the model to capture dynamic transformations where the factors driving prediction differ from those characterizing the response, contrasting with standard factor models that enforce a restrictive symmetry.

Stationarity and Identification

The MARCF model defined in (ref) belongs to the general class of bilinear MAR processes. Following RongChen2021MatrixAR, weak stationarity is guaranteed provided that the spectral radii of the coefficient matrices satisfy $\rho(\mathbf{A}_1)\rho(\mathbf{A}_2) < 1$, where $\mathbf{A}_i = [\mathbf{C}_i~\mathbf{R}_i]\mathbf{D}_i[\mathbf{C}_i~\mathbf{P}_i]^\top$ for $i=1,2$.

The model structure introduces identification challenges at two distinct levels. First, the coefficient matrices $\mathbf{A}_1$ and $\mathbf{A}_2$ are subject to global scaling indeterminacy: for any nonzero scalar $c$, the pairs $(\mathbf{A}_1, \mathbf{A}_2)$ and $(c \mathbf{A}_1, c^{-1} \mathbf{A}_2)$ yield identical $\mathbf{A}=\mathbf{A}_2\otimes\mathbf{A}_1$. In other words, the Kronecker product $\mathbf{A}$ is identifiable, but the row and column coefficient matrices are not. To resolve this, we adopt the constraint $\|\mathbf{A}_1\|_\mathrm{F} = \|\mathbf{A}_2\|_\mathrm{F}$. We define the identifiable true value of $\mathbf{A}$ as $\mathbf{A}^*$, and define the signal strength of the true model as $\phi=\|\mathbf{A}^*\|_\text{F}^{1/2}$ such that $\|\mathbf{A}_1^*\|_\text{F}=\|\mathbf{A}_2^*\|_\text{F} = \phi$.

Second, given a fixed $\mathbf{A}_i$, the decomposition into components $(\mathbf{C}_i,\mathbf{R}_i,\mathbf{P}_i,\mathbf{D}_i)$ is not unique. Specifically, for any invertible matrices $\mathbf{Q}_1\in\mathbb{R}^{d_i\times d_i}$ and $\mathbf{Q}_2,\mathbf{Q}_3\in\mathbb{R}^{(r_i-d_i)\times(r_i-d_i)}$, the parameter sets

equation[equation omitted — 284 chars of source]

are equivalent. To resolve this internal indeterminacy and balance these components, we impose the following constraints:

equation[equation omitted — 323 chars of source]

where $b>0$ is a balancing scalar. To ensure that $\mathbf{D}_i$ is of the same order of magnitude as $(\mathbf{C}_i,\mathbf{R}_i,\mathbf{P}_i)$, we set $b\asymp\phi^{1/3}$. This choice implies that $\|\mathbf{C}_i\|_\text{F} \asymp \|\mathbf{R}_i\|_\text{F} \asymp \|\mathbf{P}_i\|_\text{F} \asymp \|\mathbf{D}_i\|_\text{F}$ when $r_i$ and $d_i$ remain constant.

Under the balancing constraint in (ref), the components are determined uniquely up to orthogonal rotations. We treat parameterizations that differ only by such rotations as equivalent. Let $\mathbf{\Theta}$ denote the full set of component matrices. For any two sets $\mathbf{\Theta}$ and $\mathbf{\Theta}'$ satisfying the identification conditions in (ref), we define their distance as

equation[equation omitted — 613 chars of source]

where $\mathbb{O}^k$ denotes the set of $k\times k$ orthogonal matrices. This metric remains invariant under rotations and forms the basis of our theoretical analysis in Section (ref).

remark[Interpretation vs. Estimation] The parameterization with $b \asymp \phi^{1/3}$ is designed to facilitate theoretical analysis and numerical optimization. However, for model interpretation as discussed in Section (ref), it is often more intuitive to work with orthonormal bases. One can always transform the estimated parameters to an equivalent representation where $b=1$ (i.e., $\mathbf{C}_i^\top \mathbf{C}_i = \mathbf{I}_{d_i}$) by absorbing the scaling factors into $\mathbf{D}_i$. The interpretation of the subspaces and dynamic factors in Section (ref) remains invariant under this rescaling.

Estimation and Modeling Procedure

This section presents the estimation methodology for the MARCF model (ref). Given the high-dimensional and non-convex nature of the parameter space, we employ a computationally efficient strategy that combines spectral initialization with a regularized gradient descent algorithm. We first define the optimization objective, which incorporates specific regularization terms to enforce the balanced identification constraints discussed in Section (ref). Subsequently, we describe an initialization procedure that exploits the subspace structures of the auxiliary RRMAR model. Finally, we propose a data-driven procedure for sequentially selecting the model ranks $(r_1, r_2)$ and the common subspace dimensions $(d_1, d_2)$.

Regularized Gradient Descent Estimation

Consider the observed matrix-valued time series $\{\mathbf{Y}_t\}_{t=0}^T$. For a fixed set of structural dimensions $(r_1, r_2, d_1, d_2)$, the model parameters are collected in $\mathbf{\Theta}=(\{\mathbf{C}_i,\mathbf{R}_i,\mathbf{P}_i,\mathbf{D}_i\}_{i=1}^2)$. Recall that the implied coefficient matrices are given by

equation[equation omitted — 152 chars of source]

Our estimation procedure minimizes the least-squares loss function:

equation[equation omitted — 295 chars of source]

To ensure estimation stability and satisfy the identification conditions outlined in Section (ref), we incorporate two types of regularization. First, to resolve the global scaling indeterminacy between $\mathbf{A}_1$ and $\mathbf{A}_2$, we introduce a regularization term that encourages their Frobenius norms to be equal:

equation[equation omitted — 302 chars of source]

This regularization term, consistent with the related work on reduced-rank models Tu2016,Guquanquan2017, helps stabilize the estimation of these two coefficient matrices. Second, for a given balance scalar $b$, inspired by Han2022 and wang2023commonFactor, we introduce a second regularization term to enforce the identification constraints on the components in (ref)

equation[equation omitted — 287 chars of source]

The full regularized objective function is

equation[equation omitted — 230 chars of source]

where $\lambda_1$ and $\lambda_2$ are tuning parameters.

The model parameters are estimated by minimizing $\overline{\mathcal{L}}$ using a standard gradient descent algorithm, as outlined in Algorithm (ref). Here, $\nabla_{\mathbf{M}}\overline{\mathcal{L}}^{(j)}$ denotes the partial gradient of $\overline{\mathcal{L}}$ with respect to any component $\mathbf{M}\in\{\mathbf{C}_i,\mathbf{R}_i,\mathbf{P}_i,\mathbf{D}_i\}_{i=1}^2$, evaluated at the parameters $\boldsymbol\Theta^{(j)}$ after $j$ iterations. The partial gradients with respect to $\mathbf{C}_i$, $\mathbf{R}_i$, $\mathbf{P}_i$, and $\mathbf{D}_i$ are derived from the chain rule and are provided in detail in Appendix S2. Crucially, all update steps involve only matrix multiplications and additions, avoiding expensive matrix inversions or SVDs at each iteration, making the method highly scalable to high-dimensional data.

algorithm[algorithm omitted — 1,041 chars of source]

Algorithm Initialization

Since the optimization landscape of (ref) is nonconvex, the global convergence of Algorithm (ref) depends critically on high-quality initialization. We propose a spectral initialization strategy that leverages the structure of the simpler RRMAR model. The core idea is to first estimate the aggregate predictor and response subspaces, and then extract their intersection to initialize the common components.

Given the dimensions $(r_1, r_2, d_1, d_2)$ and scalar $b$, the initialization proceeds in two steps:

enumerate• First obtain a consistent estimate of the total subspaces by fitting the RRMAR model: \begin{equation} \widehat{\mathbf A}_1^{RR}, \widehat{\mathbf A}_2^{RR}:=\operatorname*{arg\,min}_{\mathrm{rank}(\mathbf A_i)\leq r_i, i=1,2}\frac{1}{2T}\sum_{i=1}^{T}\left\| \mathbf Y_t-\mathbf A_1\mathbf Y_{t-1}\mathbf A_2^\top \right\|_{\mathrm{F}}^2. \end{equation} This can be solved using the ALS method XiaoRRMAR2024. Let $\widehat{\mathbf A}_i^\mathrm{RR}$ be rescaled such that $\|\widehat{\mathbf A}_1^\mathrm{RR}\|_\text{F} = \|\widehat{\mathbf A}_2^\mathrm{RR}\|_\text{F}$. We compute the rank-$r_i$ SVD: $\widehat{\mathbf A}_i^\mathrm{RR} = \widehat{\mathbf U}_i \widehat{\mathbf S}_i \widehat{\mathbf V}_i^\top$. Here, $\widehat{\mathbf U}_i$ spans the sum of the common and response-specific subspaces, while $\widehat{\mathbf V}_i$ spans the sum of the common and predictor-specific subspaces. • Use spectral methods to decompose the subspaces of $\widehat{\mathbf A}_1^\mathrm{RR}$ and $\widehat{\mathbf A}_2^\mathrm{RR}$ to obtain $\boldsymbol\Theta^{(0)}$. We first discuss the rationale of the spectral method. Denote $\mathbf U_i:=[\mathbf C_i~\mathbf R_i]$ and $\mathbf V_i:=[\mathbf C_i~\mathbf P_i]$ with orthonormal columns. For any matrix $\mathbf{M}$ with orthonormal columns, denote $\mathcal{P}_{\mathbf{M}}:=\mathbf{M}\mathbf{M}^\top$ as the subspace projection matrix and $\mathcal{P}_{\mathbf{M}}^\perp:=\mathbf{I}_p-\mathbf{M}\mathbf{M}^\top$ as the subspace orthogonal complement projector. Straightforward algebra shows that $\mathcal{P}_{\mathbf U_i}\mathcal{P}_{\mathbf V_i}^\perp=\mathcal{P}_{\mathbf R_i}\mathcal{P}_{\mathbf P_i}^\perp$. This implies that $\mathcal{M}(\mathbf R_i)$ is spanned by the first $r_i-d_i$ left singular vectors of $\mathcal{P}_{\mathbf U_i}\mathcal{P}_{\mathbf V_i}^\perp$. A similar argument holds for $\mathcal{M}(\mathbf P_i)$. Additionally, $\mathcal{P}_{\mathbf R_i}^\perp \,\mathcal{P}_{\mathbf P_i}^\perp \bigl(\mathcal{P}_{\mathbf U_i} + \mathcal{P}_{\mathbf V_i}\bigr) \mathcal{P}_{\mathbf R_i}^\perp \,\mathcal{P}_{\mathbf P_i}^\perp \;=\;2\,\mathcal{P}_{\mathbf C_i},$ implying that $\mathcal{M}(\mathbf C_i)$ is spanned by the leading $d_i$ eigenvectors of the left-hand side. Based on these, we obtain the initial values as follows: \begin{enumerate} • For each $i=1,2$, using the SVD components in the first step, calculate the top $r_i-d_i$ left singular vectors of $\mathcal{P}_{\widehat{\mathbf U}_i}\mathcal{P}_{\widehat{\mathbf V}_i}^\perp$ and $\mathcal{P}_{\widehat{\mathbf V}_i}\mathcal{P}_{\widehat{\mathbf U}_i}^\perp$, and denote them by $\widetilde{\mathbf{R}}_i$ and $\widetilde{\mathbf{P}}_i$, respectively. Calculate the top $d_i$ eigenvectors of $\mathcal{P}_{\widetilde{\mathbf R}_i}^\perp\mathcal{P}_{\widetilde{\mathbf P}_i}^\perp\left(\mathcal{P}_{\widehat{\mathbf U}_i}+\mathcal{P}_{\widehat{\mathbf V}_i}\right)\mathcal{P}_{\widetilde{\mathbf R}_i}^\perp\mathcal{P}_{\widetilde{\mathbf P}_i}^\perp,$ and denote them by $\widetilde{\mathbf{C}}_i$. Then calculate $\widetilde{\mathbf{D}}_1=[\widetilde{\mathbf{C}}_1 ~\widetilde{\mathbf{R}}_1]^{\top}~ \widehat{\mathbf{A}}_1^\text{RR}[\widetilde{\mathbf{C}}_1~\widetilde{\mathbf{P}}_1]$ and $\widetilde{\mathbf{D}}_2=[\widetilde{\mathbf{C}}_2~ \widetilde{\mathbf{R}}_2]^{\top}~ \widehat{\mathbf{A}}_2^\text{RR}[\widetilde{\mathbf{C}}_2~\widetilde{\mathbf{P}}_2]$. • For each $i=1,2$, set $\mathbf{C}_i^{(0)}=b \widetilde{\mathbf{C}}_i, \mathbf{R}_i^{(0)}=b \widetilde{\mathbf{R}}_i, \mathbf{P}_i^{(0)}=b \widetilde{\mathbf{P}}_i$, and $\mathbf{D}_i^{(0)}=b^{-2} \widetilde{\mathbf{D}}_i$. \end{enumerate}

Selection of Ranks and Common Dimensions

Although our estimation methodology is applicable to any pre-specified $(r_1,r_2,d_1,d_2)$, the appropriate selection of these parameters is critical in practice. These structural parameters define the model's complexity and determine the interpretability of the factors. Since the true ranks are unknown in real-world applications, we propose a two-step procedure to determine $(r_1,r_2)$ and $(d_1,d_2)$ sequentially.

For rank selection, since $r_1$ and $r_2$ are typically much smaller than $p_1$ and $p_2$, we choose two upper bounds for the selection: $\bar{r}_1\ll p_1$ and $\bar{r}_2\ll p_2$. Then, we estimate the RRMAR model (ref) with ranks $\bar{r}_1$ and $\bar{r}_2$. Denote the rescaled solutions by $\widehat{\mathbf A}_1^\mathrm{RR}(\bar{r}_1)$ and $\widehat{\mathbf A}_2^\mathrm{RR}(\bar{r}_2)$, respectively. Let $\widehat{\sigma}_{i,1}\geq \widehat{\sigma}_{i,2}\geq ...\geq \widehat{\sigma}_{i,\bar{r}_i}$ be the singular values of $\widehat{\mathbf A}_i^\mathrm{RR}(\bar{r}_i)$. Motivated by Xia2015SelectR, we introduce a ridge-type ratio to determine $r_1$ and $r_2$ separately: $$ \widehat{r}_i=\underset{1 \leqslant j < \bar{r}_i}{\arg \min } \frac{\widehat{\sigma}_{i,j+1}+s(p_1,p_2, T)}{\widehat{\sigma}_{i,j}+s(p_1,p_2, T)}, $$ where we suggest $s(p_1,p_2,T)=\sqrt{(p_1+p_2)\log(T)/(20T)}$. The theoretical guarantee of this estimator is given in Section (ref), and its empirical performance is justified by the numerical experiments in Section (ref). Once $s(p_1,p_2,T)$ is properly specified and $\bar{r}_i> r_i$, the method is not sensitive to the choice to $\bar{r}_i$. Therefore, in practice, if $p_i$ is not too large and the computational cost is affordable, $\bar{r}_i$ can be chosen largely or even $\bar{r}_i=p_i$.

After $r_1$ and $r_2$ are determined, the problem reduces to selecting a model in low dimensions. We use the BIC criterion to select $d_1$ and $d_2$. Let $\widehat{\mathbf A}(\widehat{r}_1,\widehat{r}_2,d_1,d_2)$ be the estimator of $\mathbf A_2\otimes \mathbf A_1$ in (ref) under $\widehat{r}_i$ and $d_i$. Then, $$ \mathrm{BIC}(d_1,d_2) = T p_1p_2 \log \left(\sum_{t=1}^T\|\text{vec}(\mathbf{Y}_t)-\widehat{\mathbf A}(\widehat{r}_1,\widehat{r}_2,d_1,d_2) \text{vec}(\mathbf{Y}_{t-1})\|_2^2\right)+\textup{df}_{\mathrm{MARCF}}\log (T), $$ where $\textup{df}_{\mathrm{MARCF}}$ is defined in (ref) with $\widehat{r}_1$ and $\widehat{r}_2$. Then, we choose $d_1$ and $d_2$ as: $$ \widehat{d}_1,\widehat{d}_2=\min_{1\leq d_i\leq \widehat{r}_i, i=1,2}\mathrm{BIC}(d_1,d_2). $$

remarkWhile simultaneous selection of $r_i$ and $d_i$ via BIC is theoretically possible, the search space of size $O(\bar{r}_1^2 \bar{r}_2^2)$ is computationally prohibitive. The proposed sequential approach balances statistical accuracy with computational feasibility. Rolling window cross-validation may also be used for hyperparameter tuning but incurs significantly higher computational costs.

Theory

In this section, we investigate the theoretical properties of the proposed estimation methodology for the MARCF model. To facilitate a clear understanding of the results, we outline a roadmap that decouples the optimization landscape from the stochastic nature of the data.

First, in Section (ref), we establish the computational convergence of the algorithm via a deterministic analysis. Theorem (ref) demonstrates that Algorithm (ref) achieves local linear convergence, conditioned on a valid initialization and the local curvature of the loss function (characterized by Restricted Strong Convexity and Smoothness). Subsequently, in Section (ref), we provide the stochastic analysis under mild assumptions on the data-generating process. We verify that the conditions required for computational convergence hold with high probability and derive the final statistical error rates for the estimator. Finally, Section (ref) proves the consistency of the rank selection procedure.

Computational Convergence Analysis

We first analyze the convergence properties of Algorithm (ref). Due to the nonconvex nature of the objective function $\overline{\mathcal{L}}$, our analysis relies on the assumptions of restricted strong convexity (RSC) and restricted strong smoothness (RSS), which characterize the curvature of the loss function along the directions of relevant low-rank matrices. These properties are defined with respect to the Kronecker product parameter matrix $\mathbf A=\mathbf A_2\otimes \mathbf A_1$. Let $\widetilde{\mathcal{L}}(\mathbf A)$ denote the least-squares loss function with respect to $\mathbf A$:

equation[equation omitted — 188 chars of source]
definition[RSC/RSS Condition] The loss function $\widetilde{\mathcal{L}}(\mathbf{A})$ is said to be restricted strongly convex (RSC) with parameter $\alpha$ and restricted strongly smooth (RSS) with parameter $\beta$ (where $0<\alpha\leq \beta$) if, for any two parameter matrices $\mathbf{A}=\mathbf{A}_2\otimes\mathbf{A}_1$ and $\mathbf{A}'=\mathbf{A}_2'\otimes\mathbf{A}_1'$ admitting the MARCF structural decomposition presented in (ref), the following inequality holds: \begin{equation} \frac{\alpha}{2}\left\|\mathbf{A}-\mathbf{A}'\right\|_{\mathrm{F}}^2 \leqslant \widetilde{\mathcal{L}}(\mathbf{A})-\widetilde{\mathcal{L}}\left(\mathbf{A}'\right)-\left\langle\nabla \widetilde{\mathcal{L}}\left(\mathbf{A}'\right), \mathbf{A}-\mathbf{A}'\right\rangle \leqslant \frac{\beta}{2}\left\|\mathbf{A}-\mathbf{A}'\right\|_{\mathrm{F}}^2. \end{equation}

The statistical error is quantified by the following deviation bound.

definition[Deviation Bound] For given ranks $(r_1,r_2)$, common subspace dimensions $(d_1,d_2)$, and the identifiable true value $\mathbf A^*$, the deviation bound is defined as \begin{equation} \xi(r_1,r_2,d_1,d_2):=\sup_{\substack{[\mathbf C_i \mathbf R_i]\in\mathbb{O}^{p_i\times r_i},\\ [\mathbf C_i \mathbf P_i]\in\mathbb{O}^{p_i\times r_i},\\\mathbf D_i\in\mathbb{R}^{r_i},\left\| \mathbf D_i \right\|_{\mathrm{F}}=1}}\left\langle \nabla\widetilde{\mathcal{L}}(\mathbf{A}^*), [\mathbf{C}_2 \mathbf{R}_2]\mathbf{D}_2[\mathbf{C}_2 \mathbf{P}_2]^\top\otimes[\mathbf{C}_1 \mathbf{R}_1]\mathbf{D}_1[\mathbf{C}_1 \mathbf{P}_1]^\top\right\rangle. \end{equation}

To quantify the estimation error, we consider $\left\| \mathbf A_2\otimes\mathbf A_1-\mathbf A^* \right\|_{\mathrm{F}}^2$, which is invariant under scaling and rotation. For the truth $\mathbf{A}_i^*$ and its decomposition components $\mathbf{C}_i^*$, $\mathbf{R}_i^*$, and $\mathbf{P}_i^*$, we follow the identification constraint in Section (ref) to require $\|\mathbf A_1^*\|_\mathrm{F}=\|\mathbf A_2^*\|_\mathrm{F}=\phi$, $[\mathbf C_i^*~\mathbf R_i^*]^\top[\mathbf C_i^*~\mathbf R_i^*]=\phi^{2/3}\mathbf I_{r_i}$ and $[\mathbf C_i^*~\mathbf P_i^*]^\top[\mathbf C_i^*~\mathbf P_i^*]=\phi^{2/3}\mathbf I_{r_i}$, for $i=1,2$. As defined in (ref), $\mathrm{dist}(\boldsymbol \Theta^{(j)},\boldsymbol \Theta^*)^2$ measures the estimation error of the parameters after $j$ iterations, and naturally, $\mathrm{dist}(\boldsymbol \Theta^{(0)},\boldsymbol \Theta^*)^2$ represents the initialization error. Let $\underline{\sigma}:=\min(\sigma_{1,r_1},\sigma_{2,r_2})$ be the smallest sigular value among all the non-zero singular values of $\mathbf A_1^*$ and $\mathbf A_2^*$. Define $\kappa:=\phi/\underline{\sigma}$, which quantifies the ratio between the overall signal and the minimal signal. Then, we have the following sufficient conditions for local convergence of the algorithm.

theorem[Local Linear Convergence] Suppose that the RSC/RSS condition is satisfied with $\alpha$ and $\beta$. If $\lambda_1\asymp\alpha$, $\lambda_2\asymp\kappa^{-2}\phi^{8/3}\alpha$, $b\asymp\phi^{1 / 3}$, and $\eta=\eta_0\phi^{-10/3}\beta^{-1}$ with $\eta_0$ being a sufficiently small positive constant, $\xi^2\lesssim \kappa^{-6}\phi^4\alpha^3\beta^{-1}$, and the initialization error $\mathrm{dist}(\boldsymbol \Theta^{(0)},\boldsymbol \Theta^*)^2\leq c_0\kappa^{-2}\phi^{2/3}\alpha\beta^{-1}$, where $c_0$ is a small constant, then for all $j \geq 1$, \begin{equation} \begin{split} \mathrm{dist}(\boldsymbol \Theta^{(j)},\boldsymbol \Theta^*)^2 \leq & (1-C\eta_0\alpha\beta^{-1}\kappa^{-2})^j\mathrm{dist}(\boldsymbol \Theta^{(0)},\boldsymbol \Theta^*)^2 +C\kappa^4\phi^{-10/3}\alpha^{-2}\xi^2(r_1,r_2,d_1,d_2),\\ \end{split} \end{equation} \begin{equation} \begin{split} \left\| \mathbf A_1^{(j)}-\mathbf A_1^* \right\|_{\mathrm{F}}^2+\left\| \mathbf A_2^{(j)}-\mathbf A_2^* \right\|_{\mathrm{F}}^2 \lesssim & \kappa^2(1-C\eta_0\alpha\beta^{-1}\kappa^{-2})^j\left[\left\| \mathbf A_1^{(0)}-\mathbf A_1^* \right\|_{\mathrm{F}}^2+\left\| \mathbf A_2^{(0)}-\mathbf A_2^* \right\|_{\mathrm{F}}^2\right]\\ &+\kappa^4\phi^{-2}\alpha^{-2}\xi^2(r_1,r_2,d_1,d_2), \end{split} \end{equation} and \begin{equation} \begin{split} \left\| \mathbf A_2^{(j)}\otimes\mathbf A_1^{(j)}-\mathbf A_2^*\otimes\mathbf A_1^* \right\|_{\mathrm{F}}^2\lesssim & \kappa^2(1-C\eta_0\alpha\beta^{-1}\kappa^{-2})^j\left\| \mathbf A_2^{(0)}\otimes\mathbf A_1^{(0)}-\mathbf A_2^*\otimes\mathbf A_1^* \right\|_{\mathrm{F}}^2\\ &+\kappa^4\alpha^{-2}\xi^2(r_1,r_2,d_1,d_2). \end{split} \end{equation}

Theorem (ref) confirms that, provided a good initialization, Algorithm (ref) converges linearly to a neighborhood of the truth defined by the statistical noise level $\xi$. Importantly, the conditions on the hyperparameters do not explicitly depend on the ambient dimensions $p_1, p_2$ or sample size $T$, but rather on the signal structure ($\alpha, \beta, \kappa$), highlighting the scalability of the method. The RSC/RSS and initialization conditions are justified in next subsection.

Statistical Convergence Rates

We now establish the statistical guarantees by verifying the deterministic conditions of Theorem (ref) under a probabilistic framework. We impose the following standard assumptions on the data generating process.

assumption[Stationarity] The spectral radius of $\mathbf{A}^*=\mathbf A_1^*\otimes \mathbf A_2^*$ is strictly less than one.
assumption[Sub-Gaussian Noise] The white noise $\mathbf E_t$ can be represented as $\mathrm{vec}(\mathbf E_t)=\boldsymbol\Sigma_{\mathbf e}^{1/2}\boldsymbol\zeta_t$, where $\boldsymbol\Sigma_{\mathbf e}$ is a positive definite matrix, $\{\boldsymbol\zeta_t\}$ are independent and identically distributed random variables with $\mathbb{E}[\boldsymbol \zeta_t]=0$ and $\mathrm{Cov}(\boldsymbol \zeta_t)=\mathbf I_{p_1p_2}$. The entries of $\boldsymbol\zeta_t$, denoted by $\{\zeta_{it}\}_{i=1}^{p_1p_2}$, are independent $\tau^2$-sub-Gaussian, i.e., $\mathbb{E}\left[\exp \left(\mu \zeta_{i t}\right)\right] \leq \exp \left(\tau^2 \mu^2 / 2\right)$ for any $\mu \in \mathbb{R}$ and $1 \leq i \leq p_1p_2$.

Assumption (ref) ensures the existence of a unique strictly stationary solution to the proposed model. The sub-Gaussianity condition in Assumption (ref) is standard in high-dimensional time series literature, such as Zheng2020LinearRisVAR and wang2024LRtsAR.

To ensure the identifiability of the common versus specific subspaces, we require a minimal separation between the response-specific and predictor-specific components. We quantify the separation between subspaces using $\sin\theta$ distance. Let $s_{i,1} \geq \cdots \geq s_{i,r_i-d_i} \geq 0$ be the singular values of $\mathbf{R}_i^{\top} \mathbf{P}_i$. Then, the canonical angles are defined as $\theta_k(\mathbf{R}_i, \mathbf{P}_i)=\arccos \left(s_{i,k}\right)$ for $1 \leq k \leq r_i-d_i$ and $i=1,2$.

assumption[Minimal Gap between Specific Subspaces] There exists a constant $g_{\min}>0$ such that $\min_{i=1,2}\sin \theta_1\left(\mathbf{R}_i^*, \mathbf{P}_i^*\right) \geq g_{\min }$.

We quantify the dependency structure of the time series following the spectral approach as in BasuandGeorgeMichailidis2015. The characterstic matrix polynomial of the MARCF model is given by $\mathcal{A}(z)=\mathbf{I}_p-\mathbf{A}^* z$ for $z\in\mathbb{C}$. Define $\mu_{\min }(\mathcal{A})=\min _{|z|=1} \lambda_{\min }\left(\mathcal{A}^{\dagger}(z) \mathcal{A}(z)\right)$ and $\mu_{\max }(\mathcal{A})=$ $\max _{|z|=1} \lambda_{\max }\left(\mathcal{A}^{\dagger}(z) \mathcal{A}(z)\right)$, where $\mathcal{A}^{\dagger}(z)$ is the conjugate transpose of $\mathcal{A}(z)$. We define the key theoretical quantities:

equation[equation omitted — 598 chars of source]

The following proposition verifies that the RSC and RSS conditions required for optimization hold with high probability given sufficient sample size.

proposition[Condition Verification] Under Assumptions (ref)--(ref), if $T\gtrsim [M_2^{-2}\max\left(\tau,\tau^2\right) + g_{\min}^{-4}\kappa^2\underline{\sigma}^{-4}\alpha_{\textup{RSC}}^{-3}\beta_{\textup{RSS}} \tau^2M_1^2](p_1r_1+p_2r_2+4r_1^2r_2^2)$, we have that with probability at least $1-C\exp \left(-C\left(p_1r_1+p_2r_2\right)\right)$, the RSC/RSS condition in (ref) holds with $\alpha=\alpha_{\textup{RSC}}$ and $\beta=\beta_{\textup{RSS}}$ and the initial value $\mathbf{\Theta}^{(0)}$ obtained as in Section (ref) satisfies $\mathrm{dist}^2_{(0)}\lesssim \kappa^{-2}\phi^{2/3}\alpha_{\textup{RSC}}\beta_{\textup{RSS}}^{-1}$.

Combining the computational convergence with the bound on the statistical deviation $\xi$, we derive the final error rates for the MARCF estimator.

theorem[Statistical Error Rates] Under Assumptions (ref)--(ref), if the hyperparameters $\lambda_1$, $\lambda_2$, $b$, and $\eta$ are set as in Theorem (ref) with $\alpha=\alpha_\mathrm{RSC}$ and $\beta=\beta_\mathrm{RSS}$, $T\gtrsim [M_2^{-2}\max\left(\tau,\tau^2\right) + g_{\min}^{-4}\kappa^2\underline{\sigma}^{-4}\alpha_{\textup{RSC}}^{-3}\beta_{\textup{RSS}} \tau^2M_1^2](p_1r_1+p_2r_2+4r_1^2r_2^2)$, and \begin{equation} J\gtrsim\frac{\log\left(\kappa^{-2}g_{\min}^{4}\right)}{\log(1-C\eta_0\alpha_\mathrm{RSC}\beta_\mathrm{RSS}^{-1}\kappa^{-2})}, \end{equation} then, with probability at least $1-C\exp\left(-C(p_1r_1+p_2r_2)\right)$, $$ \left\| \mathbf A_1^{(J)}-\mathbf A_1^* \right\|_{\mathrm{F}}^2+\left\| \mathbf A_2^{(J)}-\mathbf A_2^* \right\|_{\mathrm{F}}^2\lesssim \kappa^4\phi^{-2}\alpha_\mathrm{RSC}^{-2}\tau^2M_1^2\frac{\mathrm{df}_\mathrm{MARCF}}{T}, $$ and $$ \left\| \mathbf A_2^{(J)}\otimes \mathbf A_1^{(J)}-\mathbf A_2^*\otimes \mathbf A_1^* \right\|_{\mathrm{F}}^2\lesssim \kappa^4\alpha_\mathrm{RSC}^{-2}\tau^2M_1^2\frac{\mathrm{df}_\mathrm{MARCF}}{T}. $$ where $\mathrm{df}_\mathrm{MARCF}$ is the effective number of free parameters defined in (ref).

Theorem (ref) establishes that the estimator achieves the parametric rate $O_p(\mathrm{df}_\mathrm{MARCF}/T)$. Since $\mathrm{df}_\mathrm{MARCF} =\mathrm{df}_\mathrm{RRMAR} - (p_1d_1 + p_2d_2)$, the MARCF model offers a demonstrable variance reduction over the standard RRMAR model by exploiting the common subspace structure. Moreover, the required number of iterations in (ref) does not explicitly depend on $p_1$, $p_2$, or $T$, implying the scalability and computational efficiency of our method even in high-dimensional regimes.

Consistency of Rank Selection

In this subsection, we establish the consistency of rank selection under the standard high-dimensional asymptotic regime lam2011estimation, LamAndYao2012, where $T$, $p_1$, and $p_2\to\infty$ with $r_1$ and $r_2$ remaining fixed.

theorem[Rank Selection Consistency] Under Assumptions (ref) and (ref), if $\bar{r}_1\geq r_1, \bar{r}_2\geq r_2$, $T\gtrsim (p_1r_1+p_2r_2+4r_1^2r_2^2)M_2^{-2}\max\left(\tau,\tau^2\right)$, $\phi^{-1}\alpha_{\textup{RSC}}^{-1}\tau M_1\sqrt{(p_1+p_2)/T}=o(s(p_1,p_2,T))$, and $s(p_1,p_2,T)=o(\sigma_{i,r_i}^{-1}\min_{1\leq j\leq r_i-1}\sigma_{i,j+1}/\sigma_{i,j})$ for $i=1,2$, then we have $\mathbb{P}(\widehat{r}_1=r_1,\widehat{r}_2=r_2)\to 1$, as $T\to \infty$ and $p_1,p_2\to\infty$.

This theorem provides sufficient conditions under which the selected ranks $\widehat{r}_1$ and $\widehat{r}_2$ converge in probability to the true ranks $r_1$ and $r_2$. In practice, if the parameters $M_1$, $M_2$, $\alpha_{\textup{RSC}}$, $\sigma_{1,r_1}$, and $\sigma_{2,r_2}$ are bounded, and $\phi$ is either bounded or diverges slowly as $p_1,p_2\to\infty$, the conditions simplify to requiring that $s(p_1,p_2,T)\to 0$ and $\sqrt{(p_1+p_2)/T}/s(p_1,p_2,T)\to 0$. These simplified conditions are more tractable in applications and help ensure the interpretability and reliability of the rank selection procedure. Particularly, the value $s(p_1,p_2,T)=\sqrt{(p_1+p_2)\log(T)/(20T)}$ suggested in Section (ref) satisfies these requirements, ensuring consistency across a wide range of settings for $p_1$, $p_2$, and $T$.

Simulation Studies

In this section, we evaluate the finite-sample performance of the proposed MARCF model through a series of simulation experiments. The objectives are three-fold: (i) to assess the accuracy of parameter estimation and model selection, (ii) to validate the model's ability to identify structured factor-driven dynamics, and (iii) to compare its performance with existing approaches.

Experiment I: Estimation Accuracy and Model Selection

In the first experiment, we generate data from the MARCF model defined in (ref). We first assess the reliability of the proposed procedures for selecting the model ranks and common dimensions. Then, we compare the MARCF model with two competing approaches: the full-rank MAR model RongChen2021MatrixAR and the RRMAR model RongChen2021MatrixAR, by evaluating their estimation accuracy of $\mathbf A_2 \otimes \mathbf A_1$.

We fix the model ranks at $r_1=2$ and $r_2=3$, and consider all combinations where $d_1\in\{0,1,2\}$ and $d_2\in\{0,1,2,3\}$. Two settings for dimensions and sample size are examined: (1) $p_1=20$, $p_2=10$, and $T=500$; and (2) $p_1=30$, $p_2=20$, and $T=1000$. For each combination of $\{p_1,p_2,d_1,d_2,T\}$, the matrices $\mathbf A_1^*$, $\mathbf A_2^*$, and the data are generated as follows:

enumerate• For $i=1, 2$, generate $\mathbf D_i$ as $\mathbf D_i=\mathbf O_{i,1}^\top\mathbf S_i\mathbf O_{i,2}$, where $\mathbf O_{i,1}$ and $\mathbf O_{i,2}$ are random orthogonal matrices, and $\mathbf S_i \in \mathbb{R}^{r_i\times r_i}$ is a diagonal matrix with entries drawn from $U(0.8,1)$ when $r_i=d_i$, and from $U(0.9,1.1)$ when $r_i\neq d_i$. • Generate two random matrices $\mathbf C_1\in\mathbb{R}^{p_1\times r_1}$ and $\mathbf C_2\in\mathbb{R}^{p_2\times r_2}$ with orthonormal columns. Throughout this section, random matrices with orthonormal columns are constructed by applying QR decomposition to matrices with i.i.d. standard normal entries. Then, generate random matrices $\mathbf M_{i,1}$ and $\mathbf M_{i,2}$, and apply QR decomposition onto $(\mathbf I-\mathbf C_i\mathbf C_i^\top)\mathbf M_{i,1}$ and $(\mathbf I-\mathbf C_i\mathbf C_i^\top)\mathbf M_{i,2}$ to obtain $\mathbf R_i$ and $\mathbf P_i$. Results with $\sin\theta_1<0.8$ are rejected. • Combine the components to form $\mathbf A_1$ and $\mathbf A_2$. We check for stationarity and repeat the above steps if $\rho(\mathbf A_1)\cdot\rho(\mathbf A_2)\geq 1$. The white noise is generated with $\mathrm{vec}(\mathbf E_t)\overset{i.i.d.}\sim N(\mathbf 0,\mathbf I)$, and $\{\mathbf Y_t\}_{t=0}^T$ is generated according to model (ref).

For each pair of $(d_1,d_2)$, the MARCF model is trained using the full procedure proposed in Section (ref), including rank selection, common dimension selection, initialization, and final gradient descent estimation (see Appendix S1 of supplementary materials for details). For rank selection, we set $\bar{r}_1=\bar{r}_2=8$ with $s\left(p_1, p_2, T\right)=\sqrt{\left(p_1+p_2\right) \log (T) /(20 T)}$. Both the RRMAR and MARCF models are estimated using our gradient descent algorithm, as RRMAR is a special case of MARCF. Based on preliminary experiments, we observed that the model performance is insensitive to variations in $\lambda_1, \lambda_2$, and $b$. Accordingly, we set these parameters to 1. We set $\lambda_1=\lambda_2=b=1$ and $\eta=0.001$. Both initialization and final estimation are performed by running Algorithm (ref) until convergence or for a maximum of 1000 iterations. No divergence was observed under these hyperparameter settings.

Table (ref) presents the accuracy of selecting $(r_1, d_1, r_2, d_2)$ across various settings. A trial is considered successful only if all four parameters are correctly identified. The results show that in both settings of $p_1$, $p_2$, and $T$, the proposed selection procedure achieves a success rate approaching 1. Even when $d_1$ and $d_2$ are close to $r_1$ and $r_2$, the procedure maintains high accuracy.

table[table omitted — 808 chars of source]

Figure (ref) displays the relative estimation errors of the three models under comparison. The error is defined as $\|\widehat{\mathbf A}_2\otimes \widehat{\mathbf A}_1-\mathbf A_2^*\otimes \mathbf A_1^*\|_\text{F}/\|\mathbf A_2^*\otimes \mathbf A_1^*\|_\text{F}$. For each model, the lower and upper lines of each error bar represent the first and third quartiles of the 500 errors, respectively and the cross markers indicate the medians. In all cases, the estimation errors of the MAR model are significantly larger than those of the other models, as expected, due to its failure to account for the low-rank structure in the data. When $d_1=d_2=0$, the MARCF model reduces to the RRMAR model, and the estimation errors of the two are identical. However, as $d_1$ and $d_2$ increase, the error bars of the MARCF model decrease and gradually diverge from those of the RRMAR model, highlighting the improved performance of MARCF model when common structures of the subspaces exist. This performance advantage becomes most pronounced when $d_1=r_1$ and $d_2=r_2$, where MARCF achieves the lowest estimation error. This trend empirically validates the efficiency gains obtained by explicitly modeling the common subspace structures.

figure[figure omitted — 752 chars of source]

Experiment II: Identification of Factor-Driven Dynamics

In the second experiment, we evaluate whether the MARCF model can correctly identify a pure factor-driven structure when the data are generated from the dynamic MFM, as defined in (ref) and (ref). The goal is to assess the proportion of trials in which the modeling precedure correctly identifies the dynamic MFM structure, and to compare the estimation errors of the factor loading projection matrices between MARCF and MFM in wang2019factor.

We set $p_1=p_2=16$, $T=800$, and consider cases where $r_1=d_1\in\{1,2,3,4\}$ and $r_2=d_2\in\{1,2,3,4\}$. The data are generated as follows:

enumerate• For $i=1,2$, generate $\mathbf \Lambda_i\in \mathbb{R}^{p_i\times r_i}$ with random orthonormal columns as factor loading matrices. Generate $\mathbf B_i\in\mathbb{R}^{r_i\times r_i}$ using the same procedure as for $\mathbf D_i$ in Experiment I. Ensure stationarity by rejecting realizations where $\rho(\mathbf B_1)\cdot\rho(\mathbf B_2)\geq 1$. • Generate two independent white noise processes $\mathbf E_t\in\mathbb{R}^{p_1\times p_2}$ and $\boldsymbol \xi_t\in\mathbb{R}^{r_1 \times r_2}$ with $\mathrm{vec}(\mathbf E_t)\stackrel{\text{i.i.d.}}{\sim} N_{p_1p_2}(\mathbf 0,\mathbf I)$ and $\mathrm{vec}(\boldsymbol \xi_t)\stackrel{\text{i.i.d.}}{\sim} N_{r_1r_2}(\mathbf 0,\mathbf I)$. Construct $\mathbf W_t=\mathbf E_t-\mathbf\Lambda_1\mathbf\Lambda_1^\top\mathbf E_t\mathbf \Lambda_2\mathbf \Lambda_2^\top$ to eliminate the moving average term in the composite noise of (ref). • Generate the factors as $\mathbf F_t=\mathbf B_1\mathbf F_{t-1}\mathbf B_2^\top+\boldsymbol\xi_t$ and the observation as $\mathbf Y_t=\mathbf\Lambda_1\mathbf F_t\mathbf\Lambda_2^\top+\mathbf W_t$.

The training procedure and tuning parameters are identical to those in Experiment I, except that the step size $\eta$ is initially set to 0.01 and reduced to 0.001 if Algorithm (ref) diverges. For the MFM, we use the TIPUP method from Chen2022TensorFM with true ranks.

Table (ref) presents the proportion of trials in which the MFM pattern is correctly identified, i.e., $\widehat{r}_1=\widehat{d}_1=r_1$ and $\widehat{r}_2=\widehat{d}_2=r_2$. The results demonstrate that the MARCF model successfully identifies the MFM with high probability, particularly when multiple factors are present. Even for $r_1=r_2=1$, the success rate is 84.8%.

table[table omitted — 570 chars of source]

Figure (ref) persents the log estimation errors of the proposed model and MFM. Only trials with correctly identified parameters are included. The error is defined as $$ \mathrm{Log ~estimation ~error}=\log\left(\left\| \widehat{\mathbf \Lambda}_i(\widehat{\mathbf \Lambda}_i^\top\widehat{\mathbf \Lambda}_i)^{-1}\widehat{\mathbf \Lambda}_i^\top-\mathbf \Lambda_i^*\mathbf \Lambda_i^{*\top} \right\|_{\mathrm{F}}\right),\quad \mathrm{for~} i=1,2. $$ Figures (ref) and (ref) show the estimation errors for the column spaces $\mathcal{M}(\boldsymbol\Lambda_1^*)$ and $\mathcal{M}(\boldsymbol\Lambda_2^*)$, respectively. In each subplot, the left (blue) boxplot represents the MFM, and the right (yellow) boxplot represents the MARCF model. For $\mathcal{M}(\mathbf{\Lambda}_1^*)$, the MARCF model shows slightly smaller errors when $r_2=1$, and significantly smaller errors when $r_2\geq2$, especially for $r_2\geq3$. Similar trends hold for $\mathcal{M}(\mathbf{\Lambda}_2^*)$. Overall, MARCF performs comparably to or better than MFM, demonstrating its effectiveness and flexibility.

figure[figure omitted — 829 chars of source]

Real Example: Quarterly Macroeconomic Data

To demonstrate the empirical utility of the MARCF framework, we apply it to a multinational macroeconomic dataset. This dataset, previously analyzed by Chen24ConstaintMFM, comprises 10 key economic indicators observed across 14 countries over 107 quarters, spanning from 1990-Q2 to 2016-Q4. The data is available in the supplementary materials for reproducibility.

The economic indicators include the Consumer Price Index (CPI) for food (CPIF), energy (CPIE), and total items (CPIT); long-term interest rates (government bond yields, IRLT) and 3-month interbank rates (IR3); total industrial production excluding construction (PTEC) and total manufacturing production (PTM); Gross Domestic Product (GDP); and the total value of exports (ITEX) and imports (ITEM). The 14 countries are Australia (AUS), Austria (AUT), Canada (CAN), Denmark (DNK), Finland (FIN), France (FRA), Germany (DEU), Ireland (IRL), the Netherlands (NLD), Norway (NOR), New Zealand (NZL), Sweden (SWE), the United Kingdom (GBR), and the United States (USA). Prior to analysis, all series were transformed to induce stationarity (via logarithmic differences), seasonally adjusted, and standardized to zero mean and unit variance.

Model Selection and Estimation

We determine the model structure based on rolling window validation. Based on the predictive performance over a 16-quarter validation set, we selected the total ranks $r_1 = 4$ and $r_2 = 5$. Subsequently, the common subspace dimensions were selected via BIC as described in Section (ref), yielding $d_1 = 3$ and $d_2 = 2$. These non-zero intersections ($d_i > 0$) and non-zero specific dimensions ($r_i > d_i$) provide strong empirical support that the data structure is hybrid: neither the purely disjoint subspaces of RRMAR nor the perfectly aligned subspaces of DMF adequately capture the complexity of global economic dynamics. For interpretation, the estimated basis matrices $[\widehat{\mathbf C}_i~\widehat{\mathbf R}_i]$ and $[\widehat{\mathbf C}_i~\widehat{\mathbf P}_i]$ are orthogonalized via QR decomposition to enforce strict orthonormality.

Interpretation of Global Economic Dynamics

We dissect the estimated dynamics by analyzing the projection matrices associated with the row (country) and column (indicator) subspaces.

\paragraph{Country-Level Dynamics (Row Space).} Figure (ref) displays the estimated projection matrices for the country dimension. The diagonal elements of these matrices reveal four distinct clusters with specific roles in the global economy: (1) an Anglo-Pacific group (Australia and UK); (2) a Core Economic Driver group (USA, Germany, Netherlands); (3) a Periphery group (Ireland, Finland); and (4) a general group including Canada and other EU members.

figure[figure omitted — 845 chars of source]

The decomposition of these groups across the subspaces offers significant economic insight:

itemize• Common Influencers: Australia and the UK (Figure (ref)) exhibit strong loadings in the common subspace. This implies a dual role: their economic conditions are highly predictive of global trends, and simultaneously, they are strongly predicted by the global state. This reflects their status as deeply integrated, open economies with significant financial sectors. • Global Drivers: The USA, Germany, and the Netherlands feature prominently in the predictor-specific subspace (Figure (ref)) but are less visible in the response-specific subspace. This asymmetry suggests a causal directionality: shocks originating in these major economies drive future global dynamics, but their own fluctuations are less determined by external factors. This aligns with the US's role as a global policy setter and Germany/Netherlands as the industrial and logistical engines of Europe. • Global Responders: In contrast, Ireland and Finland load heavily on the response-specific subspace (Figure (ref)). This characterizes them as small open economies that absorb global shocks—meaning their economic state is predicted by global factors—but their idiosyncratic movements have limited predictive power for the rest of the world.

Figure (ref) combines these components to visualize the total predictive and response footprints. The predictor footprint (Figure (ref)) is dominated by the largest economies (USA, Germany, UK, Australia), confirming that global co-movements are driven by this core group. The correlation structure (off-diagonal elements) reveals positive associations within trade-dependent clusters (Australia/NZL/Netherlands) and cultural/economic blocks (USA/UK). Interestingly, negative correlations appear between the USA and continental European powers (Germany/France), potentially reflecting divergent economic cycles or structural differences between Anglo-Saxon and Continental economic models.

figure[figure omitted — 728 chars of source]

\paragraph{Indicator-Level Dynamics (Column Space).} The analysis of economic indicators (Figures (ref) and (ref)) reveals clear lead-lag relationships.

itemize• Predictors: Real economic activity measures, including GDP, Manufacturing (PTM), Industrial Production (PTEC), and Exports (ITEX), dominate the predictor-specific subspace (Figure (ref)). This confirms their utility as leading indicators for future economic states. • Responders: Energy prices (CPIE) are dominant in the response-specific subspace (Figure (ref)) but negligible in the predictor space. This is consistent with the status of most countries in the sample as energy importers; their domestic energy prices are determined by global markets (the response) but do not themselves drive global cycles. • Commonality: Food prices (CPIF) span both subspaces, suggesting high persistence and integration into both the cost of living and production chains.
figure[figure omitted — 853 chars of source]
figure[figure omitted — 735 chars of source]

Forecasting Performance

Finally, we evaluate the out-of-sample forecasting performance. We compute one-step-ahead forecast errors using 15 rolling windows covering the period from 2013-Q1 to 2016-Q3. We compare MARCF against the RRMAR model (with identical total ranks $r_1=4, r_2=5$) and the dynamic MFM under three specifications: (1) rank selected by information criterion ($r_1=r_2=1$) Han22rkInTSFM; (2) ranks matching our common dimensions ($r_1=3, r_2=2$); and (3) a high-rank specification ($r_1=r_2=8$). For the dynamic MFM, after estimating the factor loading and process, we use a low-dimensional MAR model (ref) to fit the process for prediction.

Table (ref) reports the mean and median sums of squared out-of-sample prediction errors. The MARCF model outperforms all benchmarks. Notably, it improves upon RRMAR by exploiting the common subspace structure to reduce estimation variance, and it significantly outperforms dynamic MFM, which is constrained by the rigid assumption of symmetric predictor/response subspaces. In terms of model fit, the mean in-sample prediction errors of MFM(8,8) and RRMAR are slightly lower (239.90 and 242.31, respectively) compared to MARCF (243.25). Their inferior out-of-sample performance indicates overfitting. In contrast, MFM(1,1) and MFM(3,2) suffer from underfitting, yielding higher errors (334.30 and 276.32, respectively). The above empirical evidence supports the theoretical premise that real-world economic data occupies a middle ground, possessing both shared factors and distinct rotational dynamics, that is best captured by the hybrid MARCF framework.

table[table omitted — 464 chars of source]

Conclusion

In this paper, we have introduced the Matrix Autoregressive model with Common Factors (MARCF), a unified framework for high-dimensional matrix time series that bridges the structural gap between reduced-rank autoregression and dynamic factor models. By explicitly parameterizing the intersection between predictor and response subspaces, the MARCF framework resolves the tension between the structural rigidity of reduced-rank models (which assume distinct subspaces) and the restrictiveness of factor models (which assume identical subspaces).

Our contributions can be summarized as follows:

itemize• Methodological Unification: We proposed a flexible decomposition that separates dynamics into common (shared), predictor-specific, and response-specific components. This allows the model to adaptively capture the “subspace alignment” of the data, providing a continuum that encompasses both RRMAR and DMF as special cases while enabling the modeling of complex interaction dynamics. • Scalable Estimation: We developed a regularized gradient descent algorithm designed for the non-convex nature of the parameter space. By incorporating balanced regularization terms ($b \asymp \phi^{1/3}$), our approach resolves the scaling indeterminacies inherent in bilinear systems and avoids the computationally expensive matrix operations associated with alternating least squares, ensuring scalability to large-scale economic datasets. • Theoretical Rigor: We established the statistical consistency of the estimator and the rank selection procedure. Our analysis decouples the geometric landscape of the loss function from the stochastic properties of the data, proving that the estimator achieves the optimal convergence rate governed by the effective degrees of freedom.

The empirical application to global macroeconomic indicators highlights the practical value of this framework. The MARCF model not only improved forecasting accuracy but also provided interpretable insights into the global economy, distinguishing between “driver” countries (dominating the predictor subspace) and “responder” countries (dominating the response subspace).

Several promising directions remain for future research. First, while we focused on dense subspace structures, incorporating sparsity constraints on the basis matrices ($\mathbf{C}, \mathbf{R}, \mathbf{P}$) could further enhance interpretability and efficiency, particularly when $p_1$ and $p_2$ are very large. Second, extending this framework to tensor-valued time series with an order $K \ge 3$ would allow for the modeling of more complex structures, though this introduces additional algebraic challenges regarding tensor rank selection. Finally, developing inferential tools, such as confidence intervals for the estimated factors or impulse response functions, remains an important challenge due to the non-standard asymptotics of singular subspace identification.

{4pt}

\counterwithin{equation}{section} \counterwithin{lemma}{section} \counterwithin{figure}{section} \AtAppendix{\counterwithin{lemma}{section}} \AtAppendix{\counterwithin{proposition}{section}}