Workbook 11 — Vector Autoregression

STA 542 · Introduction to Time Series Analysis

\[\def\E{\mathbb{E}} \def\P{\mathbb{P}} \def\R{\mathbb{R}} \def\Z{\mathbb{Z}} \def\Var{\operatorname{Var}} \def\Cov{\operatorname{Cov}}\]

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\)

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))

Figure 1: 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.

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.

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))

Moment estimates from the 300 simulated observations.

Entry	Truth	Estimate
A11	0.693	0.688
A12	-0.400	-0.475
A21	0.400	0.437
A22	0.693	0.665

The estimated spectral radius is 0.816, 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\)

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)
})

Figure 2: 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.

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 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. \]

1. 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?

2. 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.

3. 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.

4. 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 11.2 (Forecasting the Rotation Model). Use the Gaussian rotation model of Example 11.3 and its simulated observations.

1. 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)\).

2. 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.

3. 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 and astsa::Lynx.

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)
})

Figure 3: 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.

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.

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.")

Rowwise least-squares estimates. Each row predicts the next year’s log record for the named species.

	Intercept	Current log hare	Current log lynx
Hare	2.002	0.740	-0.372
Lynx	0.230	0.206	0.700

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 11.3 (Assessing the hare–lynx model). Use the fitted regressions in Example 11.8, retaining their unrounded coefficients.

1. 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.

2. 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.

3. 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?

4. 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.