---
title: "Workbook 12 — Matrix Filters and VARMA"
subtitle: "STA 542 · Introduction to Time Series Analysis"
format:
  html:
    embed-resources: true
    toc: true
    toc-depth: 2
    toc-location: right
    number-sections: false
    code-copy: true
    fig-width: 8
    fig-height: 3.2
    fig-dpi: 150
execute:
  warning: false
  message: false
editor:
  markdown:
    wrap: 72
---

::: hidden
$$\def\E{\mathbb{E}}
\def\P{\mathbb{P}}
\def\R{\mathbb{R}}
\def\Z{\mathbb{Z}}
\def\Var{\operatorname{Var}}
\def\Cov{\operatorname{Cov}}$$
:::

```{=html}
<style>
:root { --wb-ink: #173753; --wb-accent: #a35f32; }
#quarto-document-content { max-width: 820px; }
body { color: #25313b; line-height: 1.62; }
h1, h2, h3 { color: var(--wb-ink); }
h1.title { font-weight: 650; letter-spacing: -.025em; }
h2 { margin-top: 2.15rem; padding-bottom: .4rem; border-bottom: 1px solid #dce4ea; }
.subtitle { color: #607080; }
a { color: #145784; }
.exercise { border-left: 3px solid var(--wb-ink); padding: .9em 1.2em; background: #f3f6f8; margin: 1.5em 0; border-radius: 0 5px 5px 0; }
.exercise p:first-child { margin-top: 0; }
figure { margin: 1.25rem 0 1.6rem; }
figure figcaption { color: #586572; font-size: .88em; line-height: 1.45; }
details.code-fold > summary { color: #586572; font-size: .88em; }
.callout-tip { --bs-callout-color: #173753; }
table { font-size: .94em; }
@media print {
  #TOC, .quarto-title-meta, .code-copy-button, details.code-fold { display: none !important; }
  body { font-size: 10pt; }
  h2 { break-after: avoid; }
  figure, table, .exercise { break-inside: avoid; }
}
</style>
```

## More Lags in the State

A VAR(1) describes dependence on the most recent vector observation.
Dependence on several lags can be handled by putting those lagged vectors
in the state. Write $Y_t=X_t-\mu$ for the centered observations.

**Definition 12.1 (VAR($p$)).** For an integer $p\ge1$, a vector
autoregression of order $p$ satisfies
$$
Y_t=A_1Y_{t-1}+\cdots+A_pY_{t-p}+Z_t,
\qquad (Z_t)\sim\mathrm{WN}_d(0,\Sigma_Z),
$$
where $A_1,\ldots,A_p$ are real $d\times d$ matrices and $\Sigma_Z$ is
positive definite.

Stack $p$ successive vectors into $S_t$ and define the matrix $J$ that
selects its first block:
$$
S_t=\begin{pmatrix}Y_t\\Y_{t-1}\\\vdots\\Y_{t-p+1}\end{pmatrix},
\qquad J=\begin{pmatrix}I_d&0&\cdots&0\end{pmatrix}.
$$
Then $Y_t=JS_t$ and
$$
S_t=FS_{t-1}+J^\top Z_t,\qquad
F=\begin{pmatrix}
A_1&A_2&\cdots&A_{p-1}&A_p\\
I_d&0&\cdots&0&0\\
0&I_d&\cdots&0&0\\
\vdots&&\ddots&&\vdots\\
0&0&\cdots&I_d&0
\end{pmatrix}.
$$
This *companion matrix* gives a VAR(1) in $\R^{dp}$; for $p=1$, $F=A_1$.
The lower block rows shift the stored observations back one position.

If $\rho(F)<1$, the geometric-series argument in Workbook 11 gives the
unique weakly stationary solution
$$
S_t=\sum_{j\ge0}F^jJ^\top Z_{t-j},\qquad
Y_t=\sum_{j\ge0}\Psi_jZ_{t-j},\qquad \Psi_j:=JF^jJ^\top.
$$
The state innovation has covariance $J^\top\Sigma_ZJ$, which is singular
when $p>1$. The convergence and covariance proof needs only a finite,
positive-semidefinite covariance, so it still applies. The condition
checks the eigenvalues of $F$, which incorporates all lag matrices.

## Moving Averages in Multiple Dimensions

A scalar moving average combines the current innovation with a fixed
number of past innovations. For vector observations, matrix coefficients
also allow an innovation in one coordinate to contribute to another.

**Definition 12.2 (VMA($q$)).** For a nonnegative integer $q$, a *vector
moving average of order $q$* has the form
$$
Y_t=Z_t+C_1Z_{t-1}+\cdots+C_qZ_{t-q},
$$
where $Y_t=X_t-\mu$, $(Z_t)\sim\mathrm{WN}_d(0,\Sigma_Z)$,
$C_1,\ldots,C_q$ are real $d\times d$ matrices, and $\Sigma_Z$ is
positive definite. For $q=0$, the sum of lagged terms is absent.

The entry $(C_j)_{ab}$ multiplies $Z_{b,t-j}$ in the equation for
$Y_{at}$. Thus row $a$ collects the contributions to coordinate $a$,
while column $b$ describes where innovation coordinate $b$ contributes.
A fixed innovation vector $Z_s$ contributes $Z_s$ to $Y_s$ and
$C_jZ_s$ to $Y_{s+j}$ for $1\le j\le q$. There is no term involving
$Z_s$ after time $s+q$.

Every choice of the coefficient matrices gives a weakly stationary
process. The sum is finite, and the white-noise assumptions give a
constant mean and covariances that depend only on the lag. In particular,
$\E[Y_t]=0$ and $\E[X_t]=\mu$. There is no stability condition to check
for this construction.

For example, take $q=1$, $\Sigma_Z=I_2$, and
$C_1=\begin{pmatrix}0&0\\1&0\end{pmatrix}$. The coordinate equations are
$$
Y_{1t}=Z_{1t},\qquad Y_{2t}=Z_{2t}+Z_{1,t-1}.
$$
The innovation $Z_{1t}$ enters the first coordinate at time $t$ and the
second coordinate at time $t+1$.

Since $\Sigma_Z=I_2$, scalar innovations with different coordinate or
time indices are uncorrelated. For each coordinate, observations at
different dates involve distinct scalar innovations, so its autocovariance
is zero at every nonzero lag. Nevertheless, $Y_{1t}$ and
$Y_{2,t+1}$ share $Z_{1t}$, giving
$$
\begin{aligned}
\Cov(Y_{2,t+1},Y_{1t})
&=\Cov(Z_{2,t+1}+Z_{1t},Z_{1t})\\
&=\Var(Z_{1t})=1.
\end{aligned}
$$
Studying the coordinates jointly reveals a lagged relationship that
their individual ACFs miss.

## Finite Linear Filters

A VMA is a finite linear filter applied to white noise. More generally,
the finite matrix filter from Workbook 10 transforms an input $V_t\in\R^d$
into
$$
W_t=\sum_{j=0}^{m}L_jV_{t-j},
$$
where $L_0,\ldots,L_m$ are fixed real $r\times d$ matrices. The output
lies in $\R^r$. If $V$ is weakly stationary, so is $W$, by Proposition 10.3.

For white-noise input $V_t=Z_t$, different innovation times are
uncorrelated.
For $0\le h\le m$, the shared innovation $Z_{t-j}$ has coefficient
$L_j$ in $W_t$ and $L_{j+h}$ in $W_{t+h}$. It appears in both sums
precisely when $0\le j\le m-h$. For $h>m$, there are no shared
innovation times. Hence
$$
\Gamma_W(h)=\sum_{j=0}^{m-h}L_{j+h}\Sigma_ZL_j^\top
\quad(0\le h\le m),\qquad
\Gamma_W(h)=0\quad(|h|>m).
$$
Negative lags follow from $\Gamma_W(-h)=\Gamma_W(h)^\top$.
For a VMA($q$), take $m=q$, $L_0=I_d$, and $L_j=C_j$ for $1\le j\le q$.

## Constructing a VARMA Process

Combining lagged observations with lagged innovations gives the vector
version of an ARMA equation.

**Definition 12.3 (VARMA($p,q$)).** A vector autoregressive moving-average
equation is
$$
Y_t=\sum_{i=1}^{p}A_iY_{t-i}+Z_t+\sum_{j=1}^{q}C_jZ_{t-j},
\qquad Y_t=X_t-\mu,
$$
where $(Z_t)\sim\mathrm{WN}_d(0,\Sigma_Z)$, $\Sigma_Z$ is positive
definite, the coefficients are real $d\times d$ matrices, and $p,q$ are
nonnegative integers. An empty sum is zero. Setting $q=0$ gives a VAR;
setting $p=0$ gives a VMA.

Consider first a VARMA(1,1),
$$
Y_t=AY_{t-1}+Z_t+CZ_{t-1}.
$$
Substituting for $Y_{t-1}$ gives
$$
Y_t=Z_t+(A+C)Z_{t-1}+ACZ_{t-2}+A^2Y_{t-2}.
$$
Continuing the substitution suggests the innovation representation
$$
Y_t=Z_t+(A+C)Z_{t-1}+A(A+C)Z_{t-2}+\cdots
    =\sum_{j\ge0}\Psi_jZ_{t-j},
$$
with $\Psi_0=I_d$ and $\Psi_j=A^{j-1}(A+C)$ for $j\ge1$. To use this series, we need a condition ensuring that it converges to a
stationary solution. For a VARMA($1,1$), the condition is $\rho(A)<1$. It ensures that $\Psi_j Z_{t-j} =  A^{j-1} \eta_{t-j+1}$ for the weakly stationary series $\eta_t = (A+C)Z_{t-1}$ does not explode.



For the general VARMA($p,q$), the math is more complicated, but similar reasoning applies. We start by rewriting the VARMA($p,q$) equation in Definition 12.3. Consider, for a formal variable $z$, the
``matrix polynomials'' given by
$$
\Phi(z):=I_d-\sum_{i=1}^{p}A_i z^i,\qquad
M(z):=I_d+\sum_{j=1}^{q}C_j z^j.
$$
Replacing each power $z^j$ by the backshift $B^j$, where
$B^jY_t=Y_{t-j}$, gives the rewritten VARMA equation
$$
\Phi(B)Y_t=M(B)Z_t.
$$
Recognize that this is just rearranging the terms in the VARMA equation of Definition 12.3. If we could write 
$$
Y_t=\Phi(B)^{-1}M(B)Z_t,
$$
and make sense of the right-hand side series, we would have a causal weak stationary solution. Inverting $\Phi(B)$ is easier said than done however, recall that $B$ is a linear operator on sequences of random variables, and so is the linear filter $\Phi(B)$. 
Some sophisticated reasoning (some spectral theory, some operator algebra)  we can prove the proposition below, which allows us to check when a causal weak stationary solution exists.


**Proposition 12.4 (A stationary VARMA representation).** If
$\det\Phi(z)\ne0$ for every complex $z$ with $|z|\le1$, the VARMA
equation has a causal weakly stationary solution with mean $\mu$ and
unique mean-square convergent representation
$$
X_t-\mu=\sum_{j\ge0}\Psi_jZ_{t-j},\qquad
\Psi(z):=\sum_{j\ge0}\Psi_jz^j=\Phi(z)^{-1}M(z),
$$
where $\Psi_0=I_d$ and the coefficients are entrywise absolutely
summable. Its matrix ACVF is
$$
\Gamma(h)=\sum_{j=0}^{\infty}\Psi_{j+h}\Sigma_Z\Psi_j^\top
\quad(h\ge0),\qquad
\Gamma(-h)=\Gamma(h)^\top.
$$

The proposition says that, if the complex matrix polynomial $\Phi(z)$ with $z \in \mathbb{C}$ has ``no roots'' inside the unit circle (meaning that $\Phi(z)$ is invertible for all $z$ with $|z|\le1$), then the linear filter $\Phi(B)$ can effectively be inverted too: $\Phi(z)^{-1}M(z)$ can be written as a unique convergent series of matrices $\Psi_j z^j$, and these $\Psi_j$ are also the coefficients of the causal weak stationary solution. Wonderful! 


**Remark 12.5 (Polynomial roots and companion eigenvalues).** For a
VAR($p$), the zeros of $\det\Phi(z)$ are the reciprocals of the
nonzero eigenvalues of its companion matrix. Thus the condition in
Proposition 12.4 agrees with $\rho(F)<1$: companion eigenvalues lie inside
the unit circle and lag-polynomial roots lie outside it.

::: {.callout-note collapse="true"}
## Convergence and covariance ♠

Terms at different innovation times are uncorrelated. Thus, for $n>m$,
$$
\begin{aligned}
\E\!\left[\left\|\sum_{j=m+1}^{n}\Psi_jZ_{t-j}\right\|_2^2\right]
&=\sum_{j=m+1}^{n}\operatorname{tr}(\Psi_j\Sigma_Z\Psi_j^\top)\\
&\le\operatorname{tr}(\Sigma_Z)
  \sum_{j=m+1}^{n}\sum_{a,b}(\Psi_j)_{ab}^2
\longrightarrow0\qquad(m\longrightarrow\infty).
\end{aligned}
$$
Absolute summability makes each entry sequence square summable, so the
partial sums converge in mean square.

Each partial sum has mean zero. Passing to the mean-square limit
preserves its mean and covariances. The innovation $Z_{t-j}$ has
coefficients $\Psi_{j+h}$ in $X_{t+h}-\mu$ and $\Psi_j$ in $X_t-\mu$.
Matching these times gives the stated covariance sum, which depends only
on $h$. The constructed process is therefore weakly stationary with mean
$\mu$. The identity $\Phi(z)\Psi(z)=M(z)$ shows that its series satisfies
$\Phi(B)(X_t-\mu)=M(B)Z_t$.
:::

## A VARMA(1,1) in Detail

**Example 12.6 (A lagged innovation in the rotation model).** Consider
$$
Y_t=AY_{t-1}+Z_t+CZ_{t-1},\qquad Y_t=X_t-\mu,
$$
with the rotation from Workbook 11 and
$$
\mu=0,\qquad A=0.8R_{\pi/6},\qquad
C=\begin{pmatrix}0.3&0.4\\0&-0.2\end{pmatrix},\qquad
\Sigma_Z=\begin{pmatrix}1&0.6\\0.6&1\end{pmatrix}.
$$
Here $\rho(A)=0.8$, so the causal solution exists.
The matrices $A$ and $C$ do not commute, and the innovation
coordinates are positively correlated.

Use the stationary causal solution. Substituting $\Psi_0=I_2$ and
$\Psi_j=A^{j-1}(A+C)$ into Proposition 12.4 gives
$$
\Gamma(0)=\Sigma_Z+
\sum_{k=0}^{\infty}A^k(A+C)\Sigma_Z(A+C)^\top(A^k)^\top.
$$
For lag one, separate the $j=0$ term and use
$\Psi_{j+1}=A\Psi_j$ for $j\ge1$:
$$
\begin{aligned}
\Gamma(1)
&=(A+C)\Sigma_Z+
  A\sum_{j=1}^{\infty}\Psi_j\Sigma_Z\Psi_j^\top\\
&=(A+C)\Sigma_Z+A\{\Gamma(0)-\Sigma_Z\}\\
&=A\Gamma(0)+C\Sigma_Z.
\end{aligned}
$$

```{r}
#| label: wb12-varma-covariances
#| code-fold: true
#| code-summary: "R code for the covariance sum"
angle_cov <- pi / 6
A_cov <- 0.8 * matrix(c(cos(angle_cov), -sin(angle_cov),
                        sin(angle_cov), cos(angle_cov)), 2, byrow = TRUE)
C_cov <- matrix(c(0.3, 0.4, 0, -0.2), 2, byrow = TRUE)
Sigma_cov <- matrix(c(1, 0.6, 0.6, 1), 2)
Gamma0_cov <- Sigma_cov
Psi_cov <- A_cov + C_cov
# Add Psi_j Sigma_Z Psi_j^T for j = 1,...,200.
for (j_cov in seq_len(200L)) {
  Gamma0_cov <- Gamma0_cov + Psi_cov %*% Sigma_cov %*% t(Psi_cov)
  Psi_cov <- A_cov %*% Psi_cov
}
Gamma1_cov <- A_cov %*% Gamma0_cov + C_cov %*% Sigma_cov
```

Summing the variance contributions through $\Psi_{200}$ gives
$$
\Gamma(0)\approx
\begin{pmatrix}
`r sprintf('%.3f', Gamma0_cov[1,1])`&`r sprintf('%.3f', Gamma0_cov[1,2])`\\
`r sprintf('%.3f', Gamma0_cov[2,1])`&`r sprintf('%.3f', Gamma0_cov[2,2])`
\end{pmatrix},\qquad
\Gamma(1)\approx
\begin{pmatrix}
`r sprintf('%.3f', Gamma1_cov[1,1])`&`r sprintf('%.3f', Gamma1_cov[1,2])`\\
`r sprintf('%.3f', Gamma1_cov[2,1])`&`r sprintf('%.3f', Gamma1_cov[2,2])`
\end{pmatrix}.
$$

To recover innovations, solve $(I_2+CB)Z_t=(I_2-AB)Y_t$. Since
$\rho(C)=0.3<1$, expanding $(I_2+Cz)^{-1}(I_2-Az)$ gives
$$
\begin{aligned}
Z_t&=\sum_{j\ge0}\Pi_jY_{t-j},\\
\Pi_0&=I_2,\qquad
\Pi_j=(-C)^j-(-C)^{j-1}A\quad(j\ge1).
\end{aligned}
$$
In particular, $\Pi_1=-(A+C)$ and $\Pi_2=C^2+CA$.
Compare $\Psi_2=A^2+AC$: construction and recovery involve different
matrix product orders. $\diamond$

With observations only from time zero onward, start from a chosen
$\widetilde Z_0$ and compute
$$
\widetilde Z_t=Y_t-AY_{t-1}-C\widetilde Z_{t-1},\qquad t\ge1.
$$
Subtracting the model equation shows that
$$
\widetilde Z_t-Z_t=(-C)^t(\widetilde Z_0-Z_0).
$$
The initialization error decays under $\rho(C)<1$. A finite initialized
recursion approximates recovery from the complete past; it need not
recover the current innovation exactly.

## Granger causality in VARMA models

Let $Y_t=X_t-\mu\in\R^d$. Compare the conditional-mean forecasts of
$Y_{i,t+1}$ from the complete joint history $\mathcal H_t$ and from
$\mathcal H_t^{(-j)}$, which omits only coordinate $j$'s history.
Coordinate $j$ Granger-causes coordinate $i$ in mean if these forecasts
differ with positive probability for some $t$.

Coordinate $i$'s own history remains available. For $d>2$, the other
$d-2$ coordinates also remain available; this is Granger causality
conditional on those coordinates. For $d=2$, it is exactly the comparison
in Definition 11.6.

For the VARMA equation $\Phi(B)Y_t=M(B)Z_t$, use the stationary causal
solution with iid mean-zero innovation vectors of finite covariance.
Suppose the model also has an inverse filter with entrywise absolutely
summable coefficients:
$$
\Pi(z):=M(z)^{-1}\Phi(z)=I_d+\sum_{k\ge1}\Pi_kz^k,
\qquad
Z_t=Y_t+\sum_{k\ge1}\Pi_kY_{t-k}.
$$
This is the general form of the innovation recovery in Example 12.6.
Moving the lagged observations to the other side gives an autoregressive
representation entirely in observations:
$$
Y_{t+1}=-\sum_{k\ge1}\Pi_kY_{t+1-k}+Z_{t+1}.
$$

Causality makes observations through $t$ depend only on innovations
through $t$. Since the innovation vectors are iid across time,
$\E[Z_{t+1}\mid\mathcal H_t]=0$. All the lagged observations are in
$\mathcal H_t$, so the forecast of coordinate $i$ is
$$
\E[Y_{i,t+1}\mid\mathcal H_t]
=-\sum_{k\ge1}\sum_{a=1}^{d}(\Pi_k)_{ia}Y_{a,t+1-k}.
$$
Row $i$ gives the forecast of coordinate $i$. Column $j$ gives the
coefficients multiplying coordinate $j$'s past.

**Proposition 12.7 (Reading Granger causality from the inverse filter).**
Consider a stationary causal VARMA with iid mean-zero innovation vectors
of finite covariance and an entrywise absolutely summable inverse filter
$\Pi$ as above. Suppose the innovation vectors have a joint density on
$\R^d$. For $i\ne j$, coordinate $j$ Granger-causes coordinate $i$ in mean
at horizon one, conditional on the remaining coordinates, if and only if
$$
(\Pi_k)_{ij}\ne0\qquad\text{for at least one }k\ge1.
$$

Thus, first form $M(z)^{-1}\Phi(z)$, then inspect row $i$, column $j$ of
its lag coefficients. A nonzero entry gives the direction $j\to i$;
zeros at every lag give no such direction. Equivalently, check whether
the $(i,j)$ entry of the matrix function $M(z)^{-1}\Phi(z)$ is identically
zero. The coefficients combine the AR and MA terms before the direction
is read.

If all these entries are zero, the displayed forecast uses only the
retained coordinates. Adding coordinate $j$'s history cannot change it.
The joint-density assumption ensures that a nonzero contribution from
the omitted coordinate cannot already be determined by the retained
histories. Gaussianity is not required.

::: {.callout-note collapse="true"}
## Why a nonzero coefficient gives a forecast gain ♠

Let $r$ be the first lag with $(\Pi_r)_{ij}\ne0$, and put
$\tau=t+1-r$. Start from $\mathcal H_t^{(-j)}$ and supply the complete
joint past through $\tau-1$ as additional information. In the joint
forecast, every term is now known except
$-(\Pi_r)_{ij}Y_{j,\tau}$: earlier observations of coordinate $j$ have
been supplied, and its later observations have zero coefficients by the
choice of $r$.

The inverse filter recovers innovations through $\tau-1$ from the supplied
past. Conditional on this past, the block $Y_\tau,\ldots,Y_t$ is a
triangular affine transformation of $Z_\tau,\ldots,Z_t$, with identity
diagonal blocks. It therefore has a joint density. After conditioning on
all retained coordinates in this block, $Y_{j,\tau}$ still has a
conditional density, so it cannot be determined exactly. The joint
forecast is consequently not determined even after the extra past has
been supplied. It cannot be determined by the smaller history
$\mathcal H_t^{(-j)}$ either, so the two conditional-mean forecasts
differ with positive probability. $\square$
:::

**Example 12.8 (Reading the directions in a bivariate model).** Take the
earlier VMA(1), with iid $N_2(0,I_2)$ innovations:
$$
Y_t=Z_t+CZ_{t-1},\qquad
C=\begin{pmatrix}0&0\\1&0\end{pmatrix}.
$$
Here $\Phi(z)=I_2$ and $C^2=0$, so
$$
\Pi(z)=(I_2+Cz)^{-1}=I_2-Cz,\qquad
\Pi_1=\begin{pmatrix}0&0\\-1&0\end{pmatrix},\qquad
\Pi_k=0\quad(k\ge2).
$$
Row 2, column 1 is nonzero, so coordinate 1 Granger-causes coordinate 2.
Row 1, column 2 is zero at every lag, so coordinate 2 does not
Granger-cause coordinate 1.

Indeed, the autoregressive representation is
$$
Y_{1,t+1}=Z_{1,t+1},\qquad
Y_{2,t+1}=Y_{1t}+Z_{2,t+1}.
$$
The joint forecasts are therefore $0$ and $Y_{1t}$. From the second
coordinate's own history, the second forecast is zero: its observations
involve $Z_{2s}$ and $Z_{1,s-1}$ for $s\le t$, all independent of
$Y_{1t}=Z_{1t}$. Adding the first coordinate's history supplies the term
$Y_{1t}$ in its forecast.

The same rule handles cancellation. In the VARMA(1,1) of Example 12.6,
retain $A$, use iid mean-zero innovations with finite covariance and a
joint density, and set $C=-A$.
Then
$$
M(z)=\Phi(z)=I_2-Az,\qquad
\Pi(z)=M(z)^{-1}\Phi(z)=I_2.
$$
Every lag coefficient $\Pi_k$ is zero, so neither coordinate Granger-causes
the other. Although the original AR matrix has $A_{12}=-0.4$, the AR and
MA terms cancel in the inverse representation. $\diamond$


Workbook 13 takes up fitting these models: estimating their coefficients
and innovation covariances from observations. The later workbook on
multivariate forecasting uses recovered innovations to calculate
conditional-mean forecasts and prediction regions.

## Exercises

::: {.exercise}
**Exercise 12.1 (Stability with Several Lags).** Consider
$$
X_t=A_1X_{t-1}+A_2X_{t-2}+Z_t,\qquad
A_1=\begin{pmatrix}0.5&0.2\\0&0.4\end{pmatrix},\quad
A_2=\begin{pmatrix}-0.1&0\\0.1&-0.2\end{pmatrix},
$$
with $(Z_t)\sim\mathrm{WN}_2(0,I_2)$.

(a) Write the $4\times4$ companion matrix. Compute its eigenvalues in R
    and report their moduli. Does the model have a stable causal solution?

(b) A colleague checks only that $\rho(A_1)<1$. Explain why that does
    not establish stability. Illustrate the failure with $A_1=0$ and
    $A_2=1.21I_2$.
:::


::: {.exercise}
**Exercise 12.2 (Matrix Order and Innovation Recovery).** Use the matrices
$A,C,\Sigma_Z$ in Example 12.6.

(a) Compute $\Psi_2$ and $\Pi_2$ numerically. Show that replacing $AC$
    by $CA$ changes $\Psi_2$. Explain which multiplication occurs first
    when an innovation contributes through $AC$.

(b) Suppose two recovery recursions begin at time zero with values that
    differ by $\delta_0$. Derive their discrepancy at time $T$ and show
    that it tends to zero as $T$ grows. Why can the recursions give close
    values without either recovering $Z_T$ exactly?
:::


::: {.exercise}
**Exercise 12.3 (Following an Innovation).** Let $(Z_t)$ be iid
$N_2(0,I_2)$ and use the same innovation sequence in the three models
$$
\begin{aligned}
\text{VAR(1):}\qquad V_t&=AV_{t-1}+Z_t,\\
\text{VMA(1):}\qquad U_t&=Z_t+AZ_{t-1},\\
\text{VARMA(1,1):}\qquad W_t&=AW_{t-1}+Z_t+AZ_{t-1},
\end{aligned}
\qquad
A=\begin{pmatrix}1/2&0\\1/2&1/2\end{pmatrix}.
$$

(a) Write both coordinate equations for each model. Identify the terms
that use a previous observation and those that use a previous innovation.
Check that the VAR and VARMA have stationary causal solutions. Does the
VMA require a stability condition on $A$?

(b) Increase $Z_{1s}$ by one, keeping every other innovation and the
initial observations unchanged. Before simulating, calculate the change
in each observation vector at times $s,s+1,s+2,s+3$. For which models does
this change vanish completely after a finite number of steps? Explain
how a direct lag-one contribution can continue to affect later
observations in the VARMA.

(c) Use the following innovations, start the recursions from zero with
$Z_0=0$, and discard the first 200 observations. Plot both coordinates of
each model over the first 80 retained times, using a common vertical scale.

```{r}
#| label: wb12-model-comparison-input
set.seed(542)
n_cmp <- 2000L # retained sample size T
burn_cmp <- 200L
z_cmp <- matrix(rnorm(2L * (n_cmp + burn_cmp)), ncol = 2)
```

Copy `z_cmp`, increase its first coordinate at retained time 20 by one,
and rerun all three models. Keep every other innovation unchanged.
Subtract the original paths from the new paths and plot the changes in
both coordinates over lags $h=0,\ldots,10$. Compare with (b). Why does
using the same innovations after the change isolate its effect?

(d) For each model, use all $T=2000$ retained observations to estimate
$$
g(h):=\Cov(X_{2,t+h},X_{1t}),\qquad h=0,\ldots,6,
$$
with that model's process in place of $X$. Use
$$
\widehat g(h)
=\frac1T\sum_{t=1}^{T-h}
 (x_{2,t+h}-\bar x_2)(x_{1t}-\bar x_1).
$$
Plot the estimates against $h$. Which model's population covariances are
zero at every $h\ge2$? Explain this from its innovation contributions.
Why need its sample estimates not be exactly zero?

(e) From the VAR and VARMA paths, calculate
$$
R_t^V=V_t-AV_{t-1},\qquad R_t^W=W_t-AW_{t-1}.
$$
Write these vectors in terms of innovations and estimate
$\Cov(R_{2,t+1},R_{1t})$ for each. Why can the extra MA term in the VARMA
not simply be absorbed into a white-noise error after removing $AW_{t-1}$?
Use the results to explain what order one means in the AR and MA parts.
:::