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.
72,420 characters · 14 sections · 29 citation commands
The Cointegrated Matrix Autoregressive Model
Multivariate time-series data have traditionally been treated as vectors despite the fact that the relationships among variables can give rise to more complex structures. Recent studies have focused on observations arranged in matrix form, extending the vector autoregressive (hereafter VAR) model to its matrix counterpart, commonly referred to as the matrix autoregressive (hereafter MAR) model.
Matrix-valued time series naturally arise when data are collected at the intersection of two classifications. In economics, a common example is the analysis of macroeconomic indicators across countries or regions. For example, wang2009bayesian examines the simultaneous evolution of employment statistics in US states and industrial sectors, while zhang2024additive investigates relationships among the main macroeconomic indicators in different countries. In such cases, it is natural to structure the data as matrices, with the economic indicators on the rows and the corresponding regional entities on the columns (or vice versa).
Another prominent application involves the analysis of bilateral trade flows, where the matrix representation effectively captures import-export relationships between countries. This framework facilitates the modeling of trade dynamics, as illustrated in chen2019modeling.
In finance, stock market trends exhibit dependencies not only over time and across economic variables but also across different sectors, such as manufacturing and transportation. Furthermore, measures of financial connectedness are often derived from sequences of adjacency matrices that represent networks of financial assets; billio2021matrix. These applications have spurred the development of high-dimensional methods for matrix-valued time series, including autoregressive models, billio2023bayesian, matrix panel regression models, kapetanios2021estimation, and matrix factor models, chen2023statistical.
The utility of matrix-valued time series extends beyond economics. In spatio-temporal analysis, where classifications often involve geographical and temporal dimensions, the matrix framework provides a natural representation (e.g., sun2023matrix and yu2024dynamic). Neuroscience is another domain where matrix-valued data is prevalent. Functional magnetic resonance (fMRI) data, for example, capture interactions across different brain regions, with dependencies often visualized through matrices representing connectivity patterns; samadi2025matrix.
The autoregressive representations are of particular interest in economics. The MAR model, as proposed in chen2021autoregressive, is expressed as:
where $X_t$ is an $m \times n$ matrix, $\Lambda$ and $\Psi$ are square matrices of dimensions $m$ and $n$, respectively, and $E_t$ is an $m \times n$ white noise matrix. The vectorization of (ref) yields the following:
which corresponds to the classical VAR form under specific (Kronecker-based) parameter restrictions.
The matrix representation offers two main advantages. First, it provides a significant reduction in dimensionality, decreasing the number of coefficients from $n^2m^2$ to $n^2 + m^2$. Second, it preserves the inherent structure of the data, facilitating the interpretation of relationships between rows and columns. This paper focuses on the latter advantage, examining the interpretability of MAR coefficients in the context of non-stationary and cointegrated systems.
Multivariate cointegrated systems, such as the vector error correction model (hereafter VECM), are widely valued for their economic interpretability. The VECM is represented as:
where $Y_t$ is a $k-$dimensional vector, $\alpha$ and $\beta$ are $k \times r$ matrices, and $r<k$ is the cointegration rank. Economically, the adjustment matrix $\alpha$ captures the behavior of agents responding to equilibrium errors by driving the economic variables back to the steady-state position defined by $\beta' Y = 0$. Specifically, $\beta$ characterizes the relations to which the variables are attracted, while $\alpha$ quantifies the direction and speed of their adjustment in response to deviations from equilibrium.
The connection between the coefficients of the VECM and the economic dynamics allows us to test numerous relevant economic hypotheses (see johansen1991estimation). In particular, hypotheses about the cointegration matrix $\beta$ enable the assessment of the equilibrium relationships suggested in the macroeconomic literature. Furthermore, we can interpret hypotheses on the adjustment matrix $\alpha$ as assertions about adjustment mechanisms, providing insights into which variables remain unaffected by disequilibrium conditions, possibly establishing non-causality relationships (see, for example, juselius2007taking and juselius2017real).\footnote{A null row in the adjustment matrix $\alpha$ implies weak exogeneity for the long-run coefficients $\beta$, as defined in engle1983exogeneity (see johansen1992testing for a formal proof).} Moreover, testing restrictions on $\alpha$ allows the identification of the variables serving as driving trends.
This paper extends the analytical framework of cointegrated systems to matrix-valued time series. This topic has recently received attention, with a proposal put forward by li2024cointegrated and the work of hecq2024detecting. In their modeling strategy, the cointegrated matrix structure is obtained by imposing restrictions on the coefficients of a vector error-correction model and is driven by a bilinear equilibrium condition. We take a different perspective and derive the error-correction representation directly from an underlying matrix autoregressive (MAR) specification. This distinction is important. In our framework, the matrix structure is preserved both in the autoregressive representation in levels and in the corresponding error-correction form, so that the matrix nature of the model is intrinsic to the data-generating process, rather than being a consequence of a particular parametrization. Moreover, our model is driven by two sets of linear equilibrium conditions acting separately but jointly, over rows and over columns. Further, we show that existing cointegrated models for matrix-valued time series do not generally share this property, as their matrix form depends on whether the model is written in autoregressive or error-correction form.
The second contribution of the paper is interpretative. In our formulation, the left and right cointegration matrices admit separate economic interpretations, as they characterize long-run equilibrium relations across the rows and columns of the observed matrix process, respectively. This feature is not available in existing alternatives, where the equilibrium condition is intrinsically bilinear and therefore combines row and column effects into a single object. As a result, those models do not allow one to disentangle the equilibrium structure along the two dimensions of the data. Our approach, instead, makes it possible to study row and column steady states separately, thereby providing a more flexible and informative framework for testing economically meaningful hypotheses, in close analogy with the role played by cointegration vectors in the standard VECM literature.
We provide further contributions related to the estimation of the model. We discuss both the identification of the two cointegration ranks that characterized the cointegrated MAR model and the estimation of the cointegration matrices. In addition, we also tackle an inferential perspective and introduce approaches for testing hypotheses on the cointegration and adjustment matrices. Finally, we present an empirical analysis that examines the relationship between industrial production and inflation expectations in a selected group of European countries, demonstrating the practical application of the proposed methodology. The goodness of fit of the model is assessed by verifying the stationarity of the estimated cointegration relationships and by comparing the matrix-based estimation with the standard vector representation. Furthermore, we illustrate that employing a conventional Vector Error Correction Model does not facilitate an interpretation in terms of parallel structures. This highlights the improvement in interpretability achieved through our proposed model.
The paper is structured as follows. Section 2 proposes the model. Section 3 clarifies our contribution to the literature by comparing our model with existing formulations and highlighting the strengths of our approach. Section 4 presents an estimation algorithm and outlines procedures for testing restrictions on model coefficients and rank identification, while Section 5 evaluates the performance of these methods through Monte Carlo experiments. Section 6 provides an empirical illustration, and Section 7 concludes. Proofs of the lemmas and theorems are provided in Appendix A.
Consider the matrix autoregressive process of order one
where $\Lambda\in\mathbb{R}^{m\times m}$, $\Psi\in\mathbb{R}^{n\times n}$, $E_t$ is a $m\times n$ white-noise process, and $X_t$ is a $m\times n$ observation of a matrix-valued time series.
Theorem (ref) shows that the vectorized process $x_t=\operatorname{vec}(X_t)$ is cointegrated. Therefore, the matrix autoregressive representation in (ref), together with Assumption (ref), gives rise to a cointegrated model with MAR representation, which we call a Cointegrated Matrix Autoregressive model (C-MAR).
We now study the possibility of a matrix error-correction representation. From (ref), one obtains the standard vector error-correction form
Consider the following lemma.
It follows from Lemma (ref) that the impact matrix of the VECM associated with a C-MAR model cannot be written in matrix form.
To find an alternative matrix representation, note that under Assumption (ref), \[ \Pi_1:=\Lambda-I_m, \qquad \Pi_2:=\Psi-I_n \] have ranks $r_1$ and $r_2$, respectively. Hence, rank factorizations of the form \[ \Pi_1=\tau\gamma', \qquad \Pi_2=\varphi\theta', \] are admissible, with $\tau,\gamma\in\mathbb{R}^{m\times r_1}$ and $\varphi,\theta\in\mathbb{R}^{n\times r_2}$ of full column rank. It follows that
Solving for $\Delta X_t=X_t-X_{t-1}$ yields the error-correction representation
We refer to model (ref) as the error-correction representation of the C-MAR model (hereafter, ECC-MAR).
As emphasized in the introduction, cointegrated models allow for an interpretation in terms of adjustment toward a long-run equilibrium. In a standard VECM, the equilibrium condition is given by $\beta'x_t=0$. In fact, although $x_t$ is $I(1)$, the linear combination $\beta'x_t$ is $I(0)$. This means that the system reverts toward the position $\beta'x_t=0$ despite the variables being non-stationary.
To understand whether ECC-MAR carries the same interpretation, we study the existence of mean-reverting linear combinations associated with the ECC-MAR coefficients $\gamma$ and $\theta$.
It follows from Theorem (ref) that processes $\gamma'X_t$ and $X_t\theta$ are mean reverting. Consequently, $\gamma'X_t$ and $X_t\theta$ can be interpreted as two distinct equilibrium errors: one pertaining to the relations among the rows of $X_t$, and the other pertaining to the relations among the columns of $X_t$.
If both $\gamma'X_{t-1}=0$ and $X_{t-1}\theta=0$, then the deterministic error-correction component in (ref) vanishes, and therefore $\mathbb{E}(\Delta X_t\mid \mathcal{F}_{t-1})=0.$ Thus, if both $\gamma'X_t=0$ and $X_t\theta=0$, there is no systematic tendency for the process to move away from its current position. In this case, the deterministic error-correction component vanishes, so the system is in equilibrium. Consequently, $\gamma$ and $\theta$ play an analogous role to the cointegration matrix $\beta$ in the standard VECM, and may be interpreted as row and column cointegration matrices, respectively.
Similarly, matrices $\tau$ and $\varphi$ play a role analogous to the adjustment matrix in a VECM. They capture the speed of adjustment toward the long-run equilibrium and indicate which variables do not respond to disequilibrium conditions.
To illustrate the adjustment mechanism implied by ECC-MAR, rewrite (ref) as \[ \Delta X_t = \tau \big(\gamma' X_{t-1} + \gamma' X_{t-1} \theta \varphi'\big) + X_{t-1} \theta \varphi' + E_t. \] The first term in parentheses reflects the row disequilibrium at time $t-1$. The second component, $\gamma' X_{t-1} \theta \varphi'$, captures the row disequilibrium induced by the adjustment mechanism acting on the columns of $X_t$. Indeed, $X_{t-1}\theta\varphi'$ represents the system's response to deviations from the column equilibrium condition $X_{t-1}\theta=0$, while $\gamma'X_{t-1}\theta\varphi'$ measures the row disequilibrium generated by that column correction. Therefore, the adjustment toward row equilibrium is driven by the composite disequilibrium term $\gamma' X_{t-1} + \gamma' X_{t-1}\theta\varphi'.$ The matrix $\tau$ determines how such deviations are corrected over time.
A symmetric argument applies to the column disequilibrium. In particular, the adjustment coefficient $\varphi$ governs the response to the composite disequilibrium term $X_{t-1}\theta + \tau\gamma'X_{t-1}\theta.$
In conclusion, the third term in (ref), namely $\tau\gamma'X_{t-1}\theta\varphi',$ reconciles the two parallel adjustment mechanisms: one targeting $\gamma'X_t=0$ and the other targeting $X_t\theta=0$. These two adjustment mechanisms may otherwise point in different directions, and the interaction term ensures their consistency within a unified matrix error-correction structure.
Finally, we study the relationship between ECC-MAR and its vector counterpart. Vectorizing (ref), we obtain
The following theorem proves that the impact matrix of the VECM generated by the ECC-MAR representation has a reduced rank and provides a representation of the cointegration and adjustment matrices.
The natural extension of the ECC-MAR representation in equation (ref) is
The correspondence between the ECC-MAR above and the C-MAR($p$) is not straightforward as for the case where the autoregressive form is of order 1, and requires assumptions on the relations among the MAR coefficients. To illustrate, consider the MAR($p$) model
and its alternative representation
Vectorizing, one obtains
To obtain the ECC-MAR representation in (ref), it is necessary that $\sum_{i=j}^p \Psi_i \otimes \Lambda_i=A_j \otimes B_j$ $\forall j=\{1,2,\dots,p\}$. In other words, it is necessary that the sum of these Kronecker products is itself a Kronecker product.
This is not generally the case, as the rank condition illustrated in Lemma (ref) is not ensured for the sum of Kronecker products. The following Theorem illustrates a condition ensuring that the sum of Kronecker products is itself a Kronecker product. We first introduce two Lemmas.
To ensure the duality between MAR and ECC-MAR representations, the condition that $\sum_{i=j}^p \Psi_i \otimes \Lambda_i$ is itself a Kronecker product is not sufficient. In fact, it is required that the resulting matrices can be decomposed in such a way to obtain the terms governing the adjustment mechanism, namely, it is required that $A \otimes B=\sum_{i=1}^p \Psi_i \otimes \Lambda_i$ have eigenvalues lower and equal to 1 in absolute value, so that they can be decomposed as $B=I+\Pi_1$ and $A=I+\Pi_2$, where $\Pi_{1,2}$ are rank deficient.
The stability conditions hold in the usual way, by analyzing the eigenvalues of the companion matrix.
From Lemmas (ref) and (ref), the following Theorem immediately results:
Concepts related to cointegration and reduced-rank modeling have already been incorporated into the literature on matrix-valued time series. For example, xiao2022reduced introduce the Reduced Rank Matrix Autoregressive model (RRMAR), which extends the MAR model in (ref) by imposing reduced-rank restrictions on the autoregressive coefficient matrices. However, the first contribution to explicitly introduce cointegration into matrix-valued autoregressive models is li2024cointegrated; see also the extension in hecq2024detecting. In its simplest form, their model is
where $\Pi_1=\tau\gamma'$ and $\Pi_2=\varphi\theta'$ are reduced-rank matrices of ranks $r_1$ and $r_2$, respectively, with full-column-ranks $\tau,\gamma\in\mathbb{R}^{m\times r_1}$ and $\varphi,\theta\in\mathbb{R}^{n\times r_2}$.
The corresponding vectorized representation is a standard VECM with coefficient matrix constrained to have Kronecker-product form:
To ensure that the vectorized process $x_t=\operatorname{vec}(X_t)$ is $I(1)$ and cointegrated, the following Johansen-type condition is imposed on the characteristic polynomial of the vectorized system.
This construction differs fundamentally from the C-MAR proposed in the present paper for two main reasons. First, the model in (ref) is postulated directly in error-correction form and is not derived from a MAR representation in levels. Second, its equilibrium structure is intrinsically bilinear and does not disentangle the equilibrium relations of the rows and columns.
To illustrate the first point, vectorizing (ref) yields the VECM representation in (ref), where $\varphi\otimes\tau$ and $\theta'\otimes\gamma'$ are the adjustment and cointegration matrices, respectively. However, (ref) does not imply a MAR representation in levels of the form (ref). Indeed, the MECM can be written in VAR form as
Consider the following lemma.
Lemma (ref) shows that the VAR representation associated with (ref) does not admit a MAR representation in levels. Therefore, the term cointegrated MAR would be inappropriate. For this reason, we refer to (ref) as a Matrix Error Correction Model (hereafter MECM), in order to distinguish it from the C-MAR developed in the previous section, which instead originates from a MAR.
The absence of a MAR representation is not merely a mathematical detail; it has important consequences for the interpretation of the equilibrium structure in the MECM. In this model, the stationary long-run component is the bilinear term $\gamma'X_t\theta$. Hence, the equilibrium condition is given by $\gamma'X_t\theta = 0,$, which in general cannot be decomposed into two separate one-sided equilibrium relations.
To further illustrate the point, consider the following theorem:
Theorem (ref) further clarifies why the equilibrium structure of the MECM is fundamentally different from that of the C-MAR. Since the one-sided processes $\gamma'X_t$ and $X_t\theta$ are themselves $I(1)$ rather than $I(0)$, they cannot be interpreted as equilibrium errors or mean-reverting deviations from a steady state. In a cointegrated system, only stationary combinations can play this role, because an equilibrium relation must remain bounded over time and pull the system back toward a long-run position. Here, instead, the only stationary object is the bilinear combination $\gamma'X_t\theta$. Therefore, the long-run equilibrium in the MECM is intrinsically joint, involving the row and column structures simultaneously, and cannot be decomposed into two separate equilibrium relations.
In light of the previous discussion, the C-MAR proposed in this paper should be considered conceptually distinct from the matrix cointegration models already available in the literature. Although both the MECM and the C-MAR are designed to model cointegrated matrix-valued time series, they differ in both structure and interpretation. The MECM is formulated directly in error-correction form and delivers an intrinsically bilinear equilibrium condition. By contrast, the C-MAR originates from a genuine MAR representation and gives rise to two separate equilibrium errors, which admit an independent interpretation in terms of row- and column long-run relations.
We estimate the ECC--MAR model by maximum likelihood and implement the estimation through an alternating algorithm, in the spirit of the iterative likelihood-based procedures proposed by li2024cointegrated.
We begin from the general ECC--MAR$(p)$ in equation (ref). To obtain a likelihood-based estimator, we impose the additional assumption that the innovation matrix is Gaussian with separable covariance as in chen2021autoregressive:
where $\Sigma_r \in \mathbb{R}^{m\times m}$ and $\Sigma_c \in \mathbb{R}^{n\times n}$ are positive definite. Equivalently, $\mathrm{vec}(E_t)\sim N\!\left(0,\Sigma_c\otimes\Sigma_r\right).$ Let
where $\vartheta = \big( \tau,\gamma,\phi,\theta,\Gamma_{1,1},\Gamma_{2,1},\dots,\Gamma_{1,p-1},\Gamma_{2,p-1} \big),$ and define the structural residual
Conditional on the initial values and under (ref), we have $\mathcal{U}_t(\vartheta)=E_t \sim MN_{m,n}(0,\Sigma_r,\Sigma_c).$
It follows that the conditional likelihood of the sample is
where $T_0 := T-p+1$. Taking logarithms yields the Gaussian log-likelihood
To estimate (ref), we employ an alternating maximum likelihood algorithm. The idea is to update the coefficients of the row structure while holding fixed the coefficients of the column structure, and then to proceed symmetrically in the opposite direction.
Note that post-multiplying (ref) by $\varphi_\perp$ eliminates the right-hand error-correction terms and yields
The projected errors satisfy $E_t\varphi_\perp \sim MN\!\left(0,\Sigma_r,\varphi_\perp'\Sigma_c\varphi_\perp\right)$. Hence, defining $C_\varphi:=\varphi_\perp'\Sigma_c\varphi_\perp$ and right-whitened variables $\widetilde Y_t^{L}:=\Delta X_t\varphi_\perp C_\varphi^{-1/2},$ $\widetilde X_t^{L}:=X_{t-1}\varphi_\perp C_\varphi^{-1/2},$ $ \widetilde Z_{i,t}^{L}:=\Delta X_{t-i}\Gamma_{2,i}'\varphi_\perp C_\varphi^{-1/2},$ equation (ref) becomes
with $\widetilde E_t^{L}\sim MN(0,\Sigma_r,I)$.
Similarly, pre-multiplying (ref) by $\tau_\perp'$ eliminates the left-hand error-correction terms and gives
Again, the projected errors satisfy $\tau_\perp'E_t \sim MN\!\left(0,\tau_\perp'\Sigma_r\tau_\perp,\Sigma_c\right).$ Therefore, defining $C_\tau:=\tau_\perp'\Sigma_r\tau_\perp$ and $\widetilde Y_t^{R}:=\Delta X_t'\tau_\perp C_\tau^{-1/2},$ $\widetilde X_t^{R}:=X_{t-1}'\tau_\perp C_\tau^{-1/2},$ $\widetilde Z_{i,t}^{R}:=\Delta X_{t-i}'\Gamma_{1,i}'\tau_\perp C_\tau^{-1/2},$ we obtain
with $\widetilde E_t^{R}\sim MN(0,\Sigma_c,I)$.
Equations (ref) and (ref) can then be estimated using standard Gaussian reduced-rank regression methods; see anderson1951estimating and johansen1995likelihood. After whitening, the auxiliary matrix equations can be read column by column as standard multivariate reduced-rank regressions with common coefficient matrices. The corresponding likelihood therefore decomposes as a sum over all time-column pairs $(t,j)$, so that the estimation can be based on pooled cross-products of the transformed variables. Appendix B illustrates the Johansen estimation procedure step by step in the present setting.
The factorization of the long-run component into $(\tau,\gamma)$ and $(\varphi,\theta)$ is not unique, as in standard reduced-rank representations, so the economically meaningful objects are the associated row and column cointegration spaces, with specific matrix representatives selected by normalization. Likewise, under the separable Gaussian assumption, the covariance factors $\Sigma_r$ and $\Sigma_c$ are identified only up to a common scale factor, since only the Kronecker product $\Sigma_c\otimes\Sigma_r$ is uniquely determined.
The preceding results motivate the following alternating maximum likelihood procedure.
Two identification remarks are in order. First, the cointegration and adjustment matrices are identified only up to the usual rank-preserving normalizations, so inference is formulated in terms of the associated cointegration spaces. Second, under the separable Gaussian assumption $\mathrm{vec}(E_t)\sim N(0,\Sigma_c\otimes\Sigma_r)$, the covariance factors $\Sigma_r$ and $\Sigma_c$ are identified only up to a common scale normalization.
The estimation procedure requires prior knowledge of the cointegration ranks $r_1$ and $r_2$. By Theorem (ref), these ranks are linked to the cointegration rank $r$ of the vectorized system through
where $m$ and $n$ are known. Since $r_1$ and $r_2$ are ranks, they must be non-negative integers. The first step is therefore to estimate the overall cointegration rank $r$ of the vectorized model $\mathrm{vec}(X_t)$ using the standard Johansen tests, such as the trace test or the maximum eigenvalue test. The asymptotic properties of these tests are well established; see Chapter 11 of johansen1995likelihood, in particular Theorem 11.1.
In many cases, solving (ref) over the set of admissible integer values yields a unique pair $(r_1,r_2)$. For example, if $m=4$, $n=3$, and $\hat r=8$, then (ref) implies the unique admissible solution $(r_1,r_2)=(2,1)$. In this case, the values of $(r_1,r_2)$ are completely determined by $\hat r$, so the only possible source of error is an incorrect estimate of the vectorized cointegration rank $r$.
When more than one admissible pair satisfies (ref), the ambiguity can be resolved using the stationarity implications of the ECC--MAR model. This is formalized in the following proposition.
Proposition (ref) applies only when more than one admissible pair remains. In that case, each candidate model is estimated and one checks whether the transformed processes $\hat\gamma'X_t$ and $X_t\hat\theta$ are stationary in all components. The correct pair $(r_1,r_2)$ is identified as the one for which both sets of transformed series are stationary, while any candidate that overstates one of the two ranks must generate at least one non-stationary component. Stationarity can be assessed by standard unit root tests, such as the Dickey--Fuller test; see dickey1979distribution.
In summary, the determination of $(r_1,r_2)$ proceeds in two steps. First, estimate the vectorized rank $r$ and enumerate all admissible integer pairs that satisfy (ref). Second, if this yields more than one admissible pair, estimate the competing specifications and select the one for which $\hat\gamma'X_t$ and $X_t\hat\theta$ are stationary in all components.
From an economic point of view, the cointegration matrices $\gamma$ and $\theta$ describe the long-run equilibrium structure of the row and column spaces, respectively, while the adjustment matrices $\tau$ and $\varphi$ describe the speed with which deviations from these equilibrium relations are corrected. This interpretation of the ECC--MAR coefficients makes it possible to associate parameter restrictions with economically meaningful hypotheses.
A first class of economically relevant hypotheses concerns uniform restrictions on the cointegration matrices:
where $H_\gamma$ and $H_\theta$ are design matrices. These restrictions are imposed on the entire cointegration space and are useful, for example, to test whether a variable is long-run excludable, through a zero-row restriction, or whether a long-run homogeneity relation holds across all equilibrium relations. Economically, these tests assess whether a given variable contributes to the long-run structure at all or whether certain transformed variables, such as spreads, ratios, or real quantities, provide a more meaningful description of equilibrium behavior.
A second class of hypotheses concerns whether a known vector belongs to the cointegration space, namely, whether
for known vectors $g$ and $t$. These restrictions are useful when economic theory suggests a specific equilibrium relation, for example, a spread, a unit-elasticity relation, or the stationarity of a single variable. In this case, the test asks whether the proposed vector belongs to the stationarity space spanned by the estimated cointegration relations.
A third class of hypotheses concerns restrictions on the adjustment matrices:
Of particular interest is the case of zero-row restrictions, which corresponds to weak exogeneity. If a row of $\tau$ or $\varphi$ is zero, the corresponding variable does not respond to disequilibria in the long run. Economically, such variables can be interpreted as drivers of the common stochastic trends, whereas variables with non-zero adjustment coefficients are those that bear the burden of restoring equilibrium. Hence, tests on the adjustment matrices help distinguish between the variables that push the system and those that pull it back toward its long-run path.
These three classes of restrictions cover the hypotheses that are most relevant in empirical applications of the ECC--MAR model. Many additional restrictions can be handled. For a detailed treatment of such cases in the standard C-VAR model, see Chapters 10 and 11 of juselius2006cointegrated.
The ECC--MAR model admits a natural likelihood-based testing framework through the auxiliary regressions (ref) and (ref). As shown above, after whitening these auxiliary systems are standard Gaussian reduced-rank regressions, so that inference can be conducted by likelihood ratio tests, exactly as in the classical cointegrated VAR model. The only difference is that, in the present matrix setting, the relevant sample moment matrices are pooled across the transformed column-wise observations, exactly as in the estimation step. Under the null, the corresponding likelihood ratio statistics are asymptotically chi-square, with degrees of freedom determined by the number of imposed restrictions. Appendix C reports the explicit construction of each test statistic.
We conducted a Monte Carlo study to evaluate three aspects of the proposed methodology: (i) finite-sample estimation accuracy, (ii) rank selection, and (iii) empirical performance of the restriction tests.
Data are generated from the ECC-MAR model in (ref). The matrices $\tau$ and $\varphi$ are constructed with the first $m-r_1$ and $n-r_2$ rows set to zero, respectively, while the remaining entries are drawn from independent standard normal distributions. The matrices $\gamma$ and $\theta$ are generated independently of standard normal distributions. This construction ensures $\mathrm{rank}(\tau\gamma')=r_1$ and $\mathrm{rank}(\varphi\theta')=r_2$, thus inducing $r_1$ row-wise and $r_2$ column-wise cointegration relations. To obtain a $I(1)$ process, we retain only parameter draws satisfying the stability conditions on the non-unit eigenvalues of $I_{r_1}+\gamma'\tau$ and $I_{r_2}+\theta'\phi$. Innovations are generated as i.i.d. matrix-valued white noise with $E[E_t]=0$ and $\mathrm{Var}(\mathrm{vec}(E_t))=I_{m\times n}$. From $X_0=0$, the process is generated recursively from (ref). A burn-in of 100 observations is used throughout.
We vary the sample size $T$, the matrix dimensions $(m,n)$, and the cointegration ranks $(r_1,r_2)$. Following li2024cointegrated, we consider $(m,n)\in\{(4,3),(6,5),(8,7)\}$. For $(m,n)=(4,3)$ we set $(r_1,r_2)\in\{(1,1),(2,2),(3,2)\}$, while for the other two dimensions we consider $(r_1,r_2)\in\{(1,1),(2,2),(3,3)\}$.
We compare the finite-sample accuracy of the proposed ECC-MAR estimator with that of a standard CVAR estimated on the vectorized process $\mathrm{vec}(X_t)$. In this subsection, the ranks $r_1$, $r_2$, and $r$ are assumed known and fixed at their true values, so the analysis focuses exclusively on the estimation of the long-run structure.
Since individual coefficient matrices are not uniquely identified, we assess performance through the cointegration space. Let $\beta_T$ denote the true cointegration matrix implied by the DGP, constructed from the true $\gamma$ and $\theta$ according to (ref). Let $\hat\beta_{\mathrm{ECC}}$ and $\hat\beta_{\mathrm{CVAR}}$ denote the corresponding estimators obtained from the ECC-MAR and CVAR models. For each estimator, we consider the projection matrix $P_\beta=\beta(\beta'\beta)^{-1}\beta'.$ For each Monte Carlo replication, the accuracy of the estimation is measured by the spectral norm distance $\|\hat P_{\beta}-P_T\|$, where $P_T$ is the projector associated with the true cointegration space.
Figure (ref) summarizes the results. As expected, both estimators improve with the sample size, with a clear reduction in dispersion from $T=100$ to $T=1000$. In general, the ECC-MAR estimator consistently outperforms the unrestricted CVAR in recovering the cointegration space. The difference between the two estimators is larger when the cointegration ranks are low, and it tends to shrink as $r_1$ and $r_2$ increase. Equivalently, the gain from exploiting the matrix structure is more pronounced when the number of common stochastic trends, $mn-r=(m-r_1)(n-r_2)$, is large, whereas it becomes smaller when $r$ approaches $mn$.
This pattern is clearly visible within the same matrix dimension. For example, in the $(m,n)=(4,3)$ designs, the gap between ECC-MAR and CVAR is evident for $(r_1,r_2)=(1,1)$, whereas the performance of the two estimators are comparable for $(r_1,r_2)=(3,2)$. A similar pattern emerges in the other configurations. For instance, for $(m,n)=(8,7)$, the advantage of ECC-MAR is more pronounced when $(r_1,r_2)=(1,1)$ than when $(r_1,r_2)=(3,3)$.
The first step is to determine the total cointegration rank $r$ of the vectorized process $\mathrm{vec}(X_t)$. This can be done using standard CVAR procedures, such as the Johansen trace test. Since our focus is on the identification of the structural ranks $(r_1,r_2)$ in the ECC-MAR model, we do not further investigate the determination of $r$. Moreover, misspecification of $r$ affects both the ECC-MAR and the unrestricted CVAR benchmark in a similar way. Consequently, in this subsection we treat $r$ as known.
Given $r$, the problem reduces to identifying the admissible pair $(r_1,r_2)$ consistent with (ref). In four of the nine designs considered in the simulation study, this mapping is one-to-one, so the structural ranks are uniquely determined and no additional step is required. This happens for $(m,n,r)=(4,3,6)$, $(4,3,11)$, $(6,5,10)$, and $(8,7,14)$, corresponding respectively to $(r_1,r_2)=(1,1),(3,2),(1,1)$ and $(1,1)$. However, in the remaining five cases, the same total rank is compatible with two different pairs $(r_1,r_2)$.
Following Section 4.1, for each admissible pair $(r_1,r_2)$ we estimate the ECC-MAR model under the corresponding rank restrictions and obtain $\hat\gamma$ and $\hat\theta$. We then test the stationarity of the implied cointegrating components by applying Augmented Dickey--Fuller tests to the estimated series $\{\hat\gamma'X_t\}$ and $\{X_t\hat\theta\}$. The admissible pair is selected when these transformed processes are found to be stationary for one candidate but not for the other.
Table (ref) reports the results for the five ambiguous configurations. In all cases, the correct rank pair is the outcome most frequently selected, and its frequency increases markedly with the sample size. By contrast, the incorrect admissible pair is only chosen rarely, often with frequencies close to zero even at $T=100$.
The main difficulty is the occurrence of undefined outcomes. Two cases may arise. In Und.\,1, both admissible pairs pass the stationarity checks; in Und.\,2, neither pair does. The second type is dominant in small samples, but decreases sharply as $T$ increases, reflecting the low power of ADF tests in short samples. For example, in configuration $(8,7,36)$ its frequency drops from 41.6% with $T=100$ to 7.8% with $T=1000$.
The behavior of Und.\,1 is more subtle. Its frequency tends to increase with the sample size, most visible in the case $(4,3,10)$, where it increases from 14.2% to 19.8%. This can occur when the model is estimated under a rank larger than the true one. In that case, the additional estimated cointegrating vector cannot belong exactly to the true cointegration space, but as the sample grows, it may become nearly collinear with it. Equivalently, the estimated cointegrating space may become increasingly close to the true one while preserving the full column rank. As a result, the extra linear combination can appear progressively more stationary, increasing the probability that the ADF test rejects the unit-root null.
We examine this mechanism in the configuration $(m,n,r)=(4,3,10)$, where the increase in Und.\,1 is strongest. Let $\gamma$ be the true cointegration matrix, $\gamma_\perp$ its orthogonal complement, and $\gamma_T$ the misspecified estimate obtained under incorrect rank $r_1=3$. Define $Q_{\gamma_\perp}=\mathrm{orth}(\gamma_\perp),$ $ Q_{\gamma,T}=\mathrm{orth}(\gamma_T),$ and $M_T = Q_{\gamma,T}'Q_{\gamma_\perp},$ $ d_T=\|M_T\|_F.$ Since both matrices have orthonormal columns, $d_T$ measures the extent to which the estimated cointegrating space is not orthogonal to the true common-trend space. Smaller values of $d_T$ therefore indicate that the misspecified estimated space is closer to the true cointegration space. In the Monte Carlo experiment, the event $d_{1000}<d_{100}$ occurs in 89.8% of replications, supporting the interpretation above.
In general, four stylized facts emerge. First, the correct rank pair is the outcome that is the most frequently selected even for small sample sizes, and its frequency increases with the sample size. Second, the incorrect admissible pair is only rarely selected. Third, Und.\,1 tends to become more frequent as $T$ grows. Fourth, Und.\,2 is mainly a small-sample phenomenon and decreases substantially with larger samples. Hence, the proposed procedure performs satisfactorily in general, although undefined cases remain its main limitation. In such cases, one may rely on an auxiliary criterion, such as an information criterion or prior structural knowledge, to select among the remaining admissible pairs.
Next, we study the finite-sample behavior of the proposed hypothesis tests for the ECC-MAR model. Unlike the previous Monte Carlo experiments, where the DGP parameters were randomly generated, here we consider a fixed design in order to assess size and power more directly.
We consider an ECC-MAR model with $(m,n)=(4,3)$ and $(r_1,r_2)=(2,2)$, with parameters \[ \tau =
, \qquad \gamma =
, \qquad \phi =
, \qquad \theta =
. \]
For the size analysis, we consider three true null hypotheses: (i) the first row of $\tau$ is zero, corresponding to a weak exogeneity restriction; (ii) $[\,1\;\;1\;\;0\,]\theta = 0$; and (iii) $[\,0\;\;-0.5\;\;-0.5\;\;1\,]' \in \mathrm{span}(\gamma)$.
Table (ref) shows that the tests have a satisfactory finite-sample size. For restrictions on cointegration matrices, rejection frequencies move toward the nominal level as $T$ increases. The weak exogeneity test is already well sized even in the smallest sample with an empirical rejection frequency that fluctuates around the nominal level and do not exhibit a systematic trend with \(T\).
We then examine power under false null hypotheses. Specifically, we test $\tau_{(3)}=0$, the linear restrictions $[\,1\;\;1.1\;\;0\,]\theta=0$, $[\,1\;\;1.3\;\;0\,]\theta=0$, and $[\,1\;\;1.5\;\;0\,]\theta=0$, and the membership restrictions $[\,0\;\;-0.5\;\;-0.6\;\;1\,]'$, $[\,0\;\;-0.5\;\;-0.8\;\;1\,]'$, and $[\,0\;\;-0.5\;\;-1\;\;1\,]' \in \mathrm{span}(\gamma)$.
The power results in Table (ref) are very strong. For $T\geq 500$, empirical power is one for all alternatives. At $T=100$, power is still 1 for the weak exogeneity test and above 97% for the membership tests on $\mathrm{span}(\gamma)$. The only case in which power is more moderate is the restriction on $\theta$ under the closest alternative, $[\,1\;\;1.1\;\;0\,]\theta=0$, where the rejection frequency is 69.8%. As expected, power increases monotonically as the alternative moves farther from the null, reaching 99.6% for hypothesis $[\,1\;\;1.5\;\;0\,]\theta = 0$.
Overall, the simulations indicate that the proposed tests combine satisfactory size with very strong power, especially in moderate and large samples.
To illustrate the empirical usefulness of the proposed C-MAR model, we consider a macroeconomic application with $m=2$ variables observed for $n=3$ European countries. Specifically, we study the joint dynamics of consumer inflation expectations ($\pi_e$) and industrial production ($\mathrm{In}$) for Belgium (B), Germany (D), and Austria (A).
From an economic perspective, stronger industrial production is expected to signal robust demand conditions, inducing firms to raise prices and consumers to revise inflation expectations upward. We therefore expect a long-run relation of the form $\pi_e=\varphi\,\mathrm{In}$ with $\varphi>0$. Across countries, the common monetary policy and the harmonized institutional environment suggest similar long-run dynamics in both inflation expectations and industrial production.
The analysis is based on monthly data from January 2000 to December 2019, thus excluding the Covid-19 period and the Russia--Ukraine war. Inflation expectations are proxied by consumer survey data from FRED, while industrial production indices are taken from Eurostat; all series are seasonally adjusted. Figures (ref) and (ref) plot the data.
A preliminary analysis indicates that all six variables are $I(1)$, that the appropriate lag order is $p=1$, and that the total cointegration rank of the vectorized system is $r=4$. For brevity, the Dickey--Fuller tests, the BIC-based lag-order selection, and the Johansen trace test results are reported in Appendix D.
By Theorem (ref), the vectorized cointegration rank satisfies $r = nr_1 + mr_2 - r_1r_2.$ Since here $m=2$ and $r=4$, it follows that $r_1=1$ and $r_2=1$. Estimating the model by the procedure of Section 4 yields the unrestricted estimates
These estimates imply a long-run relation of the form $\pi_e = 0.294\,\mathrm{In},$ which is consistent with the prior expectation of a positive long-run association between inflation expectations and industrial production.
We next test weak exogeneity by imposing row restrictions on the adjustment matrices $\tau$ and $\varphi$. The results are reported in Table (ref).
Consistent with our hypotheses, the results reveal a long-term relationship between price expectations and industrial production of the form $\pi_e = 0.294 \text{In}$. The positive coefficient confirms the anticipated association.
Next, we test for restrictions on the adjustment matrices $\varphi$ and $\tau$ to determine which variable adjusts in response to a disequilibrium. In other words, we test whether the $i^{\text{th}}$ row of $\varphi$ and $\tau$ is null. The results of the weak exogeneity test are recorded in table (ref).
The results indicate that industrial production is weakly exogenous with respect to $\gamma$, while inflation expectations adjust to restore the equilibrium relation. This is consistent with the interpretation of industrial production as the driving stochastic trend and inflation expectations as the adjusting variable.
Turning to the cross-country dimension, Germany appears to be weakly exogenous, suggesting that it acts as the driving trend within the country equilibrium relation, while Belgium and Austria adjust to restore the long-run balance. This is economically plausible given Germany’s dominant role in the euro-area economy.
We next assess the adequacy of the estimated C-MAR specification by examining whether the estimated long-run relations are indeed stationary. If $\hat\gamma$ and $\hat\theta$ correctly identify the row-wise and column-wise cointegration spaces, then the transformed processes $\hat\gamma'X_t$ and $X_t\hat\theta$ should be $I(0)$, even though the entries of $X_t$ are individually $I(1)$. We test stationarity of the five resulting linear combinations with a Dickey--Fuller test. The results, reported in Table (ref), reject the unit-root null in all five cases, supporting the interpretation of $\hat\gamma'X_t=0$ and $X_t\hat\theta=0$ as equilibrium relations.
Overall, the C-MAR specification provides a coherent description of the data by jointly modeling row and column cointegration. The empirical results support a positive long-run relationship between inflation expectations and industrial production, with industrial production acting as the common driving trend. At the country level, the estimates point to a long-run equilibrium across countries, with Germany behaving as the dominant stochastic trend and Belgium and Austria adjusting toward equilibrium.
For comparison, we estimate a standard VECM for the six vectorized series. The estimated model is
Unlike the C-MAR model, the unrestricted VECM does not separate the row and column cointegration mechanisms, making economic interpretation substantially more difficult. In particular, the signs and magnitudes linking industrial production and inflation expectations vary markedly across cointegrating vectors and across normalizations, so that the implied long-run relationships are not stable across countries. For example, one normalization may suggest a positive association between industrial production and inflation expectations in one country, while another normalization implies a negative relation for another country. This lack of robustness makes the economic interpretation of the unrestricted VECM considerably less transparent than that of the matrix model.
This instability is not surprising. As shown in Section 3, the unrestricted vector representation mixes the two underlying structures, whereas the C-MAR model isolates them explicitly through separate row-wise and column-wise cointegration components. As a result, the VECM coefficients do not admit a decomposition with the same structural interpretation as in the matrix specification.
The same issue arises when considering the adjustment mechanism. While the C-MAR model identifies industrial production as the driving force behind inflation expectations and Germany as the driving force in the cross-country relation, the unrestricted VECM does not yield a similarly clear pattern. In fact, the null of a zero row in the adjustment matrix is rejected for all six variables, as shown in Table (ref).
A further distinction concerns the interpretation of weak exogeneity itself. In the C-MAR framework, weak exogeneity is formulated at the structural level, that is, for entire row or column components such as industrial production, inflation expectations, or country effects. By contrast, in the standard VECM it can only be tested series by series. Hence, one cannot test whether inflation expectations as a whole are weakly exogenous, but only whether each of $\pi_{e,B}$, $\pi_{e,D}$, and $\pi_{e,A}$ is weakly exogenous separately.
Overall, the comparison highlights the main empirical advantage of the C-MAR approach: by preserving the matrix structure of the data, it yields a much more coherent and economically interpretable description of both the long-run relations and the adjustment dynamics than the unrestricted VECM.
This paper extends the analysis of matrix autoregressive models to cointegrated systems, thereby contributing to the growing literature on models for matrix-valued observations.
We argue that the use of a matrix-based model, rather than a standard vector autoregressive framework, is motivated by two main advantages: a reduction in estimation complexity and a clearer interpretation of the dynamics through the separation of row-wise and column-wise effects. Although previous contributions on cointegration in MAR settings addressed the first objective, they did not provide cointegration coefficients with a clear and economically meaningful interpretation within the MAR framework.
To fill this gap, we propose a novel cointegrated MAR model featuring parameters with direct economic content: cointegration matrices, which describe the long-run equilibrium relations toward which the system is attracted, and adjustment matrices, which measure the speed at which deviations from equilibrium are corrected. The presence of interpretable coefficients makes it possible to formulate economically meaningful restrictions and to test them using standard statistical procedures.
In contrast to earlier work, we do not start from an error-correction specification whose coefficients admit a Kronecker-product decomposition. Instead, we begin with a MAR representation and then derive the corresponding error-correction form. A key strength of our approach is that the matrix structure is preserved in both the autoregressive and the error-correction representations. By contrast, when the error-correction representations proposed in earlier contributions are vectorized, the resulting autoregressive coefficients generally do not admit a Kronecker-product decomposition.