Modewise Additive Factor Model for Matrix Time Series

Modewise Additive Factor Model for Matrix Time Series

Introduction

Factor analysis is a fundamental tool for understanding common dependence among high-dimensional observations. Classical vector factor models provide a parsimonious representation of co-movements through a small number of latent factors, treating each observation as a long vector driven by shared components. These models have been rigorously developed for vector-valued time series in statistics, economics and finance, with seminal works connor1986performance, bai2002determining, bai2003inferential, forni2004generalized, forni2005generalized,lam2011estimation, lam2012factor,fan2016projected establishing consistency and inferential theory as both sample size and dimension grow. However, when observations possess natural matrix or tensor structure, such as panels of units and variables, spatiotemporal measurements, or import-export volumes across product categories, vector factor models fail to capture important multi-way dependencies wang2019factor,chen2023statistical. Vectorizing a matrix or tensor severs inherent row and column relationships, discarding structural information and often yielding less efficient and less interpretable estimates.

This limitation has motivated substantial recent interest in matrix- and tensor-valued factor models. By preserving multi-way structure, these models capture within- and across-mode dependencies that would be lost through flattening, often yielding more interpretable factors and improved estimation efficiency. A prominent class of such models employs a multiplicative architecture, where a latent factor tensor interacts with loading matrices along each mode. For a matrix time series $\{\bfm X_t\}_{t=1}^n$, a prototypical example is the Tucker factor model wang2019factor:

where $\bfm F_t\in\mathbb{R}^{r_1 \times r_2}$ is a matrix of latent factors and $\bfm A\in\mathbb{R}^{d_1\times r_1},\bfm B\in\mathbb{R}^{d_2\times r_2}$ are row and column loading matrices. Entry-wise, $X_{t,ij} = \bfm a_i^\top \bfm F_t \bfm b_j + E_{t,ij}$, so the latent factors act globally across all rows and columns. Recent research has extended this model and its estimation to more refined methodologies, including iterative refinement, formal inference, and more realistic matrix time-series settings such as time-varying and threshold models chen2020constrained,chen2023statistical,yu2022projected,yu2024dynamic,chen2025factor,chen2024time,liu2022identification. Further, this Tucker factor model framework extends naturally to tensor time series chen2022factor, han2022rank,han2024tensor, babii2025tensor, chen2024semi,chen2024rank, while a related but more constrained variant is the CP factor model chang2023modelling, chang2024identification, chen2026estimation, han2024cp,chen2025diffusion, bolivar2025threshold, which assumes a superdiagonal latent factor tensor.

Despite their advantages over vector factor models, these multiplicative frameworks share a fundamental structural limitation: co-movement across the entire data array is driven by a single global factor $\bfm F_t$. Row-specific and column-specific effects are not independent sources of variation but rather different projections of the same underlying factors. For instance, consider a retail chain analyzing sales data across stores (rows) and product categories (columns). Tucker and CP factor models implicitly assume that store-level dynamics (e.g., local competition, management quality) and product-level dynamics (e.g., supply chain disruptions, marketing campaigns) are manifestations of the same core factors, an assumption that may be overly restrictive in practice.

To overcome this limitation, lam2025matrix introduced a main-effect Tucker factor model with spectral decomposition-based estimation, while yuan2023twoway proposed a two-way dynamic factor model using quasi-maximum likelihood estimation. In this paper, we introduce and analyze a Modewise Additive Factor Model (MAFM) for matrix time series that separately captures row-specific and column-specific latent effects. Each observation $\bfm X_t \in \mathbb{R}^{d_1 \times d_2}$ decomposes into a row-factor component, a column-factor component, and noise:

where $\bfm F_t \in \mathbb{R}^{d_1 \times r_1}$ and $\bfm G_t \in \mathbb{R}^{d_2 \times r_2}$ are latent factor matrices for rows and columns, $\bfm A \in \mathbb{R}^{d_2 \times r_1}$ and $\bfm B \in \mathbb{R}^{d_1 \times r_2}$ are the corresponding loading matrices, and $\bfm E_t \in \mathbb{R}^{d_1 \times d_2}$ is a noise matrix. Write $\bfm F_t = (\mathbf{f}_{t,1},\mathbf{f}_{t,2},\dots,\mathbf{f}_{t,d_1})^\top$ with $\mathbf{f}_{t,i}\in\mathbb{R}^{r_1}$ being the latent factor for row $i$, and $\bfm G_t = (\bfm g_{t,1},\bfm g_{t,2},\dots,\bfm g_{t,d_2})^\top$ with $\bfm g_{t,j}\in\mathbb{R}^{r_2}$ being the latent factor for column $j$. Each entry of $\bfm X_t$ decomposes as

where the first term captures row-specific effects and the second captures column-specific effects. By disentangling these mode-specific sources of variation, our model offers greater flexibility than multiplicative frameworks.

To illustrate, revisit the retail example. Let $\bfm X_t$ contain sales figures with rows indexing stores and columns indexing product categories. The factor $\mathbf{f}_{t, i}$ captures store-specific drivers (e.g., location, management) at time $t$, while $\bfm g_{t,j}$ captures product-specific drivers (e.g., seasonality, branding). The loading $\bfm a_j \in \mathbb{R}^{r_1}$ reflects how store factors influence product $j$, and $\bfm b_i \in \mathbb{R}^{r_2}$ reflects how product factors influence store $i$. The additive structure $\mathbf{f}_{t,i}^\top \bfm a_j + \bfm b_i^\top \bfm g_{t,j}$ naturally separates these distinct sources of variation, offering a more realistic representation when store-level and product-level dynamics are driven by fundamentally different factors.

In this paper, we develop a computationally efficient two-stage procedure for MAFM. The first stage, Modewise Inner-product Eigendecomposition (MINE), computes eigendecompositions of the modewise sample covariance matrices of $\bfm X_t$ to obtain initial estimates. While the signal from one mode induces structural bias when estimating the other, this bias has weaker signal strength than the target component and can be controlled. The second stage, COMplement-Projected Alternating Subspace Estimation (COMPAS), iteratively refines these estimates through orthogonal complement projections, a key methodological innovation. Unlike the orthogonal projections in Tucker and CP factor models han2024tensor, han2024cp, our orthogonal complement projections completely eliminate cross-modal interference: when estimating the loading matrix for one mode, projecting onto the orthogonal complement perfectly removes the signal from the other mode. Our theoretical analysis provides details of these improvements.

In the theoretical analysis, we establish statistical upper bounds on the estimation errors of the factor loading matrices for the proposed algorithms. MINE yields good initial estimates despite bias from cross-modal contamination. COMPAS eliminates this bias through iterative refinement. Notably, unlike orthogonal projection-based algorithms in Tucker and CP factor models, COMPAS does not amplify the signal-to-noise ratio through projection; nevertheless, it achieves comparable convergence rates (see Remark (ref)). We also prove that using only a partial complement is statistically suboptimal, justifying the use of full orthogonal complement projections in COMPAS.

Beyond estimation, we establish asymptotic normality of the estimated factor loading matrices and develop consistent estimators for variance-covariance components under proper conditions, yielding a data-driven statistical inference framework which is an aspect largely underexplored in the factor model literature. These results enable valid confidence interval construction and hypothesis testing, tools that are essential for empirical applications. Our inferential approach, based on spectral representations of the estimators, differs fundamentally from techniques used in vector factor models bai2003inferential and Tucker factor models chen2023statistical, yu2022projected. Our approach can be extended naturally to other matrix and tensor factor model problems, such as inference of factor loading matrices using lagged autocovariance matrices han2024tensor, chang2023modelling. Additionally, we derive sharp tail probability inequalities for quadratic forms of matrix time series under exponential-type tail conditions and temporal dependence, which may be of independent interest for the analysis of matrix and tensor time series models.

Notation

We use lowercase letters (e.g., $a, b, c$) for scalars, bold lowercase letters (e.g., $\bfm a, \bfm b$) for vectors, and bold uppercase letters (e.g., $\bfm A, \bfm B$) for matrices. For two sequences of real numbers $\{a_n\}$ and $\{b_n\}$, write $a_n = O(b_n)$ if $|a_n| \le C|b_n|$ for some constant $C > 0$ and all sufficiently large $n$, and $a_n = o(b_n)$ if $a_n/b_n \to 0$. Write $a_n \lesssim b_n$ (resp. $a_n \gtrsim b_n$) if $a_n \le Cb_n$ (resp. $a_n \ge Cb_n$) for some $C > 0$, and $a_n \asymp b_n$ if both $a_n \lesssim b_n$ and $a_n \gtrsim b_n$ hold. The notation $O_{\mathbb{P}}(\cdot)$ and $o_{\mathbb{P}}(\cdot)$ is used in their standard probabilistic senses. Write $a \wedge b = \min\{a, b\}$ and $a \vee b = \max\{a, b\}$.

For a matrix $\bfm M$, $\Vec(\bfm M)$ denotes the vectorization obtained by stacking columns. The $p \times p$ identity matrix is $\bfm I_p$. The Stiefel manifold $\mathbb{O}^{p \times r} = \{\bfm U \in \mathbb{R}^{p \times r} : \bfm U^{\top}\bfm U = \bfm I_{r}\}$ is the set of $p\times r$ matrices with orthonormal columns. For any $\bfm U \in \mathbb{O}^{p\times r}$, its orthogonal complement is denoted $\bfm U_\perp \in \mathbb{O}^{p \times (p-r)}$. For symmetric matrices $\bfm A, \bfm B \in \mathbb{R}^{p \times p}$, write $\bfm A \preceq \bfm B$ (resp. $\bfm A \prec \bfm B$) if $\bfm B - \bfm A$ is positive semidefinite (resp. positive definite), and similarly for $\succeq$ and $\succ$.

For a matrix $\bfm A \in \mathbb{R}^{p_1 \times p_2}$ of rank $r$, let its singular value decomposition (SVD) be $\bfm A = \bfm U_1 \ensuremath{\boldsymbol{\Lambda}} \bfm U_2^{\top}$, where $\bfm U_1 \in \mathbb{O}^{p_1 \times r}$ and $\bfm U_2 \in \mathbb{O}^{p_2 \times r}$ are the left and right singular matrices, and $\ensuremath{\boldsymbol{\Lambda}} = {\rm diag}(\sigma_1(\bfm A), \ldots, \sigma_r(\bfm A))$ contains the singular values in descending order: $\sigma_1(\bfm A) \ge \cdots \ge \sigma_r(\bfm A) > 0$. We write $\sigma_{\max}(\bfm A) = \sigma_1(\bfm A)$ and $\sigma_{\min}(\bfm A) = \sigma_r(\bfm A)$. The projection matrices onto the left and right singular subspaces are $\bfm P_{\bfm A} = \bfm U_1\bfm U_1^{\top}$ and $\bfm P_{\bfm A^{\top}} = \bfm U_2\bfm U_2^{\top}$. The spectral and Frobenius norms are denoted as $\|\bfm A\|_2=\|\bfm A\|=\sigma_1(\bfm A)$ and $\|\bfm A\|_{\mathrm{F}}$, respectively. For two orthonormal matrices $\bfm U, \widehat \bfm U \in \mathbb{O}^{p \times r}$, let $\sigma_1\ge \sigma_2 \ge \cdots \ge \sigma_r\ge 0$ be the singular values of $\bfm U^\top \widehat \bfm U$. A natural measure of distance between their column spaces is $\|\widehat \bfm U\widehat \bfm U^\top - \bfm U\bfm U^\top\|_2=\sqrt{1-\sigma_r^2},$ which equals the sine of the largest principal angle between the column spaces.

Organization

The remainder of the paper is organized as follows. Section (ref) discusses identifiability of the modewise additive factor model. Section (ref) presents the two-stage estimation procedure: MINE initialization and COMPAS refinement. Section (ref) establishes statistical convergence rates of our estimators, asymptotic normality of factor loading matrices, and consistent estimation of variance-covariance components. Section (ref) develops sharp Bernstein-type inequalities for matrix time series. Section (ref) reports simulation results examining estimation accuracy and inferential validity under various settings. Section (ref) applies the proposed method to real-world matrix time series data. Section (ref) concludes the paper. Additional simulations, all proofs and technical lemmas are deferred to the Appendix.

Identifiability of Modewise Additive Factor Model

Identifiability is a foundational prerequisite for the estimation and interpretation of factor models. Like existing vector and Tucker factor models, our Modewise Additive Factor Model (MAFM) in (ref) is invariant under the transformation $(\bfm A\bfm H_{\bfm A}, \bfm B\bfm H_{\bfm B}, \bfm F_t\bfm H_{\bfm A}^{-1}, \bfm G_t\bfm H_{\bfm B}^{-1})$ for any invertible matrices $\bfm H_{\bfm A} \in \mathbb{R}^{r_1 \times r_1}$ and $\bfm H_{\bfm B} \in \mathbb{R}^{r_2 \times r_2}$. Nevertheless, the column spaces of $\bfm A$ and $\bfm B$, i.e., the factor loading spaces, are uniquely defined. The primary objective of this work is to identify and consistently estimate these subspaces.

To resolve the rotational and scaling ambiguities, we introduce a canonical representation of the model. Let $\bfm A = \bfm U_{\bfm A}\ensuremath{\boldsymbol{\Lambda}}_{\bfm A}\bfm W_{\bfm A}^\top$ and $\bfm B = \bfm U_{\bfm B}\ensuremath{\boldsymbol{\Lambda}}_{\bfm B}\bfm W_{\bfm B}^\top$ denote the SVD of the loading matrices. Absorbing the right singular matrices and singular values into the factors yields the canonical form:

where $\cano F_t = \bfm F_t \bfm W_{\bfm A}\ensuremath{\boldsymbol{\Lambda}}_{\bfm A} \in \mathbb{R}^{d_1 \times r_1}$ and $\cano G_t = \bfm G_t \bfm W_{\bfm B}\ensuremath{\boldsymbol{\Lambda}}_{\bfm B} \in \mathbb{R}^{d_2 \times r_2}$. This representation forms the basis of our estimation procedure.

To ensure that the loading spaces spanned by $\bfm U_{\bfm A}$ and $\bfm U_{\bfm B}$ are uniquely recoverable from the observed data $\{\bfm X_t\}_{t=1}^n$, we impose the following mild structural assumptions:

(C1) The factor processes $\{\bfm F_t\}_{t=1}^n$ and $\{\bfm G_t\}_{t=1}^n$ are weakly stationary with mean zero and mutually independent.

(C2) For any nonzero vectors $\bfm u \in \mathbb{R}^{d_1}$ and $\bfm v \in \mathbb{R}^{d_2}$, the covariance matrices of $\bfm A\bfm F_t^\top\bfm u$ (row-factor part) and $\bfm B\bfm G_t^\top\bfm v$ (column-factor part) have ranks $r_1$ and $r_2$, respectively.

Condition (C1) ensures that the row and column factors are uncorrelated with mean zero, allowing their contributions to the signal covariance to be separated. The weak stationarity further guarantees a consistent covariance structure across all time points, providing a well-defined population target for estimation. Specifically, the covariance tensor of the signal component $\bfm M_t$ admits the additive decomposition

where $\otimes$ denotes the tensor outer product, and the mode-$k$ product ${\cal X}\times_k \bfm A$ between $\mathcal{X}\in\mathbb{R}^{d_1\times\cdots\times d_M}$ and $\bfm A\in\mathbb{R}^{\widetilde d_k\times d_k}$ produces a tensor of size $d_1\times\cdots\times d_{k-1}\times\widetilde d_k\times d_{k+1}\times\cdots\times d_M$ with entries $({\cal X}\times_k \bfm A)_{i_1\cdots i_{k-1} j i_{k+1}\cdots i_M} = \sum_{i_k=1}^{d_k}\mathcal{X}_{i_1\cdots i_k\cdots i_M}A_{j i_k}$. This decomposition, a direct consequence of the additive model structure and the uncorrelatedness in (C1), is fundamental to separating the two factor loading spaces. Condition (C2) ensures that $r_1$ and $r_2$ are the intrinsic dimensions of the row-factor and column-factor components, which is crucial for identifiability.

Proposition (ref) establishes the identifiability of the factor loading spaces.

Estimation

In this section, we develop estimators for the factor loading matrices $\bfm U_{\bfm A}$ and $\bfm U_{\bfm B}$. The additive structure of the model poses a unique statistical challenge: when estimating the loading matrix $\bfm U_{\bfm A}$, the column-factor signal $\bfm B\bfm G_t^\top$ acts as structured contamination in $\bfm X_t^\top \bfm X_t$, and vice versa. This cross-modal interference can introduce substantial bias into standard spectral decomposition based estimators. To address this issue, we propose a computationally efficient two-stage spectral procedure that decouples the row and column factors. Our approach first extracts initial estimates and then iteratively refines them through orthogonal complement projections.

\noindentInitialization. We construct initial estimators for the factor loading matrices via a two-step PCA. Specifically,

where Eigen$_m$ stands for the matrix composed of the top $m$ eigenvectors corresponding to the largest $m$ eigenvalues. This procedure, detailed in Algorithm (ref) and named Modewise INner-product Eigendecomposition (MINE), provides initializations for our iterative refinement.

The rationale follows from the modewise covariance structure of the observations. Assuming $\cano\bfm F_t,\cano\bfm G_t,\bfm E_t$ are mutually independent, the population column-wise covariance of $\bfm X_t$ is

The dominant low-rank component is the first term, $\bfm U_{\bfm A}\mathbb{E}[\cano\bfm F_t^\top \cano\bfm F_t]\bfm U_{\bfm A}^\top$, which has rank $r_1$. Although the second term $\mathbb{E}[\cano\bfm G_t\bfm U_{\bfm B}^{\top} \bfm U_{\bfm B}\cano\bfm G_t^\top]$ also contributes spiked eigenvalues, these grow at a slower rate. Consequently, the top $r_1$ eigenvectors of this modewise covariance primarily capture the loading matrix $\bfm U_{\bfm A}$. A symmetric argument applies to the loading matrix $\bfm U_{\bfm B}$.

\noindentIterative Refinement. After obtaining initial estimates via MINE (Algorithm (ref)), we apply an orthogonal complement projection algorithm to refine them. The key idea is to isolate the signal from one mode by annihilating the signal from the other. It is motivated by the following observation. Let $\bfm U_{\bfm A\perp},\bfm U_{\bfm B\perp}$ be the orthogonal complements of $\bfm U_{\bfm A},\bfm U_{\bfm B}$, respectively. Multiplying $\bfm X_t$ by the orthogonal complement $\bfm U_{\bfm A\perp}$, model (ref) yields

The orthogonal complement projection eliminates the row-factor component $\widetilde{\bfm F}_t \bfm U_{\bfm A}^\top$, while preserving the column-factor component. Under proper conditions on the projected noise matrix $\bfm E_t \bfm U_{\bfm A\perp}$, estimation of the loading matrix $\bfm U_{\bfm B}$ based on $\bfm X_t \bfm U_{\bfm A\perp}$ can be made significantly more accurate. A symmetric argument applies to estimating $\bfm U_{\bfm A}$ via multiplying $\bfm X_t^\top$ by $\bfm U_{\bfm B\perp}$. We note that since $r_1\ll d_2 $ and $ r_2\ll d_1$, unlike orthogonal projections in Tucker and CP factor models han2024tensor,han2024cp, the orthogonal complement projection does not reduce the noise variation.

In practice, we do not know $\bfm U_{\bfm A\perp},\bfm U_{\bfm B\perp}$. Analogous to back-fitting algorithms, we iteratively estimate the loading matrices $\bfm U_{\bfm B}$ and $\bfm U_{\bfm A}$ at iteration $i$ based on $\bfm X_t \widehat\bfm U_{\bfm A\perp}^{(i-1)} ,\bfm X_t^\top \widehat\bfm U_{\bfm B\perp}^{(i)}$, using the estimate $\widehat\bfm U_{\bfm A\perp}^{(i-1)}$ obtained in the previous iteration and the estimate $\widehat\bfm U_{\bfm B\perp}^{(i)}$ obtained in the current iteration. As we shall show in the next section, such an iterative procedure leads to a much improved statistical rate in the high dimensional matrix factor model scenarios, as if $\bfm U_{\bfm A\perp},\bfm U_{\bfm B\perp}$ are known and we indeed observe $\bfm X_t \bfm U_{\bfm A\perp},\bfm X_t^\top \bfm U_{\bfm B\perp}$ that follows model (ref).

The complete procedure is presented in Algorithm (ref), named COMplement-Projected Alternating Subspace Estimation (COMPAS). Although projection onto any subset of columns from $\bfm U_{\bfm A\perp}$ (or $\bfm U_{\bfm B\perp}$) eliminates the corresponding row-factor component (or column-factor component), our theoretical analysis (Theorem (ref)) establishes that projection onto the full orthogonal complement yields optimal estimation efficiency. Hence, the use of the complete orthogonal complement is a critical design choice.

\noindentFactor Estimation. Once the loading matrices have been estimated via COMPAS (Algorithm (ref)), the final step is to recover the latent factors. We construct estimators by projecting the observed data onto the estimated signal subspaces, thereby isolating the contribution of each mode. Let $\widehat{\bfm U}_{\bfm A}$ and $\widehat{\bfm U}_{\bfm B}$ be the final estimates of the loading matrices from COMPAS, with $\widehat{\bfm U}_{\bfm A\perp}$ and $\widehat{\bfm U}_{\bfm B\perp}$ being their respective orthogonal complements. The estimators for the canonical factors $\cano F_t$ and $\cano G_t$ are

or equivalently, $\widehat \bfm F_t = \bfm X_t \widehat \bfm U_{\bfm A}, \widehat \bfm G_t^\top = \widehat \bfm U_{\bfm B}^\top \bfm X_t \widehat \bfm U_{\bfm A \perp} \widehat \bfm U_{\bfm A \perp}^\top.$ With the estimated loading matrices and factors, we construct the fitted value $\widehat{\bfm X}_t$ for the observation at time $t$ by removing the component orthogonal to both factor loading spaces:

Theoretical Properties

In this section, we shall investigate the statistical properties of the proposed algorithms described in the last section. Our theories provide consistency guarantees and characterize statistical error rates for estimating the factor loading matrices $\bfm U_{\bfm A}$ and $\bfm U_{\bfm B}$, under proper regularity conditions. We also establish the asymptotic distributions for each row of the estimated factor loading matrices.

Convergence Rates for Factor Loading Estimation

To present theoretical properties of the proposed procedures, we impose the following assumptions.

Assumption (ref) parallels noise conditions in lam2011estimation,lam2012factor, han2024tensor, han2024cp, chen2022factor. It accommodates general dependence patterns among individual time series while maintaining analytical tractability. The Gaussian assumption, imposed for technical convenience, ensures fast convergence rates in our analysis. In principle, this could be relaxed to sub-Gaussian or heavier-tailed distributions. However, such generalizations would substantially complicate the formulas, conditions, and statistical outcomes in our time series framework without yielding additional insights. To maintain focus on the core contributions, we adopt the Gaussian error assumption, which simplifies exposition without compromising the fundamental methodology or findings.

Assumption (ref) accommodates a broad class of time series models, including causal ARMA processes with continuously distributed innovations tsay2010analysis,tsay2018nonlinear,fan2003nonlinear. For notational simplicity, we impose a common $\theta_1$ for both factor processes, though different values could be accommodated without altering the proofs. The independence between factor processes $\bfm F_t$ and $\bfm G_t$ is a structural assumption that ensures separation: row factors drive across-column variation, while column factors drive across-row variation, without overlap.

Assumption (ref) imposes key structural properties on the latent factors. Part (i) is technically important: the independence of rows within each factor matrix enables us to leverage concentration inequalities for independent exponential-tail random variables when analyzing the cross-sectional dimension of factor rows. Part (ii) is a standard condition ensuring that the eigenvalues of the covariance of each factor row are uniformly bounded above and below. As detailed in the remark below, this assumption translates directly into bounds on the eigenvalues of key population covariance matrices, thereby quantifying the signal strength of the factor processes.

Assumption (ref) requires that the tail probabilities of factor rows decay exponentially, with $\theta_2=2$ corresponding to sub-Gaussian distributions. Combining Assumptions (ref) and (ref), we develop novel matrix Bernstein type inequalities (e.g., Lemma (ref)) for temporally dependent random matrices with exponential tails, extending classical results for i.i.d. random matrices. These concentration bounds may be of independent interest.

Assumption (ref) is similar to the signal strength condition in lam2012factor,han2024tensor,chen2026estimation, and aligns with the pervasive condition on factor loadings stock2002forecasting, bai2003inferential. It plays a key role in identifying the common factors and idiosyncratic noises in (ref). The indices $\delta_0,\delta_1$ measure factor loading strength, or the rate of signal strength growth as the dimension $d_k$ increases: $\delta_0$ characterizes the strength of the strongest factors, while $\delta_1$ corresponds to the weakest. When $\delta_0=\delta_1=0$, the factors are called strong factors; otherwise, the factors are called weak factors. For notational simplicity, we adopt common values of $\delta_0$ and $\delta_1$ across both dimensions, $d_1$ and $d_2$, although different values can be accommodated without affecting the analysis. The gap $\delta_1-\delta_0$ further characterizes the heterogeneity in loading strengths.

Under these assumptions, we now derive non-asymptotic error bounds for our two-stage estimation procedure. The first result establishes the accuracy of the initial factor loading estimator $\widehat{\bfm U}_{\bfm A}^{(0)}$ and $\widehat{\bfm U}_{\bfm B}^{(0)}$ from Algorithm (ref).

The error bounds in Theorem (ref) comprise three distinct components. The first term corresponds to the average of the noise-by-noise cross-products. If all factors are strong factors with $\delta_0=\delta_1=0$, this term reduces to order $1/d_2$ or $1/d_1$, consistent with non-iterative Tucker factor models chen2023statistical. This rate illustrates the “blessing of dimensionality” barigozzi2023quasi: larger cross-sectional dimensions help average out noise. The second term arises from the average of the signal-by-noise cross-products. Under strong factors, the order of rates $(nd_1)^{-1/2}$ and $(nd_2)^{-1/2}$ also matches Tucker factor models han2024tensor. The third term captures cross-modal interference, reflecting the interaction between row and column factors as well as contamination from the irrelevant factor mode.

Next, we analyze the statistical performance of the iterative refinement (Algorithm (ref)). The algorithm eliminates cross-modal interference by projecting onto the orthogonal complement of the estimated factor loading spaces. For instance, when estimating $\bfm U_{\bfm B}$, projection via $\bfm U_{\bfm A\perp}$ removes the row-factor component $\widetilde{\bfm F}_t \bfm U_{\bfm A}^\top$ in (ref). Ideally (i.e., when the true complement projection matrices $\bfm U_{\bfm A\perp},\bfm U_{\bfm B\perp}$ are used), this would eliminate the third terms in (ref) and (ref), reducing the error bounds to

In addition, define the contraction rate of each iteration,

The following theorem provides conditions under which the iterative refinement achieves the ideal rates.

The key insight in analyzing COMPAS (Algorithm (ref)) is that each iteration achieves error contraction, with contraction rates $\rho_1$ and $\rho_2$ for estimating $\bfm U_{\bfm A}$ and $\bfm U_{\bfm B}$, respectively. Unlike Tucker and CP factor models han2024tensor, han2024cp, COMPAS does not reduce noise variation; instead, it eliminates only the third term in the initialization bounds (ref) and (ref), which represents cross-modal interference. The contraction factors $\eta_{\bfm A}$ and $\eta_{\bfm B}$ govern the vanishing rate of this term, rather than affecting the oracle accuracy of the loading space estimators as in han2024tensor. This reflects the role of alternating orthogonal complement projection: each mode benefits from refined estimates of the other, so that the error $\| \widehat{\bfm U}_{\bfm A} \widehat{\bfm U}_{\bfm A}^\top - \bfm U_{\bfm A} \bfm U_{\bfm A}^\top \|$ includes the leading error term of $\widehat{\bfm U}_{\bfm B}$ scaled by $\eta_{\bfm A}$, and vice versa.

The first two error terms from initialization persist in the refined bound. The first term corresponds to the average of the noise-by-noise cross-products, while the second captures the average of the signal-by-noise cross-products. This decomposition clarifies how factor loading strength, noise variance, sample size, and dimensionality jointly determine estimation accuracy, paralleling the error structure in Tucker and CP factor models.

A key design choice in Algorithm (ref) is the use of the full orthogonal complement (e.g., $\bfm U_{\bfm A\perp} \in \mathbb{R}^{d_2 \times (d_2-r_1)}$) in the projection step. The matrix $\bfm U_{\bfm A\perp}$ spans all $d_2-r_1$ directions orthogonal to $\bfm U_{\bfm A} \in \mathbb{R}^{d_2 \times r_1}$, thereby leveraging the maximum information unrelated to $\bfm U_{\bfm A}$ when estimating $\bfm U_{\bfm B}$. One might consider a computationally cheaper variant that projects onto a smaller subspace of the complement. Specifically, suppose that we replace $\bfm U_{\bfm A\perp} \in \mathbb{O}^{d_2 \times (d_2-r_1)}$ and $\bfm U_{\bfm B\perp} \in \mathbb{O}^{d_1 \times (d_1-r_2)}$ with partial complements $\bfm V_{\bfm A} \in \mathbb{O}^{d_2 \times s_2}$ and $\bfm V_{\bfm B} \in \mathbb{O}^{d_1 \times s_1}$, respectively, where $s_2<d_2-r_1$ and $s_1<d_1-r_2$. Let $\cano \bfm V_{\bfm A} \in \mathbb{O}^{d_2 \times s_2}$ and $\cano \bfm V_{\bfm B} \in \mathbb{O}^{d_1 \times s_1}$ be the corresponding estimates of $\bfm V_{\bfm A}$ and $\bfm V_{\bfm B}$ from Algorithm (ref), and let $\cano \bfm U_{\bfm A}$ and $\cano \bfm U_{\bfm B}$ denote the resulting loading matrix estimates. Theorem (ref) formally confirms that using fewer than the full complement directions incurs a larger error bound, reflecting a loss of statistical efficiency.

Compared with Theorem (ref), which uses the full complement, Theorem (ref) shows that restricting to $s_1 < d_1-r_2$ and $s_2 < d_2-r_1$ directions inflates the signal-by-noise terms by replacing $\sqrt{n d_1}$ and $\sqrt{n d_2}$ with $\sqrt{n s_1}$ and $\sqrt{n s_2}$. This efficiency loss has a clear mechanistic explanation. The update for $\bfm U_{\bfm B}$ relies on the row-wise sample covariance of the projected data $\bfm X_t \bfm U_{\bfm A\perp}$. The signal strength in this projection, which is essential for separating signal from noise, is determined by the smallest nonzero eigenvalue of the signal covariance matrix. As established in Remark (ref), this eigenvalue is proportional to the rank of the projection matrix. Using the full complement projection ensures maximum signal variation for updating each mode, leading to tighter error contraction. In contrast, a partial complement of rank $s_2 < d_2-r_1$ discards signal information, inflating the final estimation error. While the algorithm remains valid with partial complement projections, the convergence rates and error bounds are strictly worse.

Statistical Inference for Factor Loadings

In factor models, each row of the loading matrix represents how a particular unit (e.g. a country or an asset) loads onto the latent factors. For example, in finance, asset pricing models posit that each stock’s return depends linearly on a set of common factors, with the factor loadings capturing the asset’s risk exposures or beta coefficients, to systematic market factors ROSS1976341, connor1986performance,lettau2022high,giglio2022factor,conlon2025asset. In macroeconomics, a loading row may measure a country's sensitivity to global business cycle shocks bai2003inferential,houssa2008monetary,forni2021policy. Inference on these row-specific loadings is crucial for testing hypotheses about unit heterogeneity (e.g., which countries significantly load on a global factor) or about the stability of the factor structure over time (e.g., whether certain assets cease to load on a market factor).

In this section, we establish the asymptotic distributions of the estimated factor loading matrices. Let

denote the error bounds for estimating loading matrices $\bfm U_{\bfm A}$ and $\bfm U_{\bfm B}$ after applying Algorithms (ref) and (ref), as established in (ref) and (ref). By Theorem (ref), $\| \widehat{\bfm U}_{\bfm A}\widehat{\bfm U}_{\bfm A}^{\top} - \bfm U_{\bfm A}\bfm U_{\bfm A}^{\top} \| \le \overline{\varepsilon}_{\bfm A}$ and $\| \widehat{\bfm U}_{\bfm B}\widehat{\bfm U}_{\bfm B}^{\top} - \bfm U_{\bfm B}\bfm U_{\bfm B}^{\top} \| \le \overline{\varepsilon}_{\bfm B}$ with probability at least $1- \exp(-cn^{\varphi}(d_1 \wedge d_2)^{\varphi/2}) - \exp(-c (d_1 \wedge d_2))$. Recall the noise covariance matrices $\bfsym \Sigma_{\bfm E} = \operatorname{Cov}\left(\Vec\left(\bfm E_t\right)\right)\in \mathbb{R}^{d_1d_2 \times d_1d_2}$ and $\bfsym \Sigma_{\bfm E^{\top}} = \operatorname{Cov}\left(\Vec\left(\bfm E_t^{\top}\right)\right) \in \mathbb{R}^{d_1d_2 \times d_1d_2}$. Let

be the population row-wise covariance matrices of the projected factors $\bfm U_{\bfm B\perp}\bfm U_{\bfm B\perp}^{\top}\cano \bfm F_t$ and $\bfm U_{\bfm A\perp}\bfm U_{\bfm A\perp}^{\top}\cano \bfm G_t $, respectively.

Theorem (ref) establishes the asymptotic normality of each row of the estimated factor loading matrices, with means aligned to the true values via rotations $\bfm R_{\bfm A},\bfm R_{\bfm B}$. The asymptotic variances involve the column-wise covariance matrices of the projected factors ($\bfsym \Sigma_{\bfm F}$ and $\bfsym \Sigma_{\bfm G}$), and the covariance matrices of the signal-by-noise cross-products ($\bfsym \Sigma_{\bfm A,i}$ and $\bfsym \Sigma_{\bfm B,j}$), analogous to statistical inference in conventional factor models. However, our inferential framework relies on a fine-grained analysis of the covariance matrix of the projected data. Specifically, we employ spectral representation tools from xia2021normal to achieve a second-order accurate characterization of matrix SVD, enabling precise control of the leading term and sharp bounds on higher-order terms. This differs fundamentally from techniques used in vector factor models bai2003inferential and Tucker factor models chen2023statistical, yu2022projected, which adopt the approach of bai2003inferential. Our framework offers two advantages: first, we establish explicit convergence rates for asymptotic normality, which previous methods cannot provide; second, the technique in bai2003inferential relies on a manually designed rotation matrix for error decomposition, which does not extend to our setting due to the modewise additive structure and high dimensionality of the factor processes. Our approach can be generalized to inference in other matrix and tensor factor models, such as CP factor models chang2023modelling,chang2024identification.

The convergence rates of asymptotic normality $\Omega_{\bfm A}$ and $\Omega_{\bfm B}$ comprise three terms. Taking $\Omega_{\bfm A}$ as an example: (i) the first term $r_1\sqrt{1/(nd_1)}d_2^{\delta_1-\delta_0}$ is the Berry Esseen type bound for signal-by-noise cross-product of the first-order term; (ii) the second term, $ d_2^{(\delta_1-\delta_0)/2} \sqrt{r_1\log(nd_1)} \overline{\varepsilon}_{\bfm B}$, arises from estimation error in $\bfm U_{\bfm B}$ propagating to the first-order term; (iii) the third term $d_2^{(\delta_1-\delta_0)/2} \cdot \sqrt{n} d_2^{(1-\delta_0)/2}d_1^{1/2} (\overline{\varepsilon}_{\bfm A}^2 + \overline{\varepsilon}_{\bfm A}\overline{\varepsilon}_{\bfm B} )$ reflects the noise-by-noise crossproducts of the first order term, while the normal approximation is derived based on the first-order signal-by-noise cross-products induced by the perturbation $\widehat{\bfm U}_{\bfm A} - \bfm U_{\bfm A}\bfm R_{\bfm A}$.

Data-Driven Statistical Inference

Although statistical inference for loading matrices has been studied in vector bai2003inferential, Tucker chen2023statistical, and CP chen2026estimation factor models, practical procedures for constructing confidence intervals remain underexplored, where covariance matrix estimation is crucial, particularly under weak factor settings. In this subsection, we develop a data-driven inference approach based on plug-in estimators. We first estimate the projected factor covariance matrices, which capture the row-wise covariances of the projected factors $\bfm U_{\bfm B\perp} \bfm U_{\bfm B\perp}^{\top}\cano \bfm F_t$ and $\bfm U_{\bfm A\perp} \bfm U_{\bfm A\perp}^{\top}\cano \bfm G_t$. Using the estimated factors in (ref), the plug-in estimators for the projected factor covariance matrices in (ref) are

A direct error metric $\|\widehat{\bfsym \Sigma}_{\bfm F} - \bfsym \Sigma_{\bfm F}\|$ and $\|\widehat{\bfsym \Sigma}_{\bfm G} - \bfsym \Sigma_{\bfm G}\|$ is uninformative, since the covariance matrices scale with dimensions, $\|\bfsym \Sigma_{\bfm F}\| = O_{\mathbb{P}}(d_1d_2^{1-\delta_0})$ and $\|\bfsym \Sigma_{\bfm G}\| = O_{\mathbb{P}}(d_2d_1^{1-\delta_0})$, and are subject to rotational ambiguity. Instead, we adopt standardized error metrics $\|\bfm R_{\bfm A}^{\top}\bfsym \Sigma_{\bfm F}^{-1/2}\bfm R_{\bfm A}\widehat{\bfsym \Sigma}_{\bfm F}\bfm R_{\bfm A}^{\top}\bfsym \Sigma_{\bfm F}^{-1/2}\bfm R_{\bfm A} - \bfm I_{r_1}\|$ and $\|\bfm R_{\bfm B}^{\top}\bfsym \Sigma_{\bfm G}^{-1/2}\bfm R_{\bfm B}\widehat{\bfsym \Sigma}_{\bfm G}\bfm R_{\bfm B}^{\top}\bfsym \Sigma_{\bfm G}^{-1/2}\bfm R_{\bfm B} - \bfm I_{r_2}\|$, which accounts for both the scale of $\cano \bfm F_t,\cano \bfm G_t$ and rotation.

Theorem (ref) shows that the convergence rates are controlled by the loading matrix estimation errors $\overline{\varepsilon}_{\bfm A},\overline{\varepsilon}_{\bfm B}$. It shows that accurate loading matrix estimation is sufficient for accurate estimation of the factor covariance matrices. The factors $d_1^{\delta_1-\delta_0}$ and $d_2^{\delta_1-\delta_0}$ indicate that heterogeneity in factor loading strengths (larger gap $\delta_1-\delta_0$) degrades the precision of the covariance estimate.

Next, we construct plug-in estimators for the covariance matrices $\bfsym \Sigma_{\bfm A,i}$ and $\bfsym \Sigma_{\bfm B,j}$. Define the noise covariance estimators

for $\bfsym \Sigma_{\bfm E}$ and $\bfsym \Sigma_{\bfm E^{\top}}$ (defined in Assumption (ref)). The plug-in estimators for $\bfsym \Sigma_{\bfm A,i}$ and $\bfsym \Sigma_{\bfm B,j}$ are

respectively. As in Theorem (ref), since $\|\bfsym \Sigma_{\bfm A,i}\| = O_{\mathbb{P}}(d_2^{1-\delta_0})$ and $\|\bfsym \Sigma_{\bfm B,j}\| = O_{\mathbb{P}}(d_1^{1-\delta_0})$ with rotational ambiguity, we adopt standardized error metrics $\|\bfm R_{\bfm A}^{\top}\bfsym \Sigma_{\bfm F}^{-1/2}\bfm R_{\bfm A}\widehat{\bfsym \Sigma}_{\bfm A,i}\bfm R_{\bfm A}^{\top}\bfsym \Sigma_{\bfm F}^{-1/2}\bfm R_{\bfm A} - \bfm I_{r_1}\|$ and $\|\bfm R_{\bfm B}^{\top}\bfsym \Sigma_{\bfm G}^{-1/2}\bfm R_{\bfm B}\widehat{\bfsym \Sigma}_{\bfm B,j}\bfm R_{\bfm B}^{\top}\bfsym \Sigma_{\bfm G}^{-1/2}\bfm R_{\bfm B} - \bfm I_{r_2}\|$.

Theorem (ref) shows that the convergence rates are governed by the loading matrix estimation errors and a factor $(d_1+d_2)^{1/\varphi}/n$ reflecting temporal dependence. The terms $d_1^{\delta_1-\delta_0}$ and $d_2^{\delta_1-\delta_0}$ capture the effect of heterogeneity in factor loading strengths.

Combining the covariance matrices estimators from Theorems (ref) and (ref), we establish a feasible data-driven inference procedure.

In Theorem (ref), the rotation matrices $\bfm R_{\bfm A}$ and $\bfm R_{\bfm B}$ naturally cancel in the plug-in procedure. Compared with Theorem (ref), the error bounds are inflated by factors $d_1^{(\delta_1-\delta_0)/2}$ and $d_2^{(\delta_1-\delta_0)/2}$ due to repeated use of the estimated factors $\widehat{\bfm F}_t$ and $\widehat{\bfm G}_t$, respectively. In addition, accurate loading matrix estimation $\widehat{\bfm U}_{\bfm A}$ and $\widehat{\bfm U}_{\bfm B}$ is required to counteract the factor strength contributions $d_2^{1-\delta_0}$ from $\bfm A$ and $d_1^{1-\delta_0}$ from $\bfm B$ when estimating $\bfsym \Sigma_{\bfm E}$ and $\bfsym \Sigma_{\bfm E^{\top}}$. The term $(d_1+d_2)^{1/\varphi}/n$ accounts for dependence between loading matrix estimation and factors/noises.

Matrix Bernstein Inequality for Quadratic Forms of Matrix Time Series

Tail probability inequalities play an important role in statistical estimation and inference. In this section, we derive powerful tail probability inequalities for quadratic forms of matrix time series under exponential-type tail conditions. The following lemma, proved in Appendix C.3, extends classical matrix Bernstein inequalities for i.i.d. random matrices to quadratic forms of dependent matrix time series and is of independent interest. We refer readers to Appendices C.2 and C.4 for tail probability inequalities involving Gaussian random matrices and matrix time series with exponential-type tails.

Simulation

In this section, we evaluate the empirical performance of the two-stage estimation procedure, MINE and COMPAS, and compare it with a computationally efficient variant of COMPAS that uses partial (subset) orthogonal complement projections, validating the empirical efficiency loss of this practical alternative. We also examine the asymptotic normality of rows of the estimated loading matrices and verify that data-driven statistical inference procedures achieve nominal coverage rates, demonstrating accurate uncertainty quantification in finite samples.

Data generation. The matrix time series $\{\bfm X_t\}_{t=1}^n \in \mathbb{R}^{d_1 \times d_2}$ are generated according to MAFM model (ref). Each row of the latent row-factor matrix $\mathbf{f}_{t,i} \in \mathbb{R}^{r_1}$ follows an independent VAR(1) process $\mathbf{f}_{t,i} = \ensuremath{\boldsymbol{\Phi}}_{F,i} \mathbf{f}_{t-1,i} + \boldsymbol{\varepsilon}_{t,i}$ with $\boldsymbol{\varepsilon}_{t,i} \stackrel{\text{i.i.d.}}{\sim} {\cal N}(0, \bfm I_{r_1})$, and similarly for each row of the latent column-factor matrix $\bfm g_{t,j} \in \mathbb{R}^{r_2}$. To induce heterogeneous temporal dynamics, each coefficient matrix $\ensuremath{\boldsymbol{\Phi}}_{F,i}$ (or $\ensuremath{\boldsymbol{\Phi}}_{G,j}$) is constructed as $\bfm Q_i \bfm D_i(\phi_1, \phi_2) \bfm Q_i^\top$, where $\bfm D_i(\phi_1, \phi_2) = {\rm diag}(\underbrace{\phi_1, \ldots, \phi_1}_{\lceil r/2\rceil}, \underbrace{\phi_2, \ldots, \phi_2}_{r-\lceil r/2\rceil})$ with $r = r_1$ for row factors and $r = r_2$ for column factors, and $\bfm Q_i$ is an orthogonal matrix from QR decomposition of a standard Gaussian random matrix. Each coordinate is uniformly assigned to one of three representative eigenvalue pairs $(\phi_1, \phi_2) \in \{(0.90, 0.70),(0.50, -0.50),(-0.90, -0.70)\}$, representing highly persistent, moderately persistent, and oscillatory dynamics, respectively.

The loading matrices $\bfm A \in \mathbb{R}^{d_2 \times r_1}$ and $\bfm B \in \mathbb{R}^{d_1 \times r_2}$ are constructed via SVD to satisfy the loading strength conditions in Assumption (ref). Specifically, singular values are geometrically spaced between $d_k^{(1-\delta_0)/2}$ and $d_k^{(1-\delta_1)/2}$ for $k \in \{1,2\}$. This ensures $\|\bfm A\|_2 \asymp d_2^{(1-\delta_0)/2} $ and $\sigma_{r_1}(\bfm A) \asymp d_2^{(1-\delta_1)/2}$, with analogous conditions for $\bfm B$. We consider two factor loading strength regimes. In the strong factor loading regime, $(\delta_0,\delta_1)=(0,0)$, so that all row and column factors are pervasive and the loading matrix singular values grow at the same rate; in the weak factor loading regime, $(\delta_0,\delta_1)=(0.3,0.5)$, which increases the spread between the largest and smallest singular values for the loading matrices. The noise matrix $\bfm E_t$ has i.i.d. $N(0,\sigma^2_\varepsilon)$ entries with $\sigma^2_\varepsilon=1$, independent of the row and column factor processes.

Finite-sample validation of estimation error bounds. Figures (ref) and (ref) display the loading space estimation accuracy across varying dimensions $d_1 = d_2 \in \{50, 100, 200\}$ and sample sizes $T \in \{100, 200, 400, 800\}$ under both strong and weak factor loading strength regimes. For each configuration $(d, T)$, we fix one realization of $(\bfm A, \bfm B)$ and generate $500$ independent Monte Carlo replicates with new factor and noise processes. The figures display the average logarithm of estimation error $\log({\cal D}(\bfm U_{\bfm A}, \widehat \bfm U_{\bfm A}))$, where ${\cal D}(\bfm U_{\bfm A}, \widehat \bfm U_{\bfm A})=\|\widehat \bfm U_{\bfm A} \widehat \bfm U_{\bfm A}^\top - \bfm U_{\bfm A}\bfm U_{\bfm A}^\top\|_2$, for three methods: MINE (Algorithm (ref), initialization only), P-COMPAS (Algorithm (ref) with partial complement projections and MINE initialization), and COMPAS (Algorithm (ref) with full complement projections and MINE initialization).

The results reveal several key findings. First, the proposed COMPAS consistently achieves the lowest loading space estimation error across all configurations, substantially improving over the MINE initialization. Second, P-COMPAS exhibits degraded performance compared to COMPAS, confirming that partial (subset) orthogonal complement projections sacrifice statistical accuracy for computational efficiency, a tradeoff formalized in Theorem (ref). Third, a striking pattern emerges when comparing the two loading strength regimes: despite similarly poor MINE initialization in both cases, the strong factor loading regime (Figure (ref)) enables substantially better final recovery than the weak loading regime (Figure (ref)). This occurs because stronger loadings provide larger signal strength, allowing iterative refinement to effectively leverage the modewise additive structure and achieve sharper error bounds, whereas weaker factors fundamentally limit recovery accuracy. Analogous results for $\widehat \bfm U_{\bfm B}$ appear in Appendix (ref).

Asymptotic normality for rows of the loading matrices. Figures (ref) and (ref) assess the asymptotic normality for the estimated loading matrix $\widehat \bfm U_{\bfm A}$, evaluating the oracle asymptotic distribution from Theorem (ref) and the feasible data-driven inference procedure from Theorem (ref). We focus on the first row of the loading matrix $\widehat \bfm U_{\bfm A}$ under configurations $d_1=d_2\in\{20,50,100\} $ and $T = 200$ in the strong factor loading regime across $500$ Monte Carlo replicates. Corresponding results for $\widehat \bfm U_{\bfm B}$ appear in Appendix (ref). For each replicate, we construct both oracle standardized rows of the loading matrix using true population covariance matrices and feasible standardized rows of the loading matrix using plug-in covariance estimators from the estimated factors and residual-based noise estimates. We assess the quality of normal approximation using the first coordinate of the standardized first row of $\widehat\bfm U_{\bfm A}-\bfm U_{\bfm A}\bfm R_{\bfm A}$ through Q-Q plots (Figure (ref)) and histograms overlaid with the standard normal density (Figure (ref)). The plots for other components are similar.

Both diagnostics demonstrate strong agreement with ${\cal N}(0,1)$: Q-Q plots align closely with the reference line, and histograms match the theoretical density well. These results validate the asymptotic normality theory and confirm that the feasible data-driven inference procedure accurately captures the asymptotic distribution of the estimated loading matrices without requiring knowledge of population covariance matrices.

Real Data Analysis

We evaluate the proposed MAFM on a quarterly OECD macroeconomic panel in which each observation is naturally represented as a country-by-indicator matrix observed over time.

Data and preprocessing. The data are compiled from OECD releases that provide consistent quarterly coverage for a common set of advanced economies and macroeconomic indicators. After intersecting country availability across sources, we obtain a balanced panel of $d_1 = 18$ countries observed over $n=131$ quarters (1991Q2–2023Q4) with $d_2 = 9$ indicators capturing prices, real activity, financial conditions, and labor markets. Following standard macroeconometric practice, price indices, production indices, and GDP are transformed via log first differences, while interest rates and unemployment are first differenced and scaled by $1/100$. At each quarter $t$, the transformed indicators form the raw matrix $\bfm Y_t \in \mathbb{R}^{d_1 \times d_2}$, with rows corresponding to countries and columns to indicators.

To remove scale differences across indicators, we apply pooled standardization: for each indicator $j\in[d_2]$, we compute the pooled mean $\mu_j$ and standard deviation $\sigma_j$ across all country-quarter observations and set $X_{t,ij} = (Y_{t,ij}-\mu_j)/\sigma_j$. This follows standard practice in macroeconomic panel studies ottinger2025history. All estimation and forecasting use the standardized series $\{\bfm X_t\}_{t=1}^n$.

Competing methods. We compare MAFM using MINE initialization and COMPAS refinement against four alternative approaches: (1) TwDFM yuan2023twoway, which decomposes matrix time series into separate row and column components with autoregressive dynamics via two-step quasi-maximum likelihood estimation; (2) Tucker factor model han2024tensor with cross-covariance based estimation; (3) CP factor model chen2026estimation with cross-covariance based estimation; and (4) VecFactor bai2003inferential, a vectorized baseline that flattens matrix observations into vectors before applying standard factor analysis with cross-covariance estimation. For fair comparison, all methods use AR models to forecast factors, which are then used to reconstruct predicted observations.

Factor rank selection. Figure (ref) displays eigenvalue proportion diagnostics for factor rank selection. For MAFM, the A-side (row-factor component) scree plot suggests $r_1=4$ factors (clear elbow at component 4), while the B-side (column-factor component) plot suggests $r_2 = 2$ (sharp decay after component 2), yielding $(r_1, r_2) = (4, 2)$. For robustness, we also examine $(4, 1)$, $(5, 1)$, and $(5, 2)$. Competing methods are evaluated across multiple factor ranks as shown in Table (ref).

In-sample fit. We assess in-sample performance using all observations for estimation. For a given method with fitted values $\widehat{\bfm X}_t$, define the fitting error as $\text{Fit-err} = (d_1d_2n)^{-1} \sum_{t=1}^{n} \| \bfm X_t - \widehat\bfm X_t \|_{\rm F}^2$, and the coefficient of determination as $R^2 = 1 - \sum_{t=1}^{n} \| \bfm X_t - \widehat \bfm X_t \|_{\rm F}^2/\sum_{t=1}^{n} \| \bfm X_t - \bar{x} \mathbf{1}_{d_1} \mathbf{1}_{d_2}^\top \|_{\rm F}^2$, where $\bar{x} = (d_1d_2n)^{-1} \sum_{t=1}^{n} \mathbf{1}_{d_1}^\top \bfm X_t \mathbf{1}_{d_2} $ is the grand mean. Table (ref) reports Fit-err and in-sample $R^2$ for MAFM and the four competing methods across factor rank specifications. Table (ref) reports the total computation times for all configurations. MAFM consistently outperforms alternatives across all specifications while maintaining computational efficiency. Notably, MAFM is substantially faster than TwDFM, which achieves only slightly lower in-sample $R^2$ values. This demonstrates the computational advantage of COMPAS over quasi-maximum likelihood estimation.

Out-of-sample forecasting. We evaluate predictive performance using an expanding-window design with one-step-ahead forecasts. At each forecast origin $w$, we re-estimate the model using $\{\bfm X_t\}_{t=1}^{w}$ through the following steps: (i) estimate loading matrices $\widehat{\bfm U}_{\bfm A}, \widehat{\bfm U}_{\bfm B}$ and factor sequences $\{\widehat{\bfm F}_t\}_{t=1}^{w}, \{\widehat{\bfm G}_t\}_{t=1}^{w}$; (ii) fit independent $AR(p)$ models to each factor component with $p$ selected by AIC; and (iii) generate factor forecasts $\widehat{\bfm F}_{w+1}, \widehat{\bfm G}_{w+1}$ and reconstruct $\widehat{\bfm X}_{w+1|w} = \widehat{\bfm F}_{w+1}\widehat{\bfm U}_{\bfm A}^{\top} + \widehat{\bfm U}_{\bfm B}\widehat{\bfm G}_{w+1}^{\top}$. The one-step-ahead forecast error is ${\rm FE}(w+1) = (d_1d_2)^{-1} \|\bfm X_{w+1}-\widehat \bfm X_{w+1|w}\|_{\rm F}^2 $. To assess recent forecasting performance, we report the averaged forecast errors over the last $h$ forecasts: $\overline{{\rm FE}}_{h} = h^{-1}\sum_{w=n-h}^{n-1} {\rm FE}(w+1)$ for $h \in \{5, 10, 15, 20,25,30\}$. Table (ref) reports results for MAFM and competing methods. MAFM consistently achieves the lowest forecast errors across all horizons, confirming that its modewise additive structure yields genuine predictive gains beyond in-sample fit.

\noindentInterpretation of estimated loadings. Using the feasible data-driven inference procedure in our MAFM framework (Theorem (ref)), we construct 95% confidence intervals for the factor loadings, displayed in Figures (ref) and (ref). The uniformly tight confidence intervals indicate that the main patterns are well-identified and precisely estimated.

For the row-factor loadings (indicators), Factor 1 exhibits a clear global macro comovement pattern: broadly positive loadings on prices, output, production, and interest rates, paired with a strongly negative loading on unemployment. The remaining factors capture more specialized contrasts among specific indicator subgroups. For the column-factor loadings (countries), Factor 1 loads positively and significantly on essentially every country, indicating a pervasive common component affecting all economies. In contrast, Factor 2 is geographically concentrated, Norway and the Netherlands exhibit large opposite-signed loadings while most other countries remain near zero. Overall, the decomposition reveals one dominant common factor alongside higher-order factors that capture selective cross-sectional variation.

Conclusion

In this paper, we propose a modewise additive factor model for matrix time series that extends Tucker and CP factor models by separately capturing row-specific and column-specific effects, rather than attributing all variation to global factors. To address the cross-modal interference inherent in this structure, we develop a two-stage estimation procedure, MINE for initialization and COMPAS for iterative refinement via orthogonal complement projections. We establish convergence rates and asymptotic distributions under high-dimensional settings where each mode dimension is comparable to or exceeds the sample size. As a technical contribution of independent interest, we derive sharp Bernstein-type inequalities for quadratic forms of dependent matrix time series. Simulations demonstrate the advantages of our method, and an application to OECD macroeconomic data illustrates its practical relevance.

\setcounter{page}{1}