---
title: "Workbook 13 — Fitting Time Series Models"
subtitle: "Regression, Innovations, and Error Correction"
format:
  html:
    embed-resources: true
    toc: true
    toc-depth: 2
    number-sections: false
    code-copy: true
    fig-width: 9
    fig-height: 3.6
    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\N{\mathbb{N}}
\def\cH{\mathcal{H}}
\def\Var{\operatorname{Var}}
\def\Cov{\operatorname{Cov}}
\def\Corr{\operatorname{Corr}}$$
:::

```{=html}
<style>
.exercise { border-left: 3px solid #00539B; padding: .7em 1.1em; background: #f6f8fa; margin: 1.4em 0; border-radius: 0 4px 4px 0; }
figure figcaption { color: #555; font-size: .9em; }
</style>
```

```{r}
#| label: setup
#| include: false
options(digits = 4)
blue <- "#00539B"
orange <- "#C65D21"
```

## Fitting a VAR with Ordinary Least Squares

To fit a VAR($p$) to $d$ series, we fit one multiple regression for each
series. Its current value is the response; the preceding $p$ values of
all $d$ series are the $pd$ predictors. Each regression uses the same
predictors, with its own intercept and coefficients. Collecting the
equations gives the vector model
$$
X_t=c+A_1X_{t-1}+\cdots+A_pX_{t-p}+Z_t,
$$
where $c$ collects the intercepts, $A_k$ collects the coefficients on
lag $k$, and $Z_t$ collects the current errors.

Given observed vectors $x_1,\ldots,x_T$, we pair the response $x_{i,t}$
from series $i$ with every coordinate of $x_{t-1},\ldots,x_{t-p}$.
Times $t=p+1,\ldots,T$ give $T-p$ regression observations; the first
$p$ vectors supply the initial predictors. The equation being fitted is
$$
X_{i,t}=c_i+\sum_{k=1}^p\sum_{j=1}^d
(A_k)_{ij}X_{j,t-k}+Z_{i,t}.
$$
The coefficient on series $j$ at lag $k$ is $(A_k)_{ij}$. Thus the
coefficients for response series $i$ form row $i$ of each matrix $A_k$.

The errors in different equations may be correlated at the same time.
We assume i.i.d. innovation vectors with mean zero and covariance
$$
\Cov(Z_t)=\Sigma_Z,
$$
where $\Sigma_Z$ is an unknown, finite positive-definite matrix that is
the same at every time. Its diagonal entries are innovation variances;
its off-diagonal entries are covariances between errors in different
equations at the same time. Those entries need not be zero. The
independence assumption concerns innovation vectors at different times.

Ordinary least squares minimizes the sum of squared residuals in each
response equation. Together, these regressions minimize
$$
\sum_{t=p+1}^T
\left\|x_t-c-\sum_{k=1}^p A_kx_{t-k}\right\|_2^2.
$$
The squared vector norm adds the squared residuals across the $d$
equations. Each equation has its own coefficients, so minimizing its
sum does not constrain the other equations. If the predictor columns,
including the intercept, are linearly independent, the $d$ OLS fits
uniquely determine $\widehat c$ and $\widehat A_1,\ldots,\widehat A_p$.

Each equation has its own unrestricted coefficients and uses the same
predictors. Consequently, using $\Sigma_Z^{-1}$ to weight the residual
vectors gives the same OLS coefficient estimates, even when the errors
are correlated across equations. $\Sigma_Z$ does enter the covariances
between coefficient estimates from different equations. We estimate it
from residual squares and cross-products after fitting.

To use these regressions for prediction, we assume that the current error
has mean zero given the past:
$\E[Z_t\mid\mathcal H_{t-1}]=0$, where $\mathcal H_{t-1}$ contains the
observations through time $t-1$. The lagged predictors are already known
given that history, so the regression's linear expression is the
conditional mean of $X_t$. After fitting, we insert the most recent $p$
observations to estimate the next conditional mean:
$$
\widehat x_{T+1\mid T}
=\widehat c+\widehat A_1x_T+\cdots+\widehat A_px_{T+1-p}.
$$
The coefficient estimates and the observed lags together determine this
one-step forecast.

## Output Growth and Unemployment

**Example 13.1 (Predicting unemployment changes).** When production grows,
firms often need more labor. The negative empirical association between
output growth and changes in unemployment is known as
[Okun's law](https://www.stlouisfed.org/publications/regional-economist/october-2013/output-and-unemployment-how-do-they-relate-today).
Does recent output growth also help predict the next change in unemployment,
after allowing for unemployment's own past?

The dataset
[`astsa::econ5`](https://search.r-project.org/CRAN/refmans/astsa/html/econ5.html)
contains quarterly U.S. real gross national product (GNP), in billions of
1982 dollars, and unemployment as a percentage of the labor force. We use
1948 Q3 through 1988 Q2. Write $Q_t$ for real GNP and $R_t$ for the
unemployment rate, and model
$$
X_t=\begin{pmatrix}G_t\\D_t\end{pmatrix},
\qquad G_t=100\{\log Q_t-\log Q_{t-1}\},
\qquad D_t=R_t-R_{t-1}.
$$
Here $G_t$ approximates quarterly percentage growth, without annualizing;
$D_t$ is a change in percentage points. A move in unemployment from 5% to
6% gives $D_t=1$.

```{r}
#| label: macro-data
data(econ5, package = "astsa")
macro <- window(econ5, end = c(1988, 2))
x <- cbind(growth = 100 * diff(log(macro[, "gnp"])),
           du = diff(macro[, "unemp"]))
```

```{r}
#| label: fig-macro-changes
#| echo: false
#| fig-cap: "Real data: quarterly GNP growth and unemployment changes, 1948 Q4–1988 Q2. Their units differ; each panel has its own vertical scale. Large declines in output often accompany increases in unemployment."
par(mfrow = c(1, 2), mar = c(4, 4.3, 2, 1), bg = "white")
plot(time(x), x[, "growth"], type = "l", col = blue, lwd = 1.5,
     xlab = "Year", ylab = "Quarterly growth (approx. %)", main = "Output growth")
abline(h = 0, col = "gray70", lty = 3)
plot(time(x), x[, "du"], type = "l", col = orange, lwd = 1.5,
     xlab = "Year", ylab = "Change (percentage points)",
     main = "Unemployment changes")
abline(h = 0, col = "gray70", lty = 3)
```

### Fit the Regression Equations

We predict this quarter's unemployment change from growth and
unemployment changes in each of the preceding two quarters. These are
the response and four predictors in the unemployment equation:
$$
\begin{aligned}
D_t={}&c_2+(A_1)_{21}G_{t-1}+(A_1)_{22}D_{t-1}\\
 &\quad +(A_2)_{21}G_{t-2}+(A_2)_{22}D_{t-2}+Z_{2,t}.
\end{aligned}
$$
Predicting growth from the same four predictors gives the other
regression. Together they form a VAR(2), with the growth coefficients
in the first rows of $A_1,A_2$ and the unemployment coefficients in the
second rows. The following code pairs responses with their predictors
and displays the first three observations for the unemployment regression.

```{r}
#| label: var-design
# embed orders columns as current vector, lag 1, then lag 2.
paired <- as.data.frame(embed(x, 3))
names(paired) <- c("growth", "du", "growth_l1", "du_l1",
                   "growth_l2", "du_l2")
knitr::kable(head(paired[, -1], 3), digits = 3, row.names = FALSE)
```

`lm` fits the unemployment equation from these columns. `vars::VAR`
constructs the lagged predictors and fits both equations by OLS:

```{r}
#| label: var-ols
ols_u <- lm(du ~ growth_l1 + du_l1 + growth_l2 + du_l2, data = paired)
fit_var <- vars::VAR(x, p = 2, type = "const")
coef(summary(fit_var$varresult$du))
```

The `lm` call and `fit_var$varresult$du` give the same coefficients.
The coefficient-table row `growth.l1` estimates $(A_1)_{21}$: the
coefficient on the first growth lag in the unemployment equation. Its estimate is
`r sprintf("%.3f", coef(fit_var$varresult$du)["growth.l1"])`.
Holding the other three predictors fixed, one percentage point more
output growth in the preceding quarter lowers the fitted unemployment
change by about
`r sprintf("%.3f", -coef(fit_var$varresult$du)["growth.l1"])`
percentage points.

The standard error measures uncertainty in this coefficient estimate.
It combines the residual variance with the spread and correlation of
the lagged predictors, using the usual OLS covariance formula. The
regression observations overlap in time, however: a response at one time
becomes a predictor at a later time. We therefore need a time-series
justification for using that formula.

For these standard errors, we use the stable causal VAR with the
i.i.d. innovations specified above. Causality and independence across
time make the current innovation independent of past predictors, so
$$
\Cov(Z_t\mid\mathcal H_{t-1})=\Sigma_Z.
$$
Thus its covariance remains constant even after conditioning on the
predictors. Stability ensures that past innovations have diminishing
effects; with the finite
moments, this allows the regression averages to converge. Together these
conditions give consistent OLS estimates and justify the usual
coefficient covariance formula in large samples. Gaussianity is not needed.

`confint` uses the estimate plus or minus a Student $t$
critical value times its standard error, with regression residual degrees
of freedom; the coefficient p-values use the same $t$ reference. These
are large-sample approximations for the time-series regression, even
under Gaussian innovations.

```{r}
#| label: var-interval
ci_growth <- confint(fit_var$varresult$du, "growth.l1")
ci_growth
```

The approximate 95% interval is
[`r sprintf("%.3f", ci_growth[1])`, `r sprintf("%.3f", ci_growth[2])`]
percentage points of unemployment change per percentage point of growth.
Both endpoints are negative. This interval concerns the first growth
lag; the predictive contribution of growth also includes the second lag.

### Does Growth Add Predictive Information?

The question about growth as a whole compares two conditional-mean
forecasts of unemployment change: one using the histories of both series,
and one using only unemployment's history. Growth *Granger-causes*
unemployment change in mean if adding growth's history changes that
forecast with positive probability. The direction matters: predicting
unemployment from growth is a different question from predicting growth
from unemployment.

For this stable VAR with i.i.d. innovations having a joint density,
no Granger causality from growth to unemployment change is equivalent to
$$
(A_1)_{21}=(A_2)_{21}=0.
$$
With both coefficients zero, the unemployment conditional mean uses only
its own lagged values. Conversely, the joint-density condition rules out
exact relations that would make a nonzero growth contribution already
determined by unemployment's history. Thus a nonzero growth coefficient
adds predictive information. The joint test assesses the two growth lags
together.

`causality` tests these restrictions together, accounting for the
estimated variances and covariance of the two growth coefficients. It
scales their Wald statistic and compares it with an approximate $F$
reference distribution under the joint null, using the conventional VAR
coefficient covariance estimate.

```{r}
#| label: var-granger
growth_test <- vars::causality(fit_var, cause = "growth")$Granger
growth_test
```

The p-value is `r format.pval(growth_test$p.value, digits = 2)`, giving
strong evidence against the two zero restrictions. In the fitted VAR(2),
recent output growth adds predictive information beyond unemployment's
own history. The negative growth coefficients associate stronger recent
growth with a smaller subsequent rise, or a larger fall, in unemployment.
This is a predictive conclusion. An effect of a policy intervention
would require assumptions about how that intervention changes the system.

The test uses the innovation assumptions of the fitted VAR. One
implication of i.i.d. innovations is zero serial correlation. The
following check combines residual autocorrelations and cross-correlations
through lag 12 and compares its statistic with an approximate chi-squared
reference distribution. Workbook 14 develops this joint test:

```{r}
#| label: var-residual-check
macro_check <- vars::serial.test(fit_var, lags.pt = 12,
                                 type = "PT.adjusted")$serial
macro_check
```

Its p-value is `r sprintf("%.2f", macro_check$p.value)`, so this check
finds no evidence of remaining autocorrelation through lag 12 at the
5% level.

### Recover the Mean and Innovation Covariance

The `const` coefficients are equation intercepts. The lagged terms also
contribute to the stationary mean: taking expectations in the VAR gives
$\mu=c+A_1\mu+A_2\mu$. Hence $c=(I_2-A_1-A_2)\mu$, and we solve this
equation using the fitted coefficients to estimate $\mu$. In each
coefficient matrix, rows correspond to responses and columns to lagged
predictors, in the order of `x`:

```{r}
#| label: var-mean
A_hat <- vars::Acoef(fit_var)
c_hat <- vapply(fit_var$varresult, function(eq) coef(eq)["const"], numeric(1))
mu_hat <- solve(diag(2) - A_hat[[1]] - A_hat[[2]], c_hat)
names(mu_hat) <- colnames(x)
mu_hat
vars::roots(fit_var)
```

The reported companion-matrix eigenvalue moduli are all below one, so
the fitted model is stable. The entries of `mu_hat` are its mean quarterly
growth and mean unemployment change, in the units defined above.

The remaining one-step errors estimate the innovation covariance:
$$
\widehat z_t=x_t-\widehat c-\widehat A_1x_{t-1}
                         -\widehat A_2x_{t-2}.
$$
`summary(fit_var)$covres` averages their squares and cross-products with a
regression degrees-of-freedom adjustment. Its diagonal entries estimate
the two innovation variances; the off-diagonal entry estimates their
same-quarter covariance.

```{r}
#| label: var-covariance
Sigma_hat <- summary(fit_var)$covres
Sigma_hat
```

This covariance describes variation left in the two response equations.
The earlier interval instead used the standard error of an estimated
coefficient, whose squared value is stored in
`vcov(fit_var$varresult$du)["growth.l1", "growth.l1"]`.
$\diamond$

## Fitting a VARMA in Two Stages {#when-the-regressors-include-unobserved-innovations}

A VARMA adds past innovations to the predictors used by a VAR:
$$
X_t=c+\sum_{i=1}^p A_iX_{t-i}+Z_t+\sum_{j=1}^q C_jZ_{t-j}.
$$
For a $d$-dimensional series, each response equation would have
$pd+qd$ predictors and an intercept. The observation lags are available,
but the innovation lags are not: the data contain $x_t$, not $z_t$.
Computing those errors requires the coefficients we are trying to
estimate. We first need a way to approximate the missing predictors.

We first fit a VAR with many observation lags and use its prediction
errors to approximate the innovations. Invertibility explains why this
can work. For a stationary causal, invertible VARMA with mean zero, the
inverse filter from Workbook 12 gives an infinite autoregression,
$$
X_t=\sum_{k=1}^{\infty}F_kX_{t-k}+Z_t.
$$
With i.i.d. mean-zero innovations of finite positive-definite covariance,
the remaining error $Z_t$ has conditional mean zero given past
observations. Thus the sum is the conditional-mean prediction. The
matrices $F_k$ account for both AR and MA effects and need not vanish
after lag $p$: approximating this prediction can require many observation
lags even when the original VARMA has only a few lags. This argument
does not require Gaussian innovations.

The first stage fits a finite VAR with $m$ observation lags. Its residuals
are
$$
\widetilde z_t=x_t-\sum_{k=1}^m\widehat F_kx_{t-k}.
$$
If these fitted predictions approximate the infinite-lag predictions,
then $\widetilde z_t$ approximates the innovation $z_t$. The approximation
has two sources of error: omitted lags beyond $m$, and estimation of the
retained coefficients. More lags reduce truncation but require estimating
more coefficients. For an unknown nonzero mean, include an intercept in
this preliminary regression.

The second stage treats these residuals as supplied predictors. Regress
$x_t$ on $x_{t-1},\ldots,x_{t-p}$ and
$\widetilde z_{t-1},\ldots,\widetilde z_{t-q}$, including an intercept
unless it is known to be zero. The observation-lag coefficients estimate
$A_1,\ldots,A_p$; the residual-lag coefficients estimate
$C_1,\ldots,C_q$. Both stages are ordinary least-squares calculations.

Even knowing the innovations would not always let us estimate separate
AR and MA coefficients. Consider a mean-zero VARMA(1,1) with $C=-A$:
$$
X_t=AX_{t-1}+Z_t-AZ_{t-1},
\qquad X_t-Z_t=A(X_{t-1}-Z_{t-1}).
$$
If all eigenvalues of $A$ lie inside the unit circle, the stationary
causal solution is $X_t=Z_t$. Thus $A=0$, $A=0.5I_d$, and many other
choices, each paired with $C=-A$, describe exactly the same white-noise
process. These representations are all stable and invertible. The
observations cannot distinguish their separate AR and MA coefficients;
this is a failure of *identification*.

In this case, the second-stage predictor blocks $X_{t-1}$ and $Z_{t-1}$
are identical. The regression can estimate their combined coefficient
$A+C$, whose true value is zero, but cannot estimate $A$ and $C$
separately. More observations do not resolve the duplication. Replacing
the innovations by estimated residuals may produce numerically distinct
predictor columns, but that estimation error supplies no information
about the separate matrices.

Near cancellation can also cause trouble in a finite sample: the two
predictor blocks can be nearly redundant, so substantially different
AR and MA estimates can give very similar fitted predictions. Accurate
prediction alone therefore does not show that individual coefficients
have been estimated precisely.

For the exact cancellation above, we use the white-noise model
$X_t=Z_t$ and omit the redundant AR and MA terms. More general VARMA
models can have redundant combinations of matrix coefficients. A suitable
set of coefficient restrictions selects one representation from those
that describe the same process. Merely choosing $p,q$ and checking
stability and invertibility does not guarantee unique coefficients;
see [Wilms et al. (2019, Section 2)](https://arxiv.org/html/1707.09208v3#S2).

**Example 13.2 (Estimating the hidden regressors).** Take $d=2$, $p=q=1$,
and $c=0$, writing $A=A_1$ and $C=C_1$. Reuse the matrices and innovation
covariance from Example 12.6, taking iid Gaussian innovations for this simulation:
$$
A=0.8\begin{pmatrix}\cos(\pi/6)&-\sin(\pi/6)\\
\sin(\pi/6)&\cos(\pi/6)\end{pmatrix},\qquad
C=\begin{pmatrix}0.3&0.4\\0&-0.2\end{pmatrix},\qquad
Z_t\overset{\mathrm{iid}}{\sim}N_2\!\left(0,
\begin{pmatrix}1&0.6\\0.6&1\end{pmatrix}\right).
$$

Here $A+C$ is lower triangular with positive diagonal entries, so it is
nonsingular. This is sufficient to identify $A$ and $C$ within the
stable, invertible VARMA(1,1) family. To see why, the first two
coefficients in its infinite autoregression satisfy
$$
F_1=A+C,\qquad F_2=-CF_1.
$$
The observation process determines these conditional-mean prediction
coefficients. When $F_1$ is nonsingular, they in turn determine
$$
C=-F_2F_1^{-1},\qquad A=F_1-C.
$$
Thus the separate matrices are meaningful estimation targets in this
example.

Both the stationary mean and the regression intercept are zero. We generate
500 observations after discarding 500 startup observations. Saving the
true innovations lets us assess how well the preliminary residuals
recover them.

```{r}
#| label: simulate-varma
#| code-fold: true
#| code-summary: "Generate the VARMA sample (seed 54213)"
set.seed(54213)
n <- 500L                         # retained sample size T
burn <- 500L
A <- 0.8 * matrix(c(cos(pi/6), -sin(pi/6), sin(pi/6), cos(pi/6)),
                  2, byrow = TRUE)
C <- matrix(c(.3, .4, 0, -.2), 2, byrow = TRUE)
Sigma <- matrix(c(1, .6, .6, 1), 2)
z <- matrix(rnorm(2 * (n + burn)), ncol = 2) %*% chol(Sigma)
y <- matrix(0, n + burn, 2)
for (tt in 2:nrow(y))
  y[tt, ] <- A %*% y[tt - 1, ] + z[tt, ] + C %*% z[tt - 1, ]
y <- y[burn + seq_len(n), ]
z <- z[burn + seq_len(n), ]
colnames(y) <- c("first", "second")
```

Take $m=8$ for the preliminary VAR. Each equation then uses 16
observation-lag coefficients and 492 regression rows. Exercise 13.2
compares the resulting estimate with preliminary orders 4 and 12.

```{r}
#| label: preliminary-var
preliminary <- vars::VAR(y, p = 8, type = "none")
z_tilde <- matrix(NA_real_, n, 2)
z_tilde[9:n, ] <- residuals(preliminary)
```

```{r}
#| label: fig-recovered-innovations
#| echo: false
#| fig-cap: "Simulated data: each point compares a true innovation with the corresponding residual from the preliminary VAR(8). The diagonal marks exact recovery. True innovations would be unavailable for this comparison in real data."
par(mfrow = c(1, 2), mar = c(4, 4, 2, 1), bg = "#F1F3F5")
for (j in 1:2) {
  plot(z[9:n, j], z_tilde[9:n, j], pch = 16, cex = .45,
       col = adjustcolor(blue, .5), asp = 1,
       xlab = "True innovation", ylab = "Preliminary residual",
       main = paste("Coordinate", j))
  abline(0, 1, col = orange, lwd = 1.5)
}
```

The points lie close to the diagonal in both coordinates. The correlations
between true innovations and preliminary residuals are
`r sprintf("%.3f", cor(z[9:n, 1], z_tilde[9:n, 1]))` and
`r sprintf("%.3f", cor(z[9:n, 2], z_tilde[9:n, 2]))`.
Here the finite VAR recovers much of the innovation variation, making
its lagged residuals useful approximations to the missing predictors.

The second regression uses times 10 through 500: time 9 is the first
available preliminary residual, so its lag first enters at time 10.
For each response coordinate, the four predictors at time $t$ are
$x_{1,t-1},x_{2,t-1},\widetilde z_{1,t-1},\widetilde z_{2,t-1}$.
The first two coefficients form a row of $A$, and the last two form a row
of $C$. The `lm` coefficient matrix puts predictors in rows and responses
in columns, so the code transposes each two-row block to recover our
matrix convention.

```{r}
#| label: two-stage-fit
ii <- 10:n
regressors <- cbind(y[ii - 1, ], z_tilde[ii - 1, ])
stage2 <- lm(y[ii, ] ~ regressors - 1)
A_two <- t(coef(stage2)[1:2, ])
C_two <- t(coef(stage2)[3:4, ])
A_two
C_two
```

Fitting only a VAR(1) to these observations omits $CZ_{t-1}$.
That omitted term is generally correlated with $X_{t-1}$, which also
contains $Z_{t-1}$. The resulting OLS coefficients need not estimate $A$.

```{r}
#| label: omitted-ma
A_omit <- vars::Acoef(vars::VAR(y, p = 1, type = "none"))[[1]]
knitr::kable(data.frame(
  Fit = c("True value", "VAR(1): MA term omitted", "Two regression stages"),
  A12 = c(A[1, 2], A_omit[1, 2], A_two[1, 2])), digits = 3)
```

The VAR(1) estimate of $A_{12}$ is
`r sprintf("%.3f", A_omit[1, 2])`, compared with the true value $-0.400$.
Including the estimated innovation lag gives
`r sprintf("%.3f", A_two[1, 2])`, much closer to the truth in this sample.
The comparison illustrates the consequence of leaving the correlated
MA term out of the regression.

The two stages give preliminary estimates of the VARMA coefficients
without specifying an innovation density. The second-stage `lm`
standard errors treat the estimated innovation lags as observed
predictors; they do not account for the first-stage estimation. We use
this calculation for coefficient estimates, rather than interpreting its
ordinary regression p-values as VARMA tests.
$\diamond$

## Reading a VARMA Fit

The two-stage fit holds preliminary innovation estimates fixed during
the second regression. A likelihood fit instead evaluates each proposed
set of VARMA coefficients using the predictions and errors that those
coefficients imply. To assign a density to the observations, we now add
a distributional assumption to the stationary causal, invertible model:
$$
X_t=c+\sum_{i=1}^p A_iX_{t-i}+Z_t+\sum_{j=1}^q C_jZ_{t-j}.
$$
Assume
$$
Z_t\overset{\mathrm{iid}}{\sim}N_d(0,\Sigma_Z),
\qquad \Sigma_Z\text{ positive definite}.
$$
The orders $p,q$ are fixed; $c$, the matrices $A_i,C_j$, and $\Sigma_Z$
are unknown. Independence is across time. Coordinates of the same
innovation vector may be correlated, as described by $\Sigma_Z$.

As in Gaussian linear regression, maximizing the likelihood favors small
squared prediction errors. For a VARMA, we must reconstruct those errors
from each proposed set of coefficients before we can evaluate the fit.

### Evaluating a Proposed Fit

To calculate a prediction error, we need preceding observations and
innovations. Some of that history is missing at the start of the sample.
One approach treats the first $s=\max(p,q)$ observed vectors as given
and chooses $q$ starting residuals, for example zeros. For each proposed
set of coefficients, we then calculate the subsequent errors in time order:
$$
e_t=x_t-c-\sum_{i=1}^p A_ix_{t-i}-\sum_{j=1}^q C_je_{t-j},
\qquad t=s+1,\ldots,T.
$$
Each error becomes an estimated innovation for later predictions.
Changing the coefficients therefore changes both the predictions and the
innovation history used to make them. Evaluating a new proposal requires
repeating the recursion.

Evaluating the Gaussian innovation density at each error and adding its
logarithm gives the initialized Gaussian log likelihood:
$$
\ell_{\mathrm{init}}
=-\frac{(T-s)d}{2}\log(2\pi)
 -\frac{T-s}{2}\log\det\Sigma_Z
 -\frac12\sum_{t=s+1}^T e_t^\top\Sigma_Z^{-1}e_t.
$$
For fixed $\Sigma_Z$, only the last term depends on the coefficients.
Maximizing this log likelihood therefore minimizes a sum of squared
residual lengths, measured relative to the innovation covariance:
$$
\sum_{t=s+1}^T e_t^\top\Sigma_Z^{-1}e_t
=\sum_{t=s+1}^T\left\|\Sigma_Z^{-1/2}e_t\right\|_2^2.
$$
Here $\|v\|_2^2$ adds the squared coordinates of $v$, and
$\Sigma_Z^{-1/2}$ is the symmetric inverse square root of $\Sigma_Z$.
It standardizes the errors for their variances and same-time correlations.
If $\Sigma_Z=\sigma^2 I_d$, this criterion is simply
$\sigma^{-2}\sum_t\|e_t\|_2^2$, a constant multiple of the usual residual
sum of squares.

For a VAR, conditioning on the first $p$ observations gives residuals
that are linear in the unknown coefficients. As we saw at the start of the
workbook, the shared predictors and unrestricted coefficients let us
minimize this criterion by fitting each equation with OLS, even when
$\Sigma_Z$ is not diagonal. For a VARMA, the lagged errors depend on
the same coefficients we are fitting. The residuals are then nonlinear
functions of the coefficients, so a single OLS fit does not minimize
the criterion. Numerical optimization repeatedly changes the
coefficients, reconstructs the errors, and compares the resulting fit.

We also need to estimate $\Sigma_Z$. Its log-determinant term in the
log likelihood can no longer be treated as constant. For each proposed
set of coefficients, the likelihood-maximizing covariance is
$$
\widehat\Sigma_Z=\frac{1}{T-s}\sum_{t=s+1}^T e_te_t^\top,
$$
provided this matrix is positive definite. After this substitution,
maximizing the likelihood amounts to minimizing
$\log\det\widehat\Sigma_Z$ over the coefficients. Thus the fit uses
the residual variances and cross-products together, rather than just
adding unweighted squared errors. A two-stage fit can supply starting
coefficients for this optimization.

Choosing starting residuals makes this an approximation to the Gaussian
likelihood based only on the observed data. An exact calculation accounts
for uncertainty about the missing initial history, including its effect
on prediction-error covariances. `MTS::VARMA` uses conditional Gaussian
fitting, with starting residuals rather than an exact treatment of that
initial uncertainty.

### Read the R Fit

[`MTS::VARMA`](https://search.r-project.org/CRAN/refmans/MTS/html/VARMA.html)
computes its own preliminary VAR and residual regression estimates,
then jointly optimizes the AR and MA coefficients using the conditional
Gaussian criterion. The call below does not use our eight-lag estimates
as its starting values.

The arguments `p` and `q` specify the numbers of observation and
innovation lags. For the simulated example, take $p=q=1$ and omit the
intercept because the model has known mean zero:

```{r}
#| label: package-varma
#| results: hide
fit_ma <- MTS::VARMA(y, p = 1, q = 1, include.mean = FALSE)
```

For this VARMA(1,1), `Phi` is $\widehat A$. The package uses a minus
sign for its MA term, so $\widehat C$ is `-Theta`:

```{r}
#| label: package-varma-matrices
A_ml <- fit_ma$Phi
C_ml <- -fit_ma$Theta
```

The standard-error matrix `secoef` stores predictors in rows and
responses in columns. Transposing its first two rows gives the standard
errors for $A$; transposing its next two gives those for $C$.
The following table compares both estimates with the simulation truth.

```{r}
#| label: varma-comparison
SE_A <- t(fit_ma$secoef[1:2, ])
SE_C <- t(fit_ma$secoef[3:4, ])
entry <- c("11", "12", "21", "22")
comparison <- data.frame(
  Entry = c(paste0("A", entry), paste0("C", entry)),
  Truth = c(t(A), t(C)),
  Two_stage = c(t(A_two), t(C_two)),
  Package = c(t(A_ml), t(C_ml)),
  SE = c(t(SE_A), t(SE_C)))
knitr::kable(comparison, digits = 3)
```

For $A_{12}$, the two-stage and package estimates both round to
`r sprintf("%.3f", A_ml[1, 2])`; likelihood fitting changes little here.
For $C_{12}$, the package estimate
`r sprintf("%.3f", C_ml[1, 2])` is closer to the true $0.400$ than the
two-stage estimate `r sprintf("%.3f", C_two[1, 2])`.
The $C_{21}$ estimate moves farther from its true value of zero.
Likelihood fitting optimizes the Gaussian criterion; in a finite sample,
this can improve some coefficient estimates while worsening others.

The coefficient $C_{12}$ describes the contribution of the lagged innovation
in coordinate 2, $Z_{2,t-1}$, through the MA term, holding the lagged observations and the
other innovation coordinate fixed. An innovation also enters $X_{t-1}$,
so its full contribution one period later includes the AR term. Substituting
the equation for $X_{t-1}$ into $X_t=AX_{t-1}+Z_t+CZ_{t-1}$ gives the
coefficient $A+C$ on $Z_{t-1}$. In this simulation,
$A_{12}+C_{12}=-0.4+0.4=0$, even though $C_{12}$ is positive.

### What Supports the Reported Uncertainty?

The uncertainty calculation assumes stable and invertible dynamics.
For this VARMA(1,1), the corresponding checks on the fitted matrices
are the AR and MA spectral radii:

```{r}
#| label: varma-roots
c(AR = max(Mod(eigen(A_ml)$values)),
  MA = max(Mod(eigen(C_ml)$values)))
```

Both radii are below one, so the fitted model is stable and invertible.

The package estimates coefficient uncertainty from the curvature of the
Gaussian log likelihood near the fitted parameters, allowing the
coefficients to vary jointly. A likelihood that falls slowly as a
coefficient changes indicates less precise estimation. Inverting the
negative log likelihood's curvature matrix gives an estimated covariance
matrix for the coefficient estimates; the square roots of its diagonal
entries give the standard errors.

Identification is essential for this calculation. In the cancellation
case $C=-A$, different coefficient pairs give the same observation
distribution, so the exact likelihood cannot distinguish them. Its
curvature cannot provide separate standard errors for $A$ and $C$.
Near cancellation can leave very little curvature in some directions,
producing large or numerically unstable standard errors. A numerical
fit returning coefficient estimates does not by itself establish that
those coefficients are identifiable.

For a correctly specified, identifiable model with stable and invertible
dynamics, i.i.d. Gaussian innovations, and sufficient likelihood curvature
in every parameter direction, large-sample theory approximates a
coefficient estimate by a normal distribution. This gives the interval
$\widehat A_{12}\pm1.96\,\operatorname{SE}(\widehat A_{12})$.
For testing $A_{12}=0$, divide the estimate by its standard error and
compare with a standard normal distribution:

```{r}
#| label: varma-interval
A_ml[1, 2] + c(-1, 1) * 1.96 * SE_A[1, 2]
2 * pnorm(-abs(A_ml[1, 2] / SE_A[1, 2]))
```

`MTS` uses this normal approximation for its printed coefficient
p-values, even though it labels the statistic a `t value`. This tests one
coefficient within the chosen model and lag orders. It is not a general
Granger-causality test: in a VARMA, forecasts from observed history combine
the AR and MA matrices through the inverse filter in Workbook 12.

We can minimize the same Gaussian criterion when the innovation
distribution is unknown or non-Gaussian. This is called *Gaussian
quasi-likelihood*: it uses the model's predictions and error covariance
to choose coefficients without asserting a Gaussian innovation density.

The coefficient standard errors and normal tests above remain
asymptotically valid if we retain the model and curvature conditions
and replace Gaussian innovations by i.i.d. mean-zero innovations with
the same positive-definite covariance and finite fourth moments. This
result still assumes independence across time and constant conditional
covariance. Time-varying conditional variances can require different
standard errors; serial dependence can also indicate missing dynamics.

The normal approximation describes the sampling variation of the
coefficient estimates. The innovations themselves may have non-Gaussian
tails, so these coefficient tests do not justify Gaussian probabilities
for unusually large errors. Exercise 13.4 examines this distinction.

$\diamond$

## Fitting a Cointegrated VARMA {#fitting-a-cointegrated-var}

Cointegration allows each series to be $I(1)$ while some linear
combinations of their levels are $I(0)$. To see how these combinations
enter a VARMA fit, start with a model for the levels. Let $V_t$ collect
$d$ individually $I(1)$ series and write
$$
V_t=c+\sum_{i=1}^p A_iV_{t-i}
    +Z_t+\sum_{j=1}^q C_jZ_{t-j},
\qquad p\geq1,\quad q\geq0.
$$
Suppose there are exactly $r$ independent cointegrating relations,
where $0<r<d$. We will rewrite the model to show how these relations
can help predict changes.

Subtract $V_{t-1}$ from both sides and collect the lagged differences:
$$
\nabla V_t=c+\Pi V_{t-1}
 +\sum_{k=1}^{p-1}B_k\nabla V_{t-k}
 +Z_t+\sum_{j=1}^q C_jZ_{t-j},
$$
where
$$
\Pi=\sum_{i=1}^p A_i-I_d,
\qquad B_k=-\sum_{i=k+1}^p A_i.
$$
This is an algebraic rewriting of the levels equation. The term
$\Pi V_{t-1}$ is the only one that still involves levels rather than
changes or innovations.

Cointegration supplies stationary combinations of those levels.
Collect $r$ independent cointegrating vectors in the columns of a
$d\times r$ matrix $\beta$. We fit a model in which levels enter the
change equations through $\beta^\top V_{t-1}$. Each equation uses its
own linear combination of these $r$ quantities. Collecting those
coefficients in a $d\times r$ adjustment matrix $\alpha$ gives
$$
\Pi V_{t-1}=\alpha\beta^\top V_{t-1},
\qquad \Pi=\alpha\beta^\top.
$$
We require $\alpha$ to have $r$ independent columns, so
$\operatorname{rank}(\Pi)=r$. This is the rank restriction retained
when fitting the cointegrated model.

We also choose to restrict the intercept to $c=\alpha\rho$, where
$\rho$ is an $r$-vector of constants. Then the intercept and level
term combine as
$$
c+\Pi V_{t-1}=\alpha(\beta^\top V_{t-1}+\rho).
$$
Define the departures from the long-run relations by
$U_t=\beta^\top V_t+\rho$. Each coordinate of $U_t$ is $I(0)$;
the constants allow the relations to fluctuate around nonzero levels.
Substitution gives the *error-correction form*
$$
\nabla V_t=\alpha U_{t-1}
 +\sum_{k=1}^{p-1}B_k\nabla V_{t-k}
 +Z_t+\sum_{j=1}^q C_jZ_{t-j}.
$$
Each column of $\alpha$ describes how the series respond to one
departure. The $d\times d$ matrices $B_k$ multiply recent changes,
and the $d\times d$ matrices $C_j$ multiply past innovations. The
restriction on $c$ puts the constants inside the relations, with no
separate intercept in the change equations.

A relation's scale is arbitrary: multiplying it by a nonzero constant
and dividing its adjustment coefficients by that constant leaves every
prediction unchanged. We therefore fix a normalization when fitting
and reporting their coefficients. Normalization fixes how the relations
are written without changing the model's predictions.

### Fitting the Relations and Dynamics

If $\beta$ and $\rho$ were known, we could calculate the departures
from observed levels. The lagged changes are observed too. For $q>0$,
however, the innovation lags are still missing. Regressing changes only
on departures and lagged changes leaves out the MA contribution, which
can be correlated with those predictors. Knowing the long-run relations
therefore does not by itself make this an ordinary OLS fit.

For Gaussian likelihood fitting, assume
$$
Z_t\overset{\mathrm{iid}}{\sim}N_d(0,\Sigma_Z),
\qquad \Sigma_Z\text{ positive definite},
$$
with each innovation independent of preceding observations. The unknown
$\Sigma_Z$ is constant over time and may have nonzero off-diagonal
entries. We use an invertible MA representation, so past innovations
can be recovered from the full observed history. The conditional-mean
prediction includes all three terms: the level departures, lagged
changes, and lagged innovations. The current $Z_t$ has conditional mean
zero.

In a finite sample, reconstruct the innovations for each proposed set
of coefficients, as in the preceding VARMA fit. Treat the first
$s=\max(p,q)$ observations as given and choose starting residuals.
Then calculate
$$
\begin{aligned}
e_t={}&\nabla v_t-\alpha(\beta^\top v_{t-1}+\rho)\\
 &-\sum_{k=1}^{p-1}B_k\nabla v_{t-k}
  -\sum_{j=1}^q C_je_{t-j},
\qquad t=s+1,\ldots,T.
\end{aligned}
$$
Changing the relations, adjustment coefficients, or short-run dynamics
changes this entire residual sequence. The starting residuals approximate
unobserved initial innovations; invertibility makes their effect on later
errors decay. Up to a parameter-independent
constant, the initialized Gaussian log likelihood is
$$
\ell_{\mathrm{init}}
=-\frac{T-s}{2}\log\det\Sigma_Z
 -\frac12\sum_{t=s+1}^T\left\|\Sigma_Z^{-1/2}e_t\right\|_2^2.
$$
Thus the squared-error connection remains, but the errors now depend
on the unknown long-run relations as well as the AR and MA dynamics.

With the orders $p,q$ and rank $r$ fixed, numerical optimization changes
the coefficients, reconstructs the errors, and maximizes this criterion.
The fit must preserve the rank and constant restriction. It also uses
an identifiable VARMA representation, as discussed above.

### The Cointegrated VAR Case

When $q=0$, the innovation lags disappear. This special case is a
cointegrated VAR, usually fitted in its *vector error-correction model*
(VECM) form. With $\beta,\rho$ supplied, each change equation is now
an ordinary regression: its response is one coordinate of $\nabla v_t$,
and its predictors are the $r$ departures in $u_{t-1}$ and all coordinates
of $\nabla v_{t-1},\ldots,\nabla v_{t-p+1}$. Times $p+1,\ldots,T$
supply $T-p$ observations. OLS estimates $\alpha$ and the $B_k$.

When the relations are unknown, an unrestricted regression on
$v_{t-1}$ and lagged changes generally estimates a full-rank level
matrix. It need not preserve $\Pi=\alpha\beta^\top$ with rank $r$,
or give an intercept of the required form $\alpha\rho$. The joint
fit must estimate the relations while retaining these restrictions.

For any proposed $\beta,\rho$, calculate the departures and fit all
the change equations by OLS. These regressions give the best
$\alpha,B_1,\ldots,B_{p-1}$ for that proposed relation. As before,
common predictors and unrestricted coefficients in these regressions
make OLS appropriate even with a nondiagonal $\Sigma_Z$. Let $e_t$ be
their vector of residuals. The Gaussian log likelihood developed above
uses their squared lengths $e_t^\top\Sigma_Z^{-1}e_t$; estimating the
covariance gives
$$
\widehat\Sigma_Z=\frac{1}{T-p}\sum_{t=p+1}^T e_te_t^\top.
$$
The joint fit chooses the relations that minimize
$\log\det\widehat\Sigma_Z$, with the change equations refitted for
each proposed set of relations. Thus the relations are chosen by how
well they help predict changes in the whole system.

For $q=0$, the Johansen procedure solves this problem by
*reduced-rank regression*.
In [`urca::ca.jo`](https://search.r-project.org/CRAN/refmans/urca/html/ca.jo.html),
the calculation first regresses both the current changes and the
lagged levels on the lagged changes and keeps their residuals. The
levels also include a constant column when the relations have constants.
This removes the part explained by recent changes. An eigenvalue
calculation then finds the $r$ combinations of the remaining level
variables most strongly associated with the remaining changes. These
combinations solve the Gaussian likelihood problem with rank $r$.

After we specify the rank, `urca::cajorls` takes the corresponding
relations from `ca.jo`, normalizes their coefficients, and fits the
change equations by OLS using their lagged departures. This last OLS
step completes the joint fit: it uses the relations already selected by
the system likelihood.

## Treasury Rates and Adjustment

We illustrate the $q=0$ case with Treasury rates, using `ca.jo` and
`cajorls`. These functions fit a cointegrated VAR; they do not estimate
the MA coefficients in the general model above.

A longer-term Treasury yield reflects expected future short rates and a
term premium for holding the longer bond
([Federal Reserve Bank of New York](https://www.newyorkfed.org/research/data_indicators/term-premia-tabs)).
This motivates looking for co-movement between short and long yields.
A relation between their levels may help predict their changes, even when
a model of recent changes alone misses that information.

**Example 13.3 (Predicting rate changes from a level departure).** The dataset
[`urca::Raotbl1`](https://rdrr.io/cran/urca/man/Raotbl1.html) contains
natural logarithms of U.S. three-month and ten-year Treasury rates, quoted
in percent per year. Keep the 143 jointly observed quarters,
1953 Q2–1988 Q4. Let $S_t$ and $L_t$ denote these log rates and write
$V_t=(S_t,L_t)^\top$. The supplied columns are already logged.

```{r}
#| label: rates-data
data(Raotbl1, package = "urca")
observed <- complete.cases(Raotbl1[, c("r3m", "r10y")])
rates <- ts(as.matrix(Raotbl1[observed, c("r3m", "r10y")]),
            start = c(1953, 2), frequency = 4)
colnames(rates) <- c("short", "long")
```

```{r}
#| label: fig-rate-levels
#| echo: false
#| fig-cap: "Real data: three-month and ten-year U.S. Treasury rates, 1953 Q2–1988 Q4. The left panel undoes the supplied logarithms to show percent-per-year units; the right shows the logarithms used for fitting. Shared persistent movement motivates a candidate long-run relation."
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1), bg = "white")
matplot(as.numeric(time(rates)), exp(rates), type = "l", lty = c(1, 2), lwd = 1.5,
        col = c(blue, orange), xlab = "Year", ylab = "Rate (% per year)",
        main = "Rates on their original scale")
legend("topleft", c("Three-month", "Ten-year"), col = c(blue, orange),
       lty = c(1, 2), lwd = 1.5, bty = "n", cex = .8)
matplot(as.numeric(time(rates)), rates, type = "l", lty = c(1, 2), lwd = 1.5,
        col = c(blue, orange), xlab = "Year", ylab = "Log rate",
        main = "The variables in the model")
```

To measure how far the short log rate is from a proposed relation
$S_t=a+bL_t$, subtract the level implied by the long log rate:
$$
U_t=S_t-a-bL_t.
$$
For this to be a cointegrating relation under the course definition,
$S_t$ and $L_t$ must each be $I(1)$ and $U_t$ must be $I(0)$.
The time-series plots alone cannot determine these properties. As an
initial check, apply an ADF test to each log rate, using an intercept and
an AIC-selected lag order up to four. The unit-root null is rejected when
the statistic is more negative than its critical value.

```{r}
#| label: rate-integration-check
#| code-fold: true
#| code-summary: "Check the individual log rates"
adf_checks <- sapply(1:2, function(j) {
  test <- urca::ur.df(rates[, j], type = "drift", lags = 4, selectlags = "AIC")
  c(Statistic = test@teststat[1], Critical_5pct = test@cval[1, "5pct"])
})
colnames(adf_checks) <- c("Short log rate", "Long log rate")
adf_checks  # Reject the unit-root null when the statistic is more negative.
```

The short-rate statistic is `r sprintf("%.2f", adf_checks[1, 1])`, below
the 5% critical value `r sprintf("%.2f", adf_checks[2, 1])`; the long-rate
test does not reject. Thus these tests do not support treating both rates
as $I(1)$. We nevertheless fit a model in which this single level
departure helps predict both rate changes. We impose that shared
structure and estimate its coefficients, including $a$ and $b$.
The calculation shows how to fit and read such a model; the evidence
above leaves its cointegration interpretation in doubt.

For each rate, predict its current change using the preceding level
departure and the two preceding changes in both rates. A positive
departure means that the short log rate is above $a+bL_t$. The two
change equations together are
$$
\nabla V_t=
\begin{pmatrix}\alpha_S\\\alpha_L\end{pmatrix}U_{t-1}
+B_1\nabla V_{t-1}+B_2\nabla V_{t-2}+Z_t.
$$
This is the rank-one, $q=0$ model above, with $\beta=(1,-b)^\top$,
$\rho=-a$, and $p=3$. The same departure enters both equations,
with adjustment coefficients $\alpha_S$ and $\alpha_L$. A negative
$\alpha_S$, for example, makes a positive departure lower the predicted
short-rate change.

### Fit the Relation, Then the Change Equation

There are two sets of coefficients to estimate: $a,b$ determine the
long-run relation, while $\alpha_S,\alpha_L,B_1,B_2$ describe the subsequent
changes. A direct two-step approach first regresses $S_t$ on an intercept
and $L_t$. OLS chooses $\widehat a,\widehat b$ to minimize
$\sum_t(S_t-a-bL_t)^2$. From the observations, this gives the estimated
departures
$$
\widehat u_t=s_t-\widehat a-\widehat b\,l_t.
$$
These are the residuals from the levels regression.

Now hold the estimated departures fixed. To fit the short-rate equation,
regress the current change $\nabla s_t$ on the five predictors
$$
\widehat u_{t-1},\quad
\nabla s_{t-1},\quad\nabla l_{t-1},\quad
\nabla s_{t-2},\quad\nabla l_{t-2}.
$$
The first coefficient estimates $\alpha_S$: it measures how the short
rate adjusts to the preceding departure after allowing for recent
changes. The other four coefficients estimate the first rows of $B_1$
and $B_2$. Using $\nabla l_t$ as the response with the same predictors
fits the long-rate equation.

The constant is already inside $\widehat u_{t-1}$, so the change
regression uses `-1` to omit a separate intercept. The code below fits
the levels relation and then the short-rate change equation.

```{r}
#| label: ecm-ols
long_run <- lm(short ~ long, data = as.data.frame(rates))
u_hat <- residuals(long_run)
dr <- diff(rates)
ii <- 3:nrow(dr)
ecm_data <- data.frame(
  ds = dr[ii, "short"], dl = dr[ii, "long"],
  u_l1 = u_hat[ii],
  ds_l1 = dr[ii - 1, "short"], dl_l1 = dr[ii - 1, "long"],
  ds_l2 = dr[ii - 2, "short"], dl_l2 = dr[ii - 2, "long"])
ecm_ols <- lm(ds ~ u_l1 + ds_l1 + dl_l1 + ds_l2 + dl_l2 - 1,
              data = ecm_data)
coef(long_run)
coef(ecm_ols)
```

Because `dr[ii, ]` is the change from level row `ii` to row `ii + 1`,
`u_hat[ii]` is the preceding level departure. Its fitted coefficient is
`r sprintf("%.3f", coef(ecm_ols)["u_l1"])`.
Its negative sign means that a positive departure lowers the fitted next
short-rate change, holding recent changes fixed.

### Read the System Fit

The levels regression chose $a,b$ to make $s_t-a-b l_t$ small.
We now use the joint Gaussian fit to choose $a,b$ by how well the
relation helps predict both rate changes. Its estimates need not equal
the levels-regression estimates: the two methods minimize different
criteria.

In `ca.jo`, `ecdet = "const"` puts the constant inside the relation,
and `spec = "transitory"` puts the level departure at time $t-1$, as in
our equation. Since $\nabla V_{t-2}=V_{t-2}-V_{t-3}$, rewriting the
model in levels uses lags through three. The argument `K` counts these
levels lags, so we set `K = 3` to include two lagged differences.

The same call computes trace tests for the rank of the lagged-level
matrix. Each compares the best Gaussian likelihood under a rank
restriction with a less restricted fit; large statistics argue against
the restriction. The two rows below test rank zero and rank at most one.
Their critical values come from the cointegration model's large-sample
theory and depend on the specification of constants and trends.

```{r}
#| label: johansen-fit
jo <- urca::ca.jo(rates, type = "trace", ecdet = "const",
                  K = 3, spec = "transitory")
knitr::kable(data.frame(
  Null = c("Rank zero", "Rank at most one"),
  Statistic = rev(jo@teststat), Critical_5pct = rev(jo@cval[, "5pct"])),
  digits = 2)
```

At 5%, the trace tests reject rank zero and do not reject rank at most
one. These results favor a rank-one level matrix for the specified lags
and constant. They do not resolve the earlier short-rate result: a
stationary short rate alongside an $I(1)$ long rate can itself produce
rank one, without a cointegrating relation between two $I(1)$ rates.
We retain the imposed rank for the calculation below.

[`cajorls`](https://search.r-project.org/CRAN/refmans/urca/html/cajorls.html)
with rank one uses the first relation estimated by `ca.jo` and normalizes
the short-rate coefficient to one. It then fits both change regressions
using this relation and the four lagged-change predictors. It returns
the normalized relation in `beta` and these OLS fits in `rlm`.

```{r}
#| label: vecm-output
fit_ec <- urca::cajorls(jo, r = 1)
fit_ec$beta
coef(fit_ec$rlm)
```

The normalization fixes the coefficient of $S_t$ at one. The fitted
relation is
$$
\widehat u_t=s_t
`r sprintf("%+.3f", fit_ec$beta[2, 1])`\,l_t
`r sprintf("%+.3f", fit_ec$beta[3, 1])`.
$$
In `coef(fit_ec$rlm)`, each column is a response equation. The row `ect1`
contains $\widehat\alpha_S$ and $\widehat\alpha_L$; `short.dl1`
multiplies $\nabla S_{t-1}$. The remaining lag rows are read in the same
way. The constant in the displayed relation is $-\widehat a$.

The table gives approximate 95% intervals for the adjustment coefficients
under the specified cointegrated model and its i.i.d. Gaussian
innovations. `summary(fit_ec$rlm)` computes ordinary regression standard
errors with the fitted level relation held fixed. In a correctly
specified cointegrated system, estimating that relation does not change
the leading large-sample distribution of the adjustment estimates. This
justifies using those regression standard errors here; it is a result
specific to cointegration, rather than a general rule for estimated
predictors.

Approximate normality gives the intervals
$\widehat\alpha\pm1.96\,\operatorname{SE}(\widehat\alpha)$.
Gaussianity is sufficient here; the same coefficient approximation also
holds with i.i.d. mean-zero innovations of constant positive-definite
covariance and finite fourth moments.
The calculation treats the rank, lag orders, and placement of the constant
as fixed. It describes coefficient uncertainty under the assumed model;
it does not resolve the conflicting evidence about whether both rates
are $I(1)$.

```{r}
#| label: vecm-adjustment
alpha_hat <- coef(fit_ec$rlm)["ect1", ]
ec_summaries <- summary(fit_ec$rlm)
alpha_se <- vapply(ec_summaries,
                   function(s) coef(s)["ect1", "Std. Error"], numeric(1))
knitr::kable(data.frame(Equation = c("Short log rate", "Long log rate"),
                        Adjustment = alpha_hat, SE = alpha_se,
                        Lower = alpha_hat - 1.96 * alpha_se,
                        Upper = alpha_hat + 1.96 * alpha_se),
              row.names = FALSE, digits = 3)
```

A departure $\widehat u_{t-1}=0.01$ means the short rate is about 1% above
the level implied by the fitted relation. Holding lagged changes fixed,
the adjustment term contributes approximately
`r sprintf("%.4f", .01 * alpha_hat[1])` to the next short log-rate change:
a fall of about `r sprintf("%.3f", -alpha_hat[1])` percent in the short
rate. In this fitted model, most adjustment occurs through the short rate.
The long-rate coefficient is small relative to its standard error, and
its interval includes zero.

```{r}
#| label: fig-rate-relation
#| echo: false
#| fig-cap: "Real data and fitted model: the estimated log-rate departure and the ACF of short-rate residuals after including that departure and two lags of changes. The residual ACF shows the serial dependence left by this fit."
u_ec <- drop(cbind(rates, 1) %*% fit_ec$beta)
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1), bg = "white")
plot(as.numeric(time(rates)), u_ec, type = "l", col = blue, lwd = 1.5,
     xlab = "Year", ylab = "Departure (log-rate units)",
     main = "Estimated long-run departure")
abline(h = 0, col = "gray60", lty = 3)
acf(residuals(fit_ec$rlm)[, 1], lag.max = 12,
    main = "Short-rate residuals", xlab = "Lag (quarters)")
```

`vec2var` converts the fitted VECM to its equivalent levels VAR for methods
such as residual checking. It does not fit a new model.

```{r}
#| label: vecm-check
fit_levels <- vars::vec2var(jo, r = 1)
rates_check <- vars::serial.test(fit_levels, lags.pt = 12,
                                 type = "PT.adjusted")$serial
rates_check
```

The joint residual check has p-value
`r sprintf("%.2f", rates_check$p.value)`, so it finds no evidence of
remaining autocorrelation through lag 12 at 5%. The fitted rank-one
model predicts that departures are followed mainly by short-rate
adjustment. Deciding whether to interpret this as adjustment between
cointegrated rates still requires resolving the conflicting integration
evidence; a residual autocorrelation test does not address that question.
$\diamond$

## Exercises

::: {.exercise}
**Exercise 13.1 (Read the unemployment equation).** Use Example 13.1.

(a) Identify the matrix entry estimated by `growth.l2`. Interpret its
estimate and a 95% interval in the units of the two series.

(b) Is `const` the estimated mean unemployment change? Use the fitted
matrix equation to obtain that mean.

(c) The p-value for `du.l1` is large. Does that settle whether unemployment
changes Granger-cause output growth? Name the restrictions for that question.
:::


::: {.exercise}
**Exercise 13.2 (Read the innovation coefficients).** Use Example 13.2.

(a) Explain why the preliminary residuals are needed. Which regressors
would be available if the true innovations were supplied?

(b) Read $C_{12}$ and its standard error from `fit_ma`. Interpret its
contribution through the MA term, respecting the package's sign convention.
Why does $C_{12}$ alone not give the full one-period contribution of the
innovation in coordinate 2 to coordinate 1?

(c) Repeat both regression stages, using 4 and then 12 lags in the
preliminary VAR. Compare the resulting second-stage $A_{12}$ estimates
with the eight-lag fit. Does closer agreement with the
truth in this one simulation identify a generally best preliminary order?
:::


::: {.exercise}
**Exercise 13.3 (Read the adjustment).** Use Example 13.3.

(a) Translate the three entries of `fit_ec$beta` into
$U_t=S_t-a-bL_t$. Which entry gives $-a$?

(b) Set $u_{t-1}=0.02$. Calculate the contribution of the adjustment term
to each next log-rate change, holding lagged changes fixed. Explain the
signs and distinguish a percent change in a yield from a percentage-point
change in a yield.

(c) Why does a VAR fitted only to `diff(rates)` leave out part of this
model? What does `K = 3` mean in the Johansen call?
:::


::: {.exercise}
**Exercise 13.4 (A good mean fit need not give the right tails).** Consider
the stationary AR(1)
$$
X_t=0.6X_{t-1}+Z_t,\qquad
Z_t\overset{\mathrm{iid}}{\sim}0.95N(0,1)+0.05N(0,25).
$$
Generate the sample below and fit `arima(x_mix, order = c(1, 0, 0),
method = "ML")`. Here `order = c(1, 0, 0)` specifies one AR lag, no
differencing, and no MA lags. For this undifferenced model,
[`arima`](https://search.r-project.org/R/refmans/stats/html/arima.html)
reports the fitted mean $\mu$ as `intercept`, so its equation is
$X_t-\mu=\phi(X_{t-1}-\mu)+Z_t$. The returned `sigma2` estimates the
innovation variance. The option `method = "ML"` uses Gaussian densities,
so for these mixture innovations it is Gaussian quasi-likelihood. For
this stationary AR(1), it consistently estimates the AR coefficient,
mean, and innovation variance despite the non-Gaussian mixture density.

```{r}
#| label: mixture-sample
set.seed(542134)
z_mix <- rnorm(800, sd = ifelse(runif(800) < .05, 5, 1))
x_mix <- as.numeric(filter(z_mix, filter = .6, method = "recursive"))[301:800]
```

(a) Report the fitted AR coefficient, mean, and innovation variance. What
are their values in the generating model?

(b) Does an AR estimate near $0.6$ establish that the innovations are
Gaussian? Compute $\P(|Z_t|>3)$ under the supplied mixture and under a
Gaussian distribution with the same variance.
:::


*Further references:* [vars documentation](https://cran.r-project.org/package=vars),
[urca documentation](https://cran.r-project.org/package=urca), and
[MTS documentation](https://cran.r-project.org/package=MTS); Tsay (2014),
*Multivariate Time Series Analysis with R and Financial Applications*,
Chapters 2–4, for regression, VARMA estimation and error correction.
The code uses `vars` 1.6-1, `MTS` 1.2.1, `urca` 1.3-4 and `astsa` 2.5.