Workbook 14 — Assessing Fitted Vector Models

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

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.

Figure 1: 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")

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

Figure 2: 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")

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.

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

Q12	Reference df	p-value
44.323	40	0.294

For this fit, the p-value is 0.294, 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 and astsa::rec 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?

fish <- wb14_fish_data()

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

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.

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

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

Figure 5: 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(2)"]], H = 24)
wb14_plot_correlations(fish_fits[["VAR(6)"]], H = 24)

Figure 6: 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")

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.

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

Model	Q24	Reference df	p-value
VAR(1)	293.0	92	<1e-04
VAR(2)	246.3	88	<1e-04
VAR(6)	63.1	72	0.763
VAR(12)	34.4	48	0.930

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.

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

Model	Response months	Parameters	AIC	BIC
VAR(1)	290	31	2152.0	2265.8
VAR(2)	290	35	2128.6	2257.0
VAR(6)	290	51	1947.1	2134.2
VAR(12)	290	75	1956.9	2232.2

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.

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

Model	Response quarters	Parameters	AIC	BIC
VAR1	103	9	395.5	419.2
VAR2	103	13	384.8	419.1
VARMA11	103	13	363.2	397.4

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.

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)

Horizon	SOI	Recruitment
1	0.11	52.14
6	-0.15	40.64
12	0.38	79.78

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.

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

Figure 7: 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")

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)

Procedure	1 month RMSE	6 month RMSE	12 month RMSE
VAR-BIC	9.04	18.12	20.78
AR-BIC	9.58	24.39	27.39
VAR-AIC	8.79	18.55	23.17
Seasonal mean	26.76	26.16	25.76

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.

Figure 8: 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"))

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 and Keokuk, 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 cover 2010–2024; source details and quality flags 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.

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

Figure 9: 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)

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.

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)

Figure 10: 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)

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.

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

Model	Q28	Reference df	p-value
VAR(1)	5156.4	108	< 1e-06
VAR(2)	468.3	104	< 1e-06
VAR(7)	140.3	84	0.000114

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.

Figure 11: 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)

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.

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

Figure 12: 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)

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)

Model	1 day RMSE	3 day RMSE	7 day RMSE
VAR(2)	5.29	9.64	17.31
VAR(7)	5.28	9.49	17.20
AR(7)	5.84	11.62	19.68
Persistence	5.74	11.98	20.86
Previous year	60.64	60.63	60.65

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.

Figure 13: 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)

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

Find its mean and covariance matrix at lag zero. Explain why the vector process is weakly stationary.
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.
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 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.

Count the free parameters in each model, including the covariance. Determine which model AIC prefers and which BIC prefers.
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.
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?
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 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. \]

Show that $(Z_t)$ is weakly stationary with mean zero and covariance $I_2$, and that every nonzero-lag covariance matrix is zero.
Calculate $\Cov(Z_{1,t}^2,Z_{1,t-1}^2)$. Use the result to decide whether the vectors are independent across time.
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 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.

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

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.
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.
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 14.6 (More memory or an upstream measurement?). Use river_table and river_eval from Example 14.8.

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