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.
81,852 characters · 15 sections · 29 citation commands
Decomposing Co-Movements in Matrix-Valued Time Series: A Pseudo-Structural Reduced-Rank Approach
\newtheorem{remark}{Remark}
Keywords: Co-movements, common features, matrix-valued time series, reduced rank, structured parameterization
JEL: C32, C55, F20
In recent decades, macroeconomic and financial time series have expanded both in number and complexity. This complexity is inherently multi-dimensional -- researchers routinely observe entities (e.g., countries, firms) across multiple indicators (e.g., GDP growth, inflation, unemployment) over time. The time series thus oftentimes display a matrix structure: A series of matrix data are observed over time. Traditional methods, such as vector autoregressions (VARs; e.g., lutkepohl2005new) would typically vectorize matrix-valued data into a single high-dimensional vector which loses the row- and column-specific information of the data. Panel VARs holtz1988estimating, koop_forecasting_2019 can handle multiple cross-sections but do not naturally model the two-way dependencies that are intrinsic to matrix-valued data. The methods for matrix-valued time series (MVTS; chen_autoregressive_2021, tsay_matrix-variate_2024, zhang2024additive, samadi_matrix-valued_2025) address this gap by preserving and exploiting the matrix structure of the data. However, MVTS models typically still encounter difficulties with the large number of parameters that need to be estimated. Recently, however, reduced-rank matrix autoregressive models (RRMAR, xiao_reduced_forth) have shown their promise to this end, by exploiting low-rank decompositions of autoregressive coefficient matrices.
In this paper, we use the RRMAR to decompose contemporaneous co-movements in matrix-valued time series into three interpretable components, in line with economic intuition: row-specific, column-specific, and joint (row and column) co-movements, thereby isolating co-movement relations within and across the matrix dimensions. We obtain this decomposition by casting the RRMAR in a pseudo-structural form, imposing identification restrictions (analogous to those in structural VARs) that rotate factor matrices so that co-movements partition into the three components. We interpret these co-movements through the lens of serial correlation common features engle_testing_1993, where a common serial correlation occurs if a linear combination of a series is free of serial correlation even though each individual series displays it. The pseudo-structural specification then permits standard asymptotic inference on row-, column-, and joint co-movement parameters.
Furthermore, we propose a carefully-designed algorithm to identify the co-movements and offer practical guidance, through a BIC-type criterion, on how to select the reduced ranks. The finite-sample properties of estimation, inference, and rank selection procedures are evaluated through Monte Carlo simulations, including experiments with rank under- and over-specification. The simulations demonstrate good accuracy in parameter recovery in cases of rank under- and over-specification, while valid coverage of confidence intervals is only present in cases of rank over-specification. Moreover, rank selection procedures using BIC often identify the true rank, even in cases of rank over-specification for one of the matrix-valued time series dimensions. Finally, we illustrate the pseudo-structural framework in two empirical examples that differ markedly in their co-movement structure. The first examines macroeconomic indicators across several North American and Eurozone economies, where we uncover distinct row-, column-, and joint co-movement patterns that reflect both within-country relations and cross-country linkages. The second considers coincident and leading indexes across multiple U.S. states, a setting in which our rank selection procedure indicates an absence of low-rank contemporaneous relations. This underscores the ability of our rank selection procedure to adapt to both low-rank and full-rank environments.
Our decomposition builds on existing work of matrix-valued time series (MVTS) and reduced-rank regression but differs in both emphasis and methodology. xiao_reduced_forth introduce the reduced-rank matrix autoregressive model (RRMAR), provide estimation procedures, establish parameter consistency, propose an extended BIC for rank selection, and evaluate model performance in forecasting exercises. In contrast, we show that the RRMAR can be used to characterize contemporaneous co-movements in matrix-valued time series and to decompose those co-movements into three interpretable components (row, column, and joint interactions). Our simulation study examines the kernel densities of the decomposed parameters under three rank-specification scenarios — overestimation, underestimation, and correct specification — and complements xiao_reduced_forth by reporting results for the traditional BIC and for models with longer lag lengths.
Beyond the RRMAR literature, MVTS factor models wang_factor_2019, chen_factor_2022, chen2023inference have been used to extract common row/column factors for dimension reduction, and reduced-rank regression methods cubadda_reduced_2022, cubadda_vector_2017, cubadda_representation_2019, cubadda_dimension_2022 have been applied to forecasting and co-movement detection escribano1994cointegration, cubadda_studying_2009. Moreover, this dimension reduction can be extended to the nonstationary case, where cointegration may occur over the row and column dimensions li_cointegrated_2024, chen2025inference, hecq2025cointegrated, lopetuso_cointegrated_2025. By contrast, our pseudo-structural reduced-rank formulation explicitly separates contemporaneous interactions of stationary matrix-valued time series into row, column, and joint components — a decomposition not provided by existing high-dimensional methods that separate lagged from contemporaneous predictors wang_high-dimensional_2022, wang_high-dimensional_2023. Our inferential procedure targets these distinct contemporaneous components directly, yielding interpretable measures of how rows and columns interact contemporaneously in matrix-valued time series.
The remainder of this paper is structured as follows. Section (ref) starts by building the foundation for the pseudo-structural representation of the RRMAR. Section (ref) details the estimation of the pseudo-structural model and how we select the ranks of the coefficient matrices. Section (ref) contains a simulation study in which we evaluate estimation and inference of pseudo-structural parameters and our proposed rank selection procedure. Section (ref) gives two examples of our method applied to coincident and leading indicators of U.S.\ states and various economic indicators of different countries. Section (ref) concludes and discusses future research avenues.
A word on the notation. Throughout this paper, we denote scalars by small letters $x$, vectors by boldface small letters $\mathbf{x}$, and matrices by boldface capital letters $\mathbf{X}$. For a generic matrix $\mathbf{X}$, we call $\mathbf{X}^\top$, $\|\mathbf{X}\|_F$, $\text{vec}(\mathbf{X})$, respectively, the transpose, the Frobenius norm, and column-wise vectorization. Finally, denote the nullspace and column space of a matrix by $\mathcal{N}(\cdot)$ and $\mathcal{C}(\cdot)$, respectively.
We begin by reviewing the reduced rank matrix-valued time series models under study in Section (ref), and then define contemporaneous co-movements within the framework of serial correlation common features. In Section (ref), we present the equivalent pseudo-structural form for these reduced rank matrix-valued time series models. Section (ref) intuitively discusses the link between the serial correlation common features and the reduced rank matrix autoregressive model through an example of countries and economic indicators.
We review the reduced-rank matrix autoregressive (RRMAR) xiao_reduced_forth model that forms the basis of our analysis of contemporaneous co-movements in matrix-valued time series. The model, with one autoregressive lag, is given by
where $\mathbf{Y}_{t} \in \mathbb{R}^{N_{1} \times N_{2}}$ is the response matrix, $\mathbf{U}_{1}, \mathbf{U}_{3} \in \mathbb{R}^{N_{1} \times r_{1}}$, and $\mathbf{U}_{2}, \mathbf{U}_{4} \in \mathbb{R}^{N_{2} \times r_{2}}$ are the coefficient matrices with (reduced) row rank $1 \leq r_{1} \leq N_{1}$ and column rank $1 \leq r_{2} \leq N_{2}$. We assume that the errors follow a matrix-valued normal distribution dawid1981matrix, namely
where $\boldsymbol{\Sigma}_{1} \in \mathbb{R}^{N_1 \times N_1}$ and $\boldsymbol{\Sigma}_{2} \in \mathbb{R}^{N_2 \times N_2}$ are positive definite row- and column-covariance matrices, and the notation $MVN(\cdot,\cdot,\cdot)$ denotes the matrix-valued normal distribution and $N(\cdot,\cdot)$ denotes the multivariate normal distribution. To simplify the analysis, we assume that each series has been demeaned over time to remove the constant term. Defining $\mathbf{y}_{t} = \operatorname{vec}(\mathbf{Y}_{t})$, the equivalent vectorized form is
For stationarity of the model, we require that the spectral radius of $\mathbf{A} \in \mathbb{R}^{N_1N_2 \times N_1 N_2}$ is strictly less than one. The Kronecker structure of $\mathbf{A}$ suggests the presence of shared dynamic behavior across the different series, which we formalize using the concept of common contemporaneous co-movements introduced by engle_testing_1993.
We introduce two left null space matrices, $\boldsymbol{\delta} \in \mathbb{R}^{N_{1} \times (N_{1} - r_{1})}$ and $\boldsymbol{\gamma} \in \mathbb{R}^{N_{2} \times (N_{2} - r_{2})}$, which annihilate the row and column dynamics, respectively -- that is, they satisfy $\boldsymbol{\delta}^{\top} \mathbf{U}_{1} = \mathbf{0}$ and $\boldsymbol{\gamma}^{\top} \mathbf{U}_{2} = \boldsymbol{0}$, where $\boldsymbol{0}$ denotes a conformable matrix of zeros. In economic terms, the null space matrices $\boldsymbol{\delta}$ and $\boldsymbol{\gamma}$ identify linear combinations of rows and columns of $\mathbf{Y}_t$ that remove the serial correlation generated by $\mathbf{U}_1$ and $\mathbf{U}_2$, and these linear combinations reveal the presence of common contemporaneous co-movements. However, due to the Kronecker structure of the coefficient matrix, we have three matrices in total that annihilate the serially correlated component $(\mathbf{U}_{2} \otimes \mathbf{U}_{1}) (\mathbf{U}_{4} \otimes \mathbf{U}_{3})^{\top} \mathbf{y}_{t-1}$:
This raises the question: How do the individual null spaces $\boldsymbol{\delta}$ and $\boldsymbol{\gamma}$ interact to annihilate the Kronecker product dynamics $\mathbf{U}_{2} \otimes \mathbf{U}_{1}$? Proposition (ref) resolves this by providing an explicit form for the null space of $(\mathbf{U}_{2} \otimes \mathbf{U}_{1})^{\top}$ which decomposes into three orthogonal components; its proof is given in Appendix (ref).
The three orthogonal components in the decomposition of Proposition (ref) correspond to (i) column-specific co-movements $\left(\mathcal{N}(\mathbf{U}_{2}^{\top}) \otimes \mathcal{C}(\mathbf{U}_{1})\right)$, (ii) row-specific co-movements $\left(\mathcal{C}(\mathbf{U}_{2}) \otimes \mathcal{N}(\mathbf{U}_{1}^{\top})\right)$, and (iii) joint co-movements $\left(\mathcal{N}(\mathbf{U}_{2}^{\top}) \otimes \mathcal{N}(\mathbf{U}_{1}^{\top})\right)$. However, the interpretation of the joint co-movement component is not immediately obvious. In what follows, we first show how each of the three components can be formalized in a pseudo-structural model; which is an equivalent representation of the RRMAR model in (ref). We end with a toy example to illustrate how each component can be intuitively interpreted.
vahid_common_1993 construct a pseudo-structural form that is algebraically equivalent to a reduced-rank VAR and decomposes dynamics into contemporaneous and lagged components. Here, to make the paper self-contained, we first review their approach and then extend it to matrix-valued time series by deriving an analogous pseudo-structural representation for the RRMAR. To keep notation (relatively) compact, we give the pseudo-structural model corresponding to the RRMAR in (ref) with a single lag. The correspondence also holds for multi-lag RRMAR models, albeit with additional notation (see Remark (ref)).
Let the reduced-rank VAR($1$) for an \(N\)-dimensional stationary vector process with rank \(r\) be given by \[ \mathbf{y}_{t} = \mathbf{A}\mathbf{B}^\top \mathbf{y}_{t-1} + \mathbf{e}_{t}, \] where \(\mathbf{y}_t\in\mathbb{R}^N\), \(\mathbf{A},\mathbf{B}\in\mathbb{R}^{N\times r}\), and \(\operatorname{rank}(\mathbf{A}\mathbf{B}^\top)= 1 \leq r<N\). Hence, there exists a left null space basis \(\boldsymbol{\psi}\in\mathbb{R}^{N\times(N-r)}\) with \(\boldsymbol{\psi}^\top\mathbf{A}=\mathbf{0}\). Rotating \(\boldsymbol{\psi}\) so its top block is the identity, we may write \[ \boldsymbol{\psi}=
, \] with $\boldsymbol{\psi}^*\in\mathbb{R}^{r\times (N-r)}$. Then \(\boldsymbol{\psi}^\top\mathbf{y}_t\) yields \((N-r)\) pseudo-structural equations, while the remaining \(r\) equations form reduced-form regressions that use unrestricted lags as instruments. Stacking these relations gives
where \(\boldsymbol{\Omega} \in \mathbb{R}^{N \times N}\) encodes contemporaneous relations and \(\boldsymbol{\Pi} \in \mathbb{R}^{N \times N}\) collects the lagged dynamics.
We now apply the same logic to the RR-MAR in (ref). Partition \(\mathbf{Y}_t\) into four (unbalanced) blocks:
where $\mathbf{Y}_{11,t} \in \mathbb{R}^{(N_{1} - r_{1}) \times (N_{2} - r_{2})}$, $\mathbf{Y}_{12,t} \in \mathbb{R}^{(N_{1} - r_{1}) \times r_{2}}$, $\mathbf{Y}_{21,t}\in \mathbb{R}^{r_{1} \times (N_{2} - r_{2})}$, and $\mathbf{Y}_{22,t} \in \mathbb{R}^{r_{1} \times r_{2}}$. We rotate the parameters $\boldsymbol{\delta} \in \mathbb{R}^{N_1 \times (N_1 - r_1)}$ and $\boldsymbol{\gamma} \in \mathbb{R}^{N_2 \times (N_2 - r_2)}$ to have the $(N_{1} - r_{1})$- and $(N_{2} - r_{2})$-dimensional identity sub-matrix respectively
where $\boldsymbol{\delta}^* \in \mathbb{R}^{r_1 \times (N_1-r_1)}$ and $\boldsymbol{\gamma}^* \in \mathbb{R}^{r_2 \times (N_2-r_2)}$.
We then have three pseudo-structural equations corresponding to the three annihilation conditions given in Proposition (ref). More specifically,
The first two equations capture the row-specific and column-specific co-movement restrictions; the third corresponds to the joint row- and column-wise co-movement structures. Although parts of the first and second relations are algebraically implied by the third, they also contain components that capture distinct aspects of the co-movement structure. We thus retain these non-redundant components when constructing the pseudo-structural form.
To obtain the pseudo-structural form, we vectorize the partitioned matrix time series in equation (ref)
where $\operatorname{vecb}(\cdot)$ denotes the block vectorization, and define $\boldsymbol{\Omega} \in \mathbb{R}^{N_1 N_2 \times N_1 N_2}$ and $\boldsymbol{\Pi} \in \mathbb{R}^{N_1 N_2 \times N_1 N_2}$ as {
} where $\boldsymbol{\Omega}$ contains the contemporaneous row- and column-specific relations in $\boldsymbol{\delta}^{*}$ and $\boldsymbol{\gamma}^{*}$, while $\boldsymbol{\Pi}$ holds the lags with a matrix autoregressive restriction (see e.g., chen_autoregressive_2021) as instruments. This results in the pseudo-structural form
This pseudo-structural form has exactly the same number of parameters as the RRMAR, as given by the formula $2r_{1}N_{1} - r_{1}^2 + 2r_{2}N_{2} - r_{2}^2$.
To illustrate the reduced-rank restriction and provide some intuition for the row- and column-specific co-movements as well as the joint co-movements, we consider a toy example with $N_1=3$ economic indicators corresponding to $N_2 = 4$ countries, in line with the $N_1 \times N_2 = 3 \times 4$ matrix-valued set-up of our empirical application in Section (ref). We denote the countries as the United States (USA), Canada (CAN), Germany (DEU), and France (FRA), and the economic indicators as the interest rate (IR), gross domestic product (GDP), and manufacturing production (PROD). Vectorizing the matrix-valued times series yields a $12 \times 12$ coefficient matrix $\mathbf{A}$, which is illustrated in Figure (ref). We discuss three cases of reduced rank matrix autoregressive models, (i) partially reduced only among the rows, (ii) partially reduced only among the columns, and (iii) reduced in both rows and columns.
Firstly, our focus is on a toy example where we reduce the rank of the economic indicator dimension (rows of the matrix-valued time series) of the coefficient matrix to two, hence $r_1 = 2$, and leave the country dimension at full rank with $r_2=4$. Thus, only $\boldsymbol{\delta}$ annihilates the dynamics of the system, and the rank of two implies one SCCF co-movement relation for all three economic indicators. In our toy example, visualized in the left matrix of Figure (ref), for simplicity, we let the GDP and PROD of each country co-move (as highlighted through the same coloring for each of the four countries in Figure (ref), left panel), but the IR is unrestricted (i.e.\ unrestricted elements are neither displayed nor colored in Figure (ref), left panel). We can quantify this relationship with $a_{1, i} = -\delta_{1}^* a_{2, i}$, $a_{4, i} = -\delta_{1}^* a_{5, i}$, $a_{7, i} = -\delta_{1}^* a_{8, i}$, and $a_{10, i} = - \delta_{1}^* a_{11, i}$ where $\boldsymbol\delta^* = (\delta_{1}^*, \delta_{2}^*)^\top$ with $\delta_{1}^*$ an arbitrary constant representing the scale with respect to the GDP (and the negative sign will become clear once we discussed the link with the pseudo-structural form below), and $\delta_{2}^* = 0$ since IR is left unrestricted (for simplicity). This toy example illustrates we can capture within-country co-movements, where (certain) indicators co-move for each country, but the countries themselves do not. Put differently, we can remove the serially correlated component through a linear combination of any country's GDP and the same country's PROD, indicative of a serial correlation common feature in the model. These within-country co-movements, captured by the $\boldsymbol{\delta} = (1, \delta_1^*, \delta_2^*)^\top$, can similarly be represented by the pseudo-structural system of equations. This form yields four relationships corresponding to each of the four countries in our toy example. We obtain
This normalization provides the contemporaneous co-movement relations with respect to the first economic indicator (GDP), and if $\delta_{2}^{*}$ is zero, we obtain the toy example given above. For instance, the co-movement relation for the economic indicators of the USA given in equation (ref) would be $y_{11t} = -\delta_1^* y_{21t}$ plus a white noise process represented by the first element of $\boldsymbol{\delta}^\top \mathbf{E}_t$. The connection with the toy example can now be directly made since the scaling term for the any country's PROD to recover the same country's GDP series is given by $-\delta_1^*$.
Secondly, we can have a low rank in the country dimension (columns of the matrix-valued time series), as shown on the right matrix of Figure (ref). Suppose the country dimension has reduced rank $r_2=3$, while the indicator dimension remains full rank with $r_1=3$. In this case, there exists a vector $\boldsymbol{\gamma}$ that annihilates the system’s dynamics. In the toy example, only the USA and FRA move together, as highlighted through the same coloring for each of the four indicators in Figure (ref), right panel); CAN and DEU remain unrestricted (and are therefore neither colored nor displayed). We can quantify this relationship by saying $a_{1, i} = - \gamma_{3}^* a_{10, i}$, $a_{2, i} = - \gamma_{3}^* a_{11, i}$, and $a_{3, i} = - \gamma_{3}^* a_{12, i}$ for $i = 1, \dots, 12$ and where $\boldsymbol\gamma^*=(\gamma_{1}^*, \gamma_{2}^*, \gamma_{3}^*)$, with $\gamma_{3}^*$ an arbitrary constant representing the scale of the last country (FRA) with respect to the first country (USA), and the other two elements are zero since the countries CAN and DEU are left unrestricted (for simplicity). Thus, USA and FRA co-move with one another by the scale of $\gamma_{3}^*$, given the normalization $\boldsymbol{\gamma} =(1, \gamma_{1}^{*}, \gamma_{2}^{*}, \gamma_{3}^{*})^\top$, while CAN and DEU, on the other hand, are not bound to move in tandem with any other country in our example. These dynamics thus represent instances of across-country co-movements, they can be represented by the pseudo-structural system of equations in an analogous manner as the within-country co-movements discussed above.
Finally, consider the case in which the rank is reduced along both dimensions. In line with the previous two examples, we have $N_1 = 3$ economic indicators and $N_2 = 4$ countries, but with a reduced rank in both dimensions, namely, $r_1 = 2$ and $r_2 = 3$. As in our previous examples, the matrices $\boldsymbol{\delta}$ and $\boldsymbol{\gamma}$ are defined as {4pt}
with the restriction $\delta_2^* = \gamma_1^* = \gamma_2^* = 0$ to align with our previous examples where GDP co-moves with industrial production and the USA co-moves with France. The difference here being that the system has a rank reduction in both dimensions, thus an interaction between the two co-movement relations will occur. Figure (ref) displays the full left null space structure for such a configuration. What we observe are the following three groupings of the null space matrix:
Here, we see that the row-specific and the column-specific co-movements display a particular structure in the full left null space, as visualized by the blue and red coloring: GDP co-moves with the industrial production (row-specific in blue) and the USA co-moves with France (column-specific in red). The joint co-movements form the overlap between the row-specific and the column-specific structures, as can be seen from the first column of the null space matrix that highlights (in purple) the interaction of the two null spaces.
This implies that the USA GDP co-moves not only with the USA industrial production, but also with the French GDP and the French industrial production.
As a final illustration, consider the special case where both $r_{1} = r_{2} = 1$. Then the left null space becomes large (dimension $11$ in the $12$-series toy example). We define $\boldsymbol{\delta}$ and $\boldsymbol{\gamma}$ as
One convenient construction of $\boldsymbol{\delta}$ and $\boldsymbol{\gamma}$ that yields an explicit null space matrix is given in Figure (ref). As can be seen, the terms associated with $\boldsymbol{\delta}$ form a block structure in the null space (as indicated in blue). The structure associated with $\boldsymbol{\gamma}$ is indicated in red. The terms associated with $\boldsymbol{\gamma} \otimes \boldsymbol{\delta}$ form the overlap between the two structures as visible from the purple shading in the figure. The figures thus show how the joint, row-specific, and column-specific null space components combine, which provides intuition for modeling the orthogonal subspaces explicitly, as in Proposition (ref).
This section first describes the full-information maximum likelihood (FIML) estimator for the pseudo-structural model (Section (ref)) and then discusses practical rank selection (Section (ref)).
From the pseudo-structural model with one lag in equation (ref), we define a $p$ lag pseudo-structural model as
where $\boldsymbol{\Pi}_j$ has the same structure as $\boldsymbol{\Pi}$ in equation (ref), but now in terms of $\mathbf{U}_{3,j}$ and $\mathbf{U}_{4,j}$. For given ranks $r_{1}$ and $r_{2}$, we estimate $\boldsymbol{\Omega}$, the lag coefficients $\{\boldsymbol{\Pi}_{j}\}_{j=1}^{p}$, and the covariance matrices $\boldsymbol{\Sigma}_{1}$, and $\boldsymbol{\Sigma}_{2}$ by maximizing the full-information likelihood $\mathcal{L}(\boldsymbol{\theta})$ defined by
where $\boldsymbol{\theta}$ collects the parameters $\boldsymbol{\delta}^*$, $\boldsymbol{\gamma}^*$, $\boldsymbol{\Sigma}_{1}$, $\boldsymbol{\Sigma}_{2}$, $\mathbf{U}_{3,j}$, and $\mathbf{U}_{4,j}$. Note that $\boldsymbol{\Omega}$ is omitted from the $\log | \cdot |$ terms because the determinant of $\boldsymbol{\Omega}$ is always one by construction (see Section (ref)).
Because the optimization problem is non-convex, we recommend initializing the estimator both from the RRMAR solution (given from xiao_reduced_forth) and from several randomized starts to reduce the risk of convergence to local optima. The RRMAR initialization exploits the multidimensional structure of the coefficient matrices and yields a consistent starting point close to a global optimum. We perform optimization with gradient-based methods (BFGS, nocedal2006numerical) and monitor convergence across starts to avoid local maxima and saddle points. Further algorithmic details appear in Appendix (ref).
In practice, the true ranks $r_1$ and $r_2$ are unknown and must be estimated. We use standard information criteria-- Akaike Information Criteria (AIC, akaike_new_1974) and Bayesian Information Criteria (BIC, schwarz_estimating_1978) --to jointly obtain an estimate for the ranks for the RRMAR. We therefore label this the rank selection criteria in the remainder of the paper. Define
where $\mathcal{L} (\widehat{\boldsymbol{\theta}})$ log likelihood value obtained from estimating the pseudo-structural form and $c_{T}$ is a penalty term such that $c_{T} = 2$ for the AIC and $c_{T} = \ln(T)$ for the BIC, and $\phi(r_1, r_2)$ is the number of parameters in a model with $p$ lags and reduced ranks $r_1$ and $r_2$, given by
Then the selected ranks (for fixed $p$) can be obtained based on the minimum value derived from either AIC or BIC.
We conduct two Monte Carlo experiments to assess (i) the inferential performance of our pseudo-structural parameter estimates and (ii) the effectiveness of the rank-selection procedure described in Section (ref). The first experiment examines the sampling densities and coverage probabilities of $\widehat{\boldsymbol{\delta}}$ and $\widehat{\boldsymbol{\gamma}}$. The second evaluates how well the procedure recovers the true ranks and the autoregressive lag order.
For each experiment, we consider the matrix dimensions: $N_1\times N_2 = 3\times 4$. Data are generated from the pseudo-structural model (ref) using specified matrices $\boldsymbol{\Omega}$ and $\boldsymbol{\Pi}$. We draw the true parameter vectors $\boldsymbol{\delta}^*$ and $\boldsymbol{\gamma}^*$ independently from standard normal distributions and construct $\boldsymbol{\Omega}$ as in (ref). The matrix $\boldsymbol{\Pi}$ is formed by sampling each column of $\mathbf{U}_{3,i}$ and $\mathbf{U}_{4,i}$ for $i = 1, \dots, p$ from independent standard normal distributions. The error matrices $\mathbf{E}_t$ are i.i.d.\ from a standard matrix-normal with mean $\mathbf{0}$ and row/column covariances $\boldsymbol{\Sigma}_1$ and $\boldsymbol{\Sigma}_2$. For each configuration we simulate $T+50$ observations, discard the first 50 as burn-in, and report results for $T=100$ and $T=250$. Throughout we fix the signal-to-noise ratio (SNR)\footnote{Defined as the ratio of the largest eigenvalue of the coefficient to the largest eigenvalue of the error covariance matrix.} at $0.7$.
To evaluate the estimation and inferential performance of our pseudo-structural estimates, we conduct $1000$ simulation runs for a matrix-valued time series with dimension $N_{1} \times N_{2} = 3 \times 4$ and ranks $r_1\times r_2 = 2\times 2$. Results for the $N_{1} \times N_{2} = 3\times 6$ case with ranks $r_1\times r_2 = 2\times 2$ are similar and available in Appendix (ref).
First, we estimate $\boldsymbol{\delta}^*$ under rank settings: (i) correctly estimated ranks $\widehat{r}_1\times \widehat{r}_2= 2\times 2$, (ii) underestimated $2\times 1$, and (iii) overestimated rank $2\times 3$ in the second dimension. For each simulation run, we
Figure (ref) displays kernel densities for the two components of $\widehat{\boldsymbol{\delta}}^*$. When the rank in the second dimension is specified correctly or overestimated, the densities are closely aligned and centered on the true values for both sample sizes. Under rank underestimation, the densities remain well-behaved but show modestly increased variance. Turning to Figure (ref), we plot the empirical coverage rates of the 95% confidence intervals for $\delta^*_1$ and $\delta^*_2$. Both correct and overestimated ranks in the second dimension achieve coverage close to the nominal 95% level. When the rank is underestimated, however, coverage deteriorates to around 92%.
Second, we estimate $\boldsymbol{\gamma}^*$ under rank settings: (i) correctly estimated ranks $2\times 3$, (ii) underestimated $1\times 3$, and (iii) overestimated $3\times 3$ in the first dimension. Figures (ref) and (ref) present the corresponding density and coverage plots. As with $\boldsymbol{\delta}$, correct specification or overestimation yields densities centered on the true values. Underestimating the rank of the first dimension increases the standard deviation of $\widehat{\boldsymbol{\gamma}}$ and reduces coverage to about 89%.
To evaluate our rank-selection procedures we estimate the pseudo-structural model, per simulation design, for all possible rank combinations and use the information criteria from Section (ref) to select the ranks $r_1$ and $r_2$. We hereby first set the autoregressive lag order to its true value $p=1$. We report the average selected rank (“Average Rank”), its standard deviation (“Std.\ Rank”), and the frequency of correct selection (“Freq.\ Correct”) over 100 simulation runs, where we fix $p$ to the true value. Results are summarized in Tables (ref) for the simulation design and different true rank settings.
For $N_1=3$, $N_2=4$ (Table (ref)), BIC selects the true rank at a high rate and outperforms AIC. For example, with true rank $(1,1)$ AIC selects the correct ranks about $68$% and $83$% of the time for the first and second dimensions when $T=100$, increasing to $73$% and $87$% at $T=250$. BIC selects the correct ranks at rates of $83$% and $98$% for $T=100$, rising to $89$% and $99$% for $T=250$. For partially reduced ranks (reduced over only one dimension), AIC selects the true rank at over $90$% for $T=100$, improving with larger $T$, while BIC is accurate across all settings. When the true rank is full, both AIC and BIC select the correct rank in all reported simulations.
In the setting where the true rank is $(2,1)$, AIC correctly selects the ranks at a rate of at least $76$%; similar or slightly higher rates can be observed when increasing the sample size. BIC selects the correct rank at a rate of at least $92$%, going up to at least $98$% for $T=250$. In this setting, we also explored how well information criteria perform in case we underestimate the true rank of the first dimension ($\widehat{r}_1 = 1$), when estimating the rank of the second dimension. In this case, both AIC and BIC still maintain $100$% accuracy when selecting the second rank, even when $T=100$. Overestimating the true rank of the first dimension ($\widehat{r}_1 = 3$) leads to the same conclusion.
We also investigate the joint performance of rank and lag selection using information criteria. To this end, we conduct an additional simulation study evaluating the ability of AIC and BIC to recover the true rank-lag specification, restricting the maximum lag order to two. The experimental setup is identical to the previous simulations, including the fully reduced case, partially reduced cases and the case with no rank reduction.
Table (ref) reports the results when the true lag is one. Overall, BIC identifies the correct rank and lag with high accuracy and consistently outperforms AIC. For example, when the true ranks are $(1,1)$, AIC selects the correct ranks only about $47$% and $43$% of the time and the true lag just $6$% of the time. By contrast, BIC selects the correct ranks $94$% and $100$% of the time, with the correct lag chosen approximately $89$% of the time. In the partially reduced cases, AIC struggles to identify the reduced rank dimension, correctly selecting it only about $40$% of the time, although it consistently chooses the full rank dimension ($100$%). This missclassification extends to lag selection: AIC tends to favor higher lag orders, while BIC always selects the true lag (one). A similar pattern emerges in the full rank case, where AIC again favors higher lags, while BIC selects the true lag with much greater reliability.
We repeat the joint rank and lag selection experiment when the true lag is two. The rank configurations (fully reduced, partially reduced, and no reduction) match the previous experiment, with the only change being that the data generating process now has lag order $p=2$. Table (ref) reports the average rank selected, standard deviation of the rank, and the frequency of correct selection. The main findings are, overall, qualitatively similar to the $p=1$ experiment, namely, BIC again outperforms AIC in jointly recovering both rank and lag, while AIC continues to overselect the reduced rank dimensions. However, unlike the $p=1$ case, both AIC and BIC now select the lag $100$% of the time.
For example, when the true ranks are $(1,1)$, AIC correctly identifies the ranks in approximately $75$% and $81$% of cases and recovers the true lag ($p=2$) $100$% of the time. By contrast, BIC correctly recovers the ranks about $87$% and $100$% of the time and selects the correct lag in $100$% of replications. In the partially reduced case, AIC continues to struggle with the reduced rank dimension (about $70$% in the worst case) while always selecting the full rank dimension ($100$%) and the lag ($100$%). BIC, however, selects the reduced rank dimension and lag correctly at higher rates (approximately $78$% in the worst case). For the full rank case, both AIC and BIC select the full rank dimensions and the lag order at near $100$% of the times.
We present two applications of the proposed pseudo-structural method. Section (ref) studies macroeconomic indicators across North American and Eurozone countries; Section (ref) examines coincident and leading indexes across U.S. States.
We consider an application with data from $N_1 = 3$ macroeconomic indicators in $N_2 = 4$ countries. The sample is quarterly from 1991Q1 to 2019Q4, totaling $T = 116$ observations. The indicators are real GDP, manufacturing production (PROD), and ten-year government bond yields (IR). The countries are the United States (USA), Canada (CAN), Germany (DEU), and France (FRA). Data are obtained from the Organization for Economic Cooperation and Development (OECD) at \url{https://data-explorer.oecd.org/}. For data transformations, we take the first differences of the interest rates, while GDP and manufacturing production are taken in log-differences. The transformed series for the three indicators and four countries are visualized in Figure (ref).
Our objective is to identify co-movement structures across indicators and countries. To this end, we compute the rank selection criteria with both AIC and BIC and with a maximum of four lags. The results point toward a reduced-rank structure along both economic indicator and country dimensions: AIC selects ranks $(\widehat{r}_{1}, \widehat{r}_{2}) = (2,3)$ with one lag, indicating reduced-rank structure along both dimensions. BIC, on the other hand, selects ranks $(\widehat{r}_{1}, \widehat{r}_{2}) = (2,1)$ with one lag, thereby signaling an even lighter reduced-rank structure along the indicator and country dimensions. Based on our simulation evidence favoring BIC, we proceed with the BIC-selected ranks $(2,1)$. For notational purposes, we write $y_{i,j,t}$ with $i$ indexing the indicator (GDP, PROD, IR), $j$ indexing the country (USA, CAN, DEU, FRA), and $t$ indexing time.
{\bf Row-Specific Co-movements.} Under the estimation strategy in Section (ref) and with estimated $\widehat{r}_1 = 2$, we obtain
where standard errors (in parenthesis) are from the observed information (Hessian at the MLE). The implied row-specific relation is
Thus, manufacturing production loads positively on GDP with coefficient $0.323$ (statistically significant, $p < 0.01$), while interest rate coefficient is not significant. These coefficients apply uniformly across the four countries.
{\bf Column-Specific Co-movements.} In addition to the row-specific co-movements, we obtain column-specific co-movement relations. Given $\widehat{r}_2 = 1$, we have
which yields three column-specific relations
for $i \in \{\text{GDP, PROD, IR}\}$. Thus, for all economic indicators, France loads positively on i) the USA with a coefficient of $1.190$ (statistically significant, $p < 0.01$), ii) Canada with a coefficient of $1.305$ (statistically significant, $p < 0.01$), and iii) Germany with a coefficient of $1.370$ (statistically significant, $p < 0.01$).
{\bf Joint Co-movements.} Combining the row and column components produces three joint co-movement equations: {
}
Each equation integrates i) row-specific co-movements (within country indicator relations, e.g., GDP-PROD coefficient $0.323$), ii) column-specific co-movements (cross-country relations, e.g., USA-FRA coefficient $1.190$), and iii) joint co-movements (within and across country relations, given by $\delta^*_i \gamma^*_j$ for $i \in \{\text{PROD, IR}\}$ and $j \in \{\text{USA, CAN, DEU}\}$)\footnote{Standard errors of the product can be obtained via the Delta method given the standard errors of the row-specific and column-specific parameters.}. The first three terms of each equation give us the row-specific and column-specific co-movement relations already given in the prior subsections. However, we see additionally that French industrial production co-moves with the GDP of the USA, Canada, and Germany with coefficients $0.384$, $0.421$, and $0.442$ respectively (all statistically significant, $p < 0.01$). The interest rate coefficients are again not significant across the countries.
We analyze two indicators ($N_1 = 2)$ --the monthly coincident and leading indices (CI and LI)-- for $N_2 = 9$ North Central U.S. States: Illinois (IL), Indiana (IN), Iowa (IA), Michigan (MI), Minnesota (MN), North Dakota (ND), South Dakota (SD), Ohio (OH), and Wisconsin (WI). The data for these series are sourced from the Federal Reserve Bank of Philadelphia at \url{https://www.philadelphiafed.org/surveys-and-data/regional-economic-analysis/}.
The coincident indexes combine four state-level indicators to summarize prevailing economic conditions. This includes indicators such as non-farm payroll employment, average hours worked in manufacturing by production workers, the unemployment rate, and wage and salary deflated by the consumer price index. These four variables are combined through a dynamic single-factor model, as outlined by stock_new_1989 into a single coincident index. The coincident index is seasonally adjusted on the sample spanning from January 1982 to February 2020, excluding the pandemic period, resulting in a sample size of $T=458$ monthly observations. Furthermore, we compute monthly growth rates for each coincident indicator as visualized in the top panel of Figure (ref).\footnote{For Michigan, we adjust the level value of July 1998 by averaging the values of June and August. This adjustment is made to mitigate the impact of a drop in levels, which otherwise generates two outliers in the first differences of the data.}
The leading indices are provided as monthly growth rates and are constructed to forecast six-month growth in the coincident index. Their components include the coincident index itself, state-level housing permits, state initial unemployment claims, delivery times from the Institute for Supply Management manufacturing survey, and the interest rate spread between the 10-year Treasury bond and the three-month Treasury bill. These leading indexes are plotted in the bottom panel of Figure (ref).
Because leading indices forecast future coincident performance, we do not necessarily expect contemporaneous co-movement between CI and LI. Indeed, each state's leading index is designed to forecast the six-month growth rate of its coincident index. The nature of the two indexes is thus different; coincident indicators reflect current economic conditions while leading indicators aim to predict future trends in the coincident index. Similarly, co-movement between the states may or may not occur since neighboring economies are often interconnected and interact with one another.
To investigate the possible presence of co-movement within and between the dimensions of this $2\times 9$ matrix-valued time series, we compute our rank-lag selection criteria with a maximum of three lags. AIC selects a full-rank model with three lags, while BIC selects a full-rank model with two lags. These results thus suggest no evidence for co-movement among the indicators or the states and highlight the capability of our selection procedure to handle full-rank coefficients and multiple lags when supported by the data.
This paper introduces a pseudo-structural framework for reduced-rank matrix-valued time series models, enabling the decomposition of contemporaneous co-movements into row-specific, column-specific, and joint components. By leveraging the Kronecker structure of the reduced-rank matrix autoregressive model, we derive interpretable linear combinations that annihilate serial correlation, akin to common feature analysis in vector autoregressions. Our approach provides explicit inference on these co-movement structures, supported by an estimation procedure to navigate the non-convex landscape and a rank-lag selection criterion validated in simulations. Empirical applications reveal distinct co-movement patterns: Macroeconomic indicators across countries exhibit strong row- and column-specific linkages, while coincident and leading indexes across U.S. states show no evidence of low-rank structures, underscoring the adaptability of our method to both reduced- and full-rank settings. All code and replication material for this paper can be found in the Julia repository PseudoStructuralComovements on the second author's GitHub page \url{https://github.com/ivanuricardo/pseudostructuralcomovements}.
The proposed methodology can be extended in several directions. One could consider a tensor-valued time series instead of a matrix-valued time series. The left null space would then require additional terms to further disaggregate the co-movement structure into three co-movement dimensions instead of simply row- and column-wise co-movement relations. Another interesting future research direction is the investigation of non-contemporaneous co-movements in matrix-valued time series (e.g., cubadda2001noncontemporaneous). This would allow for adjustment delays in the co-movement relations across the different rows and columns of the matrix-valued time series.
{\bf Acknowledgements.} The last author was financially supported by the Dutch Research Council (NWO) under grant number VI.Vidi.211.032.