---
title: "Workbook 11 — Vector Autoregression"
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>
```

## From Scalar to Vector Autoregression

A scalar AR(1) uses the previous observation to model the current one.
For a vector observation, each coordinate can depend on every coordinate
at the previous time. The coefficients therefore form a matrix. We retain
the stationarity convention from the workbook on multivariate series:
finite second moments, a constant mean $\E[X_t]=\mu$, and an ACVF
$\Gamma(h)=\Cov(X_{t+h},X_t)$ that depends only on the lag, with
$\Gamma(-h)=\Gamma(h)^\top$.

**Definition 11.1 (Vector white noise and VAR(1)).** A sequence
$(Z_t)_{t\in\Z}$ is *vector white noise*, written
$\mathrm{WN}_d(0,\Sigma_Z)$, if it has finite second moments and
$$
\E[Z_t]=0,\qquad \Cov(Z_t)=\Sigma_Z,\qquad
\Cov(Z_t,Z_s)=0\quad(t\ne s).
$$
A *vector autoregression of order one*, or *VAR(1)*, satisfies
$$
X_t-\mu=A(X_{t-1}-\mu)+Z_t,
\qquad A\in\R^{d\times d},\quad \mu\in\R^d.
$$
We take $\Sigma_Z$ to be positive definite. A solution is *causal* if
$X_t$ is a function of the current and past innovations
$(Z_t,Z_{t-1},\ldots)$.

Row $i$ of $A$ gives the equation for coordinate $i$: $A_{ij}$ multiplies
$X_{j,t-1}-\mu_j$. For example, the first equation of a bivariate VAR(1) is
$$
X_{1t}-\mu_1=A_{11}(X_{1,t-1}-\mu_1)
             +A_{12}(X_{2,t-1}-\mu_2)+Z_{1t}.
$$
The entry $A_{12}$ concerns a lagged value of the second coordinate.
An off-diagonal entry of $\Sigma_Z$ concerns innovation coordinates at the
same time. 

## Matrix Powers and Stationarity

Iterating the recursion $m$ times gives
$$
X_t-\mu=A^m(X_{t-m}-\mu)+\sum_{j=0}^{m-1}A^jZ_{t-j}.
$$
Here $A^0=I_d$, the identity matrix. The matrix $A^j$ describes how an
innovation contributes $j$ steps later. A stationary causal construction
requires these contributions to be small enough for their variances to sum.

An *eigenvalue* $\lambda$ of $A$ has a nonzero, possibly complex vector
$v$ with $Av=\lambda v$, so $A^jv=\lambda^jv$. The *spectral radius*
$\rho(A)$ is the largest modulus (absolute value) of an eigenvalue. We use this matrix fact
without proof: $\rho(A)<1$ if and only if the powers decay geometrically,
in the sense that some $C<\infty$ and $0<q<1$ satisfy
$$
\|A^jv\|_2\le Cq^j\|v\|_2\quad\text{for all real }v\text{ and }j\ge0.
$$
The Euclidean norm is $\|v\|_2=(v^\top v)^{1/2}$. This bound allows an
initial increase in length; it requires eventual decay. The criterion
concerns the whole matrix, and individual entries need not be below one.

**Proposition 11.2 (Stable VAR(1)).** If $\rho(A)<1$, Definition 11.1 has
a unique weakly stationary solution. It is causal and has the mean-square
convergent representation
$$
X_t=\mu+\sum_{j=0}^{\infty}A^jZ_{t-j}.
$$
Its mean is $\mu$, and its covariance matrices satisfy
$$
\begin{aligned}
\Gamma(0)&=\sum_{j=0}^{\infty}A^j\Sigma_Z(A^j)^\top
 =A\Gamma(0)A^\top+\Sigma_Z,\\
\Gamma(h)&=A^h\Gamma(0)\quad(h\ge0).
\end{aligned}
$$

*Proof.* In a weakly stationary solution the mean square of
$A^m(X_{t-m}-\mu)$ tends to zero by the geometric bound. The remaining
terms are uncorrelated, with summable second moments, so their series
converges in mean square. The resulting series solves the recursion and
has a constant mean and lag-dependent covariances. Iteration forces every
weakly stationary solution to equal this series.

The current $Z_t$ is uncorrelated with every past term in the series.
Taking covariance in the recursion gives
$\Gamma(0)=A\Gamma(0)A^\top+\Sigma_Z$ and
$\Gamma(1)=A\Gamma(0)$. Iteration gives $\Gamma(h)=A^h\Gamma(0)$;
opposite lags follow by transposition. $\square$

The inequality $\rho(A)<1$ is the causal stability condition. As with a
scalar AR(1), equations outside this regime can sometimes have stationary
solutions that depend on future noise. Such solutions are outside the
causal construction used here.

## Rotation and Decay

**Example 11.3 (Rotating innovations).** Consider the bivariate VAR(1)
$$
X_t=AX_{t-1}+Z_t,\qquad A=rR_\theta,\qquad
R_\theta=\begin{pmatrix}
\cos\theta&-\sin\theta\\
\sin\theta&\cos\theta
\end{pmatrix}.
$$
Here $X_t=(X_{1t},X_{2t})^\top$, $0<r<1$, $\theta\in\R$, and
$(Z_t)\sim\mathrm{WN}_2(0,\sigma^2I_2)$ with $\sigma^2>0$.
We use the stationary causal solution from Proposition 11.2, with mean zero.

Multiplication by $R_\theta$ rotates a vector counterclockwise through
$\theta$ without changing its length. Therefore $A^j=r^jR_{j\theta}$:
an innovation's contribution after $j$ periods is rotated through $j\theta$
and scaled by $r^j$. The eigenvalues $re^{i\theta}$ and $re^{-i\theta}$
have modulus $r$.

Since $R_{j\theta}R_{j\theta}^\top=I_2$, Proposition 11.2 gives
$$
\Gamma(0)=\frac{\sigma^2}{1-r^2}I_2,\qquad
\Gamma(h)=\frac{\sigma^2r^h}{1-r^2}R_{h\theta}\quad(h\ge0).
$$
The current coordinates are uncorrelated, but cross-lag covariances can be
nonzero. Decay makes the infinite covariance sum finite.

At the boundary $r=1$, start instead at $X_0=0$ and take iid
$Z_t\sim N_2(0,\sigma^2I_2)$. Past innovations retain their magnitude,
and $\Cov(X_t)=t\sigma^2I_2$. If $R_\theta=I_2$, this is an ordinary
bivariate random walk. Otherwise,
$\nabla X_t=(R_\theta-I_2)X_{t-1}+Z_t$ has growing variance,
so $X_t$ is not $I(1)$ under our definition.

The time-varying matrix filter $Y_t:=R_{-t\theta}X_t$ removes the
rotation:
$$
Y_t-Y_{t-1}=R_{-t\theta}Z_t.
$$
These rotated innovations remain iid $N_2(0,\sigma^2I_2)$, so the filtered
process $(Y_t)$ is an ordinary bivariate random walk.

For the simulation and subsequent fit, take $r=0.8$, $\theta=\pi/6$,
and $\sigma^2=1$. Use iid $Z_t\sim N_2(0,I_2)$ and generate $T=300$
observations, initializing
$X_0\sim N_2(0,I_2/(1-r^2))$ independently of future innovations.
This gives the stationary Gaussian law from the first observation.
The figure displays the first 60 observations. $\diamond$

```{r}
#| label: fig-wb11-rotation
#| code-fold: true
#| code-summary: "R code for the rotation and stationary sample"
#| fig-cap: "Left: calculated contributions of the innovation (1, 0) at ages 0 through 10; arrows follow increasing age. Right: the first 60 vector observations from the stationary rotation model (simulated data, seed 54211). Each coordinate is measured in units of the innovation standard deviation."
set.seed(54211)
r <- 0.8
theta <- pi / 6
rotation <- function(angle) {
  matrix(c(cos(angle), sin(angle), -sin(angle), cos(angle)), 2, 2)
}
A <- r * rotation(theta)
ages <- 0:10
shock <- t(vapply(ages, function(j) {
  as.numeric(r^j * rotation(j * theta) %*% c(1, 0))
}, numeric(2)))
n_obs <- 300L
x_all <- matrix(0, n_obs + 1L, 2)
x_all[1, ] <- rnorm(2, sd = 1 / sqrt(1 - r^2))
z <- matrix(rnorm(2 * n_obs), ncol = 2)
for (t in seq_len(n_obs)) x_all[t + 1L, ] <- A %*% x_all[t, ] + z[t, ]
x <- x_all[-1, , drop = FALSE] # rows are times; columns are coordinates
par(mfrow = c(1, 2), mar = c(3.6, 3.8, 2.3, 0.8),
    mgp = c(2.35, 0.65, 0), las = 1, bg = "white")
plot(shock, type = "n", asp = 1, xlim = c(-1.05, 1.15),
     ylim = c(-1.05, 1.15), xlab = "First coordinate",
     ylab = "Second coordinate", main = "A past innovation")
abline(h = 0, v = 0, col = "grey85")
arrows(shock[-nrow(shock), 1], shock[-nrow(shock), 2],
       shock[-1, 1], shock[-1, 2], length = 0.065, col = "#1C5276")
points(shock, pch = 19, cex = 0.65, col = "#1C5276")
text(shock[c(1, 4, 7, 11), ], labels = ages[c(1, 4, 7, 11)],
     pos = c(3, 3, 2, 4), cex = 0.8)
matplot(1:60, x[1:60, ], type = "n", xlab = "Time index",
        ylab = "Value", main = "Stationary sample")
usr <- par("usr")
rect(usr[1], usr[3], usr[2], usr[4], col = "grey94", border = NA)
abline(h = 0, col = "white", lwd = 2)
matlines(1:60, x[1:60, ], lty = c(1, 2), lwd = 1.2,
         col = c("#1C5276", "#AA5F32"))
legend("topright", c("Coordinate 1", "Coordinate 2"), bty = "n",
       lty = c(1, 2), col = c("#1C5276", "#AA5F32"), cex = 0.72)
par(mfrow = c(1, 1))
```

## Fitting from Covariances

Suppose we observe $X_1,\ldots,X_T$ from the VAR(1) model
$$
X_t-\mu=A(X_{t-1}-\mu)+Z_t,\qquad
(Z_t)\sim\mathrm{WN}_d(0,\Sigma_Z).
$$
We assume $\rho(A)<1$ and use the stationary causal solution. The unknown
parameters are the mean vector $\mu$, the coefficient matrix $A$, and the
positive-definite innovation covariance matrix $\Sigma_Z$.

We want to estimate how each current coordinate depends on the preceding
vector. The covariance between consecutive vectors gives a way to recover
$A$.

Proposition 11.2 gives
$$
\Gamma(1)=A\Gamma(0),\qquad
\Gamma(0)=A\Gamma(0)A^\top+\Sigma_Z.
$$
The series for $\Gamma(0)$ includes the positive-definite term
$\Sigma_Z$, so $\Gamma(0)$ is invertible. Multiplying the first identity
on the right by $\Gamma(0)^{-1}$ and rearranging the second gives
$$
A=\Gamma(1)\Gamma(0)^{-1},\qquad
\Sigma_Z=\Gamma(0)-A\Gamma(0)A^\top.
$$

Let $\bar X=T^{-1}\sum_{t=1}^T X_t$ and define the sample covariance
matrices
$$
\widehat\Gamma(h)=\frac1T\sum_{t=1}^{T-h}
 (X_{t+h}-\bar X)(X_t-\bar X)^\top,\qquad h=0,1.
$$
The *moment estimators* replace population moments by these averages:
$$
\widehat\mu=\bar X,\qquad
\widehat A=\widehat\Gamma(1)\widehat\Gamma(0)^{-1},\qquad
\widehat\Sigma_Z=\widehat\Gamma(0)
 -\widehat A\widehat\Gamma(0)\widehat A^\top,
$$
provided $\widehat\Gamma(0)$ is invertible. If the sample mean and these
sample covariance matrices converge in probability to their population
counterparts, these estimators are consistent. Iid innovations with finite second moments
in the stable causal model provide one sufficient setting: the resulting
process is stationary and ergodic, so its integrable coordinate products
satisfy the time-average law.

For a realized record $x_1,\ldots,x_T$, substituting $x_t$ for $X_t$ in
these formulas gives the corresponding estimates. The simulated record is
stored as rows of `x`. The following code fits an unrestricted matrix $A$,
without imposing the rotation form.

```{r}
#| label: fit-var1
n_obs <- nrow(x)
mu_hat <- colMeans(x)
u <- sweep(x, 2, mu_hat)
G0_hat <- crossprod(u) / n_obs
G1_hat <- t(u[-1, , drop = FALSE]) %*%
  u[-n_obs, , drop = FALSE] / n_obs
A_hat <- G1_hat %*% solve(G0_hat)
Sigma_hat <- G0_hat - A_hat %*% G0_hat %*% t(A_hat)
rho_hat <- max(Mod(eigen(A_hat)$values))
```

```{r}
#| label: tab-wb11-fit
#| echo: false
fit_table <- data.frame(
  Entry = c("A11", "A12", "A21", "A22"),
  Truth = as.vector(t(A)),
  Estimate = as.vector(t(A_hat))
)
knitr::kable(fit_table, digits = 3,
             caption = "Moment estimates from the 300 simulated observations.")
```

The estimated spectral radius is `r sprintf('%.3f', rho_hat)`, compared
with the true value $0.8$. The fitted recurrence is therefore stable for
this sample. Stability is assessed from the eigenvalues of the fitted
matrix, rather than from its individual entries.

A regression fit estimates the equivalent equation
$X_t=c+AX_{t-1}+Z_t$, where $c=(I_d-A)\mu$. Regress each current
coordinate on an intercept and all lagged coordinates, using $t=2,\ldots,T$;
the slopes form the corresponding row of $A$. These least-squares estimates
need not equal the moment estimates in a finite sample: the regressions use
paired lagged observations and their own sample means. Correlation between
the current innovation coordinates is allowed in either approach.

## Vector and Scalar Forecasts

Return to the stationary causal model with known parameters and iid
innovations of mean zero and covariance $\Sigma_Z$. Let
$\mathcal H_T=(X_1,\ldots,X_T)$ be the observed vector history.
The causal solution uses only current and past innovations. The iid
assumption therefore makes future innovations independent of this history.
The squared Euclidean error adds the squared coordinate errors. Its
optimal forecast is the conditional mean
$m_h=\E[X_{T+h}\mid\mathcal H_T]$, because any forecast $a$ computed from
this history satisfies
$$
\E[\|X_{T+h}-a\|_2^2\mid\mathcal H_T]
=\E[\|X_{T+h}-m_h\|_2^2\mid\mathcal H_T]+\|m_h-a\|_2^2.
$$
This is the scalar squared-error decomposition summed over coordinates.

**Proposition 11.4 (VAR(1) forecasts).** Under these assumptions, for
integers $h\ge1$,
$$
\widehat X_{T+h\mid T}=\mu+A^h(X_T-\mu),\qquad
V_h:=\Cov(X_{T+h}-\widehat X_{T+h\mid T}\mid\mathcal H_T)
 =\sum_{j=0}^{h-1}A^j\Sigma_Z(A^j)^\top.
$$
*Proof.* Forward iteration writes $X_{T+h}$ as the displayed forecast
plus $\sum_{j=0}^{h-1}A^jZ_{T+h-j}$. These future innovations are
independent of the observed history and each other, so the sum has
conditional mean zero and covariance $V_h$. $\square$

The vector forecast has conditional mean squared Euclidean error
$\operatorname{tr}(V_h)$, the sum of the diagonal entries. For a fixed
vector $b\in\R^d$, the scalar target $b^\top X_{T+h}$ has forecast
$b^\top\widehat X_{T+h\mid T}$ and conditional mean squared error
$b^\top V_hb$. Taking $b=e_i$, the vector with entry one at position $i$
and zero elsewhere, gives the forecast of coordinate $i$ and error
variance $(V_h)_{ii}$.

In Example 11.3 the one-step forecast of the first coordinate is
$$
\widehat X_{1,T+1\mid T}
=r\cos\theta\,X_{1T}-r\sin\theta\,X_{2T}.
$$
For $\theta=\pi/6$, both current coordinates enter even though the
target is scalar. The
choice of target specifies what is predicted; the available history
specifies what information can be used.

With fitted parameters, the corresponding plug-in forecast is
$\widehat\mu+\widehat A^h(x_T-\widehat\mu)$. Substitution in the formula
for $V_h$ estimates innovation uncertainty; it does not account for
uncertainty in the estimated parameters. Exercise 11.2 computes both the
known-parameter and fitted forecasts in the rotation example.

## Granger Causality

A signal observed now may help predict a response recorded one period
later. Does the response's history also help predict the next signal
value? The two directions need not give the same answer.

**Example 11.5 (A delayed noisy response).** Let $(Z_t)$ be iid
$N_2(0,I_2)$ and, for $t\in\Z$, define
$$
X_{2t}=Z_{2t},\qquad
X_{1t}=cX_{2,t-1}+Z_{1t},\qquad c\ne0.
$$
The second coordinate is the signal; the first is its delayed response
with independent noise added. Each coordinate separately is an iid
sequence: $X_{2t}\sim N(0,1)$ and $X_{1t}\sim N(0,1+c^2)$. The pair
is a stationary VAR(1), with
$$
A=\begin{pmatrix}0&c\\0&0\end{pmatrix},\qquad \Sigma_Z=I_2.
$$

Does observing the second series improve the forecast of $X_{1,t+1}$
obtained from the first series's own past? Write the first series's
complete history and the complete joint history through time $t$ as
$$
\mathcal H_t^{(1)}=(\ldots,X_{1,t-1},X_{1t}),\qquad
\mathcal H_t=(\ldots,X_{t-1},X_t).
$$

The next response is $X_{1,t+1}=cZ_{2t}+Z_{1,t+1}$. Both shocks
are independent of $\mathcal H_t^{(1)}$, whose observations use signal
shocks only through time $t-1$. Consequently,
$$
\E[X_{1,t+1}\mid\mathcal H_t^{(1)}]=0,\qquad
\E[(X_{1,t+1}-0)^2]=1+c^2.
$$

The joint history also reveals $X_{2t}=Z_{2t}$. Only the response noise
$Z_{1,t+1}$ remains unknown, so
$$
\E[X_{1,t+1}\mid\mathcal H_t]=cX_{2t},\qquad
\E[(X_{1,t+1}-cX_{2t})^2]=1.
$$
Adding the signal history reduces the mean squared prediction error by
$c^2$, although the response has no serial dependence of its own.

In the reverse direction, $X_{2,t+1}=Z_{2,t+1}$ is independent of the
entire joint history through time $t$. Its best forecast is zero, with
mean squared error one, whether the response history is supplied or
not. The predictive improvement therefore goes in only one direction.

For the figure, take $c=2$ and simulate 30 responses together with the
signal values needed to generate them. The response forecasts use the
known model: zero from its own history and $2X_{2,t-1}$ from the joint
history when predicting $X_{1t}$. $\diamond$

```{r}
#| label: fig-wb11-granger
#| fig-height: 6
#| code-fold: true
#| code-summary: "R code for the signal, response, and forecasts"
#| fig-cap: "Top: the signal. Bottom: the delayed noisy response and its one-step forecasts. A signal at time t contributes to the response at time t + 1; each plotted forecast of the response at time t uses observations through t − 1. The joint-history forecast follows the delayed signal, while the own-history forecast is zero. Values are in innovation standard deviations. Simulated data, c = 2, seed 54213."
local({
  set.seed(54213)
  n_gc <- 30L
  c_gc <- 2
  signal_gc <- rnorm(n_gc + 1L) # signal times 0,...,30
  noise_gc <- rnorm(n_gc)
  joint_gc <- c_gc * signal_gc[seq_len(n_gc)]
  response_gc <- joint_gc + noise_gc # response times 1,...,30
  old_par <- par(no.readonly = TRUE)
  on.exit(par(old_par))
  par(mfrow = c(2, 1), mar = c(4.2, 4.2, 2.1, 0.8),
      mgp = c(2.6, 0.7, 0), las = 1, bg = "white")
  panel_gc <- function(values, label, title) {
    plot(0, 0, type = "n", xlim = c(0, n_gc),
         ylim = extendrange(values, f = 0.25),
         xlab = "Time", ylab = label, main = title, xaxt = "n")
    axis(1, at = seq(0, n_gc, by = 5))
    bounds_gc <- par("usr")
    rect(bounds_gc[1], bounds_gc[3], bounds_gc[2], bounds_gc[4],
         col = "grey94", border = NA)
    abline(v = seq(0, n_gc, by = 5), col = "white")
  }
  panel_gc(signal_gc, expression(X[2*t]), "Signal")
  abline(h = 0, col = "grey70")
  lines(0:n_gc, signal_gc, type = "o", pch = 16, cex = 0.45,
        col = "#1C5276", lwd = 1.2)
  panel_gc(c(response_gc, joint_gc), expression(X[1*t]),
           "Response and one-step forecasts")
  abline(h = 0, col = "grey35", lty = 3, lwd = 1.4)
  lines(seq_len(n_gc), joint_gc, col = "#1C5276", lty = 2, lwd = 1.6)
  lines(seq_len(n_gc), response_gc, type = "o", pch = 16,
        cex = 0.45, col = "#AA5F32", lwd = 1.1)
  legend("topright",
         c("Observed response", "Joint-history forecast", "Own-history forecast"),
         col = c("#AA5F32", "#1C5276", "grey35"),
         lty = c(1, 2, 3), lwd = c(1.1, 1.6, 1.4),
         pch = c(16, NA, NA), pt.cex = 0.6, cex = 0.75,
         bg = "grey94", box.col = NA)
})
```

For a general bivariate stationary process with finite second moments,
use $\mathcal H_t^{(1)}$ for the complete first-coordinate history and
$\mathcal H_t$ for the complete joint history through time $t$.

**Definition 11.6 (Granger causality in mean).** The second coordinate
*Granger-causes* the first at horizon one if, for some $t$,
$$
\E[X_{1,t+1}\mid\mathcal H_t]
\ne\E[X_{1,t+1}\mid\mathcal H_t^{(1)}]
\quad\text{with positive probability}.
$$
It *does not Granger-cause* the first if the forecasts agree with
probability one for every $t$. The reverse direction exchanges the
coordinate labels.

For a fixed $t$, write the joint-history and own-history conditional means
as $m_{\rm joint}$ and $m_{\rm own}$. Their unconditional mean squared
errors satisfy
$$
\E[(X_{1,t+1}-m_{\rm own})^2]
-\E[(X_{1,t+1}-m_{\rm joint})^2]
=\E[(m_{\rm joint}-m_{\rm own})^2].
$$
The right side is positive exactly when the forecasts differ with
positive probability. Granger causality in mean is therefore equivalent
to a strict reduction in mean squared prediction error at some time $t$
when the second coordinate's history is added.

In Example 11.5, the second coordinate Granger-causes the first, while
the first does not Granger-cause the second. Observing two individually
white series together can therefore reveal predictive information.

To relate this definition to the coefficient matrix, return to
$X_{t+1}-\mu=A(X_t-\mu)+Z_{t+1}$, with $\rho(A)<1$ and iid mean-zero
innovations of finite covariance $\Sigma_Z$. Use its stationary causal
solution. Observations through $t$ then depend only on innovations through
$t$, so $Z_{t+1}$ is independent of $\mathcal H_t$ and has conditional
mean zero. The two forecasts are
$$
\begin{aligned}
m_{\rm joint}
&=\mu_1+A_{11}(X_{1t}-\mu_1)+A_{12}(X_{2t}-\mu_2),\\
m_{\rm own}
&=\mu_1+A_{11}(X_{1t}-\mu_1)
  +A_{12}\bigl(\E[X_{2t}\mid\mathcal H_t^{(1)}]-\mu_2\bigr).
\end{aligned}
$$
The second equality follows by conditioning the first forecast on
$\mathcal H_t^{(1)}$.

Subtracting gives
$$
m_{\rm joint}-m_{\rm own}
=A_{12}\bigl(X_{2t}-\E[X_{2t}\mid\mathcal H_t^{(1)}]\bigr).
$$
If $A_{12}=0$, the forecasts agree, so the second coordinate does not
Granger-cause the first. No Gaussian assumption is needed for this
implication. A nonzero $A_{12}$ gives a gain precisely when the current
$X_{2t}$ is not already determined by the first coordinate's complete
history. The following result gives sufficient assumptions for this
remaining condition.

**Proposition 11.7 (The bivariate Gaussian VAR(1)).** In the stable
causal VAR(1) with iid $N_2(0,\Sigma_Z)$ innovations and positive-definite
$\Sigma_Z$, the second coordinate Granger-causes the first in mean at
horizon one if and only if $A_{12}\ne0$. The reverse causality condition
is $A_{21}\ne0$.

*Proof.* We have already shown that $A_{12}=0$ gives no forecast gain.
Suppose $A_{12}\ne0$. Even if the full vector history through $t-1$ is
supplied along with $X_{1t}$, the first equation reveals only $Z_{1t}$,
leaving conditional variance
$\Sigma_{Z,22}-\Sigma_{Z,21}^2/\Sigma_{Z,11}$ in $Z_{2t}$ and hence
in $X_{2t}$. Predicting $X_{2t}$ from only the first coordinate's
history cannot have smaller mean squared error. Therefore
$$
\E\!\left[\bigl(X_{2t}-\E[X_{2t}\mid\mathcal H_t^{(1)}]\bigr)^2\right]
\ge \Sigma_{Z,22}-\frac{\Sigma_{Z,21}^2}{\Sigma_{Z,11}}>0.
$$
Gaussianity gives the conditional-variance formula, and positive
definiteness makes it strictly positive. The forecast difference thus
has positive variance, giving Granger causality. The reverse direction
follows by exchanging indices. $\square$

Granger causality concerns the predictive information in specified
histories. An effect of intervening on a coordinate requires additional
causal assumptions. Exercise 11.1 examines the covariance structure and
longer-horizon forecasts of Example 11.5. The workbook on matrix filters
and VARMA extends these dynamics to several lags and treats recovery of
innovations and joint prediction regions.

## Exercises

::: {.exercise}
**Exercise 11.1 (Predictive Information Across Coordinates).** Let
$(Z_t)$ be iid $N_2(0,I_2)$ and consider the stationary causal model
$$
X_t=AX_{t-1}+Z_t,\qquad
A=\begin{pmatrix}0&c\\0&0\end{pmatrix},\qquad c\ne0.
$$

(a) Compute $A^2$ and the eigenvalues of $A$. Show that
    $X_t=Z_t+AZ_{t-1}$. Can $|c|$ exceed one while the model remains stable?

(b) Find $\Gamma(0)$, $\Gamma(1)$, and $\Gamma(-1)$. Show that each coordinate
    separately is white noise, but the vector process is not white noise.
    *Hint: write each coordinate in terms of the independent scalar shocks.*

(c) For an integer $h\ge2$, find the best squared-error forecast of
    $X_{1,t+h}$ from the complete first-coordinate history and then from the complete joint
    history through time $t$. Compute the two mean squared errors.
    Does the one-step predictive improvement in Example 11.5 persist
    at these longer horizons? Explain.

(d) Find the eigenvalues of
    $C=\begin{pmatrix}0.6&0.6\\0.6&0.6\end{pmatrix}$.
    Together with part (a), explain why requiring all coefficient moduli
    to be below one is neither necessary nor sufficient for causal stability.
    *Hint: apply $C$ to $(1,1)^\top$ and $(1,-1)^\top$.*
:::


::: {.exercise}
**Exercise 11.2 (Forecasting the Rotation Model).** Use the Gaussian
rotation model of Example 11.3 and its simulated observations.

(a) Show that
    $V_h=\sigma^2(1-r^{2h})(1-r^2)^{-1}I_2$.
    Find its limit as $h$ increases and explain why it equals $\Gamma(0)$.

(b) Using the 300 observations already generated, compute the fitted
    forecasts $\widehat\mu+\widehat A^h(x_T-\widehat\mu)$ for $h=1,4$.
    Report them beside the known-parameter forecasts $A^hx_T$ in a table.
    *Hint: repeated matrix multiplication computes $\widehat A^h$;
    `A_hat^h` in R raises entries separately.*

(c) For the scalar target $X_{1,T+1}-X_{2,T+1}$, calculate the
    known-parameter forecast and its conditional mean squared error.
    Explain why both coordinates can matter for a scalar target and why
    the error formula need not equal the mean squared error of the
    fitted forecast.
:::


## Hare and Lynx Revisited

**Example 11.8 (A VAR(1) for log pelt records).** Return to the annual
hare and lynx records in Example 1.14. The observations are pelts
purchased by the Hudson's Bay Company, measured in thousands, for
1845–1935: 91 paired observations. They are proxies for animal
populations. The data are available as
[`astsa::Hare`](https://search.r-project.org/CRAN/refmans/astsa/html/Hare.html)
and [`astsa::Lynx`](https://search.r-project.org/CRAN/refmans/astsa/html/Lynx.html).

```{r}
#| label: fig-wb11-harelynx
#| fig-height: 3.4
#| code-fold: true
#| code-summary: "R code for the annual log pelt records"
#| fig-cap: "Annual hare and lynx pelt records on the natural-log scale. Before taking logs, both series are measured in thousands of pelts. Real data: astsa::Hare and astsa::Lynx, 1845–1935."
hl_pelts <- cbind(Hare = as.numeric(astsa::Hare),
                  Lynx = as.numeric(astsa::Lynx))
hl_years <- as.numeric(time(astsa::Hare))
stopifnot(identical(time(astsa::Hare), time(astsa::Lynx)),
          all(hl_pelts > 0))
hl_y <- log(hl_pelts)
hl_n <- nrow(hl_y)
local({
  old_par <- par(no.readonly = TRUE)
  on.exit(par(old_par))
  par(mar = c(4.2, 4.6, 0.8, 0.8), mgp = c(2.8, 0.8, 0),
      las = 1, bg = "white")
  matplot(hl_years, hl_y, type = "l", lty = 1, lwd = 1.3,
          col = c("#1C5276", "#AA5F32"),
          xlab = "Year", ylab = "Log(pelts in thousands)")
  legend("topright", c("Hare", "Lynx"),
         col = c("#1C5276", "#AA5F32"), lty = 1, lwd = 1.3,
         bty = "n", cex = 0.85)
})
```

Let $X_{H,t}$ and $X_{L,t}$ denote the annual hare and lynx pelt records
in thousands. Write $Y_t=(Y_{H,t},Y_{L,t})^\top$, where
$Y_{H,t}=\log X_{H,t}$ and $Y_{L,t}=\log X_{L,t}$. The paired model
from Workbook 1 is
$$
Y_{t+1}=a+AY_t+Z_{t+1},\qquad
a=\begin{pmatrix}a_H\\a_L\end{pmatrix},\qquad
A=\begin{pmatrix}b_{HH}&b_{HL}\\b_{LH}&b_{LL}\end{pmatrix}.
$$
As in Example 1.14, assume
$\E[Z_{t+1}\mid\mathcal H_t]=0$ and
$\Cov(Z_{t+1}\mid\mathcal H_t)=\Sigma_Z$, where $\mathcal H_t$ is
the complete log-pelt history and $\Sigma_Z$ is fixed.

Fit each row by least squares, using the 90 transitions from years
1845–1934 to years 1846–1935. Both equations use the same current pair
as predictors.

```{r}
#| label: fit-wb11-harelynx
#| code-fold: true
#| code-summary: "R code for the two fitted log-pelt equations"
hl_lagged <- data.frame(
  H_next = hl_y[-1, 1], L_next = hl_y[-1, 2],
  H_now = hl_y[-hl_n, 1], L_now = hl_y[-hl_n, 2]
)
hl_fit_H <- lm(H_next ~ H_now + L_now, data = hl_lagged)
hl_fit_L <- lm(L_next ~ H_now + L_now, data = hl_lagged)
hl_coef <- rbind(Hare = coef(hl_fit_H), Lynx = coef(hl_fit_L))
colnames(hl_coef) <- c("Intercept", "Current log hare", "Current log lynx")
knitr::kable(hl_coef, digits = 3,
             caption = "Rowwise least-squares estimates. Each row predicts the next year's log record for the named species.")
```

The hare row estimates $a_H,b_{HH},b_{HL}$; the lynx row estimates
$a_L,b_{LH},b_{LL}$. Thus $\widehat A_{21}=\widehat b_{LH}$ is the
coefficient of the current log hare record in the next log lynx equation.
$\diamond$

::: {.exercise}
**Exercise 11.3 (Assessing the hare–lynx model).** Use the fitted
regressions in Example 11.8, retaining their unrounded coefficients.

(a) Assemble the fitted intercept vector $\widehat a$ and coefficient
    matrix $\widehat A$. Compute the eigenvalues and spectral radius of
    $\widehat A$. Is the fitted recurrence stable? Find its equilibrium
    $\widehat\mu=\widehat a+\widehat A\widehat\mu$, and explain what the
    eigenvalues imply for long-horizon fitted forecasts.

(b) Starting from the observed log records in 1935, compute the fitted
    log-pelt forecasts for 1936 and 1940. Report both species in a table.
    Does exponentiating a log forecast give the conditional mean of the
    corresponding pelt record? Explain.

(c) Plot the residual autocorrelations through lag 10 for each equation,
    and report the two lag-one values. What do the plots suggest about
    the one-lag conditional-mean restriction? Would contemporaneous
    correlation between the two residual series itself violate the model?

(d) Interpret the signs of the two fitted cross-coefficients. If the
    fitted recurrence were taken as a stationary causal VAR(1) with
    iid Gaussian innovations and positive-definite covariance, which directions of
    Granger causality in mean would it imply for the log-pelt process?
    Explain why fitting these coefficients does not establish those
    directions for the actual process, or the effects of changing an
    animal population.
:::


The workbook on matrix filters and VARMA develops additional lags for
models in which one preceding vector leaves predictive dependence
unexplained.