Impulse Response Inference for Matrix Time Series

Midwest Econometrics Group

Ivan Ricardo

Maastricht University

Alain Hecq

Maastricht University

Ines Wilms

Maastricht University

October 10, 2026

Data Structures

Univariate \((y_t)\)

\[ \begin{bmatrix} \bullet \\ \bullet \\ \bullet \end{bmatrix} \]

Multivariate \((\mathbf{y}_t)\)

\[ \begin{bmatrix} \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet \end{bmatrix} \]

Matrix-valued \((\mathbf{Y}_t)\)

Examples of matrix-valued data include

Our Setting

  • Monthly year-on-year HICP inflation, December 2001 to December 2019
  • How does an energy shock pass through into food and core inflation, and does it differ across countries?

Why Matrix-valued data?

  • Vectorizing loses dimension-specific interpretations, e.g., which country, which category
  • Allows for more parsimonious models
Coefficients (\(N_1 = 5\), \(N_2 = 3\), \(p = 2\))
VAR(\(p\)) on \(\operatorname{vec}(\mathbf{Y}_t)\) \((N_1 N_2)^2 p = 450\)
MAR(\(p\)) \((N_1^2 + N_2^2 - 1)p = 66\)

What is missing

MAR impulse responses are reported as point estimates, or with residual bootstrap bands whose validity under the Kronecker restriction has not been established.

Main Contributions

  • Asymptotic inference. Joint asymptotic distribution of the MAR(\(p\)) coefficient and covariance estimators, giving closed-form intervals
  • Bootstrap inference. ProBAB-MAR, a bias-corrected bootstrap whose bootstrap samples are always generated by a MAR, with asymptotically correct coverage
  • Structure. The rank of MAR impulse responses, and how proportional co-movement relaxes with the horizon

In simulations ProBAB intervals attain near-nominal coverage, with bands much narrower than the unrestricted VAR

From VAR to MAR

Matrix Autoregressive Model (MAR)

The MAR(\(p\)) of Chen et al. (2021) \[ \mathbf{Y}_t = \sum_{j=1}^{p} \mathbf{A}_j \mathbf{Y}_{t-j}\mathbf{B}_j^\top + \mathbf{U}_t, \qquad \operatorname{Cov}(\operatorname{vec}(\mathbf{U}_t)) = \boldsymbol{\Sigma}_2 \otimes \boldsymbol{\Sigma}_1 \]

Vectorizing gives a restricted VAR: \(\mathbf{y}_t = \sum_{j=1}^{p} \underbrace{(\mathbf{B}_j \otimes \mathbf{A}_j)}_{\mathbf{C}_j}\mathbf{y}_{t-j} + \mathbf{u}_t\)

A1 (Stability). All companion eigenvalues lie strictly inside the unit circle.

A2 (Innovations). \(\mathbf{U}_t\) is white noise with covariance \(\boldsymbol{\Sigma}_2 \otimes \boldsymbol{\Sigma}_1\), and \(\boldsymbol{\Sigma}_1, \boldsymbol{\Sigma}_2 \succ 0\).

Impulse Responses

Under A1–A2, \(\mathbf{y}_t\) has the moving average representation \[ \mathbf{y}_{t} = \sum_{h=0}^{\infty} \boldsymbol{\Theta}_h \mathbf{u}_{t-h}, \qquad \boldsymbol{\Theta}_0 = \mathbf{I}_{N_1 N_2}, \quad \boldsymbol{\Theta}_h = \sum_{j=1}^{\min(p,h)} (\mathbf{B}_j \otimes \mathbf{A}_j)\, \boldsymbol{\Theta}_{h-j} \]

A shock to entry \((r,s)\) of \(\mathbf{U}_t\) is \(\mathbf{e}_k = \mathbf{e}_s \otimes \mathbf{e}_r\), and its response at each horizon is again a matrix, \[ \operatorname{unvec}(\boldsymbol{\Theta}_h \mathbf{e}_k) \in \mathbb{R}^{N_1 \times N_2} \quad \text{(countries} \times \text{categories)} \]

Structural shocks: \(\theta_{ik,h} = \mathbf{e}_i^\top \boldsymbol{\Theta}_h \mathbf{P}\,\mathbf{e}_k\), with a separable impact matrix \(\mathbf{P} = \mathbf{P}_2 \otimes \mathbf{P}_1\) from \(\boldsymbol{\Sigma}_1 = \mathbf{P}_1\mathbf{P}_1^\top\) and \(\boldsymbol{\Sigma}_2 = \mathbf{P}_2\mathbf{P}_2^\top\)

The Response Inherits the Structure

For a MAR(\(1\)), \(\boldsymbol{\Theta}_h = \mathbf{B}_1^h \otimes \mathbf{A}_1^h\), so a shock to entry \((r,s)\) gives \[ \operatorname{unvec}(\boldsymbol{\Theta}_h \mathbf{e}_k) = \underbrace{(\mathbf{A}_1^h\mathbf{e}_r)}_{\text{row response}} \;\underbrace{(\mathbf{B}_1^h\mathbf{e}_s)^\top}_{\text{column response}} \]

A rank-one matrix where every country responds with the same pattern across categories, scaled by a country-specific factor (Chen et al. 2021).

Proposition 1: MAR(\(p\))

Under A1-A2, a Kronecker shock \(\mathbf{v} = \mathbf{v}_2 \otimes \mathbf{v}_1\) (e.g. \(\mathbf{e}_k\)), has \(\operatorname{rank}\big(\operatorname{unvec}(\boldsymbol{\Theta}_h \mathbf{v})\big) \le \min\{c_h, N_1, N_2\}\), where \(c_0 = 1\) and \(c_h = \sum_{j=1}^{\min(p,h)} c_{h-j}\).

Proportional co-movement is relaxed gradually with the horizon.

Inference for MAR Impulse Responses

Joint Asymptotic Normality

A3 (i.i.d. innovations). A2 holds, the \(\mathbf{U}_t\) are i.i.d., and \(\mathbb{E}\|\mathbf{U}_t\|_F^4 < \infty\).

Proposition 2

Let A1 and A3 hold, and let \(\widehat{\boldsymbol{\zeta}}\) be the Gaussian QMLE of \(\boldsymbol{\zeta} = (\boldsymbol{\beta}^\top, \boldsymbol{\sigma}^\top)^\top\), with coefficients \(\boldsymbol{\beta}\) and covariance parameters \(\boldsymbol{\sigma}\). Then \[ \sqrt{T}(\widehat{\boldsymbol{\zeta}} - \boldsymbol{\zeta}) \xrightarrow{d} N(\mathbf{0}, \boldsymbol{\Omega}), \qquad \boldsymbol{\Omega} = \operatorname{blkdiag}(\boldsymbol{\Omega}_\beta, \boldsymbol{\Omega}_\sigma) \]

  • The coefficient block is known (Chen et al. 2021), we provide the covariance block and the joint law
  • The covariance block depends on fourth moments of \(\mathbf{U}_t\)

From Parameters to Responses

An identified response depends on both blocks: \[ \theta_{ik,h} = g_{ik,h}(\boldsymbol{\zeta}) = \mathbf{e}_i^\top \boldsymbol{\Theta}_h(\boldsymbol{\beta})\,\mathbf{P}(\boldsymbol{\sigma})\,\mathbf{e}_k \]

A4 (Non-degeneracy). \(\nabla g_{ik,h}(\boldsymbol{\zeta}) \neq \mathbf{0}\), so the limiting variance of \(\widehat{\theta}_{ik,h}\) is positive.

Under A1, A3 and A4, the delta method gives closed-form intervals with correct asymptotic coverage at any fixed \(h\).

Closed form, but centered at a biased estimate in small samples.

Where the BAB-MAR Stops Working

Bootstrap-after-bootstrap (Kilian 1998): bootstrap to estimate the bias, correct the coefficients, then bootstrap again from the corrected model.

  • \(\widehat{\mathbf{C}}_j = \widehat{\mathbf{B}}_j \otimes \widehat{\mathbf{A}}_j\) lies on the set of Kronecker products
  • The bias estimate \(\widehat{\boldsymbol{\Psi}}_j\) need not have this form
  • So \(\widehat{\mathbf{C}}_j - \widehat{\boldsymbol{\Psi}}_j\) is not a Kronecker product

Bootstrap design and estimator disagree

Samples are drawn from a process that is not a MAR, while every replication is re-estimated as one.

Projecting Back

Nearest Kronecker product (Van Loan and Pitsianis 1993) \[ \mathcal{P}(\mathbf{C}_j) = \underset{\mathbf{V}\otimes\mathbf{U}}{\arg\min}\;\|\mathbf{C}_j - \mathbf{V}\otimes\mathbf{U}\|_F^2 \]

A rearrangement defined by \(\mathcal{R}(\cdot)\) turns this into \[ \|\mathbf{C}_j - \mathbf{V}\otimes\mathbf{U}\|_F^2 = \|\mathcal{R}(\mathbf{C}_j) - \operatorname{vec}(\mathbf{U})\operatorname{vec}(\mathbf{V})^\top\|_F^2, \] a rank-one approximation found via SVD.

Order matters

The projection says nothing about eigenvalues and can push a root outside the unit circle. Project first, then check stability.

ProBAB-MAR

  • Stage 1. Bootstrap from \(\widehat{\mathbf{C}}\), \(B_1\) times. Estimate \(\widehat{\boldsymbol{\Psi}} = \bar{\mathbf{C}}^* - \widehat{\mathbf{C}}\) and set \[\widetilde{\mathbf{C}}_j = \mathcal{P}\big(\widehat{\mathbf{C}}_j - \delta\widehat{\boldsymbol{\Psi}}_j\big) = \widetilde{\mathbf{B}}_j \otimes \widetilde{\mathbf{A}}_j\]
  • Stage 2. Bootstrap from \(\widetilde{\mathbf{C}}\), \(B_2\) times. Re-fit, correct and project again, \(\widetilde{\mathbf{C}}^* = \mathcal{P}(\widehat{\mathbf{C}}^* - \delta^*\widehat{\boldsymbol{\Psi}})\), and compute the responses
  • Interval. Percentile, from the \(B_2\) draws

Structure preserved throughout

Every matrix that generates data, and every matrix estimated from it, factors as \(\mathbf{B}_j \otimes \mathbf{A}_j\).

Full algorithm

Asymptotic Validity

Recall. A1 stability, A2 separable white noise, A3 iid with fourth moments, A4 \(\nabla g_{ik,h}(\boldsymbol{\zeta}) \neq \mathbf{0}\)

Proposition 3

Let A1 and A3 hold. The plain residual bootstrap is consistent for the MAR: \(\sqrt{T}(\widehat{\boldsymbol{\zeta}}^* - \widehat{\boldsymbol{\zeta}}) \xrightarrow{d^*} N(\mathbf{0}, \boldsymbol{\Omega})\) in probability.

Bias correction, projection and shrinkage each perturb the bootstrap DGP by \(o_p(T^{-1/2})\), negligible against the \(T^{-1/2}\) sampling error.

Corollary

Let A1–A4 hold and fix \(h \ge 0\). Then the ProBAB-MAR percentile interval satisfies \[ \lim_{T\to\infty} P\Big(\theta_{ik,h}\in \big[\widetilde{\theta}_{ik,h}^{*(\alpha/2)},\,\widetilde{\theta}_{ik,h}^{*(1-\alpha/2)}\big]\Big) = 1-\alpha \]

Assumptions

Simulation Study

Calibrated to the Application

MAR(2) at the euro-area estimates, \(T = 200\). Response of Dutch food inflation to a unit innovation in energy inflation in all five countries at once

  • Delta method falls to \(80\%\); ProBAB-MAR near nominal at every horizon
  • BAB-VAR also covers, but its bands are on average \(2.57\times\) as wide

Coverage in Small Samples

MAR(1), \(N_1 = 3\), \(N_2 = 4\), \(\rho_{\max} = 0.8\), \(T = 100\)

  • Delta method falls to \(74\%\); ProBAB at the nominal level at every horizon
  • MAR bands about two thirds the width of BAB-VAR

Near the unit root

A Common Energy Shock

  • MAR(2) and VAR(1), each at the lag order chosen by BIC and Hannan-Quinn
  • Energy ordered first; the shock raises energy inflation by 1 pp on impact in all five countries
  • An energy innovation, not an oil supply shock, so it mixes supply and demand

Same lag order

What the Comparison Shows

  • Same story, tighter bands. Point estimates track closely with ProBAB-MAR bands being about \(1.4\times\) narrower (\(1.6\times\) against a VAR(2))
  • The VAR can barely bias-correct. Shrinkage binds in \(99\%\) of VAR replications against \(23\%\) for the MAR; on average \(18\%\) of the correction is applied, against \(92\%\)
  • Heterogeneity from \(h = 2\). At \(h = 1\) the response is rank one and differences across countries emerge only later, which a MAR(1) would rule out

Kronecker structure?

Conclusion

Conclusion

  • Joint asymptotic distribution of the MAR(\(p\)) coefficient and covariance estimators, with delta-method intervals for identified responses
  • ProBAB-MAR keeps the bias correction inside the Kronecker structure and attains asymptotically correct coverage, with bands narrower than the unrestricted VAR
  • MAR(\(p\)) responses relax proportional co-movement with the horizon

Github

References

Chen, Rong, Han Xiao, and Dan Yang. 2021. “Autoregressive Models for Matrix-Valued Time Series.” Journal of Econometrics, Annals Issue: Financial Econometrics in the Age of the Digital Economy, vol. 222 (1, Part B): 539–60.
Guggenberger, Patrik, Frank Kleibergen, and Sophocles Mavroeidis. 2023. “A Test for Kronecker Product Structure Covariance Matrix.” Journal of Econometrics 233 (1): 88–112.
Inoue, Atsushi, and Lutz Kilian. 2020. “The Uniform Validity of Impulse Response Inference in Autoregressions.” Journal of Econometrics 215 (2): 450–72.
Kilian, Lutz. 1998. “Small-Sample Confidence Intervals for Impulse Response Functions.” Review of Economics and Statistics 80 (2): 218–30.
Kleibergen, Frank, and Richard Paap. 2006. “Generalized Reduced Rank Tests Using the Singular Value Decomposition.” Journal of Econometrics 133 (1): 97–126.
Montiel Olea, José Luis, Mikkel Plagborg-Møller, Eric Qian, and Christian K. Wolf. 2026. “Local Projections or Vector Autoregressions? A Primer for Macroeconomists.” NBER Macroeconomics Annual 40: 111–52.
Pesaran, M. Hashem, Til Schuermann, and Scott M. Weiner. 2004. “Modeling Regional Interdependencies Using a Global Error-Correcting Macroeconometric Model.” Journal of Business & Economic Statistics 22 (2): 129–62.
Samadi, S. Yaser, and Lynne Billard. 2024. “On a Matrix-Valued Autoregressive Model.” Journal of Time Series Analysis.
Srivastava, Muni S, Tatjana von Rosen, and Dietrich Von Rosen. 2008. “Models with a Kronecker Product Covariance Structure: Estimation and Testing.” Mathematical Methods of Statistics 17 (4): 357–70.
Sung, Bongjung, and Peter D Hoff. 2025. “Testing Separability of High-Dimensional Covariance Matrices.” arXiv Preprint arXiv:2506.17463.
Van Loan, Charles F, and Nikos Pitsianis. 1993. “Approximation with Kronecker Products.” In Linear Algebra for Large Scale and Real-Time Applications. Springer.
Wang, Dong, Xialu Liu, and Rong Chen. 2019. “Factor Models for Matrix-Valued High-Dimensional Time Series.” Journal of Econometrics 208 (1): 231–48.

Appendix

Is the Kronecker Structure Right?

Back to main

Near the Unit Root

MAR(1), \(N_1 = 3\), \(N_2 = 4\), \(T = 100\), largest root \(\rho_{\max} = 0.99\)

  • MAR bootstrap coverage falls to \(89\%\) at worst, \(94\%\) for BAB-VAR
  • The shrinkage binds in \(36\%\) of replications, switching off the correction when the bias is largest
  • Lag augmentation (Inoue and Kilian 2020) for the MAR is open

Back to main

Same Lag Order in the Application

  • VAR(2) on the vectorized system: \(450\) coefficients against \(66\) for the MAR(2)
  • ProBAB-MAR bands narrower by \(1.6\times\) in mean and \(1.4\times\) in median
  • Against \(1.4\times\) and \(1.3\times\) for the VAR(1)

Back to main

Assumptions

  • A1 (Stability). All companion eigenvalues strictly inside the unit circle; for \(p = 1\), \(\rho(\mathbf{A}_1)\rho(\mathbf{B}_1) < 1\)
  • A2 (Innovations). White noise, \(\operatorname{Cov}(\operatorname{vec}(\mathbf{U}_t)) = \boldsymbol{\Sigma}_2 \otimes \boldsymbol{\Sigma}_1\), with \(\boldsymbol{\Sigma}_1, \boldsymbol{\Sigma}_2 \succ 0\)
  • A3 (i.i.d. innovations). A2, plus i.i.d. with \(\mathbb{E}\|\mathbf{U}_t\|_F^4 < \infty\)
  • A4 (Non-degeneracy). \(\nabla g_{ik,h}(\boldsymbol{\zeta}) \neq \mathbf{0}\)
  • Proposition 1 uses A1; Propositions 2 and 3 use A1 and A3; the coverage results add A4
  • Gaussianity is not required: the estimator is quasi-ML
  • A4 excludes the reduced-form response at \(h = 0\), since \(\boldsymbol{\Theta}_0 = \mathbf{I}\) is known, and entries a Cholesky scheme sets to zero
  • The guarantee is pointwise in the parameter, not uniform up to the unit circle

Back to main

ProBAB-MAR Algorithm

  1. Fit the MAR: \(\widehat{\mathbf{C}}_j = \widehat{\mathbf{B}}_j \otimes \widehat{\mathbf{A}}_j\) and residuals \(\widehat{\mathbf{u}}_t\)
  2. Stage 1. For \(m = 1, \ldots, B_1\): resample residuals, simulate from \(\widehat{\mathbf{C}}\) with zero initial values, re-fit to get \(\widehat{\mathbf{C}}^{*(m)}\)
  3. Bias \(\widehat{\boldsymbol{\Psi}} = \bar{\mathbf{C}}^* - \widehat{\mathbf{C}}\). Take the largest \(\delta \in \{1, 0.99, \ldots, 0.01\}\) with \(\rho\big(\mathcal{P}(\widehat{\mathbf{C}} - \delta\widehat{\boldsymbol{\Psi}})\big) < 1\) and set \(\widetilde{\mathbf{C}} = \mathcal{P}(\widehat{\mathbf{C}} - \delta\widehat{\boldsymbol{\Psi}})\); no correction if none qualifies
  4. Stage 2. For \(m = 1, \ldots, B_2\): resample residuals, simulate from \(\widetilde{\mathbf{C}}\), re-fit, compute \(\delta^*\), set \(\widetilde{\mathbf{C}}^* = \mathcal{P}(\widehat{\mathbf{C}}^* - \delta^*\widehat{\boldsymbol{\Psi}})\), re-estimate \(\boldsymbol{\sigma}^*\) from its residuals, and compute \(\widetilde{\theta}^{*(m)}_{ik,h} = \mathbf{e}_i^\top \widetilde{\boldsymbol{\Theta}}^*_h \mathbf{P}^* \mathbf{e}_k\)
  5. Percentile interval from \(\{\widetilde{\theta}^{*(m)}_{ik,h}\}_{m=1}^{B_2}\)

Reusing \(\widehat{\boldsymbol{\Psi}}\) in Stage 2 needs \(B_1 + B_2\) re-estimations instead of \(B_1 B_2\).

Back to main