---
title: "Workbook 14 — Assessing Fitted Vector Models"
subtitle: "Residuals, Model Specification, and Forecast Evaluation"
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.8
    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>
: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: .91em; }
@media print {
  #TOC, .quarto-title-meta, .code-copy-button { display: none !important; }
  body { font-size: 10pt; }
  h2 { break-after: avoid; }
  figure, table { break-inside: avoid; }
}
</style>
```

```{r}
#| label: setup
#| include: false
options(digits = 6)
knitr::opts_chunk$set(comment = "#>")
wb_blue <- "#1C5276"
wb_copper <- "#AA5F32"
wb_gray <- "#657582"
```

## What Does the Fitted Model Leave Unexplained?

A fitted equation separates each observation into its prediction and an
error. In the models fitted in the workbook on regression, innovations,
and error correction, those errors approximate innovation vectors with
mean zero, constant covariance, and no correlation across time. A pattern
left in the errors can therefore reveal a feature the fitted equation
has missed.

**Definition 14.1 (Fitted innovation residuals).** For observed vectors
$x_1,\ldots,x_T\in\R^d$ and a fitted $\operatorname{VARMA}(p,q)$,
define the residuals recursively by
$$
\widehat z_t
:=x_t-\widehat c-\sum_{i=1}^p\widehat A_i x_{t-i}
       -\sum_{j=1}^q\widehat C_j\widehat z_{t-j},
\qquad t=s+1,\ldots,T,
$$
where the first $s\ge\max(p,q)$ observations and the starting residuals
supply the initial history. The MA matrices use the plus-sign convention
in the model equation. A VAR has $q=0$, so its residuals require only
observed lags.

The vectors $\widehat z_t$ are computed from the data and the fitted
coefficients. They estimate the unobserved innovations $Z_t$, whose
assumed mean and covariance structure we compare with the residuals.

**Example 14.2 (Errors in the growth–unemployment fit).** We return to the
quarterly GNP and unemployment series in
[`astsa::econ5`](https://search.r-project.org/CRAN/refmans/astsa/html/econ5.html).
As in the fitting workbook, the modeled coordinates are
$G_t=100\nabla\log Q_t$, approximately quarterly percentage GNP growth,
and $D_t=\nabla U_t$, the unemployment-rate change in percentage points,
where $Q_t$ denotes GNP and $U_t$ the unemployment rate.
Use the first 105 transformed vectors, 1948 Q4–1974 Q4, for the initial
VAR(2) fit. The later observations will be available for forecast evaluation.

The accompanying scripts
[`model-assessment.R`](R/model-assessment.R) and
[`companion-models.R`](R/companion-models.R) contain the repeated regression,
diagnostic, and forecasting calculations. The code below runs from the
workbook folder with the `astsa`, `vars`, `urca`, and `MTS` packages available.
The first fit uses the same OLS construction as the fitting workbook.

```{r}
#| label: load-and-fit-macro
source("R/model-assessment.R")
source("R/companion-models.R")
data(econ5, package = "astsa")
macro_levels <- window(econ5, end = c(1988, 2))
macro_x <- cbind(growth = 100 * diff(log(macro_levels[, "gnp"])),
                 du = diff(macro_levels[, "unemp"]))
macro_data <- list(y = as.matrix(macro_x), dates = as.numeric(time(macro_x)))
macro_fit <- wb14_fit_var(macro_data, p = 2, origin = 105,
                          response_start = 3, seasonal = FALSE)
```

To compare the two error series on one scale, divide each coordinate by
its residual root mean square. This marginal standardization does not
remove correlation between the coordinates. It makes a residual of size
two mean twice that coordinate's typical fitted error size.

```{r}
#| label: fig-macro-residuals
#| echo: false
#| fig-height: 4.4
#| fig-cap: "Real data: standardized VAR(2) residuals for GNP growth and unemployment changes through 1974 Q4. Unemployment errors are larger in the 1950s than across much of the 1960s."
wb14_plot_residuals(macro_fit, data_label = "GNP growth and unemployment changes",
                    time_label = "Year")
```

```{r}
#| label: macro-residual-plot-code
#| eval: false
wb14_plot_residuals(macro_fit, data_label = "GNP growth and unemployment changes",
                    time_label = "Year")
```

The unemployment residuals are larger through the 1950s and smaller
across much of the 1960s. This variation raises a question about the
constant innovation covariance assumed in the fit. The time plot
locates the larger errors; to examine whether earlier errors help
predict later ones, we turn to their lag correlations.

### Stationarity and White Noise

Weak stationarity permits dependence across time. For a weakly
stationary process $(Z_t)$, write
$\Gamma_Z(h)=\Cov(Z_{t+h},Z_t)$. Vector white noise additionally requires
$$
\E[Z_t]=0,\qquad \Gamma_Z(0)=\Sigma_Z,\qquad
\Gamma_Z(h)=0\quad(h\ne0).
$$
Here $\Sigma_Z$ is finite and positive definite. Its off-diagonal entries
may be nonzero: the same quarter's growth and unemployment innovations
may be correlated. The zero restriction concerns different times.

A stationary residual sequence with nonzero lag covariance can still
contain information useful for linear prediction. Even zero
*autocorrelations* in every coordinate are insufficient: an error in one
equation may predict a later error in another. Exercise 14.1 constructs
exactly this possibility. We therefore inspect the full matrix of lag
correlations.

Relabel the $T_*=T-s$ fitted residual vectors as $e_1,\ldots,e_{T_*}$.
Using zero as the model innovation mean, define
$$
\widehat\Gamma_e(h)
:=\frac{1}{T_*}\sum_{t=1}^{T_*-h}e_{t+h}e_t^\top,
\qquad 0\le h<T_*,
$$
and
$$
\widehat R_{ij}(h)
:=\frac{[\widehat\Gamma_e(h)]_{ij}}
 {\sqrt{[\widehat\Gamma_e(0)]_{ii}
              [\widehat\Gamma_e(0)]_{jj}}}.
$$
The divisor is $T_*$ at every lag. At positive $h$, row $i$ is the later
coordinate and column $j$ is the earlier coordinate. Negative lags use
$\widehat\Gamma_e(-h)=\widehat\Gamma_e(h)^\top$.

```{r}
#| label: fig-macro-correlations
#| echo: false
#| fig-height: 5.2
#| fig-cap: "Real data: residual auto- and cross-correlations from the growth–unemployment VAR(2). Positive lags compare the row coordinate with the column coordinate from an earlier quarter. Diagonal lag-zero values, which equal one, are omitted; off-diagonal lag-zero values are allowed to be nonzero. Dashed lines are the rough white-noise guides ±1.96/√103."
wb14_plot_correlations(macro_fit, H = 12, lag_unit = "quarters")
```

```{r}
#| label: macro-correlation-plot-code
#| eval: false
wb14_plot_correlations(macro_fit, H = 12, lag_unit = "quarters")
```

The strongest off-diagonal feature is the negative correlation at lag
zero: a quarter with unexpectedly high growth tends to have an
unexpectedly small unemployment increase. This contemporaneous
correlation is allowed by $\Sigma_Z$. The serial-correlation check
concerns the entries at nonzero lags.

The familiar guides $\pm1.96/\sqrt{T_*}$ give a rough scale for an
individual sample correlation under independent white noise. Estimated
parameters alter the residual distribution, and this figure displays
many correlations at once. Counting crossings is therefore not a
calibrated joint test. A joint statistic combines the evidence across
coordinates and lags.

### A Joint Check of Residual Correlation

Choose a maximum lag $H$. We ask whether any coordinate of the model
error $Z_t$ is correlated with any coordinate of $Z_{t+h}$ at one of
the lags $h=1,\ldots,H$. The hypotheses are
$$
\begin{aligned}
H_0 &: \Gamma_Z(h)=0 \quad\text{for every }h=1,\ldots,H,\\
H_1 &: \Gamma_Z(h)\ne0 \quad\text{for at least one }h=1,\ldots,H.
\end{aligned}
$$
The null rules out both auto- and cross-correlations at these nonzero
lags. We use fitted residuals to assess these restrictions on the
unobserved errors.

**Definition 14.3 (Adjusted multivariate portmanteau statistic).** For a
chosen maximum lag $H<T_*$ and nonsingular $\widehat\Gamma_e(0)$, define
$$
Q_H:=T_*^2\sum_{h=1}^H\frac{1}{T_*-h}
\operatorname{tr}\!\left\{
\widehat\Gamma_e(h)^\top\widehat\Gamma_e(0)^{-1}
\widehat\Gamma_e(h)\widehat\Gamma_e(0)^{-1}
\right\}.
$$
The trace term measures the size of the lag-$h$ covariance after
standardizing for the variances and contemporaneous correlations. Large
values of $Q_H$ give evidence against $H_0$.

The reference distribution uses stronger assumptions than $H_0$ alone:
the data follow a stable
$\operatorname{VAR}(p)$ with i.i.d., mean-zero innovations of constant
positive-definite covariance and finite fourth moments. For its OLS
residuals, we compare $Q_H$ with an approximate
$\chi^2_{d^2(H-p)}$ distribution. We take $H>p$, sufficiently
large to cover relevant lags but small relative to $T_*$, with the
order and deterministic regressors specified before this check.
The degrees of freedom account for the $d^2p$ fitted lag coefficients.

The p-value is the probability, under this reference distribution, of
a statistic at least as large as the observed $Q_H$. Reject $H_0$ when
this upper-tail probability is below the chosen significance level.
For the bivariate VAR(2) with $H=12$, the reference has
$2^2(12-2)=40$ degrees of freedom.

```{r}
#| label: macro-portmanteau
macro_q <- wb14_portmanteau(macro_fit, H = 12)
knitr::kable(row.names = FALSE, macro_q, digits = 3,
             col.names = c("Q12", "Reference df", "p-value"))
```

For this fit, the p-value is `r sprintf("%.3f", macro_q$p_value)`, so we
do not reject $H_0$ at the 5% level. The test does not detect remaining
auto- or cross-correlation through lag 12 quarters.
Dependence in squared residuals can remain when these
correlations are zero. Exercise 14.3 constructs such a process and
shows why white noise alone does not justify the i.i.d. reference
used here. $\diamond$

The same residual matrices can be plotted for VARMA and VECM fits.
Their test calibration must account for the fitted innovation recursion
or cointegration restrictions. The VAR reference above applies to the
VAR application below.

## Climate and Fish Recruitment

**Example 14.4 (Does climate history help predict recruitment?).**
The Southern Oscillation Index (SOI) summarizes a sea-level air-pressure
contrast between Tahiti and Darwin associated with tropical Pacific
climate. Fish recruitment records the entry of new fish into a population.
The series [`astsa::soi`](https://search.r-project.org/CRAN/refmans/astsa/html/soi.html)
and [`astsa::rec`](https://search.r-project.org/CRAN/refmans/astsa/html/rec.html)
contain 453 monthly index observations, January 1950–September 1987,
provided by the NOAA Pacific Fisheries Environmental Laboratory.
The recruitment observations are an index, not raw counts. Can past SOI
improve recruitment forecasts beyond what recruitment's own history supplies?

```{r}
#| label: fish-data
fish <- wb14_fish_data()
```

```{r}
#| label: fig-fish-training
#| echo: false
#| fig-height: 4.5
#| fig-cap: "Real data: the first 302 monthly SOI and recruitment observations, through February 1975. Annual variation and slower changes are both visible. The later 151 months are reserved for evaluating forecasts."
par(mfrow = c(2, 1), mar = c(2.5, 4, 2, 1), oma = c(2, 0, 0, 0),
    mgp = c(2.5, .6, 0), bg = "white")
for (j in 1:2) {
  plot(fish$dates[1:302], fish$y[1:302, j], type = "l",
       col = c(wb_blue, wb_copper)[j], lwd = 1.2,
       xlab = "", ylab = "Index", main = c("Southern Oscillation", "Fish recruitment")[j])
}
mtext("Year", side = 1, outer = TRUE, line = .5)
```

We fit the first 302 months, through February 1975. Each candidate uses
12 monthly indicators and the preceding $p$ observations of both
coordinates:
$$
X_t=b_{m(t)}+\sum_{i=1}^p A_iX_{t-i}+Z_t,
\qquad X_t=(S_t,R_t)^\top,
$$
where $S_t$ is SOI, $R_t$ is recruitment, $m(t)$ identifies the calendar
month, and each $b_m\in\R^2$ is a month-specific intercept. There is no
additional constant. Monthly intercepts and lag coefficients are
estimated jointly. Residuals subtract both the month's fitted intercept
and the fitted lag contributions; the stationarity assumption concerns
the innovations remaining after those terms.

Use the small candidate set $p\in\{1,2,6,12\}$, spanning a short memory,
half a year, and a year. Fit all four models to response months
$13,\ldots,302$. The first twelve months supply lagged predictors, leaving
$T_*=290$ common response vectors. This lets the later likelihood
comparison use exactly the same observations.

```{r}
#| label: fish-initial-fits
orders <- c(1, 2, 6, 12)
fish_fits <- lapply(orders, function(p)
  wb14_fit_var(fish, p = p, origin = 302, response_start = 13))
names(fish_fits) <- paste0("VAR(", orders, ")")
```

```{r}
#| label: fig-fish-residual-correlation-short
#| echo: false
#| fig-height: 5.2
#| fig-cap: "Real data: the seasonal VAR(2) leaves a large negative residual correlation between recruitment and SOI five months earlier. The lower-left panel places this at positive lag five; the upper-right panel shows the same pairs at negative lag five."
wb14_plot_correlations(fish_fits[["VAR(2)"]], H = 24)
```

```{r}
#| label: fig-fish-residual-correlation-six
#| echo: false
#| fig-height: 5.2
#| fig-cap: "Real data: extending the seasonal model to six lags greatly reduces the five-month cross-correlation. The regression can now use the SOI observation at that lag to predict recruitment."
wb14_plot_correlations(fish_fits[["VAR(6)"]], H = 24)
```

```{r}
#| label: fish-diagnostic-code
#| eval: false
wb14_plot_correlations(fish_fits[["VAR(2)"]], H = 24)
wb14_plot_correlations(fish_fits[["VAR(6)"]], H = 24)
```

```{r}
#| label: fig-fish-residual-time
#| echo: false
#| fig-height: 4.4
#| fig-cap: "Real data: standardized residuals of the seasonal VAR(6). A large negative recruitment residual remains near the end of the fitting period even after the lag-correlation pattern has been reduced."
wb14_plot_residuals(fish_fits[["VAR(6)"]], time_label = "Year")
```

```{r}
#| label: fish-residual-time-code
#| eval: false
wb14_plot_residuals(fish_fits[["VAR(6)"]], time_label = "Year")
```

The recruitment plot contains an error about nine residual root mean
squares below zero. An error this large is extremely unlikely under a
Gaussian innovation model with the fitted variance. The additional lags
have reduced serial correlation while leaving this distributional
mismatch visible.

```{r}
#| label: fish-portmanteau-table
fish_table <- wb14_fit_table(fish_fits, H = 24)
fish_diagnostics <- fish_table[, c("Model", "Q", "df", "p_value")]
fish_diagnostics$p_value <- format.pval(fish_diagnostics$p_value,
                                       digits = 3, eps = 1e-4)
knitr::kable(row.names = FALSE, fish_diagnostics, digits = c(0, 1, 0, 3),
             col.names = c("Model", "Q24", "Reference df", "p-value"))
```

The one- and two-lag models leave substantial serial dependence. The
lag-24 test does not reject the six- or twelve-lag model at the 5% level.
A larger p-value does not make twelve lags preferable to six: the test
is designed to detect departures, not to rank models that it fails to
reject. Additional lags also require estimating more coefficients.
$\diamond$

## Does Additional Complexity Earn Its Parameters?

Adding a lag to a VAR adds predictors to every equation. On the same
response dates, the larger model contains the smaller one by setting
the new coefficients to zero. Its maximized Gaussian likelihood therefore
cannot be smaller. Information criteria balance this improved fit against
the additional estimated parameters.

**Definition 14.5 (AIC and BIC).** Suppose a fitted model has maximized
log likelihood $\widehat\ell$, $k$ free identifiable parameters, and
$T_*$ scored response vectors. Define
$$
\operatorname{AIC}:=-2\widehat\ell+2k,
\qquad
\operatorname{BIC}:=-2\widehat\ell+k\log T_*.
$$
Smaller values are preferred within the candidate set. Here $T_*$ counts
time points, not $dT_*$ scalar entries.

Fitting and scoring on the same observations makes the fitted likelihood
look too favorable for predicting new data. Under a correctly specified,
regular likelihood model, AIC's $2k$ penalty gives a large-sample estimate
of the expected advantage on the $-2\ell$ scale. For the fish data, the large recruitment residual
calls the Gaussian assumption into question, so the Gaussian AIC remains
a penalized fit score without that predictive-risk justification.

BIC charges $\log T_*$ per parameter instead of two. Its penalty is
larger when $T_*>e^2$, so BIC requires a larger likelihood improvement
to accept the same additional coefficients. Exercise 14.2 calculates
the improvement required by each criterion.

### Comparing the Same Observations

For the initialized Gaussian criterion used in the fitting workbook,
let $e_t$ be the residual vectors on the common response dates. At the
fitted coefficients, the likelihood estimate of the innovation covariance is
$$
\widehat\Sigma_Z=\frac{1}{T_*}\sum_{t=1}^{T_*}e_te_t^\top.
$$
If this matrix is positive definite, profiling the covariance gives
$$
-2\widehat\ell
=T_*\left\{d\bigl(1+\log(2\pi)\bigr)
                  +\log\det\widehat\Sigma_Z\right\}.
$$
The likelihood divisor is $T_*$, without a regression degrees-of-freedom
correction. The parameter penalty is applied separately.

The candidates must use the same response coordinates, dates, units,
and comparable density and initialization conventions. In particular,
a joint SOI–recruitment likelihood cannot be compared numerically with
a recruitment-only likelihood. We can later compare their recruitment
forecasts because those predict the same target.

### Counting Free Parameters

For a $\operatorname{VAR}(p)$ with an intercept, the count is
$$
k=\underbrace{d^2p}_{\text{lag coefficients}}
  +\underbrace{d}_{\text{intercepts}}
  +\underbrace{d(d+1)/2}_{\text{innovation covariance}}.
$$
With twelve monthly intercept vectors, replace the intercept contribution
$d$ by $12d$. Thus the
bivariate seasonal models have $k=4p+27$. The covariance contributes
three parameters: two variances and one covariance.

```{r}
#| label: fish-ic-table
knitr::kable(row.names = FALSE, fish_table[, c("Model", "N", "k", "AIC", "BIC")],
             digits = c(0, 0, 0, 1, 1),
             col.names = c("Model", "Response months", "Parameters", "AIC", "BIC"))
```

Both criteria select six lags on this initial period. Moving from two
to six lags improves the fit enough to pay for sixteen more lag
coefficients. Moving from six to twelve improves the fitted likelihood,
but neither criterion finds the improvement sufficient for twenty-four
more coefficients. The six-lag model also removed the large five-month
residual cross-correlation.

For a VARMA, we count the free parameters in an *identified*
representation. If every AR and MA entry is free and identifiable, the
dynamic contribution is $d^2(p+q)$. The identified VARMA(1,1)
family from the fitting workbook supplies a concrete comparison:
with $d=2$ and an intercept, both VARMA(1,1) and VAR(2) have $k=13$.
One uses an innovation lag where the other uses a second observation lag.

**Example 14.6 (Observation lags or innovation lags?).** Return to the
105 growth–unemployment vectors from Example 14.2. Compare VAR(1),
VAR(2), and VARMA(1,1) on the same 103 response dates, starting at the
third vector. The VARMA is fitted with `MTS::VARMA`; its first retained
residual is initialized as $e_2=x_2-c$, following the package convention.
The table recomputes the full initialized Gaussian criteria on these
common dates instead of mixing the packages' printed scores.

```{r}
#| label: macro-information-criteria
macro_assessment <- wb14_macro_assessment()
stopifnot(nrow(macro_assessment$failures) == 0)
knitr::kable(row.names = FALSE, macro_assessment$comparison[, c("model", "n", "k", "AIC", "BIC")],
             digits = c(0, 0, 0, 1, 1),
             col.names = c("Model", "Response quarters", "Parameters", "AIC", "BIC"))
```

The VARMA has lower scores than VAR(2), with the same parameter count.
The difference comes from fitted likelihood. This fit converged, with
stable AR dynamics, an invertible MA part, and nonsingular
$\widehat A+\widehat C$, the identification condition used for this
VARMA(1,1). Exercise 14.5 compares their later forecast errors. $\diamond$

## Does the Model Predict Later Observations?

The preceding scores use observations that also determined the fitted
coefficients. To assess forecasts directly, fit through an origin $t$
and predict an observation that arrives later. A chronological split
keeps observations $1,\ldots,t$ for fitting and uses subsequent
observations for evaluation. Calendar indicators for a future month
are known at the origin; that month's SOI and recruitment are not.

### Constructing the Forecast

For a causal, invertible VARMA with independent, mean-zero innovations,
the full observed past determines the past innovations. With parameters
known, replacing future innovations by zero gives the conditional mean
given that full history. In a finite record, we instead use initialized,
reconstructed residuals and fitted coefficients. This gives the
*plug-in forecast*, an approximation that also depends on the treatment
of the missing initial history. Under merely weak white noise, the
corresponding full-history construction has a linear-prediction
interpretation; zero covariance alone does not give conditional mean zero.

At an observed origin $t$, the one-step forecast is
$$
\widehat x_{t+1\mid t}
=\widehat c+\sum_{i=1}^p\widehat A_i x_{t+1-i}
            +\sum_{j=1}^q\widehat C_j\widehat z_{t+1-j}.
$$
For $h>1$, some observation lags lie after $t$. Replace each of those
by its forecast from the same origin, and replace each innovation
with index greater than $t$ by zero. The fitted coefficients are held
fixed while computing this path.

For a VARMA(1,1), the first two steps are
$$
\begin{aligned}
\widehat x_{t+1\mid t}
 &=\widehat c+\widehat A x_t+\widehat C\widehat z_t,\\
\widehat x_{t+2\mid t}
 &=\widehat c+\widehat A\widehat x_{t+1\mid t}.
\end{aligned}
$$
The last fitted innovation affects the second forecast through the
first forecast. For the seasonal VAR, use the future month's fitted intercept vector at
each step and recursively forecast both coordinates.

```{r}
#| label: fish-first-forecast
fish_path <- wb14_forecast(fish_fits[["VAR(6)"]], fish, h = 12)
knitr::kable(row.names = FALSE, data.frame(
  Horizon = c(1, 6, 12),
  SOI = fish_path[c(1, 6, 12), 1],
  Recruitment = fish_path[c(1, 6, 12), 2]), digits = 2)
```

These predictions all use information available in February 1975.
When March is observed, a new forecast may use March's observed values.
That is a forecast from a new origin. Substituting the observed March
SOI into the February six-month prediction would instead give that
prediction information it did not have when it was issued.

### Repeating the Evaluation

An *expanding-origin evaluation* repeats this operation: fit observations
$1,\ldots,t$, forecast $t+h$, advance $t$ by one, and refit. The training
period grows as more observations become available. We use this rule
for the comparisons below.

For a set of origins $\mathcal O_h$, let
$$
u_{i,t,h}:=x_{i,t+h}-\widehat x_{i,t+h\mid t},\qquad
\operatorname{RMSE}_i(h)
:=\left\{\frac{1}{|\mathcal O_h|}
             \sum_{t\in\mathcal O_h}u_{i,t,h}^2\right\}^{1/2}.
$$
RMSE has the units of the target coordinate. Mean absolute error uses
$|u_{i,t,h}|$ in place of $u_{i,t,h}^2$ and omits the square root.
At each horizon, compare candidates on identical origins and target
observations. If coordinates have different units, report their errors
separately rather than adding them without a stated weighting rule.

There are also two different ways to handle model order. A fixed-order
procedure refits the coefficients of, say, VAR(6) at every origin. A
selection procedure refits every declared candidate and chooses the
order by that origin's AIC or BIC. In both cases, coefficients, seasonal
terms, and any order choice use only the training observations available
at that origin.

**Example 14.7 (Recruitment forecasts at several horizons).** Use
origins $t=302,\ldots,441$, giving 140 forecasts at each of the horizons
$h=1,6,12$. The common upper endpoint leaves all twelve future months
available. These origins run from February 1975 through September 1986;
the target dates depend on $h$. Compare the joint seasonal VARs with
recruitment-only seasonal AR models of the same candidate orders and
with a benchmark that predicts the historical mean for the target month.
The labels `VAR-BIC` and `VAR-AIC` denote joint VAR procedures that
reselect the order by the named criterion at each origin; `AR-BIC`
does the same within the recruitment-only family. All recruitment
forecasts are scored on the same target observations.

```{r}
#| label: rolling-fish-evaluation
fish_eval <- wb14_evaluate_fish(fish, initial = 302, orders = orders,
                               horizons = c(1, 6, 12))
stopifnot(all(fish_eval$scores$failures == 0),
          all(is.finite(fish_eval$scores$RMSE)))
```

```{r}
#| label: fig-fish-forecast-scores
#| echo: false
#| fig-height: 4.2
#| fig-cap: "Real data: recruitment RMSE from 140 expanding origins at each horizon. Including SOI improves the medium- and longer-horizon scores in this comparison. The preferred information-criterion procedure changes with horizon."
wb14_plot_rmse(fish_eval, target = "rec")
```

```{r}
#| label: fish-forecast-score-code
#| eval: false
wb14_plot_rmse(fish_eval, target = "rec")
```

```{r}
#| label: fish-forecast-score-table
fish_methods <- c("VAR-BIC", "AR-BIC", "VAR-AIC", "Seasonal mean")
fish_scores <- subset(fish_eval$scores, target == "rec" & model %in% fish_methods)
fish_score_table <- data.frame(Procedure = fish_methods)
for (h in c(1, 6, 12)) {
  a <- fish_scores[fish_scores$h == h, ]
  fish_score_table[[paste(h, "month RMSE")]] <- a$RMSE[match(fish_methods, a$model)]
}
knitr::kable(row.names = FALSE, fish_score_table, digits = 2)
```

BIC selects six VAR lags at all 140 origins. AIC selects twelve at
138 origins and six at the other two. Both had selected six on the
initial training period; the AIC choice changes as observations are
added. AIC's larger fitted models have lower one-month RMSE, while
BIC has lower twelve-month RMSE. Recruitment-only AR-BIC selects two
lags throughout. The joint model's six-month RMSE is about 26% lower
than that of this scalar procedure.

```{r}
#| label: fig-fish-six-month
#| echo: false
#| fig-height: 4.2
#| fig-cap: "Real data: six-month recruitment forecasts, plotted at their target dates. Each forecast was made six months earlier using only the observations then available. The joint model follows the low recruitment values in 1976 more closely than the scalar model."
wb14_plot_forecasts(fish_eval, target = "rec", h = 6,
                    models = c("VAR-BIC", "AR-BIC"))
```

```{r}
#| label: fish-six-month-plot-code
#| eval: false
wb14_plot_forecasts(fish_eval, target = "rec", h = 6,
                    models = c("VAR-BIC", "AR-BIC"))
```

The linear VAR-BIC procedure produces two negative six-month recruitment
forecasts. The index is nonnegative, so these predictions expose a
limitation that an average score can conceal. They are retained in the
reported errors. Predictions above the historical maximum of 100 are
extrapolations; the data documentation does not establish 100 as an
upper bound. Shumway and Stoffer (2025, Section 2.2) also describe
nonlinear SOI–recruitment relationships.

Errors from overlapping multistep forecasts can involve the same future
innovations. For forecasts made in consecutive months at a six-month
horizon, the errors involve five shared months of future innovations,
so they can be dependent.

Once this evaluation period has guided a new transformation, candidate
set, or forecasting rule, it has become development data. Assessment of
that revised rule needs another chronological evaluation period. $\diamond$

## Does Upstream Flow Improve a Downstream Forecast?

**Example 14.8 (Forecasting the Mississippi at Keokuk).** In the fish
comparison, SOI history improved forecasts beyond recruitment's own
history. A river supplies a second version of that question: does a
measurement upstream improve a forecast downstream? The USGS gauges at
[Clinton](https://waterdata.usgs.gov/monitoring-location/USGS-05420500/)
and [Keokuk](https://waterdata.usgs.gov/monitoring-location/USGS-05474500/),
Iowa, record daily mean Mississippi River discharge. Clinton is about
148 river miles upstream of Keokuk. We model flow in thousands of cubic
feet per second and forecast Keokuk.

The bundled [daily records](data/river-flow.csv) cover 2010–2024;
[source details and quality flags](data/river-flow-source.txt) accompany
them. Use 2010–2019 for the initial fit. The first plot shows those ten
years and compares daily changes at Clinton with changes at Keokuk up
to fourteen days later.

```{r}
#| label: river-data
source("R/river-assessment.R")
source("R/river-plots.R")
river <- wb14_river_data()
river_initial <- max(which(river$dates <= as.Date("2019-12-31")))
```

```{r}
#| label: fig-river-training
#| echo: false
#| fig-height: 5.8
#| fig-cap: "Real data: training flow records and the sample cross-correlations of their daily changes. Positive lag h places the Keokuk change h days after the Clinton change. The strongest displayed correlation occurs at two days."
wb14_river_plot_training(river, initial = river_initial, H = 14)
```

```{r}
#| label: river-training-plot-code
#| eval: false
wb14_river_plot_training(river, initial = river_initial, H = 14)
```

Changes at Clinton are most strongly correlated with changes at Keokuk
about two days later. This suggests that Clinton's history may contain
information that Keokuk's history alone misses. The changes make the
timing visible; the forecasts below concern flow levels. We fit the
joint process $X_t=(X_{C,t},X_{K,t})^\top$, with Clinton first and Keokuk
second.

### From Residual Patterns to Lag Choice

Use orders $p\in\{1,2,7\}$, representing one day, two days, and a week
of observed history. Annual variation is represented by harmonic
regression: an intercept and sine–cosine pairs with periods of one year
and half a year. The fitted equation is
$$
X_t=b(t)+\sum_{i=1}^p A_iX_{t-i}+Z_t,
$$
where, taking $t=1$ on January 1, 2010,
$$
b(t)=b_0+\sum_{j=1}^2\left\{
 a_j\sin\!\left(\frac{2\pi j(t-1)}{365.2425}\right)
 +b_j\cos\!\left(\frac{2\pi j(t-1)}{365.2425}\right)\right\}.
$$
The five coefficient vectors $b_0,a_1,b_1,a_2,b_2\in\R^2$ and the lag
matrices are estimated jointly. Every candidate uses the same 3,645
response days, January 8, 2010–December 31, 2019.

```{r}
#| label: river-initial-fits
river_orders <- c(1, 2, 7)
river_fits <- setNames(lapply(river_orders, function(p)
  wb14_fit_river(river, p = p, origin = river_initial)),
  paste0("VAR(", river_orders, ")"))
river_table <- wb14_fit_table(river_fits, H = 28)
```

```{r}
#| label: fig-river-residuals
#| echo: false
#| fig-height: 6.2
#| fig-cap: "Real data: residual correlations through 28 days, with common scales across orders. The top row pairs later Keokuk errors with earlier Keokuk errors; the bottom row pairs them with earlier Clinton errors. Adding lags reduces both patterns. The joint statistic also uses the other entries of the residual covariance matrices."
wb14_river_plot_residuals(river_fits, H = 28)
```

```{r}
#| label: river-residual-plot-code
#| eval: false
wb14_river_plot_residuals(river_fits, H = 28)
```

The one-lag model leaves a pronounced pattern of correlated errors.
Two lags remove much of it, and seven reduce it further. For the joint
check, use $H=28$ and reference degrees of freedom $4(28-p)$ from
Definition 14.3.

```{r}
#| label: river-diagnostic-table
river_diagnostics <- river_table[, c("Model", "Q", "df", "p_value")]
river_diagnostics$p_value <- format.pval(river_diagnostics$p_value,
                                       digits = 3, eps = 1e-6)
knitr::kable(row.names = FALSE, river_diagnostics,
             digits = c(0, 1, 0, 3),
             col.names = c("Model", "Q28", "Reference df", "p-value"))
```

All three joint tests reject the zero-covariance null through lag 28 at
the 5% level. For VAR(7), $Q_{28}$ is about 140.3 on 84 reference degrees of
freedom, with p-value about 0.000114. Serial dependence has been greatly
reduced, but some remains. The next comparison asks whether the gain in
fitted likelihood pays for the additional coefficients.

Each lag adds four coefficients. The five deterministic regressors
contribute ten, and the innovation covariance contributes three, so
$$
k=4p+10+3=4p+13.
$$
The counts are 17, 21, and 41 for one, two, and seven lags. The figure
subtracts the smallest AIC from each AIC and the smallest BIC from each
BIC, so zero identifies each criterion's preferred order.

```{r}
#| label: fig-river-information-criteria
#| echo: false
#| fig-height: 4
#| fig-cap: "Real data: information-criterion differences on the common initial response window. Both criteria select seven lags. Within each panel, the height above zero is the score lost relative to that criterion's minimum."
wb14_river_plot_ic(river_table)
```

```{r}
#| label: river-information-criterion-code
#| eval: false
wb14_river_plot_ic(river_table)
```

Moving from two to seven lags adds twenty coefficients. The likelihood
improvement is large enough to offset their cost under both AIC and
BIC; Exercise 14.6 checks the calculation. Thus the criteria favor the
same order that most reduced the residual pattern. Does it also improve
forecasts appreciably beyond the two-lag model?

### More Lags, or More Information?

Starting on December 31, 2019, expand the training period one day at a
time and forecast one, three, and seven days ahead. Refit every model at
each origin, including its seasonal coefficients. First compare the
three joint VAR orders. Then compare VAR(7) with a Keokuk-only AR(7)
using the same seasonal regressors. This second comparison holds the
lag order fixed while adding Clinton's history. The joint forecasts
recursively predict both gauges; they never use future observed
Clinton flow.

Persistence predicts the last observed Keokuk flow at every horizon.
A second benchmark uses Keokuk's flow on the same date of the preceding
year; February 29 uses February 28 of the preceding year.

The record has one missing Keokuk value, August 31, 2022. Keep that
calendar day missing. At each origin, all candidates fit the same rows
whose response and preceding seven daily vectors are observed. The
common forecast origins require the seven-day starting history, the
requested targets, and the preceding-year benchmark values to be
available. This leaves 1,808 common origins through December 24, 2024.

```{r}
#| label: river-forecast-evaluation
river_eval <- wb14_evaluate_river(river, initial = river_initial,
                                 orders = river_orders, horizons = c(1, 3, 7))
stopifnot(nrow(river_eval$failures) == 0,
          all(is.finite(river_eval$scores$RMSE)))
```

```{r}
#| label: fig-river-forecast-comparisons
#| echo: false
#| fig-height: 4.6
#| fig-cap: "Real data: Keokuk forecast RMSE on identical origins. Left: extending the joint VAR from two to seven lags gives only a small further improvement. Right: adding upstream history at a fixed seven-lag order gives a larger improvement over the Keokuk-only model. Persistence provides a common simple benchmark."
wb14_river_plot_rmse(river_eval)
```

```{r}
#| label: river-forecast-comparison-code
#| eval: false
wb14_river_plot_rmse(river_eval)
```

```{r}
#| label: river-forecast-score-table
river_methods <- c("VAR2", "VAR7", "AR7", "Persistence", "Previous year")
river_scores <- river_eval$scores
river_score_table <- data.frame(Model = c("VAR(2)", "VAR(7)", "AR(7)",
                                         "Persistence", "Previous year"))
for (h in c(1, 3, 7)) {
  a <- river_scores[river_scores$h == h, ]
  river_score_table[[paste(h, "day RMSE")]] <-
    a$RMSE[match(river_methods, a$model)]
}
knitr::kable(row.names = FALSE, river_score_table, digits = 2)
```

At three days, moving from VAR(2) to VAR(7) reduces RMSE from about
9.64 to 9.49, a 1.5% gain. Comparing VAR(7) with AR(7) instead reduces
RMSE from 11.62 to 9.49, an 18.3% gain. The same upstream information
reduces seven-day RMSE by about 12.6%. Both joint AIC and BIC select
seven lags at every origin, so their selected-model forecasts coincide
with the VAR(7) row.

```{r}
#| label: fig-river-forecast-path
#| echo: false
#| fig-height: 4.3
#| fig-cap: "Real data: three-day Keokuk forecasts during the first evaluation year, plotted at their target dates. Each prediction uses observations only through its origin, three days before the target. All three procedures track the broad spring rise; their differences are more visible around abrupt changes."
wb14_river_plot_path(river_eval, year = 2020, h = 3)
```

```{r}
#| label: river-forecast-path-code
#| eval: false
wb14_river_plot_path(river_eval, year = 2020, h = 3)
```

Clinton's history improves the Keokuk forecasts at all three horizons
in this evaluation. Most of the improvement over the scalar model is
already obtained with two joint lags; the extra five lags improve the
joint fitting likelihood much more markedly than Keokuk forecast RMSE.
Exercise 14.6
checks whether the three-day advantage of using upstream history also
appears within each evaluation year. $\diamond$

Residual diagnostics assess particular implications of a fitted model;
information criteria charge for its estimated parameters; chronological
forecast errors assess its later performance for a stated target.
A point forecast still leaves the range of plausible future
observations unspecified. The workbook on forecast uncertainty develops
prediction intervals and the additional uncertainty from estimation and
model choice.

## Exercises

::: {.exercise}
**Exercise 14.1 (White coordinates, dependent vectors).**
Let $Z_t\overset{\mathrm{iid}}{\sim}N_2(0,I_2)$ and define
$$
E_{1,t}=Z_{1,t}+aZ_{2,t-1},\qquad E_{2,t}=Z_{2,t},
\qquad a\ne0.
$$
Treat $(E_t)$ as a possible residual process.

(a) Find its mean and covariance matrix at lag zero. Explain why the
vector process is weakly stationary.

(b) Show that both coordinates have zero autocovariance at every
nonzero lag. Calculate $\Gamma_E(1)=\Cov(E_{t+1},E_t)$ and identify
which earlier error predicts which later error.

(c) Set $a=0.8$. Generate 1,500 retained vectors using seed 542141,
with an extra initial innovation vector for the lag. Plot the matrix
of sample lag correlations through lag 8. Compare the relevant
lag-one correlation with its population value. What would be missed
by plotting only the diagonal panels?
:::


::: {.exercise}
**Exercise 14.2 (How much improvement pays for a lag?).**
Two trivariate VARs, with an intercept and unrestricted innovation
covariance, are fitted to $T_*=240$ common response vectors. Moving
from order two to order three increases the maximized log likelihood
by 15.

(a) Count the free parameters in each model, including the covariance.
Determine which model AIC prefers and which BIC prefers.

(b) For two models with parameter-count difference $\Delta k>0$,
derive the improvement in maximized log likelihood needed for the
larger model to have smaller AIC, and then smaller BIC.

(c) An identified trivariate VARMA(1,1) has a free intercept and all
entries of its two coefficient matrices free. Compare its parameter
count with VAR(2). Why must identification be checked before using
this count?

(d) Consider the rank-one trivariate VECM
$$
\nabla V_t=\alpha(\beta^\top V_{t-1}+\rho)
             +B_1\nabla V_{t-1}+Z_t,
$$
where $\alpha,\beta\in\R^3$, $B_1\in\R^{3\times3}$, and
$\rho\in\R$. The first entry of $\beta$ is fixed at one; the
innovation covariance is unrestricted. Count the free parameters,
including the covariance, and compare with an unrestricted VAR(2).
:::


::: {.exercise}
**Exercise 14.3 (Dependence left in the squares).**
Let $(\varepsilon_t)$ and $(\eta_t)$ be independent sequences of
independent $N(0,1)$ variables. Define
$$
Z_t=(\varepsilon_t\varepsilon_{t-1},\eta_t)^\top.
$$

(a) Show that $(Z_t)$ is weakly stationary with mean zero and covariance
$I_2$, and that every nonzero-lag covariance matrix is zero.

(b) Calculate $\Cov(Z_{1,t}^2,Z_{1,t-1}^2)$. Use the result to decide
whether the vectors are independent across time.

(c) Explain why a residual correlation plot alone could miss this
feature and why the usual i.i.d. portmanteau reference is not justified
merely by the conclusion in part (a). Name a residual plot directed
at the feature in part (b).
:::


::: {.exercise}
**Exercise 14.4 (What does the SOI add?).**
Use `fish_eval$scores` for the RMSE summaries and
`fish_eval$forecasts` for individual forecasts, keeping the origins
and horizons of Example 14.7. The latter records the origin, target
coordinate, model, horizon `h`, observation `actual`, prediction
`forecast`, and observed-minus-predicted `error`.

(a) Compare fixed-order VAR(6) with recruitment-only AR(6) at horizons
1, 6, and 12. Report RMSE and the percentage reduction from including
SOI at each horizon. Why is this a different comparison from VAR-BIC
against AR-BIC?

(b) Repeat the fixed-order comparison for the SOI target, using a
SOI-only AR(6) as benchmark. Does the predictive improvement run in
both directions?

(c) For the first origin, compare the VAR(6) forecasts in the main
text with the observations at horizons 1, 6, and 12. Report the three
recruitment errors. Explain why these three errors alone cannot
replace the expanding-origin comparison.

(d) Suppose the six-month errors suggest adding a nonlinear SOI
term. Explain how fitting that revised model on the initial period
and selecting the nonlinear term using the displayed evaluation
scores would change the role of the evaluation period.
:::


::: {.exercise}
**Exercise 14.5 (Does the better VARMA likelihood predict better?).**
Use the macro fits in Example 14.6. The `MTS` object is available as
`macro_assessment$fits$VARMA11$object`; its fields `Ph0`, `Phi`, and
`Theta` give the intercept, AR matrix, and *minus-sign* MA matrix.
Its last residual estimates the final innovation.

(a) Compute the one- and four-quarter forecasts from the initial
origin by recursively applying the fitted VARMA equation. State
which innovation values enter each step.

(b) Use `wb14_macro_forecasts()` to evaluate VAR(1), VAR(2), and
VARMA(1,1) at 51 common expanding origins, $t=105,\ldots,155$.
Report the RMSE for each coordinate at horizons one and four, with
all models refitted at each origin. The returned `$scores` table
contains RMSE; check `$failures` before interpreting it.

(c) Compare the result with the initial AIC/BIC table. Why can the
criteria and the later errors favor different models even though
VAR(2) and VARMA(1,1) have the same parameter count?
:::


::: {.exercise}
**Exercise 14.6 (More memory or an upstream measurement?).**
Use `river_table` and `river_eval` from Example 14.8.

(a) Explain the twenty additional parameters in VAR(7) compared with
VAR(2). Use the `logLik` column of `river_table` to calculate the
improvement in fitted log likelihood. Does it exceed the thresholds
required by AIC and BIC on the 3,645 common response days?

(b) Use `river_eval$scores` at seven days. Calculate the percentage
RMSE reduction from VAR(2) to VAR(7), and from AR(7) to VAR(7).
Explain the different question answered by each comparison.

(c) The `river_eval$forecasts` table records `model`, horizon `h`,
`target_date`, and observed-minus-predicted `error`. For three-day
forecasts, calculate the RMSE of `VAR7` and `AR7` separately in each
target year, 2020–2024. Report the five pairs and their percentage
reductions. Is the pooled improvement confined to one year?
*Hint: `format(target_date, "%Y")` extracts the target year.*
:::


## Further Reading

The multivariate residual diagnostic and its adjusted form are
documented in [`vars::serial.test`](https://search.r-project.org/CRAN/refmans/vars/html/serial.html).
Hyndman and Athanasopoulos (2021),
[*Forecasting: Principles and Practice*, Section 5.10](https://otexts.com/fpp3/tscv.html),
discuss evaluation from successive forecast origins. Shumway and
Stoffer (2025), *Time Series Analysis and Its Applications*, Sections
1.1 and 2.2, introduce the SOI–recruitment data and their dependence.