Transfer Function Modelling

Libraries

library(rmarkdown)
library(knitr)
library(datasets)    
library(arfima)      
library(astsa)       
library(tidyverse)   
library(patchwork)   
library(forecast)    
library(tseries)     
library(lmtest)      
library(TSA)         

Introduction

This notebook analyses five classic data sets in which one series (the output Y_t) responds dynamically to one or more other series (the inputs X_t). Each case is modelled with a transfer function model

Y_t = \mu + \nu(B)\,X_t + \eta_t = \mu + \frac{\omega_s(B)}{\delta_r(B)}\,B^b X_t + \frac{\theta_q(B)}{\phi_p(B)}\,a_t, \qquad a_t \sim \text{WN}(0,\sigma_a^2),

where \nu(B) = \nu_0 + \nu_1 B + \nu_2 B^2 + \cdots holds the impulse response weights (\nu_k is the response of Y_t to X_{t-k}), b is the delay, \omega_s(B) = \omega_0 - \omega_1 B - \cdots - \omega_s B^s, \delta_r(B) = 1 - \delta_1 B - \cdots - \delta_r B^r, and the noise \eta_t follows an ARMA(p,q) model. In practice Y_t and X_t are the transformed or differenced series. With more than one input, one term \nu_j(B) X_{j,t} is added for every input.

Case	Data	Output Y_t	Input(s) X_t	Frequency, n	Approach
1	Box–Jenkins Series M (BJsales)	sales	leading indicator	unspecified, 150	CCF
2	Box–Jenkins Series J (gas furnace)	CO_2 in outlet gas (%)	methane feed rate	9 seconds, 296	CCF
3	SOI and Recruitment (astsa)	number of new fish	Southern Oscillation Index	monthly, 453	CCF
4	LA pollution study (astsa)	cardiovascular mortality	temperature, particulates	weekly, 508	LTF
5	Lake Shasta (astsa::climhyd)	reservoir inflow	precipitation, temperature	monthly, 454	LTF

CCF (prewhitening) approach, used for the single-input cases 1–3:

1. Stationarity: apply the same transformation and differencing to Y and X.
2. Input model: identify an ARIMA model for the input only, \phi_x(B)x_t = \theta_x(B)\alpha_t.
3. Prewhitening: filter both series with the same filter, \alpha_t = \pi(B)x_t and \beta_t = \pi(B)y_t, where \pi(B) = \phi_x(B)/\theta_x(B).
4. CCF: compute r_{\alpha\beta}(k) and the impulse response weights \hat\nu_k = r_{\alpha\beta}(k)\, s_\beta / s_\alpha.
5. Orders (b, r, s) from the pattern of \hat\nu_k.
6. Preliminary estimates \hat\omega, \hat\delta by solving the \nu-recursion.
7. Noise model: identify ARMA(p,q) for the preliminary noise \hat\eta_t.
8. Estimation of all parameters jointly by maximum likelihood.
9. Diagnostic checking: Q for the residuals \hat a_t and Q_0 for the cross-correlation between \alpha_t and \hat a_t.
10. Explicit equation Y_t = f(\text{lags of } Y, X, a), written in LaTeX and in the original variables: differencing, seasonal adjustment, centring and transformations are undone.
11. Forecasting: one-step-ahead and k-step-ahead forecasts from the explicit equation, with their accuracy (see Forecasting below).

LTF (linear transfer function) approach (Liu & Hanssens, 1982; Pankratz, 1991), used for the multiple-input cases 4–5:

1. Free-form fit: regress Y_t on lags 0, \dots, k of every input with a proxy AR(1) noise.
2. Stationarity: if the errors look non-stationary (proxy \hat\phi_1 close to 1), difference all series.
3. Orders (b, r, s) for every input from the pattern of the free-form weights and their t-ratios.
4. Noise model: identify ARMA(p,q) for the free-form errors \hat\eta_t.
5. Rational fit: estimate the transfer functions and the ARMA noise jointly by ML.
6. Diagnostics: Ljung–Box on \hat a_t and Q_0 with every prewhitened input.
7. Explicit equation, in the original variables.
8. Forecasting: one-step-ahead and k-step-ahead forecasts from the explicit equation, with their accuracy (see Forecasting below).

In every case the explicit equation is derived step by step: the operators are shown factor by factor and multiplied out, the explicit equation of the series used in the model is written down, and the definitions of these series (differences, deviations from the monthly means, centred values) are substituted, which gives the explicit equation of the original variables.

The forecasting step (Step 11 of Cases 1–3, Step 8 of Cases 4–5) has the same five parts in every case:

1. Forecast equations: the one-step-ahead and the k-step-ahead forecast equation, written out from the explicit equation.
2. Training and test periods: the model is re-estimated on a training period, and the last part of the data is kept as a test period. The benchmark is a univariate ARIMA for the output.
3. One-step-ahead forecasts of the test period: plot, and below it the accuracy (RMSE, MAE, MAPE).
4. k-step-ahead forecasts of the test period: the forecast path from the end of the training period (plot, then accuracy) and the rolling-origin evaluation (RMSE by horizon, plot, then table).
5. Forecasts beyond the end of the data: plot, and below it the forecasts with the RMSE that the test period shows for the same horizon.

The LTF approach is preferred with several inputs: the CCF approach needs a separate prewhitening filter for every input, and correlated inputs distort the individual CCFs, whereas the free-form regression estimates all lags of all inputs jointly.

Forecasting. In explicit form the model reads

Y_t = c + \sum_{i} \pi_i Y_{t-i} + \sum_{j}\sum_{k} \omega^*_{jk} X_{j,t-k} + a_t + \sum_{l} \theta^*_l a_{t-l}.

The forecast of Y_{n+h} made at the origin n replaces every unknown quantity on the right-hand side by its forecast:

\hat Y_n(h) = c + \sum_{i} \pi_i \hat Y_n(h-i) + \sum_{j}\sum_{k} \omega^*_{jk} \hat X_{j,n}(h-k) + \sum_{l} \theta^*_l \hat a_n(h-l),

with \hat Y_n(m) = Y_{n+m}, \hat X_{j,n}(m) = X_{j,n+m} and \hat a_n(m) = a_{n+m} for m \le 0 (observed values) and \hat a_n(m) = 0 for m > 0. A one-step-ahead forecast (h = 1) is made at every origin with the newly observed data. k-step-ahead forecasts (h = 1, \dots, k from one origin) are computed recursively: the forecasts \hat Y_n(1), \hat Y_n(2), \dots take the place of the unknown future values of Y.

Because the input enters as X_{t-b}, X_{t-b-1}, \dots, the forecasts up to h = b use only observed inputs. For h > b, and for every h when b = 0, the inputs themselves have to be forecast. For every input the forecasting method is selected automatically by time-series cross-validation (select_auto() below). The candidates are ARIMA (auto.arima()), exponential smoothing (ets()) and the Theta method (thetaf()), together with two simple benchmarks, the naive forecast (the last value) and the mean of the series; for the weekly, seasonal inputs of Case 4 they are ARIMA with Fourier terms for the annual cycle, STL decomposition with ETS (stlf()) and the seasonal naive forecast. From ten forecast origins spread over the last 20% of the input series, every candidate is re-estimated on the data up to the origin and forecasts the horizons that are needed later (up to K - b). The method with the smallest RMSE is re-estimated on the whole series and used.

Two versions of the k-step-ahead forecasts are compared on the test periods: ex ante forecasts use forecast inputs, as is necessary in real time, and ex post forecasts use the observed future inputs, i.e. they show what the model achieves when the inputs are known. The forecasts are evaluated in two ways. The rolling-origin evaluation repeats them from every origin of the test period and reports the RMSE separately for every horizon. The single forecast path is the set of forecasts made once, at the end of the training period, for the whole test period, as a forecaster would make them in practice. A single path depends strongly on the situation at its origin, so conclusions rest on the rolling-origin results.

Identification rules. The impulse response weights are read as follows (the approximate 95% bands of the CCF are \pm 2/\sqrt{n}):

Feature of \hat\nu_k	Order
number of initial weights that are zero	delay b
number of unpatterned weights \nu_b, \dots, \nu_{b+s-r} before the decay starts, s - r + 1	numerator order s
no decay (cut-off) / exponential decay / damped sine	r = 0 / r = 1 / r = 2

Sign conventions (read first). ccf(a, b) at lag k estimates \text{cor}(a_{t+k}, b_t), so ccf(beta, alpha) puts the response of y to past x at positive lags. R writes MA polynomials as 1 + \theta_1 B + \cdots, whereas Box and Jenkins write 1 - \theta_1 B - \cdots. In TSA::arimax(), transfer = list(c(r, s)), X-AR1 = \delta_1 in 1 - \delta_1 B, and X-MA0, X-MA1, … are the numerator coefficients with plus signs, \omega_0 + \omega_1 B + \cdots (so X-MA1 = -\omega_1 in the Box–Jenkins notation). arimax() has no delay argument, so the lagged input x_{t-b} is built by hand and the first b observations are dropped.

All tests use a significance level of \alpha = 5\%.

Helper Functions

The first helpers are the ones used in the lecture: polynomial product, explicit form of a fitted transfer function, printing of the explicit equation, and point forecasts from the explicit equation. show_explicit() is extended so that it prints long equations over several lines and also handles several inputs and a trend; num() prints small coefficients with four significant digits.

## polynomial in B = coefficient vector c(c0, c1, c2, ...)
polymul <- function(p, q) {
  out <- numeric(length(p) + length(q) - 1)
  for (i in seq_along(p)) for (j in seq_along(q))
    out[i + j - 1] <- out[i + j - 1] + p[i] * q[j]
  out
}
## (1-B)^d Y = mu + w/delta (1-B)^d X[t-b] + theta/phi a
## w: plus signs; delta, phi: c(1, -coef); theta: c(1, th1)
tf_explicit <- function(w, delta = 1, b = 0, phi = 1,
                        theta = 1, d = 0, mu = 0) {
  D <- 1                                   # (1-B)^d
  if (d > 0) for (i in 1:d) D <- polymul(D, c(1, -1))
  Ay <- polymul(polymul(delta, phi), D)    # acts on Y_t
  Cx <- polymul(polymul(w, phi), D)        # acts on X_{t-b}
  Da <- polymul(delta, theta)              # acts on a_t
  list(const = mu * sum(delta) * sum(phi),
       y = -Ay[-1], ylag = seq_along(Ay[-1]),
       x = Cx,      xlag = b + seq_along(Cx) - 1,
       a = Da[-1],  alag = seq_along(Da[-1]),
       spec = list(inputs = list(list(name = "X", w = unname(w), delta = unname(delta), b = b)),
                   phi = unname(phi), theta = unname(theta), d = d, mu = unname(mu),
                   beta = 0))                # the factors, for show_operators()
}
## number format: 4 decimals, but 4 significant digits for coefficients smaller than 0.1
num <- function(v, digits = 4) {
  d <- pmin(8, pmax(digits, digits - floor(log10(pmax(abs(v), 1e-12))) - 1))
  unname(mapply(function(x, dd) formatC(x, format = "f", digits = dd), v, d))
}
## print the explicit equation over several lines (one or several inputs, optional trend)
show_explicit <- function(e, digits = 4, per_line = 4) {
  g <- function(v) num(v, digits)
  trm <- function(v, s, lag)
    if (length(v)) paste0(ifelse(v < 0, " - ", " + "), g(abs(v)), " ", s,
                          ifelse(lag == 0, "[t]", paste0("[t-", lag, "]"))) else character(0)
  inputs <- if (is.null(e$inputs)) list(list(name = "X", x = e$x, xlag = e$xlag)) else
    e$inputs
  slope  <- if (is.null(e$slope)) 0 else e$slope
  trend  <- if (slope != 0) paste0(ifelse(slope < 0, " - ", " + "), g(abs(slope)), " t")
  groups <- c(list(c(trend, trm(e$y, "Y", e$ylag))),
              lapply(inputs, function(z) trm(z$x, z$name, z$xlag)),
              list(c(" + a[t]", trm(e$a, "a", e$alag))))
  lines <- unlist(lapply(groups, function(gr)
    vapply(split(gr, ceiling(seq_along(gr) / per_line)), paste, character(1), collapse = "")))
  cat(paste0("Y[t] = ", g(e$const), lines[1]), paste0("      ", lines[-1]), sep = "\n")
}
## point forecasts from the explicit equation (future a = 0);
## X must already contain the future values that are needed
explicit_forecast <- function(e, Y, X, a, h) {
  n <- length(Y)
  Y <- c(as.numeric(Y), rep(NA, h))
  a <- c(as.numeric(a), rep(0, h))
  for (t in (n + 1):(n + h))
    Y[t] <- e$const + sum(e$y * Y[t - e$ylag]) +
            sum(e$x * X[t - e$xlag]) + sum(e$a * a[t - e$alag])
  Y[(n + 1):(n + h)]
}

The additional helpers below automate the remaining steps: stationarity tests, the prewhitening filter, the CCF table, preliminary estimates, the Q_0 test, implied impulse response weights, the free-form lag matrix, the explicit equation for several inputs, and one-step-ahead forecasts.

## ADF and PP (H0: unit root) and KPSS (H0: level-stationary) in one table
stat_tests <- function(...) {
  s <- list(...)
  map_dfr(names(s), function(nm) {
    z   <- na.omit(as.numeric(s[[nm]]))
    adf <- suppressWarnings(adf.test(z)); pp <- suppressWarnings(pp.test(z))
    kp  <- suppressWarnings(kpss.test(z))
    tibble(series = nm, `ADF stat` = unname(adf$statistic), `ADF p` = adf$p.value,
           `PP p` = pp$p.value, `KPSS stat` = unname(kp$statistic), `KPSS p` = kp$p.value)
  }) %>% mutate(across(where(is.double), ~ round(.x, 4)))
}

## prewhitening filter pi(B) = phi(B)/theta(B), with R's MA sign (1 + theta B);
## the first length(phi) values are lost (NA)
pw_filter <- function(z, phi = numeric(0), theta = numeric(0)) {
  z <- as.numeric(z)
  w <- if (length(phi)) as.numeric(stats::filter(z, c(1, -phi), sides = 1)) else z
  if (length(theta)) {
    ok <- !is.na(w)
    w[ok] <- as.numeric(stats::filter(w[ok], -theta, method = "recursive"))
  }
  w
}

## CCF of the prewhitened series: r(k) = cor(beta[t+k], alpha[t]), nu_k = r(k) s_beta / s_alpha
ccf_table <- function(beta, alpha, lag.max = 15) {
  ok <- !is.na(alpha) & !is.na(beta)
  cc <- ccf(beta[ok], alpha[ok], lag.max = lag.max, plot = FALSE)
  tibble(k = as.numeric(cc$lag), r = as.numeric(cc$acf)) %>%
    mutate(nu = r * sd(beta[ok]) / sd(alpha[ok]),
           bound = 2 / sqrt(sum(ok)),
           significant = abs(r) > bound)
}
plot_ccf <- function(tab, title) {
  ggplot(tab, aes(k, r, fill = significant)) +
    geom_col(width = 0.6) +
    geom_hline(yintercept = c(-1, 1) * tab$bound[1], linetype = "dashed", colour = "darkorange") +
    geom_hline(yintercept = 0) +
    scale_fill_manual(values = c(`FALSE` = "grey70", `TRUE` = "steelblue"), guide = "none") +
    labs(title = title, x = "lag k", y = "cross-correlation") +
    theme_bw()
}

## preliminary estimates from the nu-recursion (Box-Jenkins signs:
## omega(B) = w0 - w1 B - ..., delta(B) = 1 - d1 B - ...); nu = weights at lags 0, 1, 2, ...
prelim_tf <- function(nu, b, r, s) {
  v <- function(j) if (j < b) 0 else nu[j + 1]
  d <- numeric(0)
  if (r > 0) {
    A   <- t(sapply(1:r, function(i) sapply(1:r, function(l) v(b + s + i - l))))
    rhs <- sapply(1:r, function(i) v(b + s + i))
    d   <- solve(matrix(A, r, r), rhs)
  }
  w <- numeric(s + 1)
  for (j in 0:s) {
    ar_part  <- if (r > 0) sum(d * sapply(1:r, function(l) v(b + j - l))) else 0
    w[j + 1] <- if (j == 0) v(b) else ar_part - v(b + j)
  }
  c(setNames(w, paste0("omega", 0:s)), if (r > 0) setNames(d, paste0("delta", 1:r)))
}

## cross-correlation check of the transfer function: residuals vs prewhitened input, M = r + s + 1
Q0_test <- function(a.hat, alpha, K = 12, M = 1) {
  ok <- !is.na(a.hat) & !is.na(alpha)
  cc <- ccf(a.hat[ok], alpha[ok], lag.max = K, plot = FALSE)
  rk <- cc$acf[cc$lag >= 0]; m <- sum(ok)
  Q0 <- m * (m + 2) * sum(rk^2 / (m - 0:K))
  tibble(K = K, Q0 = round(Q0, 4), df = K + 1 - M,
         `p-value` = round(1 - pchisq(Q0, K + 1 - M), 4))
}

## transfer-function component u_t = w(B)/delta(B) x_t (w with plus signs, zero starting values)
tf_component <- function(x, w, delta = 1) {
  x <- as.numeric(x); s <- length(w) - 1
  u <- if (s > 0) as.numeric(stats::filter(c(rep(0, s), x), w, sides = 1))[-(1:s)] else w * x
  if (length(delta) > 1) u <- as.numeric(stats::filter(u, -delta[-1], method = "recursive"))
  u
}

## impulse response weights implied by a fitted rational transfer function
nu_weights <- function(w, delta = 1, b = 0, K = 20) {
  nu <- numeric(K + 1)
  for (k in 0:K) {
    j   <- k - b
    val <- if (j >= 0 && j < length(w)) w[j + 1] else 0
    if (length(delta) > 1) for (l in 2:length(delta))
      if (k - l + 1 >= 0) val <- val - delta[l] * nu[k - l + 2]
    nu[k + 1] <- val
  }
  nu
}

## seasonal adjustment: subtract the monthly means (estimated from the first n_train values)
desea <- function(z, n_train = length(z)) {
  m  <- cycle(z)
  mm <- tapply(as.numeric(z)[1:n_train], m[1:n_train], mean)
  as.numeric(z) - as.numeric(mm[m])
}

## free-form (LTF) lag matrix: columns z_t, z_{t-1}, ..., z_{t-k}
lagmat <- function(z, k, nm) {
  L <- embed(as.numeric(z), k + 1)
  colnames(L) <- paste0(nm, "_", 0:k)
  L
}
## free-form weights with standard errors and t-ratios
ff_weights <- function(fit, inputs) {
  cf <- coef(fit); se <- sqrt(diag(fit$var.coef))
  map_dfr(inputs, function(nm) {
    id <- grep(paste0("^", nm, "_[0-9]+$"), names(cf))
    tibble(input = nm, k = as.numeric(sub(".*_", "", names(cf)[id])),
           nu = unname(cf[id]), se = unname(se[id]), t = nu / se)
  })
}
plot_ff <- function(w, title, implied = NULL) {
  p <- ggplot(w, aes(k, nu)) +
    geom_col(fill = "steelblue", width = 0.6) +
    geom_errorbar(aes(ymin = nu - 2 * se, ymax = nu + 2 * se), width = 0.2) +
    geom_hline(yintercept = 0) +
    facet_wrap(~ input, scales = "free_y") +
    labs(title = title, x = "lag k", y = "weight") +
    theme_bw()
  if (!is.null(implied))
    p <- p + geom_line(data = implied, aes(k, nu), colour = "darkorange", linetype = "dashed") +
             geom_point(data = implied, aes(k, nu), colour = "darkorange")
  p
}

## explicit equation for several inputs (d = 0, optional linear trend beta * t)
## inputs: list of list(name, w (plus signs), delta = c(1, -d1, ...), b)
tf_explicit_multi <- function(inputs, phi = 1, theta = 1, mu = 0, beta = 0) {
  dprod <- Reduce(polymul, lapply(inputs, `[[`, "delta"), 1)
  Ay <- polymul(dprod, phi)
  xs <- lapply(seq_along(inputs), function(j) {
    others <- Reduce(polymul, lapply(inputs[-j], `[[`, "delta"), 1)
    cx <- polymul(polymul(inputs[[j]]$w, others), phi)
    list(name = inputs[[j]]$name, x = cx, xlag = inputs[[j]]$b + seq_along(cx) - 1)
  })
  Da <- polymul(dprod, theta)
  list(const = mu * sum(Ay) - beta * sum((seq_along(Ay) - 1) * Ay),
       slope = beta * sum(Ay),
       y = -Ay[-1], ylag = seq_along(Ay[-1]), inputs = xs,
       a = Da[-1], alag = seq_along(Da[-1]),
       spec = list(inputs = lapply(inputs, function(z)
                     list(name = z$name, w = unname(z$w), delta = unname(z$delta), b = z$b)),
                   phi = unname(phi), theta = unname(theta), d = 0, mu = unname(mu),
                   beta = unname(beta)))     # the factors, for show_operators()
}

## one-step-ahead forecasts at the times idx from an explicit equation;
## a holds the residuals up to the forecast origin and is updated as we move on
onestep_explicit <- function(e, Y, X, a, idx) {
  Y <- as.numeric(Y); X <- as.numeric(X); a <- as.numeric(a)
  f <- numeric(length(idx))
  for (i in seq_along(idx)) {
    t <- idx[i]
    f[i] <- e$const + sum(e$y * Y[t - e$ylag]) + sum(e$x * X[t - e$xlag]) +
            sum(e$a * a[t - e$alag])
    a[t] <- Y[t] - f[i]
  }
  f
}
onestep_explicit_multi <- function(e, Y, X, a, idx) {   # X: named list of inputs
  Y <- as.numeric(Y); a <- as.numeric(a)
  f <- numeric(length(idx))
  for (i in seq_along(idx)) {
    t <- idx[i]
    xpart <- sum(sapply(e$inputs, function(z) sum(z$x * X[[z$name]][t - z$xlag])))
    f[i] <- e$const + e$slope * t + sum(e$y * Y[t - e$ylag]) + xpart +
            sum(e$a * a[t - e$alag])
    a[t] <- Y[t] - f[i]
  }
  f
}

## forecast accuracy: RMSE, MAE and MAPE (in %) of every forecast in ... against actual
accuracy_tbl <- function(actual, ...) {
  fc <- list(...)
  map_dfr(names(fc), function(nm) {
    e <- as.numeric(actual) - as.numeric(fc[[nm]])
    tibble(model = nm, RMSE = sqrt(mean(e^2)), MAE = mean(abs(e)),
           `MAPE (%)` = 100 * mean(abs(e / as.numeric(actual))))
  }) %>%
    mutate(across(where(is.double), ~ round(.x, 4)))
}
plot_holdout <- function(t, actual, ...) {
  fc <- list(...)
  d <- bind_rows(tibble(t = t, value = as.numeric(actual), series = "Actual"),
                 map_dfr(names(fc), ~ tibble(t = t, value = as.numeric(fc[[.x]]), series = .x)))
  ggplot(d, aes(t, value, colour = series, linetype = series)) + geom_line() +
    geom_point(size = 1) +
    scale_colour_manual(values = c(Actual = "black", `Transfer function` = "steelblue",
                                   `TF ex ante` = "steelblue", `TF ex post` = "forestgreen",
                                   `Univariate ARIMA` = "darkorange")) +
    scale_linetype_manual(values = c(Actual = "solid", `Transfer function` = "solid",
                                     `TF ex ante` = "solid", `TF ex post` = "dashed",
                                     `Univariate ARIMA` = "solid")) +
    labs(x = "t", y = NULL, colour = NULL, linetype = NULL) + theme_bw() +
    theme(legend.position = "bottom")
}

The function prelim_tf() solves the recursion \nu_j = \delta_1\nu_{j-1} + \cdots + \delta_r\nu_{j-r} + \omega_0 (for j = b), - \omega_{j-b} (for j = b+1, \dots, b+s), 0 (for j > b+s), with \nu_j = 0 for j < b. For (b,1,0), for instance, it returns \hat\omega_0 = \hat\nu_b and \hat\delta_1 = \hat\nu_{b+1}/\hat\nu_b.

The last group of helpers writes the explicit equation and produces and evaluates forecasts several steps ahead (see Forecasting in the introduction).

• show_operators() writes the operators A(B), C(B) and D(B) of A(B)\,Y_t = c + C(B)\,X_{t-b} + D(B)\,a_t three times: as a product of the model operators, with the estimated factors, and multiplied out; it also derives the constant c. show_equation() writes the explicit equation, and show_chain() writes several explicit equations below each other, which is used to go back to the original variables step by step. d = 1 in a step of show_chain() writes every lagged difference out, e.g. \nabla Y_{t-1} = Y_{t-1} - Y_{t-2}.
• The symbols of the variables are set with sym, either as a symbol, c(Y = "y^*") gives y^*_{t-1}, or as a template in which # stands for the time index, c(P = "\\sqrt{\\mathrm{Precip}_{#}}") gives \sqrt{\mathrm{Precip}_{t-1}}.
• monthly_const() gives the month-specific constants c_m of the explicit equation in the original variables when the series were seasonally adjusted, and dummy_terms() writes c_{m(t)} = \sum_{m=1}^{12} c_m D_{m,t} term by term with the seasonal dummies D_{m,t}.
• All of these functions write LaTeX (display math). In a chunk with results = "asis" the equations are rendered, both when the notebook is knitted and when the chunks are run in the RStudio notebook. show_explicit() (above) prints the same equation as plain text in the console.
• tf_forecast() computes \hat Y_o(1), \dots, \hat Y_o(H) from the explicit equation for one or several inputs. Without xfc it uses the observed future inputs (ex post); with xfc it uses the observed inputs up to the origin o and the forecasts returned by xfc(o, H) after it (ex ante).
• select_auto() chooses the forecasting method for an input series by time-series cross-validation; auto_methods and auto_methods_seasonal are the candidates. It returns the cross-validation table ($table), the selected model ($model) and its forecasting function ($f).
• show_forecast_eq() writes the forecast equations in LaTeX: the one-step-ahead forecast \hat Y_t(1) and the recursion for the k-step-ahead forecast \hat Y_t(k).
• rolling_eval() repeats the forecasts from every origin of a test period and reports the RMSE for every horizon h = 1, \dots, H, using all origins whose target o + h lies in the test period. The forecasts themselves are kept with the table; onestep_of() extracts those at h = 1, i.e. the one-step-ahead forecasts of the whole test period.
• accuracy_tbl() gives the RMSE, the MAE and the MAPE (in %) of one or several forecasts; acc_of() applies it to the one-step-ahead forecasts.
• plot_fc() plots forecasts of the test period with the actual values, plot_rmse() the RMSE of the rolling-origin forecasts by horizon, and plot_future() the forecasts beyond the end of the data. In the notebook every forecast plot is followed by the accuracy of the forecasts it shows.

## ---- the explicit equation in LaTeX ----

## explicit equation of a single-input model in the format of tf_explicit_multi()
as_multi <- function(e, name = "X") {
  if (!is.null(e$inputs)) return(e)
  sp <- e$spec
  if (!is.null(sp)) sp$inputs[[1]]$name <- name
  list(const = e$const, slope = 0, y = e$y, ylag = e$ylag,
       inputs = list(list(name = name, x = e$x, xlag = e$xlag)), a = e$a, alag = e$alag,
       spec = sp)
}
## polynomial c(p0, p1, ...) in B, as LaTeX or as text
poly_str <- function(p, digits = 4, latex = TRUE) {
  out <- ""
  for (i in seq_along(p)) {
    v <- p[i]; k <- i - 1
    if (abs(v) < 1e-12) next
    mag <- if (k == 0 && abs(v - 1) < 1e-10) "1" else
      if (k > 0 && abs(abs(v) - 1) < 1e-10) "" else num(abs(v), digits)
    pow <- if (k == 0) "" else if (k == 1) "B" else
      if (latex) paste0("B^{", k, "}") else paste0("B^", k)
    sgn <- if (out == "") (if (v < 0) "-" else "") else (if (v < 0) " - " else " + ")
    out <- paste0(out, sgn, mag, pow)
  }
  out
}
## one factor of an operator: (1 - 0.7264B), a constant such as 4.6942, or "" for 1
fac_str <- function(p, digits = 4) {
  p <- unname(as.numeric(p))
  if (length(p) == 1) return(if (abs(p - 1) < 1e-12) "" else num(p, digits))
  paste0("\\left(", poly_str(p, digits), "\\right)")
}
## symbol of the variable nm at time t - lag. sym holds a symbol, c(Y = "y^*") gives
## y^*_{t-1}, or a template in which # stands for the time index,
## c(P = "\\sqrt{\\mathrm{Precip}_{#}}") gives \sqrt{\mathrm{Precip}_{t-1}}
tex_at <- function(nm, ix, sym = NULL) {             # at the time index ix ("t+1", "t-2")
  s <- if (nm %in% names(sym)) sym[[nm]] else nm
  if (grepl("#", s, fixed = TRUE)) gsub("#", ix, s, fixed = TRUE) else
    paste0(s, "_{", ix, "}")
}
t_ix    <- function(s) if (s > 0) paste0("t+", s) else if (s == 0) "t" else paste0("t-", -s)
tex_var <- function(nm, lag, sym = NULL) tex_at(nm, t_ix(-lag), sym)
## forecast of nm made at the origin t, k - lag steps ahead: \hat{Y}_t(k-1)
tex_hat <- function(nm, arg, sym = NULL) {
  s <- if (nm %in% names(sym) && !grepl("#", sym[[nm]], fixed = TRUE)) sym[[nm]] else nm
  paste0(if (nchar(s) == 1) "\\hat{" else "\\widehat{", s, "}_t(", arg, ")")
}
## the d-th difference of nm at time t - lag, written out: Y_{t-1} - Y_{t-2} for d = 1
diff_str <- function(nm, lag, d, sym = NULL) {
  D <- 1
  for (i in seq_len(d)) D <- polymul(D, c(1, -1))
  out <- ""
  for (j in seq_along(D)) {
    mag <- if (abs(D[j]) == 1) "" else format(abs(D[j]))
    sgn <- if (j == 1) (if (D[j] < 0) "-" else "") else (if (D[j] < 0) " - " else " + ")
    out <- paste0(out, sgn, mag, tex_var(nm, lag + j - 1, sym))
  }
  out
}
## right-hand side of an explicit equation as LaTeX lines, grouped by variable.
## const: LaTeX terms that replace the constant, e.g. dummy_terms(); d > 0: every lag of Y and
## of the inputs is a d-th difference, written out with diff_str()
## vfun(nm, lag, dd): other symbols for the lagged terms (used for the forecast equations);
## tsym: symbol of the trend variable; a0 = FALSE leaves out the current shock a_t
eq_rhs <- function(e, sym = NULL, const = NULL, d = 0, digits = 4, per_line = 5,
                   vfun = NULL, tsym = "t", a0 = TRUE) {
  e    <- as_multi(e)
  vf   <- if (!is.null(vfun)) vfun else function(nm, lag, dd)
    if (dd == 0) tex_var(nm, lag, sym) else
      paste0("\\left(", diff_str(nm, lag, dd, sym), "\\right)")
  term <- function(v, nm, lag, dd = d) {
    if (!length(v)) return(character(0))
    vapply(seq_along(v), function(i)
      paste0(if (v[i] < 0) "- " else "+ ", num(abs(v[i]), digits), "\\,", vf(nm, lag[i], dd)),
      character(1))
  }
  wrap <- function(tt, k = per_line) {
    if (!length(tt)) return(character(0))
    unname(vapply(split(tt, ceiling(seq_along(tt) / k)), paste, character(1), collapse = " "))
  }
  ct    <- if (!is.null(const)) const else
    if (e$const < 0) paste0("{", num(e$const, digits), "}") else num(e$const, digits)
  trend <- if (e$slope != 0)
    paste0(if (e$slope < 0) "- " else "+ ", num(abs(e$slope), digits), "\\,", tsym)
  ygr   <- c(trend, term(e$y, "Y", e$ylag))
  c(if (length(ct) > 1) c(wrap(ct, 6), wrap(ygr)) else wrap(c(ct, ygr)),
    unlist(lapply(e$inputs, function(z) wrap(term(z$x, z$name, z$xlag)))),
    wrap(c(if (a0) "+ a_{t}", term(e$a, "a", e$alag, 0))))
}
## several equations below each other in one aligned display; eqs: list of list(lhs, lines)
tex_block <- function(eqs, gap = "8pt") {
  body <- vapply(eqs, function(q)
    paste0(q$lhs, " ={}& ", paste(q$lines, collapse = " \\\\\n&{}")), character(1))
  cat("\n$$\n\\begin{aligned}\n", paste(body, collapse = paste0(" \\\\[", gap, "]\n")),
      "\n\\end{aligned}\n$$\n\n", sep = "")
}
## a chain of explicit equations, e.g. the steps back to the original variables. Every step
## is a list(e, sym, const, d, lhs, per_line); lhs is by default Y_t (its d-th difference)
show_chain <- function(..., digits = 4, per_line = 5) {
  tex_block(lapply(list(...), function(q) {
    d <- if (is.null(q$d)) 0 else q$d
    list(lhs   = if (is.null(q$lhs)) diff_str("Y", 0, d, q$sym) else q$lhs,
         lines = eq_rhs(q$e, q$sym, q$const, d, digits,
                        if (is.null(q$per_line)) per_line else q$per_line))
  }))
}
## the explicit equation (sym and const as in eq_rhs())
show_equation <- function(e, sym = NULL, const = NULL, digits = 4, per_line = 5)
  show_chain(list(e = e, sym = sym, const = const), digits = digits, per_line = per_line)
## forecast equations made at the origin t: the one-step-ahead forecast \hat{Y}_t(1), with the
## observed values on the right-hand side and a_{t+1} = 0, and the k-step-ahead recursion
## \hat{Y}_t(k). const1, constk: LaTeX constants that replace the numeric one, e.g. c_{m(t+1)}
show_forecast_eq <- function(e, sym = NULL, const1 = NULL, constk = NULL, digits = 4,
                             per_line = 5) {
  one <- function(nm, lag, dd) tex_at(nm, t_ix(1 - lag), sym)
  kst <- function(nm, lag, dd) tex_hat(nm, if (lag == 0) "k" else paste0("k-", lag), sym)
  tex_block(list(
    list(lhs = tex_hat("Y", "1", sym),
         lines = eq_rhs(e, sym, const1, 0, digits, per_line, vfun = one, tsym = "(t+1)",
                        a0 = FALSE)),
    list(lhs = tex_hat("Y", "k", sym),
         lines = eq_rhs(e, sym, constk, 0, digits, per_line, vfun = kst, tsym = "(t+k)",
                        a0 = FALSE))))
}
## operators of A(B) Y_t = c + sum_j C_j(B) X_{j,t-b_j} + D(B) a_t: as a product of the model
## operators, with the estimated factors, and multiplied out; then the constant c
show_operators <- function(e, sym = NULL, digits = 4) {
  e  <- as_multi(e); sp <- e$spec
  s  <- function(nm)                                    # short symbol of an input
    if (nm %in% names(sym) && !grepl("#", sym[[nm]], fixed = TRUE)) sym[[nm]] else nm
  nt <- function(p) length(p) > 1                       # operator is not the constant 1
  dl <- lapply(sp$inputs, `[[`, "delta")
  jd <- which(vapply(dl, nt, logical(1)))               # inputs with a denominator
  sb <- function(j)
    if (length(sp$inputs) == 1) "" else paste0("_{", s(sp$inputs[[j]]$name), "}")
  dsym <- function(j, arg = "B") paste0("\\delta", sb(j), "(", arg, ")")
  Dsym <- if (sp$d == 1) "(1-B)" else if (sp$d > 1) paste0("(1-B)^{", sp$d, "}")
  Dfac <- if (sp$d == 1) "\\left(1 - B\\right)" else
    if (sp$d > 1) paste0("\\left(1 - B\\right)^{", sp$d, "}")
  vis <- function(x)                                    # rough printed length
    nchar(gsub("[{}^_ ]", "", gsub("\\\\[a-zA-Z]+", "x", gsub("\\\\(left|right|,)", "", x))))
  eqrow <- function(lhs, parts, tag = "") {             # one line if short, else one per part
    parts <- unique(parts[parts != ""])
    if (sum(vis(parts)) <= 60) return(paste0(lhs, " &= ", paste(parts, collapse = " = "), tag))
    paste0(lhs, " &= ", parts[1], tag, paste0(" \\\\\n&= ", parts[-1], collapse = ""))
  }
  row <- function(lhs, syms, facs, prod, on) {
    facs <- facs[facs != ""]
    eqrow(lhs, c(if (length(syms)) paste(syms, collapse = "\\,") else "1",
                 if (length(facs) > 1 || (length(facs) == 1 && !startsWith(facs[1], "\\left(")))
                   paste(facs, collapse = ""),
                 prod),
          paste0(" && \\text{acting on } ", on))
  }
  rows <- c(
    row("A(B)", c(vapply(jd, dsym, ""), if (nt(sp$phi)) "\\phi(B)", Dsym),
        c(vapply(jd, function(j) fac_str(dl[[j]], digits), ""), fac_str(sp$phi, digits), Dfac),
        poly_str(c(1, -e$y), digits), tex_var("Y", 0, sym)),
    vapply(seq_along(sp$inputs), function(j) {
      oth <- setdiff(jd, j)
      row(paste0("C", sb(j), "(B)"),
          c(paste0("\\omega", sb(j), "(B)"), vapply(oth, dsym, ""),
            if (nt(sp$phi)) "\\phi(B)", Dsym),
          c(fac_str(sp$inputs[[j]]$w, digits),
            vapply(oth, function(i) fac_str(dl[[i]], digits), ""), fac_str(sp$phi, digits), Dfac),
          poly_str(e$inputs[[j]]$x, digits),
          tex_var(e$inputs[[j]]$name, e$inputs[[j]]$xlag[1], sym))
    }, character(1)),
    row("D(B)", c(vapply(jd, dsym, ""), if (nt(sp$theta)) "\\theta(B)"),
        c(vapply(jd, function(j) fac_str(dl[[j]], digits), ""), fac_str(sp$theta, digits)),
        poly_str(c(1, e$a), digits), "a_{t}"))
  ## the constant: c = mu delta(1) phi(1); with a trend beta t, c = mu A(1) - beta sum_k k A_k
  A0 <- polymul(Reduce(polymul, dl, 1), sp$phi)
  f  <- c(vapply(jd, function(j) sum(dl[[j]]), 0), if (nt(sp$phi)) sum(sp$phi))
  fs <- c(vapply(jd, dsym, "", arg = "1"), if (nt(sp$phi)) "\\phi(1)")
  cst <- if (sp$beta == 0) {
    eqrow("c", c(paste0("\\hat\\mu", if (length(fs)) paste0("\\,", paste(fs, collapse = "\\,"))),
                 paste0(paste(num(c(sp$mu, f), 6), collapse = " \\times "), " = ",
                        num(e$const, 6))))
  } else {                                              # c and the coefficient beta* of t
    c(eqrow("c", c("\\hat\\mu\\,A(1) - \\hat\\beta\\textstyle\\sum_k k\\,A_k",
                   paste0(num(sp$mu, 6), " \\times ", num(sum(A0), 6), " - (", num(sp$beta, 6),
                          ") \\times (", num(sum((seq_along(A0) - 1) * A0), 6), ") = ",
                          num(e$const, 6)))),
      eqrow("\\beta^*", c("\\hat\\beta\\,A(1)", paste0(num(sp$beta, 6), " \\times ",
                                                        num(sum(A0), 6), " = ", num(e$slope, 6)))))
  }
  cat("\n$$\n\\begin{aligned}\n", paste(c(rows, cst), collapse = " \\\\[4pt]\n"),
      "\n\\end{aligned}\n$$\n\n", sep = "")
}
## the month-specific constant c_{m(t)} = sum_m c_m D_{m,t}, written out with the seasonal
## dummies D_{m,t} (1 if t falls in month m of the cycle, 0 otherwise)
dummy_terms <- function(cm, digits = 4) {
  v <- num(abs(cm), digits)
  paste0(c(if (cm[1] < 0) paste0("{-", v[1], "}") else v[1],
           paste0(ifelse(cm[-1] < 0, "- ", "+ "), v[-1])), "\\,D_{", seq_along(cm), ",t}")
}
## monthly constants c[m] of the explicit equation in the original (unadjusted) variables:
## with y*_t = Y_t - mY[m(t)] and x*_t = X_t - mX[m(t)], the explicit equation of y* becomes
## Y_t = c[m(t)] + (the same terms in the lags of Y, X and a). mX: named list of monthly means
monthly_const <- function(e, mY, mX) {
  e  <- as_multi(e, names(mX)[1])
  mo <- function(m, k) (m - k - 1) %% 12 + 1          # month of t - k when t is in month m
  unname(sapply(1:12, function(m) e$const + mY[m] - sum(e$y * mY[mo(m, e$ylag)]) -
           sum(sapply(e$inputs, function(z) sum(z$x * mX[[z$name]][mo(m, z$xlag)])))))
}

## ---- k-step-ahead forecasting ----

## forecasts of Y[o+1], ..., Y[o+H] made at the origin o (future shocks a = 0).
## X: named list of the observed inputs (names as in the explicit equation).
## xfc = NULL: the observed future inputs are used (ex post);
## xfc(o, H) = list of input forecasts for o+1, ..., o+H: only the inputs up to o are
## used (ex ante)
tf_forecast <- function(e, Y, X, a, o, H, xfc = NULL) {
  e <- as_multi(e, names(X)[1])
  if (!is.null(xfc)) {
    xf <- xfc(o, H)
    X  <- imap(X, function(z, nm) c(as.numeric(z)[1:o], xf[[nm]]))
  }
  Y <- c(as.numeric(Y)[1:o], rep(NA, H)); a <- c(as.numeric(a)[1:o], rep(0, H))
  for (t in o + 1:H) {
    xpart <- sum(sapply(e$inputs, function(z) sum(z$x * X[[z$name]][t - z$xlag])))
    Y[t]  <- e$const + e$slope * t + sum(e$y * Y[t - e$ylag]) + xpart +
             sum(e$a * a[t - e$alag])
  }
  Y[o + 1:H]
}

## first o values of a series, keeping its time-series attributes
ts_head <- function(z, o) ts(as.numeric(z)[1:o], start = start(z), frequency = frequency(z))

## point forecasts and standard errors (from the 95% interval) of a forecast object
fc_out <- function(f) {
  up <- as.matrix(f$upper)[, which(f$level == 95)]
  list(mean = as.numeric(f$mean), se = (as.numeric(up) - as.numeric(f$mean)) / qnorm(0.975))
}

## automatic forecasting methods for an input. Each one is estimated on a series z and
## returns a function(zh, h) that forecasts h steps ahead from the history zh. ARIMA, ETS
## and ARIMA + Fourier keep the parameters estimated on z; the other methods have (almost)
## no parameters to estimate and are simply re-run on zh.
auto_methods <- list(
  ARIMA = function(z) {
    m <- auto.arima(z, seasonal = FALSE)
    function(zh, h) fc_out(forecast(Arima(zh, model = m), h = h)) },
  ETS   = function(z) {
    m <- ets(z)
    function(zh, h) fc_out(forecast(ets(zh, model = m, use.initial.values = TRUE), h = h)) },
  Theta = function(z) function(zh, h) fc_out(thetaf(zh, h = h)),
  Naive = function(z) function(zh, h) fc_out(naive(zh, h = h)),    # benchmark: last value
  Mean  = function(z) function(zh, h) fc_out(meanf(zh, h = h)))    # benchmark: mean
## candidates for a seasonal series (the weekly inputs of Case 4)
auto_methods_seasonal <- list(
  `ARIMA + Fourier` = function(z) {      # dynamic harmonic regression, K chosen by AICc
    fits <- lapply(1:3, function(K) auto.arima(z, xreg = fourier(z, K), seasonal = FALSE))
    K <- which.min(sapply(fits, `[[`, "aicc")); m <- fits[[K]]
    function(zh, h) fc_out(forecast(Arima(zh, model = m, xreg = fourier(zh, K)),
                                    xreg = fourier(zh, K, h = h))) },
  `STL + ETS`      = function(z) function(zh, h) fc_out(stlf(zh, h = h)),
  `Seasonal naive` = function(z) function(zh, h) fc_out(snaive(zh, h = h)))

## choice of the forecasting method for an input series z by time-series cross-validation:
## from n_origins origins spread over the last 20% of z, every method is re-estimated on the
## data up to the origin and forecasts h = 1, ..., H; the method with the smallest RMSE is
## re-estimated on the whole series z
select_auto <- function(z, H, methods = auto_methods, v = round(0.2 * length(z)),
                        n_origins = 10) {
  m <- length(z)
  origins <- unique(round(seq(m - v, m - H, length.out = n_origins)))
  score <- imap_dfr(methods, function(make, nm) {
    e <- unlist(lapply(origins, function(o) {
      zo <- ts_head(z, o)
      as.numeric(z)[o + 1:H] - make(zo)(zo, H)$mean
    }))
    tibble(method = nm, `CV RMSE` = sqrt(mean(e^2)))
  })
  best <- score$method[which.min(score$`CV RMSE`)]
  f <- methods[[best]](z)
  list(table = mutate(score, `CV RMSE` = round(`CV RMSE`, 4), selected = method == best),
       best = best, model = spec_label(f, best), f = f)
}

## short description of the model behind a forecasting function from auto_methods
spec_label <- function(f, method) {
  m <- environment(f)$m
  if (is.null(m)) return(method)
  if (inherits(m, "ets")) return(m$method)
  arma <- paste0("ARIMA(", paste(arimaorder(m)[1:3], collapse = ","), ")")
  if (!is.null(environment(f)$K))
    return(paste0(arma, " errors + ", environment(f)$K, " Fourier pairs"))
  paste0(arma, if ("intercept" %in% names(coef(m))) " with mean" else "",
         if ("drift" %in% names(coef(m))) " with drift" else "")
}

## rolling-origin forecasts from every origin o at the horizons h = 1, ..., H that fall in the
## data (o + h <= n): transfer function with the inputs forecast (ex ante) and with the
## observed inputs (ex post) vs univariate ARIMA. Returns the RMSE for every horizon, using
## every available origin; the forecasts themselves are kept in attr(, "fc")
rolling_eval <- function(eq, Y, X, a, uni, origins, H, xfc) {
  fc <- map_dfr(origins, function(o) {
    Hc <- min(H, length(Y) - o)
    tibble(origin = o, h = 1:Hc, t = o + 1:Hc, Actual = as.numeric(Y)[o + 1:Hc],
           `TF ex ante` = tf_forecast(eq, Y, X, a, o, Hc, xfc),
           `TF ex post` = tf_forecast(eq, Y, X, a, o, Hc),
           `Univariate ARIMA` = as.numeric(forecast(Arima(ts_head(Y, o), model = uni),
                                                    h = Hc)$mean))
  })
  out <- fc %>%
    group_by(h) %>%
    summarise(origins = n(), `TF ex ante` = sqrt(mean((Actual - `TF ex ante`)^2)),
              `TF ex post` = sqrt(mean((Actual - `TF ex post`)^2)),
              `univariate ARIMA` = sqrt(mean((Actual - `Univariate ARIMA`)^2)),
              .groups = "drop") %>%
    mutate(across(where(is.double), ~ round(.x, 4)))
  attr(out, "fc") <- fc
  out
}
## the one-step-ahead forecasts of the test period (the forecasts at h = 1 of rolling_eval()),
## converted to the original scale with orig(value, t)
onestep_of <- function(ev, orig = function(v, t) v)
  filter(attr(ev, "fc"), h == 1) %>%
    mutate(across(c(Actual, `TF ex ante`, `TF ex post`, `Univariate ARIMA`), ~ orig(.x, t)))
fc_cols <- function(d) as.list(d[c("TF ex ante", "TF ex post", "Univariate ARIMA")])
acc_of  <- function(d) do.call(accuracy_tbl, c(list(d$Actual), fc_cols(d)))

## forecasts of the test period with the actual values
plot_fc <- function(t, actual, fc, title, ylab = NULL)
  do.call(plot_holdout, c(list(t, actual), fc)) + labs(title = title, y = ylab)
## RMSE of the rolling-origin forecasts by horizon
plot_rmse <- function(e, title) {
  e %>% select(h, `TF ex ante`, `TF ex post`, `Univariate ARIMA` = `univariate ARIMA`) %>%
    pivot_longer(-h, names_to = "series") %>%
    ggplot(aes(h, value, colour = series, linetype = series)) + geom_line() +
    geom_point(size = 1.5) +
    scale_colour_manual(values = c(`TF ex ante` = "steelblue", `TF ex post` = "forestgreen",
                                   `Univariate ARIMA` = "darkorange")) +
    scale_linetype_manual(values = c(`TF ex ante` = "solid", `TF ex post` = "dashed",
                                     `Univariate ARIMA` = "solid")) +
    scale_x_continuous(breaks = e$h) +
    labs(title = title, x = "horizon h", y = "RMSE", colour = NULL, linetype = NULL) +
    theme_bw() + theme(legend.position = "bottom")
}

## the last n_show observations and the forecasts beyond the end of the data
plot_future <- function(Y, fc, n_show = 40) {
  n <- length(Y); old <- (n - n_show + 1):n
  bind_rows(tibble(t = old, value = as.numeric(Y)[old], series = "Observed"),
            tibble(t = n + 0:length(fc), value = c(as.numeric(Y)[n], fc),
                   series = "Forecast")) %>%
    ggplot(aes(t, value, colour = series)) + geom_line() + geom_point(size = 1) +
    scale_colour_manual(values = c(Observed = "black", Forecast = "steelblue")) +
    labs(x = "t", y = NULL, colour = NULL) + theme_bw() + theme(legend.position = "bottom")
}

## numbers quoted in the text: fixed decimals, RMSE at horizon h, and the range (in %) by
## which column a of a rolling_eval() table is smaller than column b
fmt_num  <- function(x, digits) formatC(x, format = "f", digits = digits)
rmse_at  <- function(e, col, h, digits = 3) fmt_num(e[[col]][e$h == h], digits)
pct_smaller <- function(e, a, b, rows = seq_len(nrow(e)))
  paste(unique(round(range(100 * (1 - e[[a]][rows] / e[[b]][rows])))), collapse = " to ")

The helpers below draw the summary figure at the end of the notebook: one row per case, with the original series, the series used in the model, the impulse response weights, the residual ACF and the explicit fitted equation.

## ---- summary figure ----
fig_theme <- theme_bw(base_size = 9) +
  theme(plot.title = element_text(size = 10),
        strip.text = element_text(size = 7, margin = margin(1, 0, 1, 0)))
## several series as stacked mini-plots
fig_series <- function(lst, title, colour) {
  d <- imap_dfr(lst, ~ tibble(t = seq_along(.x), value = as.numeric(.x), series = .y)) %>%
    mutate(series = factor(series, levels = names(lst)))
  ggplot(d, aes(t, value)) + geom_line(colour = colour, linewidth = 0.25) +
    facet_wrap(~ series, ncol = 1, scales = "free_y") +
    labs(title = title, x = "t", y = NULL) + fig_theme
}
## a text panel in place of a plot
fig_note <- function(title, note)
  ggplot() + annotate("text", x = 0.5, y = 0.5, label = note, size = 3.6) +
    xlim(0, 1) + ylim(0, 1) + labs(title = title) + theme_void() +
    theme(plot.title = element_text(size = 10))
## identified (bars) and fitted (line) impulse response weights, one facet per input
fig_irf <- function(est, fit, title)
  ggplot(est, aes(k, nu)) + geom_col(fill = "grey65", width = 0.6) +
    geom_hline(yintercept = 0) +
    geom_line(data = fit, colour = "steelblue") +
    geom_point(data = fit, colour = "steelblue", size = 0.8) +
    facet_wrap(~ input, ncol = 1, scales = "free_y") +
    labs(title = title, x = "lag k", y = NULL) + fig_theme
## residual ACF with the p-value of the Ljung-Box test at lag 12
fig_acf <- function(a, fitdf, title) {
  p <- Box.test(na.omit(a), lag = 12, type = "Ljung-Box", fitdf = fitdf)$p.value
  ggAcf(na.omit(a), lag.max = 24) + fig_theme +
    labs(title = paste0(title, " (Q(12) p = ", sprintf("%.3f", p), ")"), y = NULL)
}
## the explicit equation as plotmath lines, grouped by variable as in show_equation(),
## and the noise variance; sym: templates with # for the time index, e.g. c(Y = "y#^'*'")
eq_plotmath <- function(e, sigma2, sym = NULL, digits = 4, per_line = 5) {
  e  <- as_multi(e)
  sy <- function(nm, ix)
    sub("#", ix, if (nm %in% names(sym)) sym[[nm]] else paste0(nm, "#"), fixed = TRUE)
  term <- function(v, nm, lag) {
    if (!length(v)) return(character(0))
    paste0(ifelse(v < 0, " - ", " + "), "\"", num(abs(v), digits), "\"*",
           vapply(ifelse(lag == 0, "[t]", paste0("[t-", lag, "]")), sy, character(1),
                  nm = nm))
  }
  trend  <- if (e$slope != 0)
    paste0(ifelse(e$slope < 0, " - ", " + "), "\"", num(abs(e$slope), digits), "\"*t")
  groups <- c(list(c(trend, term(e$y, "Y", e$ylag))),
              lapply(e$inputs, function(z) term(z$x, z$name, z$xlag)),
              list(c(" + a[t]", term(e$a, "a", e$alag))))
  lines <- unlist(lapply(groups, function(g)
    vapply(split(g, ceiling(seq_along(g) / per_line)), paste, character(1), collapse = "")))
  const <- paste0(if (e$const < 0) "-" else "", "\"", num(abs(e$const), digits), "\"")
  c(paste0(sy("Y", "[t]"), " == ", const, lines[1]), paste0("phantom()", lines[-1]),
    paste0("a[t] %~% N(0, \"", num(sigma2, digits), "\")"))
}
## the model and its explicit equation as a text panel (n_max lines)
fig_eq <- function(title, lines, n_max = length(lines))
  ggplot() +
    annotate("text", x = 0, y = n_max - seq_along(lines) + 1, label = lines, parse = TRUE,
             hjust = 0, size = 4) +
    scale_x_continuous(limits = c(0, 1), expand = c(0, 0)) +
    scale_y_continuous(limits = c(0.4, n_max + 0.6), expand = c(0, 0)) +
    labs(title = title) + theme_void() + theme(plot.title = element_text(size = 10))

Transfer Function Modelling

Cases 1–3 follow the eleven steps of the CCF approach; cases 4–5 follow the seven steps of the LTF approach. Every case ends with a comparison against a univariate ARIMA model of the output on a test period, with one-step-ahead and k-step-ahead forecasts, which shows what the input adds in terms of forecasting.

Case 1 — Series M: Sales and a Leading Indicator

Problem Description

Box and Jenkins’ Series M contains 150 consecutive values of the sales of a company (Y_t) together with a leading indicator (X_t). The business problem is forecasting. Management wants to know whether movements in the indicator announce movements in sales, how many periods in advance they do so, and how the effect builds up over time. If the indicator leads sales by b periods, sales can be forecast up to b periods ahead from values of the indicator that are already observed, which is exactly what makes a leading indicator valuable. The causal direction is one-way (the indicator is not driven by sales), so a single-input transfer function is the natural model. The data are BJsales and BJsales.lead in base R; this case is also the worked example of the lecture.

Data

Y <- BJsales; X <- BJsales.lead
tibble(t = seq_along(Y), `Sales Y` = as.numeric(Y), `Leading indicator X` = as.numeric(X)) %>%
  pivot_longer(-t) %>%
  ggplot(aes(t, value)) +
    geom_line(colour = "steelblue") +
    facet_wrap(~ name, ncol = 1, scales = "free_y") +
    labs(x = "t", y = NULL) +
    theme_bw()

Both series wander without a fixed mean: sales rise from about 200 to 260 with long upswings and downswings, and the indicator moves between about 10 and 14. Visually, the turning points of the indicator appear a few periods before those of sales.

Step 1 — Stationarity

ADF and PP have H_0: unit root; KPSS has H_0: level-stationary. The tests are applied to the levels and to the first differences of both series.

stat_tests(Y = Y, X = X, `diff(Y)` = diff(Y), `diff(X)` = diff(X))

In levels, ADF (p-values 0.530 and 0.692) and PP (0.888 and 0.293) do not reject a unit root, and KPSS (statistics 2.62 and 2.56, p-value \le 0.01) rejects stationarity: both series are non-stationary. After first differencing, PP rejects the unit root for both series (p-value \le 0.01) and KPSS no longer rejects stationarity (p-value \ge 0.10). The ADF test on \nabla Y_t is borderline (p-value = 0.066), but the other two tests agree, so d = 1 is used for both series. The same differencing is applied to the input and the output, y_t = \nabla Y_t and x_t = \nabla X_t.

y <- as.numeric(diff(Y)); x <- as.numeric(diff(X)); n <- length(y)

Step 2 — Input Model

The ARIMA model is identified for the input only; it defines the prewhitening filter.

(ggAcf(x, lag.max = 24) + ggtitle("ACF diff(X)")) | (ggPacf(x, lag.max = 24) + ggtitle("PACF diff(X)"))

The ACF of x_t has a single significant spike at lag 1 (-0.447) and cuts off afterwards, while the PACF decays (-0.447, -0.143, -0.116, \dots). This is the pattern of an MA(1), so the input model is ARIMA(0,1,1) for X_t.

fit_xm <- Arima(X, order = c(0, 1, 1))
fit_xm

Series: X 
ARIMA(0,1,1) 

Coefficients:
          ma1
      -0.4475
s.e.   0.0635

sigma^2 = 0.08036:  log likelihood = -23.2
AIC=50.41   AICc=50.49   BIC=56.41

coeftest(fit_xm)

z test of coefficients:

     Estimate Std. Error z value  Pr(>|z|)    
ma1 -0.447453   0.063511 -7.0453 1.851e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

\hat\theta_1 = -0.4475 (s.e. 0.0635, z = -7.05), so x_t = (1 - 0.4475B)\,\alpha_t with \hat\sigma_\alpha^2 = 0.0804.

Step 3 — Prewhitening

Inverting the input model gives the prewhitening filter \pi(B) = 1/(1 - 0.4475B), i.e. \alpha_t = x_t + 0.4475\,\alpha_{t-1}. The same filter is applied to y_t.

th <- coef(fit_xm)["ma1"]
alpha_m <- pw_filter(x, theta = th)     # = stats::filter(x, -th, method = "recursive")
beta_m  <- pw_filter(y, theta = th)
Box.test(alpha_m, lag = 12, type = "Ljung-Box", fitdf = 1)

    Box-Ljung test

data:  alpha_m
X-squared = 12.321, df = 11, p-value = 0.34

Q(12) = 12.32 with 11 degrees of freedom, p-value = 0.340: the prewhitened input \alpha_t is white noise, as required. (Without the degrees-of-freedom correction the p-value is 0.42; the conclusion is the same.)

Step 4 — Cross-Correlation and Impulse Response Weights

tab_m <- ccf_table(beta_m, alpha_m, lag.max = 12)
plot_ccf(tab_m, "CCF of prewhitened series: Series M")

tab_m %>% filter(k >= 0) %>% mutate(across(c(r, nu), ~ round(.x, 3)))

With n = 149 the bounds are \pm 2/\sqrt{149} = \pm 0.164. The cross-correlations at lags 0, 1 and 2 are small (0.062, 0.079, 0.017). At lag 3 there is a large spike, r_{\alpha\beta}(3) = 0.675 (\hat\nu_3 = 4.70), followed by values that die down gradually: 0.451, 0.339, 0.256, 0.267, 0.197 and 0.172 (\hat\nu_4, \dots, \hat\nu_9 = 3.14, 2.36, 1.78, 1.86, 1.37, 1.20). At negative lags the largest absolute cross-correlation is 0.094, inside the bounds, so there is no feedback from sales to the indicator.

Step 5 — Orders (b, r, s)

• \hat\nu_0, \hat\nu_1, \hat\nu_2 are not significant, so b = 3.
• From lag 3 the weights decay geometrically (\hat\nu_4/\hat\nu_3 = 0.67, \hat\nu_5/\hat\nu_4 = 0.75), so r = 1.
• The decay starts at \hat\nu_3 itself, with no unpatterned weights before it, so s = 0.

The tentative transfer function is \nu(B)x_t = \dfrac{\omega_0}{1 - \delta_1 B}\,x_{t-3}.

Step 6 — Preliminary Estimates

prelim_tf(tab_m$nu[tab_m$k >= 0], b = 3, r = 1, s = 0)

  omega0   delta1 
4.701306 0.667806 

\hat\omega_0 = \hat\nu_3 = 4.70 and \hat\delta_1 = \hat\nu_4/\hat\nu_3 = 0.67, which imply a steady-state gain of 4.70/(1 - 0.67) \approx 14.2. These rough values serve as a first look at the transfer function and give a first noise series.

Step 7 — Noise Model

The preliminary noise \hat\eta_t = y_t - \hat\mu - \frac{\hat\omega_0}{1-\hat\delta_1B}x_{t-3} is obtained from a quick ML fit of the transfer function with a proxy AR(1) noise.

b <- 3
pre_m <- arimax(y[(b+1):n], order = c(1, 0, 0), method = "ML",
                xtransf = data.frame(lead = x[1:(n-b)]), transfer = list(c(1, 0)))
pre_m

Call:
arimax(x = y[(b + 1):n], order = c(1, 0, 0), method = "ML", xtransf = data.frame(lead = x[1:(n - 
    b)]), transfer = list(c(1, 0)))

Coefficients:
          ar1  intercept  lead-AR1  lead-MA0
      -0.4233     0.0299    0.7241    4.7281
s.e.   0.0752     0.0147    0.0052    0.0620

sigma^2 estimated as 0.05138:  log likelihood = 9.44,  aic = -10.88

cp <- coef(pre_m)
eta_m <- y[(b+1):n] - cp["intercept"] -
  tf_component(x[1:(n-b)], cp["lead-MA0"], c(1, -cp["lead-AR1"]))
(ggAcf(eta_m, lag.max = 16) + ggtitle("ACF of preliminary noise")) |
  (ggPacf(eta_m, lag.max = 16) + ggtitle("PACF of preliminary noise"))

c(AR1 = Arima(eta_m, order = c(1, 0, 0))$aic, MA1 = Arima(eta_m, order = c(0, 0, 1))$aic)

      AR1       MA1 
-12.87755 -23.92014 

The proxy fit gives \hat\omega_0 = 4.73 and \hat\delta_1 = 0.724, close to the preliminary values. The ACF of \hat\eta_t has one significant spike at lag 1 (-0.42) and then stays inside the bands, while the PACF decays (-0.42, -0.19, -0.13, -0.21, \dots). This points to an MA(1) noise, which also has a smaller AIC than an AR(1) (-23.92 versus -12.88).

Step 8 — Joint Estimation

All parameters, (b,r,s) = (3,1,0) with MA(1) noise, are estimated jointly by maximum likelihood.

fit_m <- arimax(y[(b+1):n], order = c(0, 0, 1), method = "ML",
                xtransf = data.frame(lead = x[1:(n-b)]), transfer = list(c(1, 0)))
fit_m

Call:
arimax(x = y[(b + 1):n], order = c(0, 0, 1), method = "ML", xtransf = data.frame(lead = x[1:(n - 
    b)]), transfer = list(c(1, 0)))

Coefficients:
          ma1  intercept  lead-AR1  lead-MA0
      -0.5874     0.0305    0.7264    4.6942
s.e.   0.0712     0.0085    0.0038    0.0520

sigma^2 estimated as 0.04741:  log likelihood = 15.19,  aic = -22.38

coeftest(fit_m)

z test of coefficients:

            Estimate Std. Error  z value  Pr(>|z|)    
ma1       -0.5873571  0.0712386  -8.2449 < 2.2e-16 ***
intercept  0.0305240  0.0084742   3.6020 0.0003158 ***
lead-AR1   0.7264007  0.0037909 191.6159 < 2.2e-16 ***
lead-MA0   4.6942119  0.0519788  90.3101 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

All four coefficients are significant; the smallest test statistic is z = 3.60 for the constant. The residual variance is \hat\sigma_a^2 = 0.0474 and AIC = -22.38. The fitted model is

(1-B)\,Y_t = 0.0305 + \frac{4.6942}{1 - 0.7264B}\,(1-B)\,X_{t-3} + (1 - 0.5874B)\,a_t, \qquad a_t \sim \mathcal{N}(0,\ 0.0474).

The steady-state gain is 4.6942/(1-0.7264) = 17.16: a permanent one-unit rise of the indicator raises sales by about 17 units in the long run, starting three periods later.

Step 9 — Diagnostic Checking

Two checks are needed. (a) The residuals \hat a_t must be white noise (noise model); (b) they must be uncorrelated with the prewhitened input \alpha_t (transfer function), tested with Q_0 = m(m+2)\sum_{k=0}^{K} r^2_{\alpha\hat a}(k)/(m-k) \sim \chi^2_{K+1-M}, M = r+s+1.

a_m <- residuals(fit_m)
Box.test(a_m, lag = 12, type = "Ljung-Box", fitdf = 1)

    Box-Ljung test

data:  a_m
X-squared = 10.249, df = 11, p-value = 0.5081

Box.test(a_m, lag = 24, type = "Ljung-Box", fitdf = 1)

    Box-Ljung test

data:  a_m
X-squared = 20.377, df = 23, p-value = 0.6191

q0_m <- bind_rows(Q0_test(a_m, alpha_m[(b+1):n], K = 12, M = 2),
                  Q0_test(a_m, alpha_m[(b+1):n], K = 24, M = 2))
q0_m

jarque.bera.test(a_m)

    Jarque Bera Test

data:  a_m
X-squared = 1.2219, df = 2, p-value = 0.5428

(ggAcf(a_m, lag.max = 24) + ggtitle("ACF of residuals")) |
  plot_ccf(ccf_table(a_m, alpha_m[(b+1):n], lag.max = 16) %>% filter(k >= 0),
           "CCF of alpha and residuals (k >= 0)")

• Noise model: Q(12) = 10.25 (df = 11, p-value = 0.508) and Q(24) = 20.38 (df = 23, p-value = 0.619), so the residuals are white noise.
• Transfer function: Q_0(12) = 12.10 (df = 11, p-value = 0.356) and Q_0(24) = 20.2 (p-value = 0.631), so the residuals are uncorrelated with the prewhitened input.
• Normality: Jarque–Bera X^2 = 1.22, p-value = 0.543, so normality is not rejected.

Step 10 — Explicit Equation

The model of Step 8 was estimated for the differences y_t = \nabla Y_t = Y_t - Y_{t-1} and x_t = \nabla X_t = X_t - X_{t-1}, with \nabla = 1 - B. Multiplying it by \delta(B)\phi(B) = 1 - 0.7264B removes the denominator and gives

A(B)\,Y_t = c + C(B)\,X_{t-3} + D(B)\,a_t, \qquad A(B) = \delta(B)\phi(B)(1-B),\quad C(B) = \omega(B)\phi(B)(1-B),\quad D(B) = \delta(B)\theta(B),

with c = \mu\,\delta(1)\phi(1); the noise is MA(1), so \phi(B) = 1. tf_explicit() computes these products, and show_operators() writes every operator as a product of the model operators, with the estimated factors, and multiplied out:

cf_m <- coef(fit_m)
eq_m <- tf_explicit(w = cf_m["lead-MA0"], delta = c(1, -cf_m["lead-AR1"]), b = 3,
                    theta = c(1, cf_m["ma1"]), d = 1, mu = cf_m["intercept"])
show_operators(eq_m)

\begin{aligned} A(B) &= \delta(B)\,(1-B) = \left(1 - 0.7264B\right)\left(1 - B\right) = 1 - 1.7264B + 0.7264B^{2} && \text{acting on } Y_{t} \\[4pt] C(B) &= \omega(B)\,(1-B) = 4.6942\left(1 - B\right) = 4.6942 - 4.6942B && \text{acting on } X_{t-3} \\[4pt] D(B) &= \delta(B)\,\theta(B) = \left(1 - 0.7264B\right)\left(1 - 0.5874B\right) = 1 - 1.3138B + 0.4267B^{2} && \text{acting on } a_{t} \\[4pt] c &= \hat\mu\,\delta(1) = 0.0305240 \times 0.273599 = 0.00835133 \end{aligned}

Applying B^k Z_t = Z_{t-k} and moving the lagged values of Y to the right-hand side gives the explicit equation. Because the model was estimated for the differences, the display below goes back to the original variables in three steps:

1. the explicit equation of the differenced series \nabla Y_t and \nabla X_t, i.e. of the series to which the model was fitted (tf_explicit() with d = 0);
2. the same equation with every difference written out, \nabla Y_{t-k} = Y_{t-k} - Y_{t-k-1} and \nabla X_{t-k} = X_{t-k} - X_{t-k-1};
3. after collecting the terms, the explicit equation of sales Y_t and the leading indicator X_t themselves (tf_explicit() with d = 1).

eq_m0 <- tf_explicit(w = cf_m["lead-MA0"], delta = c(1, -cf_m["lead-AR1"]), b = 3,
                     theta = c(1, cf_m["ma1"]), d = 0, mu = cf_m["intercept"])
show_chain(list(e = eq_m0, sym = c(Y = "\\nabla Y", X = "\\nabla X")),   # 1. differences
           list(e = eq_m0, d = 1),                                       # 2. written out
           list(e = eq_m))                                               # 3. original variables

\begin{aligned} \nabla Y_{t} ={}& 0.008351 + 0.7264\,\nabla Y_{t-1} \\ &{}+ 4.6942\,\nabla X_{t-3} \\ &{}+ a_{t} - 1.3138\,a_{t-1} + 0.4267\,a_{t-2} \\[8pt] Y_{t} - Y_{t-1} ={}& 0.008351 + 0.7264\,\left(Y_{t-1} - Y_{t-2}\right) \\ &{}+ 4.6942\,\left(X_{t-3} - X_{t-4}\right) \\ &{}+ a_{t} - 1.3138\,a_{t-1} + 0.4267\,a_{t-2} \\[8pt] Y_{t} ={}& 0.008351 + 1.7264\,Y_{t-1} - 0.7264\,Y_{t-2} \\ &{}+ 4.6942\,X_{t-3} - 4.6942\,X_{t-4} \\ &{}+ a_{t} - 1.3138\,a_{t-1} + 0.4267\,a_{t-2} \end{aligned}

The last equation is the explicit equation in the original variables; it is the equation used for forecasting in Step 11. It can be checked by hand:

• Y_{t-1} and Y_{t-2} have the coefficients 1 + \hat\delta_1 = 1.7264 and -\hat\delta_1 = -0.7264. They sum to one, which is the unit root of the differencing: there is no fixed level to which sales return;
• X_{t-3} and X_{t-4} have the same coefficient \hat\omega_0 = 4.6942 with opposite signs, so sales respond to the change of the indicator three periods ago;
• the MA coefficients satisfy 0.7264 + 0.5874 = 1.3138 and 0.7264 \times 0.5874 = 0.4267;
• the constant is c = \hat\mu\,(1 - \hat\delta_1), derived in the first display.

Step 11 — Forecasting

The forecasts are computed recursively from the explicit equation of Step 10: observed values where they are available, forecasts of Y otherwise, and future shocks a_t = 0. Because b = 3, the forecasts up to h = 3 need only observed values of the indicator; from h = 4 on, the indicator itself has to be forecast.

a. Forecast equations

For the origin t, the one-step-ahead forecast and the k-step-ahead forecast of sales are

show_forecast_eq(eq_m)

\begin{aligned} \hat{Y}_t(1) ={}& 0.008351 + 1.7264\,Y_{t} - 0.7264\,Y_{t-1} \\ &{}+ 4.6942\,X_{t-2} - 4.6942\,X_{t-3} \\ &{}- 1.3138\,a_{t} + 0.4267\,a_{t-1} \\[8pt] \hat{Y}_t(k) ={}& 0.008351 + 1.7264\,\hat{Y}_t(k-1) - 0.7264\,\hat{Y}_t(k-2) \\ &{}+ 4.6942\,\hat{X}_t(k-3) - 4.6942\,\hat{X}_t(k-4) \\ &{}- 1.3138\,\hat{a}_t(k-1) + 0.4267\,\hat{a}_t(k-2) \end{aligned}

with \hat Y_t(m) = Y_{t+m}, \hat X_t(m) = X_{t+m} and \hat a_t(m) = a_{t+m} (the residuals) for m \le 0, and \hat a_t(m) = 0 for m > 0. For m > 0, \hat X_t(m) is the forecast of the indicator (ex ante) or its observed value (ex post). In the one-step-ahead equation every value on the right-hand side is observed at the origin, including X_{t-2} and X_{t-3}. In the k-step-ahead equation, \hat X_t(k-3) is an observed value as long as k \le 3 and a forecast of the indicator for k \ge 4.

b. Training and test periods

To measure the accuracy, the transfer function is re-estimated on the first 135 observations, with the same orders, and the last 15 observations form the test period. The benchmark is the model chosen by auto.arima() for sales alone. Both models keep the parameters estimated on the training data. The forecasting method for the indicator is selected by select_auto() on the training data only. From every origin t = 135, \dots, 149, forecasts are made up to K = 6 steps ahead, as far as the test period reaches.

K_m  <- 6                                                 # forecast horizons h = 1, ..., K
T0_m <- length(Y) - 15                                    # training period t = 1, ..., 135
ytr <- as.numeric(diff(Y[1:T0_m])); xtr <- as.numeric(diff(X[1:T0_m])); ntr <- length(ytr)
fit_m_tr <- arimax(ytr[(b+1):ntr], order = c(0, 0, 1), method = "ML",
                   xtransf = data.frame(lead = xtr[1:(ntr-b)]), transfer = list(c(1, 0)))
ct <- coef(fit_m_tr); round(ct, 4)

      ma1 intercept  lead-AR1  lead-MA0 
  -0.5783    0.0282    0.7256    4.7088 

eq_m_tr <- tf_explicit(w = ct["lead-MA0"], delta = c(1, -ct["lead-AR1"]), b = 3,
                       theta = c(1, ct["ma1"]), d = 1, mu = ct["intercept"])
test_m <- (T0_m + 1):length(Y)
a_tr <- c(rep(0, T0_m - length(residuals(fit_m_tr))), residuals(fit_m_tr), rep(0, length(test_m)))
f1_m <- onestep_explicit(eq_m_tr, Y, X, a_tr, test_m)    # residuals of the test period
a_tr[test_m] <- as.numeric(Y)[test_m] - f1_m
uni_m <- auto.arima(as.numeric(Y)[1:T0_m])
uni_m

Series: as.numeric(Y)[1:T0_m] 
ARIMA(1,1,1) 

Coefficients:
         ar1      ma1
      0.8864  -0.6513
s.e.  0.0635   0.1031

sigma^2 = 1.909:  log likelihood = -232.64
AIC=471.28   AICc=471.47   BIC=479.98

sel_m_tr <- select_auto(ts_head(X, T0_m), H = K_m - b)   # indicator: training data only
sel_m_tr$table

print(sel_m_tr$model)

[1] "ETS(A,N,N)"

xfc_m_tr <- function(o, H) list(X = sel_m_tr$f(ts_head(X, o), H)$mean)
eval_m <- rolling_eval(eq_m_tr, Y, list(X = X), a_tr, uni_m,
                       origins = T0_m:(length(Y) - 1), H = K_m, xfc = xfc_m_tr)

The coefficients estimated on the training data (ct) are very close to the full-sample estimates, so the model is stable. auto.arima() selects an ARIMA(1,1,1) for sales alone, and select_auto() selects ETS for the indicator.

c. One-step-ahead forecasts

The one-step-ahead forecasts of the 15 test periods are made from the origins t = 135, \dots, 149, every time with the newly observed data. With b = 3 they need no forecast of the indicator, so the ex ante and the ex post forecasts are identical (the dashed line lies on the solid one).

one_m <- onestep_of(eval_m)                               # forecasts at h = 1
plot_fc(one_m$t, one_m$Actual, fc_cols(one_m), "Series M: one-step-ahead forecasts", "sales")

acc1_m <- acc_of(one_m)                                   # accuracy of the plotted forecasts
acc1_m

The one-step-ahead forecasts of the transfer function follow the rise of sales from t = 143 on almost immediately, because the indicator rose three periods earlier; the univariate ARIMA reacts only after the rise has been observed. The transfer function reduces the RMSE from 0.9141 to 0.1733 and the MAPE from 0.3026% to 0.0558%.

d. k-step-ahead forecasts

Forecast path from the end of the training period. The forecasts made at T_0 = 135 for the whole test period (h = 1, \dots, 15) show what a forecaster would have obtained at that moment.

fk_m    <- tf_forecast(eq_m_tr, Y, list(X = X), a_tr, o = T0_m, H = length(test_m),
                       xfc = xfc_m_tr)                                   # ex ante
fk_ep_m <- tf_forecast(eq_m_tr, Y, list(X = X), a_tr, o = T0_m, H = length(test_m))  # ex post
fk_ar_m <- as.numeric(forecast(Arima(as.numeric(Y)[1:T0_m], model = uni_m),
                               h = length(test_m))$mean)
plot_fc(test_m, as.numeric(Y)[test_m],
        list(`TF ex ante` = fk_m, `TF ex post` = fk_ep_m, `Univariate ARIMA` = fk_ar_m),
        paste0("Series M: k-step-ahead forecasts from T0 = ", T0_m, " (h = 1, ..., ",
               length(test_m), ")"), "sales")

path_m <- accuracy_tbl(as.numeric(Y)[test_m], `TF ex ante` = fk_m, `TF ex post` = fk_ep_m,
                       `Univariate ARIMA` = fk_ar_m)
path_m

The first three forecasts of the path are accurate and identical in the ex ante and ex post versions. After that, the ex ante path relies on forecasts of the indicator, which cannot anticipate its later rise; it stays almost flat at about 258 and misses the increase of sales to about 263 (RMSE 3.03). With the observed indicator the path follows sales closely (RMSE 0.24). Both are better than the univariate path, which slowly declines (RMSE 4.44).

Rolling-origin evaluation by horizon. A single path depends on the situation at its origin. The rolling-origin evaluation therefore repeats the forecasts from every origin t = 135, \dots, 149, up to 6 steps ahead; the RMSE at horizon h uses every origin whose target lies in the test period, and the row h = 1 contains the one-step-ahead forecasts of part c. For h > 3 the ex ante forecasts use forecasts of the indicator, the ex post forecasts the observed indicator.

plot_rmse(eval_m, "Series M: RMSE of the k-step-ahead forecasts by horizon")

eval_m

Up to h = 3 the transfer function is far more accurate: its RMSE stays below 0.2 (0.173, 0.183 and 0.195), whereas the univariate error grows from 0.914 to 2.408, because the indicator has already announced the movement in sales. The ex ante and ex post columns are identical up to h = 3, since no forecast of the indicator is needed. From h = 4 on the picture changes. With the observed indicator (ex post) the error would stay at about 0.2, but in real time the indicator has to be forecast (here with ETS), and the ex ante RMSE jumps to 1.016, 1.866 and 2.661 at h = 4, 5, 6. It is still smaller than the univariate error (3.035 to 3.800), but the gap narrows quickly: beyond the delay, the errors of the indicator forecasts pass into the sales forecasts, multiplied by the impulse response weights of the transfer function (\hat\nu_3 = 4.69, \hat\nu_4 = 3.41, \dots).

e. Forecasts beyond the end of the data

The model estimated on all 150 observations (Step 8) forecasts sales for t = 151, \dots, 156. The forecasting method for the indicator is now selected on the whole series.

sel_m <- select_auto(X, H = K_m - b)                      # forecasting method for the indicator
sel_m$table

print(sel_m$model)

[1] "Theta"

xfc_m <- function(o, H) list(X = sel_m$f(ts_head(X, o), H)$mean)    # indicator forecasts
a_full_m <- c(rep(0, length(Y) - length(a_m)), a_m)
fc_m  <- tf_forecast(eq_m, Y, list(X = X), a_full_m, o = length(Y), H = K_m, xfc = xfc_m)

plot_future(Y, fc_m) + labs(title = "Series M: forecasts of sales for t = 151, ..., 156",
                            y = "sales")

X_used <- c(as.numeric(X), xfc_m(length(Y), K_m)$X)[length(Y) + 1:K_m - b]  # X[t-3]
tibble(h = 1:K_m, t = length(Y) + 1:K_m, `X[t-3]` = round(X_used, 3),
       `X[t-3] is` = ifelse(1:K_m <= b, "observed", "forecast"),
       `sales forecast` = round(fc_m, 2),
       `RMSE at h (test period)` = eval_m$`TF ex ante`[1:K_m])

The values after t = 150 are not yet known, so the accuracy of these forecasts cannot be measured; the last column gives the RMSE of the ex ante forecasts at the same horizon in the test period (part d), which shows the accuracy to expect. For h \le 3 only observed values of the indicator enter: the forecasts for t = 151, 152, 153 are 262.90, 264.24, 263.48, the same values as explicit_forecast(eq_m, Y, X, a_full_m, h = 3), with an expected error of about 0.17 to 0.19. To forecast the indicator, ARIMA, ETS and Theta are about equally accurate in the cross-validation (the naive forecast is a little worse, and the mean, useless for a series without a fixed level, much worse); Theta is selected. From h = 4 on, the change X_{t-3} - X_{t-4} in the equation involves these forecasts, which change only slowly (between -18.49 and -18.49). The sales forecasts for t = 154, 155, 156 (263.58, 263.72, 263.89) therefore change little: beyond the lead time of three periods, the indicator brings no new information about sales, and the expected error grows to about 2.66.

Conclusion. For Series M the transfer function (b,r,s) = (3,1,0) with MA(1) noise on the differenced series is adequate: all coefficients are significant, the residuals are white, normal and uncorrelated with the prewhitened input. This is the model reported by Box and Jenkins for Series M. The indicator leads sales by three periods. Within this lead time it reduces the forecast error by a factor of five to twelve; beyond it, the indicator has to be forecast itself, and the advantage shrinks quickly (30% at h = 6).

Case 2 — Series J: Gas Furnace

Problem Description

Box and Jenkins’ Series J comes from a gas furnace in which air and methane were combined to form a mixture of gases containing CO_2. The methane feed rate X_t (cubic feet per minute, coded) was varied, and the CO_2 concentration in the outlet gas Y_t (%) was recorded every 9 seconds, giving 296 pairs. The engineering problem is to describe the dynamic behaviour of the process: how long it takes before a change in the gas rate reaches the CO_2 measurement (the dead time b), how the response builds up and dies out, and how large the eventual change in CO_2 is (the steady-state gain). Such a model is the basis for designing a controller. The input is set by the experimenter, so there is no feedback. The data are SeriesJ (list with XJ and YJ) in the package arfima.

Data

data(SeriesJ)                          # from arfima (loaded at the top)
X <- SeriesJ$XJ; Y <- SeriesJ$YJ
tibble(t = seq_along(Y), `CO2 (%) Y` = as.numeric(Y), `Gas rate X` = as.numeric(X)) %>%
  pivot_longer(-t) %>%
  ggplot(aes(t, value)) +
    geom_line(colour = "steelblue") +
    facet_wrap(~ name, ncol = 1, scales = "free_y") +
    labs(x = "t (9-second intervals)", y = NULL) +
    theme_bw()

Both series oscillate smoothly around constant levels (about 0 for the gas rate and 53.5% for CO_2). The CO_2 series is a smoothed and inverted copy of the gas rate: when the gas rate rises, CO_2 falls a little later.

Step 1 — Stationarity

stat_tests(Y = Y, X = X)

x <- as.numeric(X); y <- as.numeric(Y); n <- length(y)

ADF (-3.82 and -4.09, p-values 0.018 and \le 0.01) and PP (p-value \le 0.01) reject a unit root for both series. KPSS rejects level stationarity (p-values \le 0.01 and 0.011), a reaction to the long, smooth swings of these strongly autocorrelated series rather than to a stochastic trend. The input model below has all its roots outside the unit circle, and Box and Jenkins model both series without differencing, so d = 0.

Step 2 — Input Model

(ggAcf(x, lag.max = 24) + ggtitle("ACF X")) | (ggPacf(x, lag.max = 24) + ggtitle("PACF X"))

fit_xj <- Arima(X, order = c(3, 0, 0))
fit_xj

Series: X 
ARIMA(3,0,0) with non-zero mean 

Coefficients:
         ar1      ar2     ar3     mean
      1.9691  -1.3651  0.3394  -0.0606
s.e.  0.0544   0.0985  0.0543   0.1898

sigma^2 = 0.03578:  log likelihood = 72.57
AIC=-135.14   AICc=-134.93   BIC=-116.69

coeftest(fit_xj)

z test of coefficients:

           Estimate Std. Error  z value  Pr(>|z|)    
ar1        1.969066   0.054385  36.2061 < 2.2e-16 ***
ar2       -1.365143   0.098538 -13.8540 < 2.2e-16 ***
ar3        0.339404   0.054328   6.2473 4.177e-10 ***
intercept -0.060643   0.189800  -0.3195    0.7493    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The ACF of X_t decays slowly in a damped wave (0.952, 0.834, 0.682, …), while the PACF has three large spikes (0.952, -0.788, 0.339) and is small afterwards. The input therefore follows an AR(3): (1 - 1.9691B + 1.3651B^2 - 0.3394B^3)(X_t - \mu_x) = \alpha_t, with all three coefficients highly significant (|z| \ge 6.25). The mean (-0.061) is not significant. Box and Jenkins obtained (1.97, -1.37, 0.34).

Step 3 — Prewhitening

phi_x <- coef(fit_xj)[c("ar1", "ar2", "ar3")]
alpha_j <- pw_filter(x - mean(x), phi = phi_x)
beta_j  <- pw_filter(y - mean(y), phi = phi_x)
Box.test(na.omit(alpha_j), lag = 12, type = "Ljung-Box", fitdf = 3)

    Box-Ljung test

data:  na.omit(alpha_j)
X-squared = 19.736, df = 9, p-value = 0.01961

Box.test(na.omit(alpha_j), lag = 24, type = "Ljung-Box", fitdf = 3)

    Box-Ljung test

data:  na.omit(alpha_j)
X-squared = 30.447, df = 21, p-value = 0.08337

The filter \pi(B) = 1 - 1.9691B + 1.3651B^2 - 0.3394B^3 is applied to the mean-corrected input and output. Q(12) = 19.74 (df = 9, p-value = 0.020) shows some remaining short-lag autocorrelation in \alpha_t, while Q(24) has p-value 0.083. The AR(3) filter is the classical one and whitens the input well enough for identification.

Step 4 — Cross-Correlation and Impulse Response Weights

tab_j <- ccf_table(beta_j, alpha_j, lag.max = 12)
plot_ccf(tab_j, "CCF of prewhitened series: Series J")

tab_j %>% filter(k >= 0) %>% mutate(across(c(r, nu), ~ round(.x, 3)))

The bounds are \pm 2/\sqrt{293} = \pm 0.117. The cross-correlations at lags 0–2 are negligible (-0.003, 0.051, -0.029). They become significant and negative from lag 3 to lag 7: -0.286, -0.336, -0.460, -0.273, -0.172, with \hat\nu_3, \dots, \hat\nu_7 = -0.552, -0.648, -0.888, -0.527, -0.332. From lag 8 onwards they are again inside the bounds. At negative lags only lag -6 is marginally outside (-0.121); one isolated value among twelve is expected by chance, so there is no feedback.

Step 5 — Orders (b, r, s)

• The first significant weight is \hat\nu_3, so b = 3.
• The weights \hat\nu_3 and \hat\nu_4 grow in magnitude up to the peak at \hat\nu_5, and only from \hat\nu_5 do they decay geometrically (\hat\nu_6/\hat\nu_5 = 0.59, \hat\nu_7/\hat\nu_6 = 0.63), so r = 1.
• There are s - r + 1 = 2 unpatterned weights (\hat\nu_3, \hat\nu_4) before the decay takes over, so s = 2.

The tentative transfer function is \dfrac{\omega_0 - \omega_1 B - \omega_2 B^2}{1 - \delta_1 B}\,x_{t-3}.

Step 6 — Preliminary Estimates

prelim_tf(tab_j$nu[tab_j$k >= 0], b = 3, r = 1, s = 2)

    omega0     omega1     omega2     delta1 
-0.5523223  0.3200399  0.5031095  0.5933936 

\hat\delta_1 = \hat\nu_6/\hat\nu_5 = 0.593, \hat\omega_0 = \hat\nu_3 = -0.552, \hat\omega_1 = \hat\delta_1\hat\nu_3 - \hat\nu_4 = 0.320 and \hat\omega_2 = \hat\delta_1\hat\nu_4 - \hat\nu_5 = 0.503. In Box–Jenkins notation this is \omega(B) = -0.552 - 0.320B - 0.503B^2 and \delta(B) = 1 - 0.593B.

Step 7 — Noise Model

b <- 3
pre_j <- arimax(y[(b+1):n], order = c(1, 0, 0), method = "ML",
                xtransf = data.frame(gas = x[1:(n-b)]), transfer = list(c(1, 2)))
round(coef(pre_j), 4)

      ar1 intercept   gas-AR1   gas-MA0   gas-MA1   gas-MA2 
   0.9693   53.6805    0.5045   -0.4889   -0.3949   -0.5815 

cp <- coef(pre_j)
eta_j <- y[(b+1):n] - cp["intercept"] -
  tf_component(x[1:(n-b)], cp[c("gas-MA0", "gas-MA1", "gas-MA2")], c(1, -cp["gas-AR1"]))
(ggAcf(eta_j, lag.max = 24) + ggtitle("ACF of preliminary noise")) |
  (ggPacf(eta_j, lag.max = 24) + ggtitle("PACF of preliminary noise"))

sapply(list(AR1 = c(1, 0, 0), AR2 = c(2, 0, 0), AR3 = c(3, 0, 0), ARMA11 = c(1, 0, 1)),
       function(o) Arima(eta_j, order = o)$aic)

       AR1        AR2        AR3     ARMA11 
128.956167   4.083638   5.180609  39.366761 

The proxy AR(1) coefficient is 0.969, so the noise is highly persistent. The ACF of \hat\eta_t decays smoothly (0.901, 0.733, 0.542, 0.368, …) and the PACF has two large values (0.901 and -0.415) followed by values close to the bounds (-0.120 at lag 3, then small). This is the pattern of an AR(2) noise, which also has the smallest AIC (4.08, against 5.18 for AR(3), 39.37 for ARMA(1,1) and 128.96 for AR(1)).

Step 8 — Joint Estimation

fit_j <- arimax(y[(b+1):n], order = c(2, 0, 0), method = "ML",
                xtransf = data.frame(gas = x[1:(n-b)]), transfer = list(c(1, 2)))
fit_j

Call:
arimax(x = y[(b + 1):n], order = c(2, 0, 0), method = "ML", xtransf = data.frame(gas = x[1:(n - 
    b)]), transfer = list(c(1, 2)))

Coefficients:
         ar1      ar2  intercept  gas-AR1  gas-MA0  gas-MA1  gas-MA2
      1.5272  -0.6289    53.3617   0.5490  -0.5310  -0.3801  -0.5180
s.e.  0.0467   0.0495     0.1375   0.0392   0.0738   0.1017   0.1086

sigma^2 estimated as 0.0571:  log likelihood = 2.08,  aic = 9.83

coeftest(fit_j)

z test of coefficients:

           Estimate Std. Error  z value  Pr(>|z|)    
ar1        1.527207   0.046724  32.6858 < 2.2e-16 ***
ar2       -0.628872   0.049471 -12.7118 < 2.2e-16 ***
intercept 53.361724   0.137495 388.0983 < 2.2e-16 ***
gas-AR1    0.549049   0.039190  14.0098 < 2.2e-16 ***
gas-MA0   -0.530985   0.073813  -7.1936 6.310e-13 ***
gas-MA1   -0.380104   0.101704  -3.7374  0.000186 ***
gas-MA2   -0.517972   0.108563  -4.7711 1.832e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

All coefficients are significant; the smallest test statistic is z = -3.74 for gas-MA1. The residual variance is \hat\sigma_a^2 = 0.0571 and AIC = 9.83. With the numerator written with plus signs, as arimax() reports it,

Y_t = 53.3617 + \frac{-0.5310 - 0.3801B - 0.5180B^2}{1 - 0.5490B}\,X_{t-3} + \frac{1}{1 - 1.5272B + 0.6289B^2}\,a_t, \qquad a_t \sim \mathcal{N}(0,\ 0.0571).

This reproduces the classical Box–Jenkins result, -(0.53 + 0.37B + 0.51B^2)/(1 - 0.57B) with AR(2) noise (1.53, -0.63). The steady-state gain is (-0.5310 - 0.3801 - 0.5180)/(1 - 0.5490) = -3.17: a sustained increase of one unit in the (coded) gas rate lowers the CO_2 concentration by about 3.2 percentage points.

Box, Jenkins and Reinsel (2008) also consider a second-order denominator. It is compared below.

fit_j22 <- arimax(y[(b+1):n], order = c(2, 0, 0), method = "ML",
                  xtransf = data.frame(gas = x[1:(n-b)]), transfer = list(c(2, 2)))
coeftest(fit_j22)

z test of coefficients:

            Estimate Std. Error  z value  Pr(>|z|)    
ar1        1.5272660  0.0467418  32.6745 < 2.2e-16 ***
ar2       -0.6288508  0.0494925 -12.7060 < 2.2e-16 ***
intercept 53.3619485  0.1376750 387.5936 < 2.2e-16 ***
gas-AR1    0.5599994  0.2025061   2.7653  0.005686 ** 
gas-AR2   -0.0079272  0.1431030  -0.0554  0.955824    
gas-MA0   -0.5313478  0.0743189  -7.1496 8.706e-13 ***
gas-MA1   -0.3744310  0.1467550  -2.5514  0.010729 *  
gas-MA2   -0.5122913  0.1502382  -3.4099  0.000650 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

c(AIC_r1 = fit_j$aic, AIC_r2 = fit_j22$aic)

   AIC_r1    AIC_r2 
 9.830622 11.826685 

With r = 2 the extra coefficient is \hat\delta_2 = -0.008 (s.e. 0.143, z = -0.06), and the standard error of \hat\delta_1 increases from 0.039 to 0.203. AIC rises from 9.83 to 11.83. Both models have the same s and hence the same effective sample, so their AICs are comparable. r = 1 is preferred.

Step 9 — Diagnostic Checking

a_j <- residuals(fit_j)
Box.test(a_j, lag = 12, type = "Ljung-Box", fitdf = 2)

    Box-Ljung test

data:  a_j
X-squared = 15.65, df = 10, p-value = 0.1101

Box.test(a_j, lag = 24, type = "Ljung-Box", fitdf = 2)

    Box-Ljung test

data:  a_j
X-squared = 27.969, df = 22, p-value = 0.1767

q0_j <- bind_rows(Q0_test(a_j, alpha_j[(b+1):n], K = 12, M = 4),
                  Q0_test(a_j, alpha_j[(b+1):n], K = 24, M = 4))
q0_j

jarque.bera.test(na.omit(a_j))

    Jarque Bera Test

data:  na.omit(a_j)
X-squared = 390.93, df = 2, p-value < 2.2e-16

(ggAcf(na.omit(a_j), lag.max = 24) + ggtitle("ACF of residuals")) |
  plot_ccf(ccf_table(a_j, alpha_j[(b+1):n], lag.max = 16) %>% filter(k >= 0),
           "CCF of alpha and residuals (k >= 0)")

The first two residuals are NA because the numerator needs x_{t-1} and x_{t-2}.

• Noise model: Q(12) = 15.65 (df = 10, p-value = 0.110) and Q(24) = 27.97 (df = 22, p-value = 0.177), so the residuals are white noise.
• Transfer function: Q_0(12) = 4.30 (df = 12 + 1 - 4 = 9, p-value = 0.891) and Q_0(24) = 20.0 (p-value = 0.518), so they are uncorrelated with the prewhitened input.
• Normality: Jarque–Bera rejects (X^2 = 390.9, p-value < 2.2\times10^{-16}). The residuals have heavy tails caused by a few large disturbances of the process. This does not bias the estimates, but prediction intervals based on normality should be read as approximate.

Step 10 — Explicit Equation

Multiplying the fitted model by \delta(B)\phi(B) = (1 - 0.5490B)(1 - 1.5272B + 0.6289B^2) gives A(B)\,Y_t = c + C(B)\,X_{t-3} + D(B)\,a_t with A(B) = \delta(B)\phi(B), C(B) = \omega(B)\phi(B), D(B) = \delta(B) (the noise has no MA part) and c = \mu\,\delta(1)\phi(1). The operators, as products of the model operators, with the estimated factors and multiplied out:

cf_j <- coef(fit_j)
eq_j <- tf_explicit(w = cf_j[c("gas-MA0", "gas-MA1", "gas-MA2")], delta = c(1, -cf_j["gas-AR1"]),
                    b = 3, phi = c(1, -cf_j[c("ar1", "ar2")]), mu = cf_j["intercept"])
show_operators(eq_j)

\begin{aligned} A(B) &= \delta(B)\,\phi(B) && \text{acting on } Y_{t} \\ &= \left(1 - 0.5490B\right)\left(1 - 1.5272B + 0.6289B^{2}\right) \\ &= 1 - 2.0763B + 1.4674B^{2} - 0.3453B^{3} \\[4pt] C(B) &= \omega(B)\,\phi(B) && \text{acting on } X_{t-3} \\ &= \left(-0.5310 - 0.3801B - 0.5180B^{2}\right)\left(1 - 1.5272B + 0.6289B^{2}\right) \\ &= -0.5310 + 0.4308B - 0.2714B^{2} + 0.5520B^{3} - 0.3257B^{4} \\[4pt] D(B) &= \delta(B) = 1 - 0.5490B && \text{acting on } a_{t} \\[4pt] c &= \hat\mu\,\delta(1)\,\phi(1) = 53.361724 \times 0.450951 \times 0.101665 = 2.446409 \end{aligned}

Applying B^k Z_t = Z_{t-k} and moving the lagged values of Y to the right-hand side gives the explicit equation. Both series were used without differencing or transformation (d = 0), so it is already the explicit equation of the original variables, the CO_2 concentration Y_t (%) and the coded gas rate X_t:

show_equation(eq_j)

\begin{aligned} Y_{t} ={}& 2.4464 + 2.0763\,Y_{t-1} - 1.4674\,Y_{t-2} + 0.3453\,Y_{t-3} \\ &{}- 0.5310\,X_{t-3} + 0.4308\,X_{t-4} - 0.2714\,X_{t-5} + 0.5520\,X_{t-6} - 0.3257\,X_{t-7} \\ &{}+ a_{t} - 0.5490\,a_{t-1} \end{aligned}

The left-hand operator has degree 3, so the explicit equation contains three lags of Y. The numerator (degree 2) multiplied by the AR(2) noise operator gives five input lags, from X_{t-3} to X_{t-7}. The noise has no MA part, and the only lagged shock, a_{t-1}, comes from the denominator \delta(B). The constant c = \hat\mu\,\hat\delta(1)\,\hat\phi(1) is derived in the first display.

Step 11 — Forecasting

The forecasts are computed recursively from the explicit equation of Step 10. With b = 3, the CO_2 forecasts for the next three sampling times (27 seconds) use only observed gas rates; from h = 4 on, the gas rate itself has to be forecast.

a. Forecast equations

For the origin t, the one-step-ahead forecast and the k-step-ahead forecast of CO_2 are

show_forecast_eq(eq_j)

\begin{aligned} \hat{Y}_t(1) ={}& 2.4464 + 2.0763\,Y_{t} - 1.4674\,Y_{t-1} + 0.3453\,Y_{t-2} \\ &{}- 0.5310\,X_{t-2} + 0.4308\,X_{t-3} - 0.2714\,X_{t-4} + 0.5520\,X_{t-5} - 0.3257\,X_{t-6} \\ &{}- 0.5490\,a_{t} \\[8pt] \hat{Y}_t(k) ={}& 2.4464 + 2.0763\,\hat{Y}_t(k-1) - 1.4674\,\hat{Y}_t(k-2) + 0.3453\,\hat{Y}_t(k-3) \\ &{}- 0.5310\,\hat{X}_t(k-3) + 0.4308\,\hat{X}_t(k-4) - 0.2714\,\hat{X}_t(k-5) + 0.5520\,\hat{X}_t(k-6) - 0.3257\,\hat{X}_t(k-7) \\ &{}- 0.5490\,\hat{a}_t(k-1) \end{aligned}

with \hat Y_t(m) = Y_{t+m}, \hat X_t(m) = X_{t+m} and \hat a_t(m) = a_{t+m} (the residuals) for m \le 0, and \hat a_t(m) = 0 for m > 0. For m > 0, \hat X_t(m) is the forecast of the gas rate (ex ante) or its observed value (ex post). The one-step-ahead forecast uses the gas rates X_{t-2}, \dots, X_{t-6}, which are all observed at the origin, and the last residual a_t, which enters through the denominator \delta(B). In the k-step-ahead equation, \hat X_t(k-3) is an observed value for k \le 3 and a forecast of the gas rate for k \ge 4.

b. Training and test periods

The model is re-estimated on the first 266 observations, and the last 30 form the test period. The benchmark is the model chosen by auto.arima() for CO_2 alone. The forecasting method for the gas rate is selected on the training data only, and from every origin t = 266, \dots, 295 forecasts are made up to K = 6 steps ahead, as far as the test period reaches.

K_j  <- 6                                                 # forecast horizons h = 1, ..., K
T0_j <- n - 30                                            # training period t = 1, ..., 266
fit_j_tr <- arimax(y[(b+1):T0_j], order = c(2, 0, 0), method = "ML",
                   xtransf = data.frame(gas = x[1:(T0_j-b)]), transfer = list(c(1, 2)))
ct <- coef(fit_j_tr)
eq_j_tr <- tf_explicit(w = ct[c("gas-MA0", "gas-MA1", "gas-MA2")], delta = c(1, -ct["gas-AR1"]),
                       b = 3, phi = c(1, -ct[c("ar1", "ar2")]), mu = ct["intercept"])
test_j <- (T0_j + 1):n
a_tr <- c(rep(0, b), residuals(fit_j_tr), rep(0, length(test_j))); a_tr[is.na(a_tr)] <- 0
f1_j <- onestep_explicit(eq_j_tr, y, x, a_tr, test_j)    # residuals of the test period
a_tr[test_j] <- y[test_j] - f1_j
uni_j <- auto.arima(y[1:T0_j])
uni_j

Series: y[1:T0_j] 
ARIMA(2,1,2) 

Coefficients:
         ar1      ar2      ma1     ma2
      1.4153  -0.6848  -0.3006  0.2986
s.e.  0.0901   0.0738   0.0957  0.0945

sigma^2 = 0.1105:  log likelihood = -83.26
AIC=176.52   AICc=176.75   BIC=194.42

c(sd_X_training = sd(x[1:T0_j]), sd_X_test = sd(x[test_j]))

sd_X_training     sd_X_test 
     1.120603      0.433584 

sel_j_tr <- select_auto(ts_head(x, T0_j), H = K_j - b)   # gas rate: training data only
sel_j_tr$table

print(sel_j_tr$model)

[1] "ARIMA(3,1,1)"

xfc_j_tr <- function(o, H) list(X = sel_j_tr$f(ts_head(x, o), H)$mean)
eval_j <- rolling_eval(eq_j_tr, y, list(X = x), a_tr, uni_j,
                       origins = T0_j:(n - 1), H = K_j, xfc = xfc_j_tr)

auto.arima() selects an ARIMA(2,1,2) for CO_2 alone, and select_auto() selects ARIMA for the gas rate. The test period is difficult for the transfer function, for two reasons. First, the gas rate moves much less than in the training period (standard deviation 0.43 against 1.12), so the input carries little information. Second, in the last ten observations the CO_2 level climbs to an unusually high level (about 58%, against a mean of 53.5%). The univariate model contains a unit root and follows such a shift immediately, whereas the stationary AR(2) noise of the transfer function pulls the forecasts back towards the mean.

c. One-step-ahead forecasts

The one-step-ahead forecasts of the 30 test periods are made from the origins t = 266, \dots, 295, every time with the newly observed data. With b = 3 they need no forecast of the gas rate, so the ex ante and the ex post forecasts are identical (the dashed line lies on the solid one).

one_j <- onestep_of(eval_j)                               # forecasts at h = 1
plot_fc(one_j$t, one_j$Actual, fc_cols(one_j), "Series J: one-step-ahead forecasts", "CO2 (%)")

acc1_j <- acc_of(one_j)                                   # accuracy of the plotted forecasts
acc1_j

Both models follow the CO_2 concentration closely one step ahead, but in this test period the univariate model is more accurate (RMSE 0.4161 against 0.5143, MAPE 0.5741% against 0.6927%). The difference arises mainly in the last ten periods, where CO_2 climbs to an unusually high level: the univariate model follows the climb at once, whereas the transfer function is pulled back towards the mean.

d. k-step-ahead forecasts

Forecast path from the end of the training period. The forecasts made at T_0 = 266 for the whole test period (h = 1, \dots, 30).

fk_j    <- tf_forecast(eq_j_tr, y, list(X = x), a_tr, o = T0_j, H = length(test_j),
                       xfc = xfc_j_tr)                                   # ex ante
fk_ep_j <- tf_forecast(eq_j_tr, y, list(X = x), a_tr, o = T0_j, H = length(test_j))  # ex post
fk_ar_j <- as.numeric(forecast(Arima(y[1:T0_j], model = uni_j), h = length(test_j))$mean)
plot_fc(test_j, y[test_j],
        list(`TF ex ante` = fk_j, `TF ex post` = fk_ep_j, `Univariate ARIMA` = fk_ar_j),
        paste0("Series J: k-step-ahead forecasts from T0 = ", T0_j, " (h = 1, ..., ",
               length(test_j), ")"), "CO2 (%)")

path_j <- accuracy_tbl(y[test_j], `TF ex ante` = fk_j, `TF ex post` = fk_ep_j,
                       `Univariate ARIMA` = fk_ar_j)
path_j

The gas rate was high (about 1) just before the origin t = 266, so the forecast path of the transfer function correctly forecasts a fall of CO_2 to about 51% in the first periods; it then settles at about 54%. The univariate model extrapolates the recent movement of CO_2 upwards, to about 58%, in the wrong direction. Neither model anticipates the climb at the end of the series, which is not announced by the gas rate. Over the 30 periods the ex ante and ex post paths are about equally accurate (RMSE 2.14 and 2.07) and clearly better than the univariate path (3.11).

Rolling-origin evaluation by horizon. From every origin t = 266, \dots, 295 both models forecast up to 6 steps ahead, as far as the test period reaches; the row h = 1 contains the one-step-ahead forecasts of part c. For h > 3 the ex ante forecasts use forecasts of the gas rate, the ex post forecasts the observed gas rate.

plot_rmse(eval_j, "Series J: RMSE of the k-step-ahead forecasts by horizon")

eval_j

At h = 1 the univariate model is more accurate (RMSE 0.416 against 0.514), and at h = 2 the two are close. From h = 3 on the transfer function wins, and its advantage grows with the horizon (2.160 against 2.956 at h = 6): the gas rates observed three to seven periods earlier determine where the CO_2 level is going, whereas the univariate model can only extrapolate. Forecasting the gas rate costs almost nothing here. The ex ante and ex post errors are practically equal, and the ex ante errors are even slightly smaller at h \ge 4: in this test period the small movements of the gas rate explain little of the movements of CO_2, in line with the low variability of the input.

e. Forecasts beyond the end of the data

The model estimated on all 296 observations (Step 8) forecasts CO_2 for t = 297, \dots, 302. The forecasting method for the gas rate is now selected on the whole series.

sel_j <- select_auto(X, H = K_j - b)                      # forecasting method for the gas rate
sel_j$table

print(sel_j$model)

[1] "ARIMA(3,1,1)"

xfc_j <- function(o, H) list(X = sel_j$f(ts_head(X, o), H)$mean)
a_full_j <- c(rep(0, n - length(a_j)), a_j); a_full_j[is.na(a_full_j)] <- 0
fc_j  <- tf_forecast(eq_j, Y, list(X = X), a_full_j, o = n, H = K_j, xfc = xfc_j)

plot_future(Y, fc_j) + labs(title = "Series J: forecasts of CO2 for t = 297, ..., 302",
                            y = "CO2 (%)")

X_used <- c(as.numeric(X), xfc_j(n, K_j)$X)[n + 1:K_j - b]          # X[t-3]
tibble(h = 1:K_j, t = n + 1:K_j, `X[t-3]` = round(X_used, 3),
       `X[t-3] is` = ifelse(1:K_j <= b, "observed", "forecast"),
       `CO2 forecast (%)` = round(fc_j, 2),
       `RMSE at h (test period)` = eval_j$`TF ex ante`[1:K_j])

The values after t = 296 are not yet known; the last column gives the RMSE of the ex ante forecasts at the same horizon in the test period (part d) as the accuracy to expect. The first three forecasts, 56.53, 56.06, 55.63%, use observed gas rates only and are the same as with explicit_forecast(). For the gas rate, ARIMA has by far the smallest cross-validation error. The selected model, ARIMA(3,1,1), is integrated, so its forecasts stay near the last observed gas rates (about -1.00) instead of returning to the mean of the series. From h = 4 on, the CO_2 forecasts use these gas rates; they continue to fall towards the mean level of about 53.5%, to 54.54% at h = 6, mainly because the stationary AR(2) noise pulls them back.

Conclusion. The gas furnace is described by (b,r,s) = (3,1,2) with AR(2) noise. All coefficients are significant, and the residuals are white and uncorrelated with the prewhitened input, although not normal. The model agrees with Box and Jenkins: a dead time of three sampling intervals (27 seconds), a response that peaks at lag 5, and a steady-state gain of -3.2. A second-order denominator is not needed. For forecasting, the transfer function pays off from h = 3 on, also when the gas rate has to be forecast, while at h = 1 the univariate model, which adapts faster to level shifts, is better in this test period.

Case 3 — SOI → Recruitment

Problem Description

The Southern Oscillation Index (SOI) measures changes in air pressure related to sea-surface temperatures in the central Pacific; its large swings are associated with El Niño and La Niña episodes. Recruitment is an index of the number of new fish entering a Pacific fish population. Both are monthly series for 1950–1987 (n = 453), analysed by Shumway and Stoffer. The ecological question is whether, and with what delay, the ocean–atmosphere conditions summarised by SOI drive recruitment. A clear lead would allow fisheries managers to anticipate recruitment several months in advance. The direction is one-way: fish cannot influence the SOI. The data are soi and rec in the package astsa.

Data

X <- soi; Y <- rec                     # from astsa (loaded at the top)
tibble(t = as.numeric(time(Y)), `Recruitment Y` = as.numeric(Y), `SOI X` = as.numeric(X)) %>%
  pivot_longer(-t) %>%
  ggplot(aes(t, value)) +
    geom_line(colour = "steelblue") +
    facet_wrap(~ name, ncol = 1, scales = "free_y") +
    labs(x = "year", y = NULL) +
    theme_bw()

tibble(month = month.abb,
       SOI = as.numeric(tapply(as.numeric(X), cycle(X), mean)),
       Recruitment = as.numeric(tapply(as.numeric(Y), cycle(Y), mean))) %>%
  mutate(across(where(is.double), ~ round(.x, 3)))

Both series fluctuate around stable levels, and both have an annual cycle. The monthly mean of SOI is positive from October to March (0.21 to 0.41) and negative from April to August (-0.28 to -0.09). Mean recruitment is about 72–74 in December–February and only 44–47 in July–August. On top of the cycle, SOI shows the irregular multi-year El Niño swings.

Step 1 — Stationarity and Seasonal Adjustment

The annual cycle is deterministic: it repeats with the same shape every year. A seasonal pattern shared by input and output produces spurious cross-correlations, so both series receive the same transformation: the monthly means are subtracted. x_t^* and y_t^* denote the adjusted series.

x <- desea(X); y <- desea(Y); n <- length(y)
stat_tests(SOI = X, REC = Y, `SOI*` = x, `REC*` = y)

ADF and PP reject a unit root for all four series (p-value \le 0.01). KPSS rejects level stationarity at 5% (p-values between \le 0.01 and 0.044), reflecting the long El Niño/La Niña swings, i.e. strong persistence rather than a unit root. In line with Shumway and Stoffer, the adjusted series are treated as stationary: d = 0.

Step 2 — Input Model

(ggAcf(x, lag.max = 36) + ggtitle("ACF SOI*")) | (ggPacf(x, lag.max = 36) + ggtitle("PACF SOI*"))

cand <- list(`AR(1)` = c(1, 0, 0), `AR(2)` = c(2, 0, 0), `AR(3)` = c(3, 0, 0),
             `ARMA(1,1)` = c(1, 0, 1), `ARMA(2,1)` = c(2, 0, 1))
imap_dfr(cand, function(o, nm) {
  f <- Arima(x, order = o, include.mean = FALSE)
  tibble(model = nm, AICc = f$aicc, BIC = f$bic,
         `Ljung-Box(12) p` = Box.test(residuals(f), 12, "Ljung-Box", fitdf = sum(o[c(1, 3)]))$p.value)
}) %>% mutate(across(where(is.double), ~ round(.x, 3)))

fit_xs <- Arima(x, order = c(1, 0, 1), include.mean = FALSE)
coeftest(fit_xs)

z test of coefficients:

     Estimate Std. Error z value  Pr(>|z|)    
ar1  0.877023   0.035165  24.941 < 2.2e-16 ***
ma1 -0.573194   0.059080  -9.702 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The ACF of x_t^* decays slowly (0.485, 0.362, 0.376, 0.358, 0.317, …) and the PACF tails off (0.485, 0.166, 0.200, 0.127, …). Neither function cuts off, which suggests a mixed model. Among the candidates, ARMA(1,1) has the smallest AICc (41.8) and BIC (54.1), and its residuals pass the Ljung–Box test at lag 12 (p-value = 0.083), whereas the pure AR models do not. The fitted input model is (1 - 0.8770B)\,x_t^* = (1 - 0.5732B)\,\alpha_t, with both coefficients highly significant.

Step 3 — Prewhitening

alpha_s <- pw_filter(x, phi = coef(fit_xs)["ar1"], theta = coef(fit_xs)["ma1"])
beta_s  <- pw_filter(y, phi = coef(fit_xs)["ar1"], theta = coef(fit_xs)["ma1"])
Box.test(na.omit(alpha_s), lag = 12, type = "Ljung-Box", fitdf = 2)

    Box-Ljung test

data:  na.omit(alpha_s)
X-squared = 16.586, df = 10, p-value = 0.08405

The filter \pi(B) = (1 - 0.8770B)/(1 - 0.5732B) is applied to both adjusted series. Q(12) = 16.59 (df = 10, p-value = 0.084), so \alpha_t behaves like white noise.

Step 4 — Cross-Correlation and Impulse Response Weights

tab_s <- ccf_table(beta_s, alpha_s, lag.max = 24)
plot_ccf(tab_s, "CCF of prewhitened series: SOI -> Recruitment")

tab_s %>% filter(k >= 0, k <= 12) %>% mutate(across(c(r, nu), ~ round(.x, 3)))

tab_s %>% filter(k < 0, significant)

The bounds are \pm 2/\sqrt{452} = \pm 0.094. Lags 0–4 are clearly insignificant (|r| \le 0.035). At lag 5 a large negative spike appears, r_{\alpha\beta}(5) = -0.361 (\hat\nu_5 = -20.4). It is followed by negative values that die down slowly: -0.267, -0.192, -0.203, -0.187, -0.182, \dots (\hat\nu_6 = -15.1, \hat\nu_7 = -10.8, …). At negative lags only lags -19 and -20 are marginally outside the bounds (0.107 and 0.099). Two isolated values among 24, far from lag 0, are compatible with chance: no feedback.

Why the seasonal adjustment matters. The same analysis on the unadjusted series, with an AR(1) prewhitening filter, gives a very different picture at negative lags.

ph <- coef(Arima(X, order = c(1, 0, 0)))["ar1"]
tab_raw <- ccf_table(pw_filter(as.numeric(Y) - mean(Y), phi = ph),
                     pw_filter(as.numeric(X) - mean(X), phi = ph), lag.max = 24)
plot_ccf(tab_raw, "CCF without seasonal adjustment (AR(1) filter)")

sum(tab_raw$significant & tab_raw$k < 0)

[1] 15

Without the adjustment, 15 of the 24 negative lags are significant, with absolute values up to 0.19. Read naively, this would suggest that recruitment “drives” the SOI. These correlations are spurious: they come from the annual cycle that both series share and that an AR(1) filter cannot remove. The positive lags show the same delay of 5 months in both versions.

Step 5 — Orders (b, r, s)

• The first significant weight is at lag 5, so b = 5.
• The weights decay from \hat\nu_5 onwards (\hat\nu_6/\hat\nu_5 = 0.74) without any unpatterned weights in front, so r = 1 and s = 0.

The tentative transfer function is \dfrac{\omega_0}{1 - \delta_1 B}\,x^*_{t-5}.

Step 6 — Preliminary Estimates

prelim_tf(tab_s$nu[tab_s$k >= 0], b = 5, r = 1, s = 0)

     omega0      delta1 
-20.3576279   0.7402023 

\hat\omega_0 = -20.36 and \hat\delta_1 = 0.740.

Step 7 — Noise Model

b <- 5
pre_s <- arimax(y[(b+1):n], order = c(1, 0, 0), method = "ML",
                xtransf = data.frame(soi = x[1:(n-b)]), transfer = list(c(1, 0)))
round(coef(pre_s), 4)

      ar1 intercept   soi-AR1   soi-MA0 
   0.8994   -0.0344    0.8270  -19.3867 

cp <- coef(pre_s)
eta_s <- y[(b+1):n] - cp["intercept"] -
  tf_component(x[1:(n-b)], cp["soi-MA0"], c(1, -cp["soi-AR1"]))
(ggAcf(eta_s, lag.max = 36) + ggtitle("ACF of preliminary noise")) |
  (ggPacf(eta_s, lag.max = 36) + ggtitle("PACF of preliminary noise"))

sapply(list(AR1 = c(1, 0, 0), AR2 = c(2, 0, 0), AR3 = c(3, 0, 0), ARMA11 = c(1, 0, 1)),
       function(o) Arima(eta_s, order = o)$aic)

     AR1      AR2      AR3   ARMA11 
3085.050 3006.955 3008.478 3019.746 

The ACF of \hat\eta_t decays (0.891, 0.714, 0.529, 0.363, …), while the PACF has two large values (0.891 and -0.392) and is small afterwards. This is the pattern of an AR(2) noise. Its AIC (3006.96) is the smallest, against 3008.48 for AR(3), 3019.75 for ARMA(1,1) and 3085.05 for AR(1).

Step 8 — Joint Estimation

fit_s <- arimax(y[(b+1):n], order = c(2, 0, 0), method = "ML",
                xtransf = data.frame(soi = x[1:(n-b)]), transfer = list(c(1, 0)))
fit_s

Call:
arimax(x = y[(b + 1):n], order = c(2, 0, 0), method = "ML", xtransf = data.frame(soi = x[1:(n - 
    b)]), transfer = list(c(1, 0)))

Coefficients:
         ar1      ar2  intercept  soi-AR1   soi-MA0
      1.2582  -0.4114    -0.1195   0.8002  -20.2993
s.e.  0.0431   0.0433     2.0994   0.0217    1.2381

sigma^2 estimated as 46.96:  log likelihood = -1498.89,  aic = 3007.77

coeftest(fit_s)

z test of coefficients:

            Estimate Std. Error  z value Pr(>|z|)    
ar1         1.258195   0.043091  29.1985   <2e-16 ***
ar2        -0.411447   0.043261  -9.5108   <2e-16 ***
intercept  -0.119469   2.099357  -0.0569   0.9546    
soi-AR1     0.800216   0.021694  36.8872   <2e-16 ***
soi-MA0   -20.299319   1.238098 -16.3956   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The two AR coefficients and both transfer function parameters are highly significant. The constant (-0.12) is not significant, as expected for mean-adjusted series. \hat\sigma_a^2 = 46.96. The fitted model is

y_t^* = -0.1195 + \frac{-20.2993}{1 - 0.8002B}\,x^*_{t-5} + \frac{1}{1 - 1.2582B + 0.4114B^2}\,a_t, \qquad a_t \sim (0,\ 46.96).

A rise of SOI is followed, five months later, by lower recruitment, and the effect decays slowly (\hat\delta_1 = 0.80). The steady-state gain is -20.30/(1 - 0.80) = -101.6: a sustained increase of 0.1 in SOI lowers recruitment by about 10 units. This is the same structure that Shumway and Stoffer reach: Recruitment depends on its own value one month earlier and on SOI five months earlier, with autocorrelated errors.

Over-fitting checks: adding a numerator term or a second denominator term does not help.

fit_s11 <- arimax(y[(b+1):n], order = c(2, 0, 0), method = "ML",
                  xtransf = data.frame(soi = x[1:(n-b)]), transfer = list(c(1, 1)))
fit_s20 <- arimax(y[(b+1):n], order = c(2, 0, 0), method = "ML",
                  xtransf = data.frame(soi = x[1:(n-b)]), transfer = list(c(2, 0)))
rbind(`(5,1,1): soi-MA1` = coeftest(fit_s11)["soi-MA1", ],
      `(5,2,0): soi-AR2` = coeftest(fit_s20)["soi-AR2", ])

                   Estimate Std. Error   z value  Pr(>|z|)
(5,1,1): soi-MA1 0.81870914 1.44873391 0.5651204 0.5719918
(5,2,0): soi-AR2 0.04151267 0.05859586 0.7084574 0.4786612

The extra numerator coefficient (z = 0.57) and the second denominator coefficient (z = 0.71) are not significant, so (b,r,s) = (5,1,0) is kept. AIC is not used here: arimax() loses the first s observations, so models with different s are fitted to different samples.

Step 9 — Diagnostic Checking

a_s <- residuals(fit_s)
Box.test(a_s, lag = 12, type = "Ljung-Box", fitdf = 2)

    Box-Ljung test

data:  a_s
X-squared = 27.641, df = 10, p-value = 0.00206

Box.test(a_s, lag = 24, type = "Ljung-Box", fitdf = 2)

    Box-Ljung test

data:  a_s
X-squared = 38.331, df = 22, p-value = 0.01682

q0_s <- bind_rows(Q0_test(a_s, alpha_s[(b+1):n], K = 12, M = 2),
                  Q0_test(a_s, alpha_s[(b+1):n], K = 24, M = 2))
q0_s

jarque.bera.test(a_s)

    Jarque Bera Test

data:  a_s
X-squared = 1728.1, df = 2, p-value < 2.2e-16

(ggAcf(a_s, lag.max = 36) + ggtitle("ACF of residuals")) |
  plot_ccf(ccf_table(a_s, alpha_s[(b+1):n], lag.max = 24) %>% filter(k >= 0),
           "CCF of alpha and residuals (k >= 0)")

The diagnostics are mixed, which is instructive.

• Transfer function: Q_0(12) = 16.63 (df = 11, p-value = 0.119), so the transfer function is adequate at K = 12. At K = 24 it is borderline (Q_0 = 35.86, p-value = 0.043).
• Noise model: The Ljung–Box test rejects (Q(12) = 27.64, df = 10, p-value = 0.002; Q(24) = 38.33, p-value = 0.017). The individual residual autocorrelations are small: the largest are -0.116 (lag 4), 0.098 (lag 6), -0.095 (lag 12) and -0.091 (lag 7), against a bound of 0.094. There is no clear pattern that an extra ARMA term could capture (see the seasonal attempt below).
• Normality: The residuals are far from normal (Jarque–Bera X^2 = 1728).

A natural attempt is to add a seasonal AR term to the noise model:

fit_s_sar <- arimax(ts(y[(b+1):n], frequency = 12), order = c(2, 0, 0),
                    seasonal = list(order = c(1, 0, 0), period = 12), method = "ML",
                    xtransf = data.frame(soi = x[1:(n-b)]), transfer = list(c(1, 0)))
coeftest(fit_s_sar)["sar1", ]

   Estimate  Std. Error     z value    Pr(>|z|) 
-0.10117293  0.04800096 -2.10772752  0.03505457 

c(AIC_AR2 = fit_s$aic, AIC_AR2_SAR1 = fit_s_sar$aic)

     AIC_AR2 AIC_AR2_SAR1 
    3007.773     3005.456 

Box.test(residuals(fit_s_sar), lag = 12, type = "Ljung-Box", fitdf = 3)

    Box-Ljung test

data:  residuals(fit_s_sar)
X-squared = 23.257, df = 9, p-value = 0.005645

The seasonal coefficient is significant (\hat\Phi_1 = -0.10, p-value = 0.035) and AIC falls slightly, from 3007.8 to 3005.5. However, the Ljung–Box test still rejects (Q(12) p-value = 0.006), so the extra term does not solve the problem.

Recruitment is known to behave non-linearly (Shumway and Stoffer fit threshold models to it), and heavy-tailed residuals inflate portmanteau statistics. The AR(2) noise is therefore kept as a parsimonious linear approximation. The point estimates are reliable, but standard errors and intervals should be read with care.

Step 10 — Explicit Equation

The model was fitted to the adjusted series, so the explicit equation is first written for y^*_t and x^*_t. Multiplying by \delta(B)\phi(B) = (1 - 0.8002B)(1 - 1.2582B + 0.4114B^2) gives A(B)\,y^*_t = c + C(B)\,x^*_{t-5} + D(B)\,a_t with the operators

cf_s <- coef(fit_s)
eq_s <- tf_explicit(w = cf_s["soi-MA0"], delta = c(1, -cf_s["soi-AR1"]), b = 5,
                    phi = c(1, -cf_s[c("ar1", "ar2")]), mu = cf_s["intercept"])
show_operators(eq_s, sym = c(Y = "y^*", X = "x^*"))

\begin{aligned} A(B) &= \delta(B)\,\phi(B) && \text{acting on } y^*_{t} \\ &= \left(1 - 0.8002B\right)\left(1 - 1.2582B + 0.4114B^{2}\right) \\ &= 1 - 2.0584B + 1.4183B^{2} - 0.3292B^{3} \\[4pt] C(B) &= \omega(B)\,\phi(B) && \text{acting on } x^*_{t-5} \\ &= -20.2993\left(1 - 1.2582B + 0.4114B^{2}\right) \\ &= -20.2993 + 25.5405B - 8.3521B^{2} \\[4pt] D(B) &= \delta(B) = 1 - 0.8002B && \text{acting on } a_{t} \\[4pt] c &= \hat\mu\,\delta(1)\,\phi(1) = -0.119469 \times 0.199784 \times 0.153252 = -0.00365780 \end{aligned}

The stars denote the series after subtracting the monthly means, y^*_t = \mathrm{Rec}_t - \overline{\mathrm{Rec}}_{m(t)} and x^*_t = \mathrm{SOI}_t - \overline{\mathrm{SOI}}_{m(t)}, where \mathrm{Rec}_t = Y_t is recruitment, \mathrm{SOI}_t = X_t the SOI, m(t) the month of t, and \overline{\mathrm{Rec}}_m, \overline{\mathrm{SOI}}_m the monthly means. To go back to recruitment and SOI themselves, these definitions are substituted. The coefficients of the lags do not change; the monthly means of all lagged terms collect in a constant that depends on the month,

c_m = c + \overline{\mathrm{Rec}}_m - \sum_{i} \pi_i\, \overline{\mathrm{Rec}}_{m-i} - \sum_{k} \omega^*_k\, \overline{\mathrm{SOI}}_{m-k} \qquad (\text{months counted modulo } 12),

where c, \pi_i and \omega^*_k are the constant and the coefficients of y^*_{t-i} and x^*_{t-k} in the explicit equation of the adjusted series. With the seasonal dummies D_{m,t}, equal to 1 if t falls in month m (m = 1 for January, …, m = 12 for December) and 0 otherwise, the month-specific constant is c_{m(t)} = \sum_{m=1}^{12} c_m D_{m,t}. The display shows the three steps: (1) the explicit equation of the adjusted series, (2) the definitions of y^* and x^* substituted, and (3) the explicit equation of recruitment and SOI, with the twelve monthly constants written out.

mm_rec <- tapply(as.numeric(Y), cycle(Y), mean)      # monthly means of recruitment
mm_soi <- tapply(as.numeric(X), cycle(X), mean)      # monthly means of SOI
c_s <- monthly_const(eq_s, mm_rec, list(X = mm_soi))
sym_s <- c(Y = "\\mathrm{Rec}", X = "\\mathrm{SOI}")  # symbols of the original variables
show_chain(
  list(e = eq_s, sym = c(Y = "y^*", X = "x^*")),                        # 1. adjusted series
  list(e = eq_s, per_line = 3,                                           # 2. substituted
       sym = c(Y = "\\left(\\mathrm{Rec}_{#} - \\overline{\\mathrm{Rec}}_{m(#)}\\right)",
               X = "\\left(\\mathrm{SOI}_{#} - \\overline{\\mathrm{SOI}}_{m(#)}\\right)")),
  list(e = eq_s, sym = sym_s, const = dummy_terms(c_s)))              # 3. original variables

\begin{aligned} y^*_{t} ={}& {-0.003658} + 2.0584\,y^*_{t-1} - 1.4183\,y^*_{t-2} + 0.3292\,y^*_{t-3} \\ &{}- 20.2993\,x^*_{t-5} + 25.5405\,x^*_{t-6} - 8.3521\,x^*_{t-7} \\ &{}+ a_{t} - 0.8002\,a_{t-1} \\[8pt] \left(\mathrm{Rec}_{t} - \overline{\mathrm{Rec}}_{m(t)}\right) ={}& {-0.003658} + 2.0584\,\left(\mathrm{Rec}_{t-1} - \overline{\mathrm{Rec}}_{m(t-1)}\right) - 1.4183\,\left(\mathrm{Rec}_{t-2} - \overline{\mathrm{Rec}}_{m(t-2)}\right) \\ &{}+ 0.3292\,\left(\mathrm{Rec}_{t-3} - \overline{\mathrm{Rec}}_{m(t-3)}\right) \\ &{}- 20.2993\,\left(\mathrm{SOI}_{t-5} - \overline{\mathrm{SOI}}_{m(t-5)}\right) + 25.5405\,\left(\mathrm{SOI}_{t-6} - \overline{\mathrm{SOI}}_{m(t-6)}\right) - 8.3521\,\left(\mathrm{SOI}_{t-7} - \overline{\mathrm{SOI}}_{m(t-7)}\right) \\ &{}+ a_{t} - 0.8002\,a_{t-1} \\[8pt] \mathrm{Rec}_{t} ={}& 1.3310\,D_{1,t} + 3.0330\,D_{2,t} + 4.6402\,D_{3,t} + 1.0426\,D_{4,t} + 1.4766\,D_{5,t} + 0.9059\,D_{6,t} \\ &{}+ 1.3209\,D_{7,t} + 2.5175\,D_{8,t} + 4.2346\,D_{9,t} + 4.6183\,D_{10,t} - 2.1622\,D_{11,t} + 2.9254\,D_{12,t} \\ &{}+ 2.0584\,\mathrm{Rec}_{t-1} - 1.4183\,\mathrm{Rec}_{t-2} + 0.3292\,\mathrm{Rec}_{t-3} \\ &{}- 20.2993\,\mathrm{SOI}_{t-5} + 25.5405\,\mathrm{SOI}_{t-6} - 8.3521\,\mathrm{SOI}_{t-7} \\ &{}+ a_{t} - 0.8002\,a_{t-1} \end{aligned}

The monthly means and the monthly constants c_m (the coefficients of D_{1,t}, \dots, D_{12,t}):

tibble(month = month.abb, `mean Rec` = round(as.numeric(mm_rec), 3),
       `mean SOI` = round(as.numeric(mm_soi), 4), `c[m]` = round(c_s, 4))

The monthly constants are the coefficients of the seasonal dummies in the explicit equation for recruitment. They range from -2.2 (Nov) to 4.6 (Mar) and are much smaller than the seasonal swing of recruitment itself (monthly means from 44 to 74), because the three lags of recruitment already carry the seasonal level: their coefficients sum to 2.0584 - 1.4183 + 0.3292 = 0.969, so c_m mainly reflects how the seasonal level of month m differs from that of the preceding months. Forecasts from this equation are identical to the forecasts of y^* plus the monthly mean of recruitment.

Step 11 — Forecasting

The forecasts are computed from the explicit equation of Step 10. In the computations the adjusted series are forecast and the monthly means of recruitment are added back afterwards, which is the same as using the equation of recruitment and SOI with the monthly constants c_m. With b = 5, recruitment can be forecast five months ahead from SOI values that are already observed; for longer horizons the adjusted SOI has to be forecast.

a. Forecast equations

For the origin t, the one-step-ahead forecast and the k-step-ahead forecast of recruitment are

show_forecast_eq(eq_s, sym = sym_s, const1 = "c_{m(t+1)}", constk = "c_{m(t+k)}")

\begin{aligned} \widehat{\mathrm{Rec}}_t(1) ={}& c_{m(t+1)} + 2.0584\,\mathrm{Rec}_{t} - 1.4183\,\mathrm{Rec}_{t-1} + 0.3292\,\mathrm{Rec}_{t-2} \\ &{}- 20.2993\,\mathrm{SOI}_{t-4} + 25.5405\,\mathrm{SOI}_{t-5} - 8.3521\,\mathrm{SOI}_{t-6} \\ &{}- 0.8002\,a_{t} \\[8pt] \widehat{\mathrm{Rec}}_t(k) ={}& c_{m(t+k)} + 2.0584\,\widehat{\mathrm{Rec}}_t(k-1) - 1.4183\,\widehat{\mathrm{Rec}}_t(k-2) + 0.3292\,\widehat{\mathrm{Rec}}_t(k-3) \\ &{}- 20.2993\,\widehat{\mathrm{SOI}}_t(k-5) + 25.5405\,\widehat{\mathrm{SOI}}_t(k-6) - 8.3521\,\widehat{\mathrm{SOI}}_t(k-7) \\ &{}- 0.8002\,\hat{a}_t(k-1) \end{aligned}

where c_{m(t+k)} is the monthly constant of Step 10 for the month of t + k, \widehat{\mathrm{Rec}}_t(m) = \mathrm{Rec}_{t+m}, \widehat{\mathrm{SOI}}_t(m) = \mathrm{SOI}_{t+m} and \hat a_t(m) = a_{t+m} (the residuals) for m \le 0, and \hat a_t(m) = 0 for m > 0. The one-step-ahead forecast needs SOI values up to t - 4 only. In the k-step-ahead equation, \widehat{\mathrm{SOI}}_t(k-5) is an observed value for k \le 5 and a forecast for k \ge 6 (ex ante; the observed value in the ex post version).

b. Training and test periods

The last 48 months are held out. The monthly means are recomputed from the training months only, and both models are fitted to the adjusted training data. The univariate benchmark is the model chosen by auto.arima() for the adjusted recruitment, and the forecasting method for the adjusted SOI is selected on the training months only. From every origin t = 405, \dots, 452 forecasts are made up to K = 12 months ahead, as far as the test period reaches. orig_s() adds the monthly means of the training months back, so that all forecasts are shown and evaluated for recruitment itself; this does not change the forecast errors.

K_s  <- 12                                                # forecast horizons h = 1, ..., K
T0_s <- n - 48                                            # training period: 405 months
xh <- desea(X, T0_s); yh <- desea(Y, T0_s)                # means of the training months
fit_s_tr <- arimax(yh[(b+1):T0_s], order = c(2, 0, 0), method = "ML",
                   xtransf = data.frame(soi = xh[1:(T0_s-b)]), transfer = list(c(1, 0)))
ct <- coef(fit_s_tr)
eq_s_tr <- tf_explicit(w = ct["soi-MA0"], delta = c(1, -ct["soi-AR1"]), b = 5,
                       phi = c(1, -ct[c("ar1", "ar2")]), mu = ct["intercept"])
test_s <- (T0_s + 1):n
a_tr <- c(rep(0, b), residuals(fit_s_tr), rep(0, length(test_s)))
f1_s <- onestep_explicit(eq_s_tr, yh, xh, a_tr, test_s)  # residuals of the test period
a_tr[test_s] <- yh[test_s] - f1_s
uni_s <- auto.arima(yh[1:T0_s])
uni_s

Series: yh[1:T0_s] 
ARIMA(2,0,1) with zero mean 

Coefficients:
         ar1      ar2      ma1
      1.5285  -0.5915  -0.3190
s.e.  0.1715   0.1592   0.2124

sigma^2 = 73.19:  log likelihood = -1443.71
AIC=2895.41   AICc=2895.51   BIC=2911.43

sel_s_tr <- select_auto(ts_head(xh, T0_s), H = K_s - b)  # SOI*: training months only
sel_s_tr$table

print(sel_s_tr$model)

[1] "ARIMA(0,1,2)"

xfc_s_tr <- function(o, H) list(X = sel_s_tr$f(ts_head(xh, o), H)$mean)
eval_s <- rolling_eval(eq_s_tr, yh, list(X = xh), a_tr, uni_s,
                       origins = T0_s:(n - 1), H = K_s, xfc = xfc_s_tr)
mm_rec_tr <- tapply(as.numeric(Y)[1:T0_s], cycle(Y)[1:T0_s], mean)   # training months
orig_s <- function(v, t) v + as.numeric(mm_rec_tr[cycle(Y)[t]])      # back to recruitment

The univariate benchmark is an ARMA(2,1) for the adjusted recruitment, and select_auto() selects ARIMA for the adjusted SOI.

c. One-step-ahead forecasts

The one-step-ahead forecasts of the 48 test months are made from the origins t = 405, \dots, 452, every time with the newly observed data. With b = 5 they need no forecast of SOI, so the ex ante and the ex post forecasts are identical (the dashed line lies on the solid one).

one_s <- onestep_of(eval_s, orig_s)                       # forecasts at h = 1, recruitment
plot_fc(one_s$t, one_s$Actual, fc_cols(one_s), "Recruitment: one-step-ahead forecasts",
        "recruitment")

acc1_s <- acc_of(one_s)                                   # accuracy of the plotted forecasts
acc1_s

Both models follow recruitment closely one month ahead, because most of its predictable variation comes from its own strong persistence. The SOI adds a smaller, though systematic, improvement: the RMSE of the transfer function is 8.2904 against 8.7985 for the univariate model, and the MAPE 8.3514% against 10.3711%.

d. k-step-ahead forecasts

Forecast path from the end of the training period. The forecasts made at T_0 = 405 (September 1983) for the 48 test months.

fk_s    <- tf_forecast(eq_s_tr, yh, list(X = xh), a_tr, o = T0_s, H = length(test_s),
                       xfc = xfc_s_tr)                                   # ex ante
fk_ep_s <- tf_forecast(eq_s_tr, yh, list(X = xh), a_tr, o = T0_s,
                       H = length(test_s))                               # ex post
fk_ar_s <- as.numeric(forecast(Arima(yh[1:T0_s], model = uni_s), h = length(test_s))$mean)
path_fc_s <- list(`TF ex ante` = orig_s(fk_s, test_s), `TF ex post` = orig_s(fk_ep_s, test_s),
                  `Univariate ARIMA` = orig_s(fk_ar_s, test_s))       # recruitment
plot_fc(test_s, as.numeric(Y)[test_s], path_fc_s,
        paste0("Recruitment: k-step-ahead forecasts from T0 = ", T0_s, " (h = 1, ..., ",
               length(test_s), ")"), "recruitment")

path_s <- do.call(accuracy_tbl, c(list(as.numeric(Y)[test_s]), path_fc_s))
path_s

At the origin the adjusted SOI was negative, and the integrated model selected for it on the training data (ARIMA(0,1,2)) forecasts it to stay at about -0.16 for the whole test period instead of returning to zero. Through the negative gain, the ex ante path therefore stays about 18 units above the seasonal level of recruitment (the monthly means), whereas the univariate path returns to the seasonal level within a few months. Recruitment did stay above its seasonal level in most of the following years, so this path happens to be the most accurate (RMSE 14.6 against 20.0 for the univariate model). It is even better than the ex post path with the observed SOI (20.8), which passes every swing of SOI on to recruitment. This is the result of one origin with a favourable starting situation, not a general property, as the rolling-origin evaluation below shows.

Rolling-origin evaluation by horizon. From every origin t = 405, \dots, 452 both models forecast up to 12 months ahead, as far as the test period reaches; the row h = 1 contains the one-step-ahead forecasts of part c. Up to h = 5 the transfer function needs only observed SOI values; for h > 5 the ex ante forecasts use forecasts of the adjusted SOI.

plot_rmse(eval_s, "Recruitment: RMSE of the k-step-ahead forecasts by horizon")

eval_s

Up to h = 5, where only observed SOI values are needed, the transfer function is more accurate at every horizon, although the gain is modest: its RMSE is 2 to 6% smaller (8.29 against 8.80 at h = 1). Beyond the delay the ranking reverses: from h = 6 on the univariate ARMA(2,1) is more accurate, by about a quarter at h = 9 (17.9 against 24.2). Part of the gap is the price of forecasting SOI (at h = 12 the ex ante RMSE is 25.0, the ex post RMSE 23.2), but even with the observed SOI the transfer function is less accurate than the univariate model beyond h = 5. At these horizons both forecasts are dominated by the noise model, and in this test period the ARMA(2,1) estimated for recruitment alone describes its medium-term persistence better than the AR(2) noise of the transfer function, which already failed the Ljung–Box test in Step 9. Over all origins, therefore, the ex post forecasts are better than the ex ante forecasts from h = 6 on, and the univariate model is the most accurate beyond the delay; conclusions about forecast accuracy rest on this evaluation rather than on the single path.

e. Forecasts beyond the end of the data

The model estimated on all 453 months (Step 8) forecasts recruitment for October 1987 to September 1988. The forecasting method for the adjusted SOI is now selected on the whole series.

sel_s <- select_auto(x, H = K_s - b)                      # forecasting method for SOI*
sel_s$table

print(sel_s$model)

[1] "ARIMA(2,1,1)"

xfc_s <- function(o, H) list(X = sel_s$f(ts_head(x, o), H)$mean)
a_full_s <- c(rep(0, n - length(a_s)), a_s)
f_adj <- tf_forecast(eq_s, y, list(X = x), a_full_s, o = n, H = K_s, xfc = xfc_s)
month_next <- ((cycle(Y)[n] + 1:K_s - 1) %% 12) + 1
fc_s <- f_adj + as.numeric(mm_rec[month_next])            # recruitment forecasts

plot_future(Y, fc_s, n_show = 60) +
  labs(title = "Recruitment: forecasts for October 1987 to September 1988", y = "recruitment")

x_used <- c(x, xfc_s(n, K_s)$X)[n + 1:K_s - b]                       # x*[t-5]
tibble(h = 1:K_s, month = paste(month.abb[month_next], c(rep(1987, 3), rep(1988, 9))),
       `SOI*[t-5]` = round(x_used, 3),
       `SOI*[t-5] is` = ifelse(1:K_s <= b, "observed", "forecast"),
       `adjusted forecast` = round(f_adj, 2), `Recruitment forecast` = round(fc_s, 2),
       `RMSE at h (test period)` = eval_s$`TF ex ante`[1:K_s])

The values after September 1987 are not yet known; the last column gives the RMSE of the ex ante forecasts at the same horizon in the test period (part d) as the accuracy to expect. The first five forecasts (October 1987 to February 1988) are 31.9, 49.9, 62.1, 55.5, 57.2; they use observed SOI values only. They lie well below the usual seasonal level for those months, for two reasons. First, recruitment itself had fallen far below its seasonal level by September 1987 (adjusted value -33.8), and it is very persistent. Second, the adjusted SOI was above its monthly average in April, May and August 1987, which lowers recruitment five months later.

From March 1988 on, the forecasts also depend on forecasts of the adjusted SOI. The selected model, ARIMA(2,1,1), is integrated: as in Step 1, the long El Niño and La Niña swings make the adjusted SOI look non-stationary. Its forecasts therefore stay at about 0.10, slightly above the monthly average, instead of returning slowly to zero as the stationary ARMA(1,1) of Step 2 would. Through the negative steady-state gain this keeps the forecasts of the adjusted recruitment about 10 units below the seasonal level until the end of the horizon (-9.5 in September 1988).

Conclusion. After removing the monthly means from both series, the prewhitened CCF shows a clean delay of five months and a slowly decaying negative response. This leads to (b,r,s) = (5,1,0) with AR(2) noise, the structure found by Shumway and Stoffer. The transfer function part passes the Q_0 check, but some small residual autocorrelation and strong non-normality remain, a sign of the non-linear behaviour of recruitment. The case also shows how a shared seasonal pattern creates spurious cross-correlations when it is not removed before prewhitening. For forecasting, SOI helps within its lead time of five months; beyond it, the univariate model was more accurate in the test period, even though a single forecast path from a favourable origin suggested the opposite.

Case 4 — LA Pollution: Temperature, Particulates and Mortality

Problem Description

The Los Angeles pollution–mortality study (Shumway, Azari and Pawitan, 1988) recorded, every week from 1970 to 1979 (n = 508), the average cardiovascular mortality in Los Angeles County, the temperature (°F) and the level of particulate pollution. The epidemiological problem is to quantify how weather and air pollution affect mortality, including effects that appear one or more weeks later, while allowing for the long-term decline of mortality over the decade. Temperature acts non-linearly, because both unusually cold and unusually hot weeks raise mortality. Following Shumway and Stoffer, it therefore enters through two inputs: the centred temperature T_t = \text{tempr}_t - 74.26 and its square T_t^2. With three inputs that are autocorrelated (all with an annual cycle), separate prewhitening filters would be awkward, so the LTF approach is used. Mortality cannot change the weather, and it plausibly does not change the pollution level, so there is no feedback. The data are cmort, tempr and part in the package astsa.

Data

M <- cmort; Te <- tempr; P <- part     # from astsa (loaded at the top)
tibble(t = as.numeric(time(M)), `Mortality Y` = as.numeric(M),
       `Temperature (F)` = as.numeric(Te), Particulates = as.numeric(P)) %>%
  pivot_longer(-t) %>%
  ggplot(aes(t, value)) +
    geom_line(colour = "steelblue") +
    facet_wrap(~ name, ncol = 1, scales = "free_y") +
    labs(x = "year", y = NULL) +
    theme_bw()

stat_tests(Mortality = M, Temperature = Te, Particulates = P)

y <- as.numeric(M); n <- length(y); tt <- 1:n
Tc <- as.numeric(Te) - mean(Te); T2 <- Tc^2; Pn <- as.numeric(P)
round(cor(cbind(T = Tc, T2 = T2, P = Pn)), 3)

        T    T2      P
T   1.000 0.076 -0.017
T2  0.076 1.000  0.086
P  -0.017 0.086  1.000

All three series have an annual cycle. Mortality peaks in winter and also declines steadily over the decade, while particulate levels are erratic with occasional high peaks. ADF and PP reject a unit root for every series (p-value \le 0.01). KPSS rejects level stationarity only for mortality (statistic 2.22), because of its downward trend, which is modelled with a deterministic trend t. The contemporaneous correlations between the inputs are small (|r| \le 0.09); they matter less than the strong autocorrelation of each input, which the free-form regression handles automatically.

Step 1 — Free-Form Model with a Proxy Noise

The free-form (linear transfer function) model includes lags 0, \dots, 6 (six weeks) of each input, the trend, and a proxy AR(1) noise:

Y_t = a + \beta t + \sum_{j=0}^{6}\nu_{T,j}T_{t-j} + \sum_{j=0}^{6}\nu_{T^2,j}T^2_{t-j} + \sum_{j=0}^{6}\nu_{P,j}P_{t-j} + \eta_t, \qquad (1-\phi_1 B)\,\eta_t = a_t .

k <- 6
XL <- cbind(trend = tt[(k+1):n], lagmat(Tc, k, "T"), lagmat(T2, k, "T2"), lagmat(Pn, k, "P"))
ff_l <- Arima(y[(k+1):n], order = c(1, 0, 0), xreg = XL)
round(coef(ff_l)[c("ar1", "intercept", "trend")], 4)

      ar1 intercept     trend 
   0.4099   78.5037   -0.0252 

w_l <- ff_weights(ff_l, c("T", "T2", "P"))
w_l %>% mutate(across(c(nu, se, t), ~ round(.x, 4)))

plot_ff(w_l, "Free-form weights with +/- 2 s.e. (proxy AR(1) noise)")

Lag k	0	1	2	3	4	5	6
T: \hat\nu_k (t)	0.032 (0.7)	-0.202 (-4.3)	-0.152 (-3.2)	-0.128 (-2.7)	-0.029 (-0.6)	-0.122 (-2.8)	-0.024 (-0.5)
T^2: \hat\nu_k (t)	0.0178 (8.0)	0.0093 (4.0)	0.0044 (1.9)	0.0018 (0.8)	0.0030 (1.3)	0.0002 (0.1)	0.0024 (1.1)
P: \hat\nu_k (t)	0.095 (3.6)	0.035 (1.4)	0.051 (1.9)	0.018 (0.6)	0.055 (1.9)	0.038 (1.4)	-0.008 (-0.3)

The trend coefficient is -0.025 per week (s.e. 0.003).

Step 2 — Stationarity of the Errors

The proxy coefficient is \hat\phi_1 = 0.41 (s.e. 0.04), far from 1. The errors are stationary, so no differencing is needed.

Step 3 — Orders (b, r, s) from the Weights

• Temperature T: \hat\nu_0 is not significant (t = 0.7), so b = 1. From lag 1 the weights are negative and die down (-0.20, -0.15, -0.13, with another significant value at lag 5), so r = 1 and s = 0: (1,1,0). A warmer week lowers mortality one week later.
• Temperature squared T^2: the weights are significant at lags 0 and 1 and roughly halve from one lag to the next (0.0178, 0.0093, 0.0044), so (0,1,0).
• Particulates P: \hat\nu_0 = 0.095 is significant. The later weights are positive but individually not significant (t between 1.4 and 1.9). Together they suggest a response that continues for several weeks, i.e. (0,1,0). The simpler alternative, a purely contemporaneous effect (0,0,0), is compared in Step 5.

Step 4 — Noise Model

eta_l <- y[(k+1):n] - coef(ff_l)["intercept"] - as.numeric(XL %*% coef(ff_l)[colnames(XL)])
(ggAcf(eta_l, lag.max = 30) + ggtitle("ACF of free-form errors")) |
  (ggPacf(eta_l, lag.max = 30) + ggtitle("PACF of free-form errors"))

sapply(list(AR1 = c(1, 0, 0), AR2 = c(2, 0, 0), ARMA11 = c(1, 0, 1)),
       function(o) Arima(eta_l, order = o)$aic)

     AR1      AR2   ARMA11 
3023.102 2976.190 2991.784 

The ACF of the free-form errors decays (0.41, 0.42, 0.23, 0.16, …) and the PACF cuts off after lag 2 (0.41, 0.30, then -0.02, -0.05, …). The noise is therefore AR(2), which also has the smallest AIC (2976.2, against 2991.8 for ARMA(1,1) and 3023.1 for AR(1)). Shumway and Stoffer also use AR(2) errors for this regression.

Step 5 — Rational Transfer Function with ARMA Noise

Temperature has b = 1, so the first observation is dropped and T_{t-1} is built by hand. The default number of optimiser iterations is not sufficient here (the fit stops with a convergence warning), so optim.control = list(maxit = 2000) is used.

idx <- 2:n
XT_l <- data.frame(Temp = Tc[idx - 1], Temp2 = T2[idx], Part = Pn[idx])
fit_l <- arimax(y[idx], order = c(2, 0, 0), method = "ML", xreg = data.frame(trend = tt[idx]),
                xtransf = XT_l, transfer = list(c(1, 0), c(1, 0), c(1, 0)),
                optim.control = list(maxit = 2000))
fit_l

Call:
arimax(x = y[idx], order = c(2, 0, 0), xreg = data.frame(trend = tt[idx]), method = "ML", 
    optim.control = list(maxit = 2000), xtransf = XT_l, transfer = list(c(1, 
        0), c(1, 0), c(1, 0)))

Coefficients:
         ar1     ar2  intercept    trend  Temp-AR1  Temp-MA0  Temp2-AR1  Temp2-MA0  Part-AR1  Part-MA0
      0.2675  0.3261    80.9972  -0.0267    0.6643   -0.2172     0.3466     0.0176    0.5632    0.1144
s.e.  0.0427  0.0433     2.3316   0.0035    0.0557    0.0312     0.1145     0.0020    0.1111    0.0216

sigma^2 estimated as 22.78:  log likelihood = -1512.01,  aic = 3044.03

coeftest(fit_l)

z test of coefficients:

            Estimate Std. Error z value  Pr(>|z|)    
ar1        0.2675246  0.0427142  6.2631 3.773e-10 ***
ar2        0.3261100  0.0433097  7.5297 5.085e-14 ***
intercept 80.9971719  2.3316285 34.7385 < 2.2e-16 ***
trend     -0.0267206  0.0035356 -7.5575 4.109e-14 ***
Temp-AR1   0.6643224  0.0556945 11.9280 < 2.2e-16 ***
Temp-MA0  -0.2171672  0.0311682 -6.9676 3.224e-12 ***
Temp2-AR1  0.3465859  0.1145277  3.0262  0.002476 ** 
Temp2-MA0  0.0175886  0.0020486  8.5856 < 2.2e-16 ***
Part-AR1   0.5631598  0.1110663  5.0705 3.968e-07 ***
Part-MA0   0.1144265  0.0215904  5.2999 1.159e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

fit_l0 <- arimax(y[idx], order = c(2, 0, 0), method = "ML", xreg = data.frame(trend = tt[idx]),
                 xtransf = XT_l, transfer = list(c(1, 0), c(1, 0), c(0, 0)),
                 optim.control = list(maxit = 2000))
c(AIC_particulates_r1 = fit_l$aic, AIC_particulates_r0 = fit_l0$aic)

AIC_particulates_r1 AIC_particulates_r0 
           3044.026            3053.628 

cf_l <- coef(fit_l)
gain_l <- c(Temp  = unname(cf_l["Temp-MA0"]  / (1 - cf_l["Temp-AR1"])),
            Temp2 = unname(cf_l["Temp2-MA0"] / (1 - cf_l["Temp2-AR1"])),
            Part  = unname(cf_l["Part-MA0"]  / (1 - cf_l["Part-AR1"])))
round(gain_l, 4)

   Temp   Temp2    Part 
-0.6470  0.0269  0.2619 

implied_l <- bind_rows(
  tibble(input = "T",  k = 0:6, nu = nu_weights(cf_l["Temp-MA0"],  c(1, -cf_l["Temp-AR1"]),  b = 1, K = 6)),
  tibble(input = "T2", k = 0:6, nu = nu_weights(cf_l["Temp2-MA0"], c(1, -cf_l["Temp2-AR1"]), b = 0, K = 6)),
  tibble(input = "P",  k = 0:6, nu = nu_weights(cf_l["Part-MA0"],  c(1, -cf_l["Part-AR1"]),  b = 0, K = 6)))
plot_ff(w_l, "Free-form weights (bars) and weights implied by the rational model (dashed)",
        implied = implied_l)

All ten coefficients are significant; the smallest test statistic is z = 3.03 for Temp2-AR1. With a constant denominator for particulates, AIC increases from 3044.0 to 3053.6, so the decaying response (0,1,0) is kept. The fitted model is

\begin{aligned} Y_t ={}& 80.997 - 0.0267\,t + \frac{-0.2172}{1 - 0.6643B}\,T_{t-1} + \frac{0.0176}{1 - 0.3466B}\,T_t^2 + \frac{0.1144}{1 - 0.5632B}\,P_t \\ &+ \frac{1}{1 - 0.2675B - 0.3261B^2}\,a_t, \qquad a_t \sim (0,\ 22.78). \end{aligned}

The implied weights (dashed lines) follow the free-form weights closely while using only two parameters per input.

Interpretation.

• Temperature. The long-run gains are -0.647 for T and 0.0269 for T^2, so the steady-state temperature effect -0.647\,T + 0.0269\,T^2 is U-shaped. Its minimum lies at T = 0.647/(2 \times 0.0269) = 12.0, i.e. about 86 °F. In Los Angeles most weeks are colder than this, so colder-than-usual weeks raise mortality, with the main effect one week later, while very hot weeks raise it again through the quadratic term.
• Particulates. Their effect starts in the same week (0.114 per unit) and decays by a factor of 0.56 per week, giving a long-run gain of 0.262. A permanent reduction of particulates by 10 units (about two-thirds of their standard deviation) would lower expected weekly mortality by about 10 \times 0.262 = 2.6 units.

Step 6 — Diagnostic Checking

For Q_0, every input is prewhitened with its own ARIMA model, chosen here by auto.arima(), and aligned with the residuals (M = r + s + 1 = 2 for each input).

a_l <- residuals(fit_l)
Box.test(a_l, lag = 12, type = "Ljung-Box", fitdf = 2)

    Box-Ljung test

data:  a_l
X-squared = 7.5298, df = 10, p-value = 0.6747

Box.test(a_l, lag = 24, type = "Ljung-Box", fitdf = 2)

    Box-Ljung test

data:  a_l
X-squared = 24.086, df = 22, p-value = 0.3428

jarque.bera.test(a_l)

    Jarque Bera Test

data:  a_l
X-squared = 15.559, df = 2, p-value = 0.0004182

inputs_l <- list(Temp = Tc, Temp2 = T2, Part = Pn); shift_l <- c(Temp = 1, Temp2 = 0, Part = 0)
q0_l <- imap_dfr(inputs_l, function(z, nm) {
  fa <- auto.arima(z, seasonal = FALSE)
  al <- as.numeric(residuals(fa))
  Q0_test(a_l, al[idx - shift_l[nm]], K = 12, M = 2) %>%
    mutate(input = nm, `input model` = paste0("ARIMA(", paste(arimaorder(fa), collapse = ","), ")"),
           .before = 1)
})
q0_l

ggAcf(a_l, lag.max = 52) + ggtitle("ACF of residuals")

• Noise model: Q(12) = 7.53 (df = 10, p-value = 0.675) and Q(24) = 24.09 (df = 22, p-value = 0.343), so the residuals are white noise, even at the seasonal lag 52.
• Transfer functions: the Q_0 statistics are not significant for any input. The p-values are 0.288 for temperature (input model ARMA(1,2)), 0.480 for temperature squared (ARMA(1,1)) and 0.109 for particulates (ARMA(3,2)).
• Normality: Jarque–Bera rejects mildly (X^2 = 15.6, p-value = 0.0004), because of a few weeks with unusually high mortality.

Step 7 — Explicit Equation

The explicit equation is obtained by multiplying the model by all denominators and by the AR operator of the noise, (1 - 0.6643B)(1 - 0.3466B)(1 - 0.5632B)(1 - 0.2675B - 0.3261B^2). Because of the trend, the result also contains a term in t (t = 1 is the first week of 1970). T_t = \text{tempr}_t - 74.26 is the centred temperature, T^2_t its square and P_t the particulates. The operators and the constant:

eq_l <- tf_explicit_multi(list(
  list(name = "T",  w = cf_l["Temp-MA0"],  delta = c(1, -cf_l["Temp-AR1"]),  b = 1),
  list(name = "T2", w = cf_l["Temp2-MA0"], delta = c(1, -cf_l["Temp2-AR1"]), b = 0),
  list(name = "P",  w = cf_l["Part-MA0"],  delta = c(1, -cf_l["Part-AR1"]),  b = 0)),
  phi = c(1, -cf_l[c("ar1", "ar2")]), mu = cf_l["intercept"], beta = cf_l["trend"])
show_operators(eq_l, sym = c(T2 = "T^2"))

\begin{aligned} A(B) &= \delta_{T}(B)\,\delta_{T^2}(B)\,\delta_{P}(B)\,\phi(B) && \text{acting on } Y_{t} \\ &= \left(1 - 0.6643B\right)\left(1 - 0.3466B\right)\left(1 - 0.5632B\right)\left(1 - 0.2675B - 0.3261B^{2}\right) \\ &= 1 - 1.8416B + 0.8945B^{2} + 0.1698B^{3} - 0.2261B^{4} + 0.04228B^{5} \\[4pt] C_{T}(B) &= \omega_{T}(B)\,\delta_{T^2}(B)\,\delta_{P}(B)\,\phi(B) && \text{acting on } T_{t-1} \\ &= -0.2172\left(1 - 0.3466B\right)\left(1 - 0.5632B\right)\left(1 - 0.2675B - 0.3261B^{2}\right) \\ &= -0.2172 + 0.2557B - 0.02442B^{2} - 0.05309B^{3} + 0.01382B^{4} \\[4pt] C_{T^2}(B) &= \omega_{T^2}(B)\,\delta_{T}(B)\,\delta_{P}(B)\,\phi(B) && \text{acting on } T^2_{t} \\ &= 0.01759\left(1 - 0.6643B\right)\left(1 - 0.5632B\right)\left(1 - 0.2675B - 0.3261B^{2}\right) \\ &= 0.01759 - 0.02630B + 0.006620B^{2} + 0.005280B^{3} - 0.002146B^{4} \\[4pt] C_{P}(B) &= \omega_{P}(B)\,\delta_{T}(B)\,\delta_{T^2}(B)\,\phi(B) && \text{acting on } P_{t} \\ &= 0.1144\left(1 - 0.6643B\right)\left(1 - 0.3466B\right)\left(1 - 0.2675B - 0.3261B^{2}\right) \\ &= 0.1144 - 0.1463B + 0.01998B^{2} + 0.03067B^{3} - 0.008592B^{4} \\[4pt] D(B) &= \delta_{T}(B)\,\delta_{T^2}(B)\,\delta_{P}(B) && \text{acting on } a_{t} \\ &= \left(1 - 0.6643B\right)\left(1 - 0.3466B\right)\left(1 - 0.5632B\right) \\ &= 1 - 1.5741B + 0.7995B^{2} - 0.1297B^{3} \\[4pt] c &= \hat\mu\,A(1) - \hat\beta\textstyle\sum_k k\,A_k \\ &= 80.997172 \times 0.0389359 - (-0.0267206) \times (-0.236029) = 3.147391 \\[4pt] \beta^* &= \hat\beta\,A(1) = -0.0267206 \times 0.0389359 = -0.00104039 \end{aligned}

The trend needs care: A(B)\,(\mu + \beta t) = \sum_k A_k\,[\mu + \beta(t - k)] = \mu A(1) + \beta A(1)\,t - \beta \sum_k k A_k, which gives the constant c and the coefficient \beta^* of t in the last two lines. Applying B^k Z_t = Z_{t-k} gives the explicit equation in the variables of the model:

show_equation(eq_l, sym = c(T2 = "T^2"))

\begin{aligned} Y_{t} ={}& 3.1474 - 0.001040\,t + 1.8416\,Y_{t-1} - 0.8945\,Y_{t-2} - 0.1698\,Y_{t-3} \\ &{}+ 0.2261\,Y_{t-4} - 0.04228\,Y_{t-5} \\ &{}- 0.2172\,T_{t-1} + 0.2557\,T_{t-2} - 0.02442\,T_{t-3} - 0.05309\,T_{t-4} + 0.01382\,T_{t-5} \\ &{}+ 0.01759\,T^2_{t} - 0.02630\,T^2_{t-1} + 0.006620\,T^2_{t-2} + 0.005280\,T^2_{t-3} - 0.002146\,T^2_{t-4} \\ &{}+ 0.1144\,P_{t} - 0.1463\,P_{t-1} + 0.01998\,P_{t-2} + 0.03067\,P_{t-3} - 0.008592\,P_{t-4} \\ &{}+ a_{t} - 1.5741\,a_{t-1} + 0.7995\,a_{t-2} - 0.1297\,a_{t-3} \end{aligned}

The left-hand operator has degree 5 (three first-order denominators and the AR(2) noise), so mortality depends on its last five values. Each input enters with five consecutive lags: temperature at lags 1–5, squared temperature and particulates at lags 0–4. The noise adds three past shocks.

No differencing and no seasonal adjustment were used; the only transformation is the centring of the temperature. Substituting Y_t = \mathrm{cmort}_t (mortality), T_t = \mathrm{tempr}_t - 74.2604 (temperature in °F), T^2_t = (\mathrm{tempr}_t - 74.2604)^2 and P_t = \mathrm{part}_t (particulates) gives the explicit equation of the original variables:

tc_l  <- num(mean(Te), 4)                                # centring constant of temperature
sym_l <- c(Y = "\\mathrm{cmort}", P = "\\mathrm{part}",  # symbols of the original variables
           T  = paste0("\\left(\\mathrm{tempr}_{#} - ", tc_l, "\\right)"),
           T2 = paste0("\\left(\\mathrm{tempr}_{#} - ", tc_l, "\\right)^{2}"))
show_equation(eq_l, sym = sym_l, per_line = 3)

\begin{aligned} \mathrm{cmort}_{t} ={}& 3.1474 - 0.001040\,t + 1.8416\,\mathrm{cmort}_{t-1} \\ &{}- 0.8945\,\mathrm{cmort}_{t-2} - 0.1698\,\mathrm{cmort}_{t-3} + 0.2261\,\mathrm{cmort}_{t-4} \\ &{}- 0.04228\,\mathrm{cmort}_{t-5} \\ &{}- 0.2172\,\left(\mathrm{tempr}_{t-1} - 74.2604\right) + 0.2557\,\left(\mathrm{tempr}_{t-2} - 74.2604\right) - 0.02442\,\left(\mathrm{tempr}_{t-3} - 74.2604\right) \\ &{}- 0.05309\,\left(\mathrm{tempr}_{t-4} - 74.2604\right) + 0.01382\,\left(\mathrm{tempr}_{t-5} - 74.2604\right) \\ &{}+ 0.01759\,\left(\mathrm{tempr}_{t} - 74.2604\right)^{2} - 0.02630\,\left(\mathrm{tempr}_{t-1} - 74.2604\right)^{2} + 0.006620\,\left(\mathrm{tempr}_{t-2} - 74.2604\right)^{2} \\ &{}+ 0.005280\,\left(\mathrm{tempr}_{t-3} - 74.2604\right)^{2} - 0.002146\,\left(\mathrm{tempr}_{t-4} - 74.2604\right)^{2} \\ &{}+ 0.1144\,\mathrm{part}_{t} - 0.1463\,\mathrm{part}_{t-1} + 0.01998\,\mathrm{part}_{t-2} \\ &{}+ 0.03067\,\mathrm{part}_{t-3} - 0.008592\,\mathrm{part}_{t-4} \\ &{}+ a_{t} - 1.5741\,a_{t-1} + 0.7995\,a_{t-2} \\ &{}- 0.1297\,a_{t-3} \end{aligned}

The centring is kept inside the brackets on purpose. Multiplying the squares out would turn the small coefficients of T^2 into large coefficients of \mathrm{tempr} and \mathrm{tempr}^2 that almost cancel each other, which makes the equation very sensitive to rounding. Because of its length, the explicit form is mainly useful for forecasting; the rational form above is easier to interpret.

Step 8 — Forecasting

The forecasts are computed recursively from the explicit equation of Step 7. Temperature enters with delay b = 1, but squared temperature and particulates enter at lag 0, so every forecast, even one step ahead, needs input values that are not yet observed. Temperature and particulates are weekly series with an annual cycle; they are forecast with the seasonal method chosen by select_auto(). Squared temperature is not forecast separately: because E(T^2) = [E(T)]^2 + \text{Var}(T), its forecast is the square of the temperature forecast plus the variance of the temperature forecast error. Simply squaring the temperature forecast would underestimate T^2 and hence mortality.

a. Forecast equations

For the origin t, the one-step-ahead forecast and the k-step-ahead forecast of mortality Y_t = \mathrm{cmort}_t are, with T_t = \mathrm{tempr}_t - 74.2604 and P_t = \mathrm{part}_t,

show_forecast_eq(eq_l, sym = c(T2 = "T^2"))

\begin{aligned} \hat{Y}_t(1) ={}& 3.1474 - 0.001040\,(t+1) + 1.8416\,Y_{t} - 0.8945\,Y_{t-1} - 0.1698\,Y_{t-2} \\ &{}+ 0.2261\,Y_{t-3} - 0.04228\,Y_{t-4} \\ &{}- 0.2172\,T_{t} + 0.2557\,T_{t-1} - 0.02442\,T_{t-2} - 0.05309\,T_{t-3} + 0.01382\,T_{t-4} \\ &{}+ 0.01759\,T^2_{t+1} - 0.02630\,T^2_{t} + 0.006620\,T^2_{t-1} + 0.005280\,T^2_{t-2} - 0.002146\,T^2_{t-3} \\ &{}+ 0.1144\,P_{t+1} - 0.1463\,P_{t} + 0.01998\,P_{t-1} + 0.03067\,P_{t-2} - 0.008592\,P_{t-3} \\ &{}- 1.5741\,a_{t} + 0.7995\,a_{t-1} - 0.1297\,a_{t-2} \\[8pt] \hat{Y}_t(k) ={}& 3.1474 - 0.001040\,(t+k) + 1.8416\,\hat{Y}_t(k-1) - 0.8945\,\hat{Y}_t(k-2) - 0.1698\,\hat{Y}_t(k-3) \\ &{}+ 0.2261\,\hat{Y}_t(k-4) - 0.04228\,\hat{Y}_t(k-5) \\ &{}- 0.2172\,\hat{T}_t(k-1) + 0.2557\,\hat{T}_t(k-2) - 0.02442\,\hat{T}_t(k-3) - 0.05309\,\hat{T}_t(k-4) + 0.01382\,\hat{T}_t(k-5) \\ &{}+ 0.01759\,\widehat{T^2}_t(k) - 0.02630\,\widehat{T^2}_t(k-1) + 0.006620\,\widehat{T^2}_t(k-2) + 0.005280\,\widehat{T^2}_t(k-3) - 0.002146\,\widehat{T^2}_t(k-4) \\ &{}+ 0.1144\,\hat{P}_t(k) - 0.1463\,\hat{P}_t(k-1) + 0.01998\,\hat{P}_t(k-2) + 0.03067\,\hat{P}_t(k-3) - 0.008592\,\hat{P}_t(k-4) \\ &{}- 1.5741\,\hat{a}_t(k-1) + 0.7995\,\hat{a}_t(k-2) - 0.1297\,\hat{a}_t(k-3) \end{aligned}

with \hat Z_t(m) = Z_{t+m} for m \le 0 (the observed values of Y, T, T^2, P and the residuals a) and \hat a_t(m) = 0 for m > 0. The trend enters with the time of the forecast, t + 1 or t + k. In the one-step-ahead equation, T^2_{t+1} and P_{t+1} belong to the forecast week: in real time they are replaced by their forecasts \widehat{T^2}_t(1) and \hat P_t(1) (ex ante), and the ex post version uses the observed values. Likewise, in the k-step-ahead equation every \hat T_t(m), \widehat{T^2}_t(m) and \hat P_t(m) with m > 0 is a forecast of the input (ex ante) or its observed value (ex post).

b. Training and test periods

To measure the accuracy, the model is re-estimated on the first 456 weeks, and the last 52 weeks form the test period. The benchmark is the seasonal ARIMA chosen by auto.arima(). Temperature and particulates are forecast with the methods selected on the training weeks only, and from every origin t = 456, \dots, 507 mortality is forecast up to K = 13 weeks ahead (one quarter), as far as the test period reaches.

K_l  <- 13                                                # forecast horizons: one quarter
T0_l <- n - 52; idt <- 2:T0_l                             # training period: 456 weeks
fit_l_tr <- arimax(y[idt], order = c(2, 0, 0), method = "ML", xreg = data.frame(trend = tt[idt]),
                   xtransf = data.frame(Temp = Tc[idt - 1], Temp2 = T2[idt], Part = Pn[idt]),
                   transfer = list(c(1, 0), c(1, 0), c(1, 0)), optim.control = list(maxit = 2000))
ct <- coef(fit_l_tr)
eq_l_tr <- tf_explicit_multi(list(
  list(name = "T",  w = ct["Temp-MA0"],  delta = c(1, -ct["Temp-AR1"]),  b = 1),
  list(name = "T2", w = ct["Temp2-MA0"], delta = c(1, -ct["Temp2-AR1"]), b = 0),
  list(name = "P",  w = ct["Part-MA0"],  delta = c(1, -ct["Part-AR1"]),  b = 0)),
  phi = c(1, -ct[c("ar1", "ar2")]), mu = ct["intercept"], beta = ct["trend"])
test_l <- (T0_l + 1):n
a_tr <- c(0, residuals(fit_l_tr), rep(0, length(test_l)))
f_tf_l <- onestep_explicit_multi(eq_l_tr, y, list(T = Tc, T2 = T2, P = Pn), a_tr, test_l)
a_tr[test_l] <- y[test_l] - f_tf_l                        # residuals of the test period
uni_l <- auto.arima(ts(y[1:T0_l], frequency = 52))
uni_l

Series: ts(y[1:T0_l], frequency = 52) 
ARIMA(2,0,2)(0,1,0)[52] with drift 

Coefficients:
         ar1     ar2      ma1     ma2    drift
      0.6084  0.0108  -0.4105  0.2747  -0.0251
s.e.  0.3489  0.3015   0.3463  0.2076   0.0163

sigma^2 = 57.28:  log likelihood = -1388.64
AIC=2789.27   AICc=2789.48   BIC=2813.28

Tc_ts <- ts(Tc, start = start(Te), frequency = 52)        # inputs as weekly time series
P_ts  <- ts(Pn, start = start(P),  frequency = 52)
sel_tempr_tr <- select_auto(ts_head(Tc_ts, T0_l), H = K_l, methods = auto_methods_seasonal)
sel_part_tr  <- select_auto(ts_head(P_ts, T0_l),  H = K_l, methods = auto_methods_seasonal)
bind_rows(temperature = sel_tempr_tr$table, particulates = sel_part_tr$table, .id = "input")

print(c(temperature = sel_tempr_tr$model, particulates = sel_part_tr$model))

                            temperature                            particulates 
"ARIMA(0,0,0) errors + 3 Fourier pairs" "ARIMA(1,1,2) errors + 2 Fourier pairs" 

## input forecasts from the origin o, given the methods selected for T and P
weather_fc <- function(sel_T, sel_P) function(o, H) {
  fT <- sel_T$f(ts_head(Tc_ts, o), H); fP <- sel_P$f(ts_head(P_ts, o), H)
  list(T = fT$mean, T2 = fT$mean^2 + fT$se^2, P = fP$mean)   # E(T^2) = (E T)^2 + Var(T)
}
xfc_l_tr <- weather_fc(sel_tempr_tr, sel_part_tr)
eval_l <- rolling_eval(eq_l_tr, ts(y, frequency = 52), list(T = Tc, T2 = T2, P = Pn), a_tr,
                       uni_l, origins = T0_l:(n - 1), H = K_l, xfc = xfc_l_tr)

On the training weeks, select_auto() selects ARIMA + Fourier for temperature and ARIMA + Fourier for particulates (ARIMA with Fourier terms is dynamic harmonic regression).

c. One-step-ahead forecasts

The one-step-ahead forecasts of the 52 test weeks are made from the origins t = 456, \dots, 507, every time with the newly observed data. Because T^2_t and P_t enter at lag 0, the ex ante forecasts use the forecast weather and pollution of the coming week, and the ex post forecasts the observed values.

one_l <- onestep_of(eval_l)                               # forecasts at h = 1
plot_fc(one_l$t, one_l$Actual, fc_cols(one_l), "LA mortality: one-step-ahead forecasts",
        "mortality")

acc1_l <- acc_of(one_l)                                   # accuracy of the plotted forecasts
acc1_l

With the observed weather and pollution of the forecast week (ex post), the transfer function reduces the RMSE from 7.3925 to 4.4449 and the MAE from 6.1605 to 3.5773: knowing the weather and the pollution of the current and previous weeks explains a substantial part of the week-to-week variation in mortality that a univariate model must treat as unpredictable. With forecast inputs (ex ante) the RMSE is 4.8424, still clearly smaller than that of the seasonal ARIMA.

d. k-step-ahead forecasts

Forecast path from the end of the training period. The forecasts made at T_0 = 456 (early October 1978) for the following 52 weeks.

fk_l    <- tf_forecast(eq_l_tr, y, list(T = Tc, T2 = T2, P = Pn), a_tr, o = T0_l,
                       H = length(test_l), xfc = xfc_l_tr)               # ex ante
fk_ep_l <- tf_forecast(eq_l_tr, y, list(T = Tc, T2 = T2, P = Pn), a_tr, o = T0_l,
                       H = length(test_l))                               # ex post
fk_ar_l <- as.numeric(forecast(Arima(ts(y[1:T0_l], frequency = 52), model = uni_l),
                               h = length(test_l))$mean)
plot_fc(test_l, y[test_l],
        list(`TF ex ante` = fk_l, `TF ex post` = fk_ep_l, `Univariate ARIMA` = fk_ar_l),
        paste0("LA mortality: k-step-ahead forecasts from T0 = ", T0_l, " (h = 1, ..., ",
               length(test_l), ")"), "mortality")

path_l <- accuracy_tbl(y[test_l], `TF ex ante` = fk_l, `TF ex post` = fk_ep_l,
                       `Univariate ARIMA` = fk_ar_l)
path_l

The forecast path from early October 1978 follows the seasonal course of mortality smoothly, with a winter peak in late November and a summer low, and misses only the sharp peak of late November 1978. Its RMSE over the year (4.90) is close to that of the ex post path (4.62) and clearly smaller than that of the seasonal ARIMA (7.47), which differences at lag 52 and therefore repeats the irregular weekly pattern of the previous year.

Rolling-origin evaluation by horizon. In real time the weather and the pollution of the coming weeks are not known when the forecast is made. From every origin t = 456, \dots, 507, temperature and particulates are therefore forecast up to 13 weeks ahead, and mortality is forecast with these input forecasts (ex ante); the ex post forecasts use the observed inputs. The row h = 1 contains the one-step-ahead forecasts of part c.

plot_rmse(eval_l, "LA mortality: RMSE of the k-step-ahead forecasts by horizon")

eval_l

The transfer function remains clearly more accurate than the seasonal ARIMA at every horizon, also when the inputs have to be forecast: its ex ante RMSE is 24 to 34% smaller (4.84 against 7.39 at h = 1, 4.50 against 6.12 at h = 13). Knowing the future weather and pollution (ex post) would lower the RMSE by another 0.3 to 0.7. The errors hardly grow with the horizon: the AR(2) noise forgets the last shocks within a few weeks, and beyond that the forecasts are driven by the trend and by the annual cycles of mortality, temperature and particulates.

e. Forecasts beyond the end of the data

The model estimated on all 508 weeks (Step 5) forecasts mortality for the next 13 weeks. The forecasting methods for temperature and particulates are now selected on the whole series.

sel_tempr <- select_auto(Tc_ts, H = K_l, methods = auto_methods_seasonal)
sel_part  <- select_auto(P_ts,  H = K_l, methods = auto_methods_seasonal)
bind_rows(temperature = sel_tempr$table, particulates = sel_part$table, .id = "input")

print(c(temperature = sel_tempr$model, particulates = sel_part$model))

                            temperature                            particulates 
"ARIMA(2,0,2) errors + 3 Fourier pairs" "ARIMA(2,1,3) errors + 2 Fourier pairs" 

xfc_l <- weather_fc(sel_tempr, sel_part)
a_full_l <- c(0, a_l)
fc_l <- tf_forecast(eq_l, y, list(T = Tc, T2 = T2, P = Pn), a_full_l, o = n, H = K_l,
                    xfc = xfc_l)
xf_l <- xfc_l(n, K_l)

plot_future(y, fc_l, n_show = 104) +
  labs(title = "LA mortality: forecasts for the next 13 weeks", y = "mortality")

tibble(h = 1:K_l, week = n + 1:K_l, `temperature (F)` = round(xf_l$T + mean(Te), 1),
       particulates = round(xf_l$P, 1), `mortality forecast` = round(fc_l, 2),
       `RMSE at h (test period)` = eval_l$`TF ex ante`[1:K_l])

The values after the first week of October 1979 are not yet known; the last column gives the RMSE of the ex ante forecasts at the same horizon in the test period (part d) as the accuracy to expect. On the whole series, select_auto() selects ARIMA + Fourier for temperature (ARIMA(2,0,2) errors + 3 Fourier pairs) and ARIMA + Fourier for particulates (ARIMA(2,1,3) errors + 2 Fourier pairs). They combine a smooth annual cycle with short-run ARMA dynamics and are clearly better than the seasonal naive forecast, which repeats the noisy values of the same week one year earlier. The forecasts cover October to December. The forecast temperature falls from 70.4 °F to about 65.2 °F in late November, and mortality is forecast to rise from 88.5 to a maximum of 93.3 at about the same time and to decline slightly afterwards. This agrees with the seasonal profile of mortality in these data, which peaks in the second half of November.

Conclusion. The LTF approach identifies a delayed, decaying effect of temperature (b = 1), a contemporaneous decaying effect of squared temperature, and a contemporaneous effect of particulates that persists for several weeks. Together with a linear trend and AR(2) noise, these form an adequate model: all coefficients are significant, the residuals are white, and every Q_0 check passes. The model reproduces the U-shaped temperature–mortality relation and the particulate effect reported by Shumway and Stoffer, and adds their timing. Because temperature and particulates follow annual cycles that can be forecast, most of the gain over the seasonal ARIMA remains when these inputs have to be forecast (ex ante).

Case 5 — Lake Shasta: Climate and Reservoir Inflow

Problem Description

The data set climhyd (package astsa) contains 454 months of climatic variables measured at Lake Shasta, a large reservoir in northern California: air temperature, dew point, cloud cover, wind speed, precipitation and the inflow to the lake. The water-resource problem is to explain and predict the inflow from the weather, which matters for reservoir operation, hydropower generation and flood control. Precipitation reaches the reservoir through surface runoff, groundwater and snowmelt, so its effect is expected to be spread over several months. Temperature acts mainly through snowmelt and evaporation. Following Shumway and Stoffer (Example 7.1), the inflow is analysed on the log scale, precipitation on the square-root scale, and temperature and precipitation are the inputs. The weather is clearly not affected by the inflow, so there is no feedback.

Data

cl <- climhyd                          # from astsa (loaded at the top)
LI <- ts(log(cl$Inflow), frequency = 12)
SP <- ts(sqrt(cl$Precip), frequency = 12)
TP <- ts(cl$Temp, frequency = 12)
tibble(month = seq_along(LI), `log Inflow Y` = as.numeric(LI),
       `sqrt Precipitation` = as.numeric(SP), Temperature = as.numeric(TP)) %>%
  pivot_longer(-month) %>%
  ggplot(aes(month, value)) +
    geom_line(colour = "steelblue") +
    facet_wrap(~ name, ncol = 1, scales = "free_y") +
    labs(x = "month", y = NULL) +
    theme_bw()

All three series are dominated by a strong annual cycle: precipitation falls mainly in the winter months, inflow is high in winter and spring, and temperature peaks in summer. Superimposed on the cycle are wet and dry years.

Seasonal Adjustment and Stationarity

As in Case 3, the deterministic annual cycle is removed from all series by subtracting the monthly means. The cycle would otherwise dominate every lagged correlation.

y <- desea(LI); p <- desea(SP); te <- desea(TP); n <- length(y)
stat_tests(`log Inflow` = LI, `sqrt Precip` = SP, Temp = TP,
           `log Inflow*` = y, `sqrt Precip*` = p, `Temp*` = te)

round(cor(cbind(`log Inflow*` = y, `sqrt Precip*` = p, `Temp*` = te)), 3)

             log Inflow* sqrt Precip*  Temp*
log Inflow*        1.000        0.643 -0.058
sqrt Precip*       0.643        1.000 -0.217
Temp*             -0.058       -0.217  1.000

ADF and PP reject a unit root for all series (p-value \le 0.01). KPSS rejects level stationarity only for the adjusted inflow (p-value \le 0.01), which reflects the persistence of wet and dry periods. Step 2 checks that this persistence does not require differencing. The adjusted inflow is strongly correlated with the adjusted precipitation (0.643). The two inputs are themselves correlated (-0.217: wet months tend to be cooler), which again calls for the LTF approach.

Step 1 — Free-Form Model with a Proxy Noise

k <- 6
XL <- cbind(lagmat(p, k, "P"), lagmat(te, k, "T"))
ff_c <- Arima(y[(k+1):n], order = c(1, 0, 0), xreg = XL)
round(coef(ff_c)[c("ar1", "intercept")], 4)

      ar1 intercept 
   0.4209    0.0000 

w_c <- ff_weights(ff_c, c("P", "T"))
w_c %>% mutate(across(c(nu, se, t), ~ round(.x, 4)))

plot_ff(w_c, "Free-form weights with +/- 2 s.e. (proxy AR(1) noise)")

Lag k	0	1	2	3	4	5	6
\sqrt{P}: \hat\nu_k (t)	0.0535 (25.1)	0.0252 (11.6)	0.0182 (8.3)	0.0130 (5.8)	0.0104 (4.8)	0.0120 (5.5)	0.0043 (2.0)
T: \hat\nu_k (t)	0.0364 (5.1)	-0.0177 (-2.4)	-0.0139 (-1.9)	-0.0279 (-3.8)	-0.0076 (-1.0)	-0.0086 (-1.2)	-0.0060 (-0.8)

Step 2 — Stationarity of the Errors

The proxy coefficient is \hat\phi_1 = 0.42, far from 1, so no differencing is needed.

Step 3 — Orders (b, r, s) from the Weights

• Precipitation: the effect starts immediately (b = 0). The first weight (0.054) is much larger than the second (0.025), after which the weights decay slowly and smoothly (0.025, 0.018, 0.013, 0.010, …). One unpatterned weight followed by a geometric decay means s - r + 1 = 1, so (0,1,1). Physically, the same-month runoff is followed by a slowly draining contribution from groundwater and snowpack.
• Temperature: the effect also starts immediately (b = 0) with a positive weight (0.036), followed by negative weights that die down. A change of sign after lag 0 needs a numerator of order 1, and the decay needs a denominator of order 1, so (0,1,1). A warm month melts snow (more inflow now) and leaves less snow for the following months (less inflow later).

Step 4 — Noise Model

eta_c <- y[(k+1):n] - coef(ff_c)["intercept"] - as.numeric(XL %*% coef(ff_c)[colnames(XL)])
(ggAcf(eta_c, lag.max = 36) + ggtitle("ACF of free-form errors")) |
  (ggPacf(eta_c, lag.max = 36) + ggtitle("PACF of free-form errors"))

sapply(list(AR1 = c(1, 0, 0), AR2 = c(2, 0, 0), ARMA11 = c(1, 0, 1), ARMA21 = c(2, 0, 1)),
       function(o) Arima(eta_c, order = o)$aic)

       AR1        AR2     ARMA11     ARMA21 
 -93.16471 -102.80583 -115.78017 -123.03125 

The ACF of the errors decays very slowly (0.42, 0.31, 0.24, 0.27, 0.27, 0.22, …) and the PACF also tails off (0.42, 0.16, 0.08, 0.15, 0.11, …). This persistent, mixed pattern points to an ARMA(1,1) noise with an AR coefficient close to 1, and ARMA(1,1) has a much smaller AIC (-115.8) than AR(1) (-93.2) or AR(2) (-102.8). ARMA(2,1) is examined in Step 5.

Step 5 — Rational Transfer Function with ARMA Noise

XT_c <- data.frame(Precip = p, Temp = te)
fit_c <- arimax(y, order = c(1, 0, 1), method = "ML", xtransf = XT_c,
                transfer = list(c(1, 1), c(1, 1)), optim.control = list(maxit = 2000))
fit_c

Call:
arimax(x = y, order = c(1, 0, 1), method = "ML", optim.control = list(maxit = 2000), 
    xtransf = XT_c, transfer = list(c(1, 1), c(1, 1)))

Coefficients:
         ar1      ma1  intercept  Precip-AR1  Precip-MA0  Precip-MA1  Temp-AR1  Temp-MA0  Temp-MA1
      0.9097  -0.6933    -0.0041      0.7570      0.0541     -0.0165    0.7513     0.040   -0.0461
s.e.  0.0509   0.0934     0.0336      0.0314      0.0020      0.0032    0.1436     0.007    0.0076

sigma^2 estimated as 0.04543:  log likelihood = 57.26,  aic = -96.51

coeftest(fit_c)

z test of coefficients:

             Estimate Std. Error z value  Pr(>|z|)    
ar1         0.9097090  0.0509401 17.8584 < 2.2e-16 ***
ma1        -0.6932611  0.0934025 -7.4223 1.151e-13 ***
intercept  -0.0041253  0.0335761 -0.1229    0.9022    
Precip-AR1  0.7569647  0.0313915 24.1137 < 2.2e-16 ***
Precip-MA0  0.0541410  0.0020434 26.4951 < 2.2e-16 ***
Precip-MA1 -0.0165275  0.0031867 -5.1864 2.144e-07 ***
Temp-AR1    0.7512519  0.1435614  5.2330 1.668e-07 ***
Temp-MA0    0.0399899  0.0069996  5.7132 1.109e-08 ***
Temp-MA1   -0.0461075  0.0076344 -6.0395 1.546e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

fit_c_t0 <- arimax(y, order = c(1, 0, 1), method = "ML", xtransf = XT_c,
                   transfer = list(c(1, 1), c(0, 1)), optim.control = list(maxit = 2000))
fit_c21 <- arimax(y, order = c(2, 0, 1), method = "ML", xtransf = XT_c,
                  transfer = list(c(1, 1), c(1, 1)), optim.control = list(maxit = 2000))
c(AIC_selected = fit_c$aic, AIC_temperature_r0 = fit_c_t0$aic, AIC_noise_ARMA21 = fit_c21$aic)

      AIC_selected AIC_temperature_r0   AIC_noise_ARMA21 
         -96.51391          -92.19665         -105.42271 

round(coef(fit_c21)[c("ar1", "ar2", "ma1")], 4)

    ar1     ar2     ma1 
 1.2260 -0.2324 -0.9367 

Mod(1 / polyroot(c(1, -coef(fit_c21)[c("ar1", "ar2")])))   # inverse AR roots of ARMA(2,1)

[1] 0.9917498 0.2342953

All transfer function and noise parameters are significant; the smallest test statistic is z = -5.19 for Precip-MA1. The constant is not, as expected for adjusted series.

• Temperature denominator. Dropping it, i.e. (0,0,1) for temperature, increases AIC from -96.5 to -92.2, so it is kept.
• ARMA(2,1) noise. It lowers AIC to -105.4, but its AR polynomial has an inverse root of 0.992, practically a unit root, while its MA coefficient (-0.937) almost cancels it. Near-cancelling factors are a sign of over-parameterisation, and ARMA(1,1) already passes all diagnostics (Step 6), so ARMA(1,1) is kept.

The selected model is

y_t^* = -0.0041 + \frac{0.0541 - 0.0165B}{1 - 0.7570B}\,P_t^* + \frac{0.0400 - 0.0461B}{1 - 0.7513B}\,T_t^* + \frac{1 - 0.6933B}{1 - 0.9097B}\,a_t, \qquad a_t \sim (0,\ 0.0454),

where y^*, P^* and T^* are the adjusted log inflow, square-root precipitation and temperature.

cf_c <- coef(fit_c)
gain_c <- c(Precip = unname(sum(cf_c[c("Precip-MA0", "Precip-MA1")]) / (1 - cf_c["Precip-AR1"])),
            Temp   = unname(sum(cf_c[c("Temp-MA0", "Temp-MA1")]) / (1 - cf_c["Temp-AR1"])))
round(gain_c, 4)

 Precip    Temp 
 0.1548 -0.0246 

implied_c <- bind_rows(
  tibble(input = "P", k = 0:6, nu = nu_weights(cf_c[c("Precip-MA0", "Precip-MA1")], c(1, -cf_c["Precip-AR1"]), b = 0, K = 6)),
  tibble(input = "T", k = 0:6, nu = nu_weights(cf_c[c("Temp-MA0", "Temp-MA1")], c(1, -cf_c["Temp-AR1"]), b = 0, K = 6)))
plot_ff(w_c, "Free-form weights (bars) and weights implied by the rational model (dashed)",
        implied = implied_c)

Interpretation.

• Precipitation. One extra unit of \sqrt{\text{precipitation}} in a month raises log inflow by 0.054 (about 5.6%) in the same month. It then adds 0.0245, 0.0185, 0.0140, \dots in the following months, decaying by about a quarter per month. A sustained increase of one unit raises log inflow by 0.155 in the long run, i.e. the inflow by about 17%.
• Temperature. A month that is 1 °C warmer than usual raises log inflow by 0.040 immediately, through snowmelt. It lowers inflow slightly in the following months (-0.016, -0.012, \dots), because the snowpack has already melted. The long-run gain is small (-0.025): temperature mainly shifts the timing of the inflow, whereas precipitation determines its amount. This agrees with Shumway and Stoffer, who find that precipitation is the dominant input.

Step 6 — Diagnostic Checking

a_c <- residuals(fit_c)
Box.test(a_c, lag = 12, type = "Ljung-Box", fitdf = 2)

    Box-Ljung test

data:  a_c
X-squared = 15.453, df = 10, p-value = 0.1164

Box.test(a_c, lag = 24, type = "Ljung-Box", fitdf = 2)

    Box-Ljung test

data:  a_c
X-squared = 25.875, df = 22, p-value = 0.2571

jarque.bera.test(na.omit(a_c))

    Jarque Bera Test

data:  na.omit(a_c)
X-squared = 45.01, df = 2, p-value = 1.684e-10

inputs_c <- list(Precip = p, Temp = te)
q0_c <- imap_dfr(inputs_c, function(z, nm) {
  fa <- auto.arima(z, seasonal = FALSE)
  Q0_test(a_c, as.numeric(residuals(fa)), K = 12, M = 3) %>%
    mutate(input = nm, `input model` = paste0("ARIMA(", paste(arimaorder(fa), collapse = ","), ")"),
           .before = 1)
})
q0_c

ggAcf(na.omit(a_c), lag.max = 36) + ggtitle("ACF of residuals")

The first residual is NA because both numerators contain B.

• Noise model: Q(12) = 15.45 (df = 10, p-value = 0.116) and Q(24) = 25.88 (df = 22, p-value = 0.257), so the residuals are white.
• Transfer functions: Q_0(12) = 10.4 for precipitation (p-value = 0.402) and 9.99 for temperature (p-value = 0.441), both with df = 12 + 1 - 3 = 10, so the residuals are uncorrelated with both prewhitened inputs.
• Normality: rejected (X^2 = 45.0), mainly because of a few months with exceptional floods.

Step 7 — Explicit Equation

The explicit equation for the adjusted series is obtained by multiplying the model by (1-0.7570B)(1-0.7513B)(1-0.9097B), the two denominators and the AR operator of the noise. The operators and the constant:

eq_c <- tf_explicit_multi(list(
  list(name = "P", w = cf_c[c("Precip-MA0", "Precip-MA1")], delta = c(1, -cf_c["Precip-AR1"]), b = 0),
  list(name = "T", w = cf_c[c("Temp-MA0", "Temp-MA1")],     delta = c(1, -cf_c["Temp-AR1"]),   b = 0)),
  phi = c(1, -cf_c["ar1"]), theta = c(1, cf_c["ma1"]), mu = cf_c["intercept"])
show_operators(eq_c, sym = c(Y = "y^*", P = "P^*", T = "T^*"))

\begin{aligned} A(B) &= \delta_{P^*}(B)\,\delta_{T^*}(B)\,\phi(B) && \text{acting on } y^*_{t} \\ &= \left(1 - 0.7570B\right)\left(1 - 0.7513B\right)\left(1 - 0.9097B\right) \\ &= 1 - 2.4179B + 1.9407B^{2} - 0.5173B^{3} \\[4pt] C_{P^*}(B) &= \omega_{P^*}(B)\,\delta_{T^*}(B)\,\phi(B) && \text{acting on } P^*_{t} \\ &= \left(0.05414 - 0.01653B\right)\left(1 - 0.7513B\right)\left(1 - 0.9097B\right) \\ &= 0.05414 - 0.1065B + 0.06445B^{2} - 0.01130B^{3} \\[4pt] C_{T^*}(B) &= \omega_{T^*}(B)\,\delta_{P^*}(B)\,\phi(B) && \text{acting on } T^*_{t} \\ &= \left(0.03999 - 0.04611B\right)\left(1 - 0.7570B\right)\left(1 - 0.9097B\right) \\ &= 0.03999 - 0.1128B + 0.1044B^{2} - 0.03175B^{3} \\[4pt] D(B) &= \delta_{P^*}(B)\,\delta_{T^*}(B)\,\theta(B) && \text{acting on } a_{t} \\ &= \left(1 - 0.7570B\right)\left(1 - 0.7513B\right)\left(1 - 0.6933B\right) \\ &= 1 - 2.2015B + 1.6143B^{2} - 0.3942B^{3} \\[4pt] c &= \hat\mu\,\delta_{P^*}(1)\,\delta_{T^*}(1)\,\phi(1) \\ &= -0.00412526 \times 0.243035 \times 0.248748 \times 0.0902910 = -0.00002252 \end{aligned}

In explicit form (first equation below), the adjusted log inflow depends on its last three values, on precipitation and temperature at lags 0–3, and on three past shocks.

To go back to the original series, write L_t = \ln \mathrm{Inflow}_t, P_t = \sqrt{\mathrm{Precip}_t} and T_t = \mathrm{Temp}_t, and substitute y^*_t = L_t - \bar L_{m(t)}, P^*_t = P_t - \bar P_{m(t)} and T^*_t = T_t - \bar T_{m(t)}. As in Case 3, the coefficients stay the same and the constant becomes a month-specific constant,

c_{m(t)} = \sum_{m=1}^{12} c_m D_{m,t}, \qquad c_m = c + \bar L_m - \sum_{i} \pi_i\, \bar L_{m-i} - \sum_{k} \omega^*_{P,k}\, \bar P_{m-k} - \sum_{k} \omega^*_{T,k}\, \bar T_{m-k},

where m is the position in the 12-month cycle (position 1 is the first month of the record) and D_{m,t} = 1 if t is at position m, 0 otherwise. The display shows the three steps: (1) the explicit equation of the adjusted series, (2) the definitions substituted, and (3) the explicit equation of the original variables, with the transformations written out:

mm_LI <- tapply(as.numeric(LI), cycle(LI), mean)     # monthly means of log inflow,
mm_SP <- tapply(as.numeric(SP), cycle(SP), mean)     # of sqrt precipitation
mm_TP <- tapply(as.numeric(TP), cycle(TP), mean)     # and of temperature
c_c <- monthly_const(eq_c, mm_LI, list(P = mm_SP, T = mm_TP))
sym_c <- c(Y = "\\ln \\mathrm{Inflow}", P = "\\sqrt{\\mathrm{Precip}_{#}}",
           T = "\\mathrm{Temp}")                     # symbols of the original variables
show_chain(
  list(e = eq_c, sym = c(Y = "y^*", P = "P^*", T = "T^*")),             # 1. adjusted series
  list(e = eq_c, per_line = 3,                                           # 2. substituted
       sym = c(Y = "\\left(L_{#} - \\bar L_{m(#)}\\right)",
               P = "\\left(P_{#} - \\bar P_{m(#)}\\right)",
               T = "\\left(T_{#} - \\bar T_{m(#)}\\right)")),
  list(e = eq_c, sym = sym_c, const = dummy_terms(c_c), per_line = 4))  # 3. original variables

\begin{aligned} y^*_{t} ={}& {-0.00002252} + 2.4179\,y^*_{t-1} - 1.9407\,y^*_{t-2} + 0.5173\,y^*_{t-3} \\ &{}+ 0.05414\,P^*_{t} - 0.1065\,P^*_{t-1} + 0.06445\,P^*_{t-2} - 0.01130\,P^*_{t-3} \\ &{}+ 0.03999\,T^*_{t} - 0.1128\,T^*_{t-1} + 0.1044\,T^*_{t-2} - 0.03175\,T^*_{t-3} \\ &{}+ a_{t} - 2.2015\,a_{t-1} + 1.6143\,a_{t-2} - 0.3942\,a_{t-3} \\[8pt] \left(L_{t} - \bar L_{m(t)}\right) ={}& {-0.00002252} + 2.4179\,\left(L_{t-1} - \bar L_{m(t-1)}\right) - 1.9407\,\left(L_{t-2} - \bar L_{m(t-2)}\right) \\ &{}+ 0.5173\,\left(L_{t-3} - \bar L_{m(t-3)}\right) \\ &{}+ 0.05414\,\left(P_{t} - \bar P_{m(t)}\right) - 0.1065\,\left(P_{t-1} - \bar P_{m(t-1)}\right) + 0.06445\,\left(P_{t-2} - \bar P_{m(t-2)}\right) \\ &{}- 0.01130\,\left(P_{t-3} - \bar P_{m(t-3)}\right) \\ &{}+ 0.03999\,\left(T_{t} - \bar T_{m(t)}\right) - 0.1128\,\left(T_{t-1} - \bar T_{m(t-1)}\right) + 0.1044\,\left(T_{t-2} - \bar T_{m(t-2)}\right) \\ &{}- 0.03175\,\left(T_{t-3} - \bar T_{m(t-3)}\right) \\ &{}+ a_{t} - 2.2015\,a_{t-1} + 1.6143\,a_{t-2} \\ &{}- 0.3942\,a_{t-3} \\[8pt] \ln \mathrm{Inflow}_{t} ={}& {-0.1811}\,D_{1,t} + 0.2193\,D_{2,t} - 0.1056\,D_{3,t} + 0.03862\,D_{4,t} - 0.003642\,D_{5,t} - 0.05818\,D_{6,t} \\ &{}+ 0.3601\,D_{7,t} + 0.05300\,D_{8,t} - 0.06760\,D_{9,t} + 0.01282\,D_{10,t} + 0.01876\,D_{11,t} - 0.007467\,D_{12,t} \\ &{}+ 2.4179\,\ln \mathrm{Inflow}_{t-1} - 1.9407\,\ln \mathrm{Inflow}_{t-2} + 0.5173\,\ln \mathrm{Inflow}_{t-3} \\ &{}+ 0.05414\,\sqrt{\mathrm{Precip}_{t}} - 0.1065\,\sqrt{\mathrm{Precip}_{t-1}} + 0.06445\,\sqrt{\mathrm{Precip}_{t-2}} - 0.01130\,\sqrt{\mathrm{Precip}_{t-3}} \\ &{}+ 0.03999\,\mathrm{Temp}_{t} - 0.1128\,\mathrm{Temp}_{t-1} + 0.1044\,\mathrm{Temp}_{t-2} - 0.03175\,\mathrm{Temp}_{t-3} \\ &{}+ a_{t} - 2.2015\,a_{t-1} + 1.6143\,a_{t-2} - 0.3942\,a_{t-3} \end{aligned}

The inflow itself follows by exponentiating, \mathrm{Inflow}_t = \exp(\ln \mathrm{Inflow}_t). The monthly means and the monthly constants c_m (the coefficients of D_{1,t}, \dots, D_{12,t}):

tibble(`month position` = 1:12, `mean log inflow` = round(as.numeric(mm_LI), 4),
       `mean sqrt precip` = round(as.numeric(mm_SP), 4),
       `mean temp` = round(as.numeric(mm_TP), 3), `c[m]` = round(c_c, 4))

As in Case 3, the monthly constants are small compared with the seasonal swing of the log inflow (its monthly means range from 4.61 to 5.89), because the three lags of the log inflow already carry most of the seasonal level.

Step 8 — Forecasting

The forecasts are computed from the explicit equation of Step 7. In the computations the adjusted series are forecast and the monthly means of the log inflow are added back afterwards, which is the same as using the equation of the original variables with the monthly constants c_m. Both inputs enter at lag 0, so even the one-step-ahead forecast needs the precipitation and temperature of the coming month. The adjusted inputs are forecast with the methods chosen by select_auto(). Exponentiating the forecasts of the log inflow gives forecasts of the inflow itself (approximately the median, not the mean, of its forecast distribution).

a. Forecast equations

For the origin t, the one-step-ahead forecast and the k-step-ahead forecast of the log inflow L_t = \ln \mathrm{Inflow}_t are, with P_t = \sqrt{\mathrm{Precip}_t} and T_t = \mathrm{Temp}_t,

show_forecast_eq(eq_c, sym = c(Y = "L"), const1 = "c_{m(t+1)}", constk = "c_{m(t+k)}")

\begin{aligned} \hat{L}_t(1) ={}& c_{m(t+1)} + 2.4179\,L_{t} - 1.9407\,L_{t-1} + 0.5173\,L_{t-2} \\ &{}+ 0.05414\,P_{t+1} - 0.1065\,P_{t} + 0.06445\,P_{t-1} - 0.01130\,P_{t-2} \\ &{}+ 0.03999\,T_{t+1} - 0.1128\,T_{t} + 0.1044\,T_{t-1} - 0.03175\,T_{t-2} \\ &{}- 2.2015\,a_{t} + 1.6143\,a_{t-1} - 0.3942\,a_{t-2} \\[8pt] \hat{L}_t(k) ={}& c_{m(t+k)} + 2.4179\,\hat{L}_t(k-1) - 1.9407\,\hat{L}_t(k-2) + 0.5173\,\hat{L}_t(k-3) \\ &{}+ 0.05414\,\hat{P}_t(k) - 0.1065\,\hat{P}_t(k-1) + 0.06445\,\hat{P}_t(k-2) - 0.01130\,\hat{P}_t(k-3) \\ &{}+ 0.03999\,\hat{T}_t(k) - 0.1128\,\hat{T}_t(k-1) + 0.1044\,\hat{T}_t(k-2) - 0.03175\,\hat{T}_t(k-3) \\ &{}- 2.2015\,\hat{a}_t(k-1) + 1.6143\,\hat{a}_t(k-2) - 0.3942\,\hat{a}_t(k-3) \end{aligned}

where c_{m(t+k)} is the monthly constant of Step 7 for the position of t + k in the 12-month cycle, \hat Z_t(m) = Z_{t+m} for m \le 0 (the observed values of L, P, T and the residuals a) and \hat a_t(m) = 0 for m > 0. In the one-step-ahead equation, P_{t+1} and T_{t+1} belong to the forecast month: in real time they are replaced by their forecasts \hat P_t(1) and \hat T_t(1) (ex ante), and the ex post version uses the observed weather. The inflow forecast is \widehat{\mathrm{Inflow}}_t(k) = \exp\big(\hat L_t(k)\big).

b. Training and test periods

The model is re-estimated on all but the last 48 months, with the monthly means also computed from the training months. The benchmark is the univariate model chosen by auto.arima() for the adjusted log inflow. The adjusted precipitation and temperature are forecast with the methods selected on the training months only, and from every origin t = 406, \dots, 453 the log inflow is forecast up to K = 12 months ahead, as far as the test period reaches. orig_c() adds the monthly means of the training months back, so that all forecasts are shown and evaluated for the log inflow itself; this does not change the forecast errors.

K_c  <- 12                                                # forecast horizons h = 1, ..., K
T0_c <- n - 48                                            # training period: 406 months
yh <- desea(LI, T0_c); ph <- desea(SP, T0_c); teh <- desea(TP, T0_c)
fit_c_tr <- arimax(yh[1:T0_c], order = c(1, 0, 1), method = "ML",
                   xtransf = data.frame(Precip = ph[1:T0_c], Temp = teh[1:T0_c]),
                   transfer = list(c(1, 1), c(1, 1)), optim.control = list(maxit = 2000))
ct <- coef(fit_c_tr)
eq_c_tr <- tf_explicit_multi(list(
  list(name = "P", w = ct[c("Precip-MA0", "Precip-MA1")], delta = c(1, -ct["Precip-AR1"]), b = 0),
  list(name = "T", w = ct[c("Temp-MA0", "Temp-MA1")],     delta = c(1, -ct["Temp-AR1"]),   b = 0)),
  phi = c(1, -ct["ar1"]), theta = c(1, ct["ma1"]), mu = ct["intercept"])
test_c <- (T0_c + 1):n
a_tr <- c(residuals(fit_c_tr), rep(0, length(test_c))); a_tr[is.na(a_tr)] <- 0
f_tf_c <- onestep_explicit_multi(eq_c_tr, yh, list(P = ph, T = teh), a_tr, test_c)
a_tr[test_c] <- yh[test_c] - f_tf_c                       # residuals of the test period
uni_c <- auto.arima(yh[1:T0_c])
uni_c

Series: yh[1:T0_c] 
ARIMA(1,0,0) with zero mean 

Coefficients:
         ar1
      0.5933
s.e.  0.0400

sigma^2 = 0.1144:  log likelihood = -135.74
AIC=275.47   AICc=275.5   BIC=283.48

sel_precip_tr <- select_auto(ts_head(ph, T0_c),  H = K_c)        # training months only
sel_temp_tr   <- select_auto(ts_head(teh, T0_c), H = K_c)
bind_rows(`sqrt precipitation*` = sel_precip_tr$table, `temperature*` = sel_temp_tr$table,
          .id = "input")

print(c(`sqrt precipitation*` = sel_precip_tr$model, `temperature*` = sel_temp_tr$model))

sqrt precipitation*        temperature* 
             "Mean"              "Mean" 

xfc_c_tr <- function(o, H) list(P = sel_precip_tr$f(ts_head(ph, o), H)$mean,
                                T = sel_temp_tr$f(ts_head(teh, o), H)$mean)
eval_c <- rolling_eval(eq_c_tr, yh, list(P = ph, T = teh), a_tr, uni_c,
                       origins = T0_c:(n - 1), H = K_c, xfc = xfc_c_tr)
mm_LI_tr <- tapply(as.numeric(LI)[1:T0_c], cycle(LI)[1:T0_c], mean)  # training months
orig_c <- function(v, t) v + as.numeric(mm_LI_tr[cycle(LI)[t]])     # back to log inflow

For the adjusted log inflow alone, auto.arima() selects an AR(1). For the weather, select_auto() selects Mean for precipitation and Mean for temperature on the training months.

c. One-step-ahead forecasts

The one-step-ahead forecasts of the 48 test months are made from the origins t = 406, \dots, 453, every time with the newly observed data. Because both inputs enter at lag 0, the ex ante forecasts use the forecast weather of the coming month, and the ex post forecasts the observed weather.

one_c <- onestep_of(eval_c, orig_c)                       # forecasts at h = 1, log inflow
plot_fc(one_c$t, one_c$Actual, fc_cols(one_c), "Lake Shasta: one-step-ahead forecasts",
        "log inflow")

acc1_c <- acc_of(one_c)                                   # accuracy of the plotted forecasts
acc1_c

With the observed weather of the month (ex post), the transfer function lowers the RMSE from 0.4079 to 0.2565 and the MAE from 0.3162 to 0.1886. The weather of the current month explains much of the deviation of the inflow from its seasonal pattern, especially in the wet months, where the univariate model reacts only after the flood has started. With forecast weather (ex ante) nothing of this remains: the RMSE is 0.4131, practically the same as that of the univariate model, because the deviations of precipitation and temperature from their monthly averages can hardly be forecast.

d. k-step-ahead forecasts

Forecast path from the end of the training period. The forecasts made at T_0 = 406 for the 48 test months.

fk_c    <- tf_forecast(eq_c_tr, yh, list(P = ph, T = teh), a_tr, o = T0_c,
                       H = length(test_c), xfc = xfc_c_tr)               # ex ante
fk_ep_c <- tf_forecast(eq_c_tr, yh, list(P = ph, T = teh), a_tr, o = T0_c,
                       H = length(test_c))                               # ex post
fk_ar_c <- as.numeric(forecast(Arima(yh[1:T0_c], model = uni_c), h = length(test_c))$mean)
path_fc_c <- list(`TF ex ante` = orig_c(fk_c, test_c), `TF ex post` = orig_c(fk_ep_c, test_c),
                  `Univariate ARIMA` = orig_c(fk_ar_c, test_c))       # log inflow
plot_fc(test_c, as.numeric(LI)[test_c], path_fc_c,
        paste0("Lake Shasta: k-step-ahead forecasts from T0 = ", T0_c, " (h = 1, ..., ",
               length(test_c), ")"), "log inflow")

path_c <- do.call(accuracy_tbl, c(list(as.numeric(LI)[test_c]), path_fc_c))
path_c

Both forecast paths return to the seasonal level of the log inflow (its monthly means) within a few months, whereas the log inflow stayed below its seasonal level during most of the test period. Neither model can anticipate such a long period of low inflow, and the two are about equally accurate (RMSE 0.541 and 0.553). With the observed weather the path would be much better (0.335).

Rolling-origin evaluation by horizon. From every origin t = 406, \dots, 453, the adjusted precipitation and temperature are forecast up to 12 months ahead, and the log inflow is forecast with these input forecasts (ex ante); the ex post forecasts use the observed weather. The row h = 1 contains the one-step-ahead forecasts of part c.

plot_rmse(eval_c, "Lake Shasta: RMSE of the k-step-ahead forecasts by horizon")

eval_c

This case shows the difference between ex post and ex ante forecasts most clearly. With the observed weather (ex post), the RMSE of the transfer function is 37 to 45% smaller than that of the univariate model, at every horizon. With forecast weather (ex ante), almost nothing of this remains. The methods selected on the training months (Mean for precipitation and Mean for temperature) forecast hardly any deviation from the monthly averages, and the transfer function is practically as accurate as the univariate model (RMSE 0.413 against 0.408 at h = 1, 0.570 against 0.575 at h = 12). The transfer function explains the inflow well once the weather is known, but in real time it forecasts the inflow no better than the univariate model.

e. Forecasts beyond the end of the data

The model estimated on all 454 months (Step 5) forecasts the next 12 months. The forecasting methods for the adjusted inputs are now selected on the whole series.

sel_precip <- select_auto(p,  H = K_c)                    # sqrt precipitation*
sel_temp   <- select_auto(te, H = K_c)                    # temperature*
bind_rows(`sqrt precipitation*` = sel_precip$table, `temperature*` = sel_temp$table,
          .id = "input")

print(c(`sqrt precipitation*` = sel_precip$model, `temperature*` = sel_temp$model))

sqrt precipitation*        temperature* 
     "ARIMA(3,0,2)"              "Mean" 

xfc_c <- function(o, H) list(P = sel_precip$f(ts_head(p, o), H)$mean,
                             T = sel_temp$f(ts_head(te, o), H)$mean)
a_full_c <- as.numeric(a_c); a_full_c[is.na(a_full_c)] <- 0
f_adj_c <- tf_forecast(eq_c, y, list(P = p, T = te), a_full_c, o = n, H = K_c, xfc = xfc_c)
m_next  <- ((cycle(LI)[n] + 1:K_c - 1) %% 12) + 1
fc_c    <- f_adj_c + as.numeric(mm_LI[m_next])            # log inflow forecasts
xf_c    <- xfc_c(n, K_c)

plot_future(LI, fc_c, n_show = 60) +
  labs(title = "Lake Shasta: forecasts of the log inflow for the next 12 months",
       y = "log inflow")

tibble(h = 1:K_c, `month position` = m_next, `P* forecast` = round(xf_c$P, 2),
       `T* forecast` = round(xf_c$T, 2), `adjusted forecast` = round(f_adj_c, 3),
       `log inflow forecast` = round(fc_c, 3), `inflow forecast` = round(exp(fc_c), 1),
       `RMSE at h (test period)` = eval_c$`TF ex ante`[1:K_c])

The values after the end of the record are not yet known; the last column gives the RMSE of the ex ante forecasts of the log inflow at the same horizon in the test period (part d) as the accuracy to expect. For the adjusted precipitation, ARIMA(3,0,2) has the smallest cross-validation error, but only just: the mean of the series, i.e. a forecast equal to the monthly average, is almost as accurate. For the adjusted temperature the selected method is Mean. The deviations of the weather from its monthly average can thus hardly be forecast from their own past, and their forecasts are close to zero. The forecasts of the adjusted log inflow are small and negative (between -0.17 and -0.04): the inflow is expected to stay slightly below its seasonal level, mainly because it was below it in the last months and is persistent. The inflow forecasts follow the seasonal cycle, from about 152 in the first month to a maximum of about 323 in month 2 of the cycle and a minimum of about 91 in month 8.

Conclusion. The inflow to Lake Shasta responds immediately to precipitation, with a long, geometrically decaying tail (0,1,1). Temperature has a positive immediate effect and a negative delayed effect (0,1,1) that nearly cancel in the long run. With ARMA(1,1) noise the model is adequate: all coefficients are significant, the residuals are white and uncorrelated with both prewhitened inputs. The LTF approach made the identification straightforward despite the correlation between the two inputs. For forecasting, the model is valuable when the weather is known (ex post); because both inputs enter at lag 0 and the weather anomalies are hardly predictable, it forecasts the inflow in real time hardly better than a univariate model.

Summary

The figure brings the five cases together. In the third column, the fitted rational transfer functions (blue) follow the identified weights (grey) closely where the response is concentrated, with only two or three parameters per input; at long lags the CCF-based estimates are noisy (for instance in Case 3), and the rational form smooths them. The residual ACFs stay inside the bands except for a few lags, in line with the Ljung–Box tests: only Case 3 rejects white noise (see its diagnostics). The explicit equations in the last column are the ones used for all forecasts.

The first table collects, for every case and input, the identified orders, the noise model, the steady-state gain, and the p-values of the Ljung–Box test (K = 12) and of the Q_0 test (K = 12).

The next two tables summarise the forecast comparisons on the test periods. The first gives the smallest delay, the largest horizon K and the forecasting method selected for the inputs on the training data. The second gives the RMSE of the transfer function with forecast inputs (ex ante), of the transfer function with observed inputs (ex post) and of the univariate ARIMA, from the rolling-origin evaluation at h = 1 and at h = K, and for the single forecast path from the end of the training period.

Forecast plots of all five cases. Every case has the same forecasting step (Step 11 of Cases 1–3, Step 8 of Cases 4–5) with four forecast plots, each followed by the accuracy of the forecasts it shows: the one-step-ahead forecasts of the test period, the k-step-ahead forecast path from the end of the training period, the RMSE of the rolling-origin forecasts by horizon, and the forecasts beyond the end of the data. The figure below collects the last kind, each in the original variable of the output (log inflow for Case 5).

Explicit equations in the original variables. For reference, the explicit equations of the five models, as derived in the cases. They are written in the original variables: the differencing of Case 1, the seasonal adjustment of Cases 3 and 5 (monthly constants with the dummies D_{m,t}) and the centring of the temperature in Case 4 are undone. The shocks satisfy a_t \sim \text{WN}(0, \sigma_a^2) with the residual variances of the fitted models.

Case 1, Series M: sales Y_t and leading indicator X_t (\hat\sigma_a^2 = 0.04741).

\begin{aligned} Y_{t} ={}& 0.008351 + 1.7264\,Y_{t-1} - 0.7264\,Y_{t-2} \\ &{}+ 4.6942\,X_{t-3} - 4.6942\,X_{t-4} \\ &{}+ a_{t} - 1.3138\,a_{t-1} + 0.4267\,a_{t-2} \end{aligned}

Case 2, Series J: CO_2 (%) Y_t and coded gas rate X_t (\hat\sigma_a^2 = 0.05710).

\begin{aligned} Y_{t} ={}& 2.4464 + 2.0763\,Y_{t-1} - 1.4674\,Y_{t-2} + 0.3453\,Y_{t-3} \\ &{}- 0.5310\,X_{t-3} + 0.4308\,X_{t-4} - 0.2714\,X_{t-5} + 0.5520\,X_{t-6} - 0.3257\,X_{t-7} \\ &{}+ a_{t} - 0.5490\,a_{t-1} \end{aligned}

Case 3, SOI and recruitment: D_{m,t} for the calendar months, m = 1 for January (\hat\sigma_a^2 = 46.9563).

\begin{aligned} \mathrm{Rec}_{t} ={}& 1.3310\,D_{1,t} + 3.0330\,D_{2,t} + 4.6402\,D_{3,t} + 1.0426\,D_{4,t} + 1.4766\,D_{5,t} + 0.9059\,D_{6,t} \\ &{}+ 1.3209\,D_{7,t} + 2.5175\,D_{8,t} + 4.2346\,D_{9,t} + 4.6183\,D_{10,t} - 2.1622\,D_{11,t} + 2.9254\,D_{12,t} \\ &{}+ 2.0584\,\mathrm{Rec}_{t-1} - 1.4183\,\mathrm{Rec}_{t-2} + 0.3292\,\mathrm{Rec}_{t-3} \\ &{}- 20.2993\,\mathrm{SOI}_{t-5} + 25.5405\,\mathrm{SOI}_{t-6} - 8.3521\,\mathrm{SOI}_{t-7} \\ &{}+ a_{t} - 0.8002\,a_{t-1} \end{aligned}

Case 4, LA mortality: t = 1 in the first week of 1970 (\hat\sigma_a^2 = 22.7801).

\begin{aligned} \mathrm{cmort}_{t} ={}& 3.1474 - 0.001040\,t + 1.8416\,\mathrm{cmort}_{t-1} \\ &{}- 0.8945\,\mathrm{cmort}_{t-2} - 0.1698\,\mathrm{cmort}_{t-3} + 0.2261\,\mathrm{cmort}_{t-4} \\ &{}- 0.04228\,\mathrm{cmort}_{t-5} \\ &{}- 0.2172\,\left(\mathrm{tempr}_{t-1} - 74.2604\right) + 0.2557\,\left(\mathrm{tempr}_{t-2} - 74.2604\right) - 0.02442\,\left(\mathrm{tempr}_{t-3} - 74.2604\right) \\ &{}- 0.05309\,\left(\mathrm{tempr}_{t-4} - 74.2604\right) + 0.01382\,\left(\mathrm{tempr}_{t-5} - 74.2604\right) \\ &{}+ 0.01759\,\left(\mathrm{tempr}_{t} - 74.2604\right)^{2} - 0.02630\,\left(\mathrm{tempr}_{t-1} - 74.2604\right)^{2} + 0.006620\,\left(\mathrm{tempr}_{t-2} - 74.2604\right)^{2} \\ &{}+ 0.005280\,\left(\mathrm{tempr}_{t-3} - 74.2604\right)^{2} - 0.002146\,\left(\mathrm{tempr}_{t-4} - 74.2604\right)^{2} \\ &{}+ 0.1144\,\mathrm{part}_{t} - 0.1463\,\mathrm{part}_{t-1} + 0.01998\,\mathrm{part}_{t-2} \\ &{}+ 0.03067\,\mathrm{part}_{t-3} - 0.008592\,\mathrm{part}_{t-4} \\ &{}+ a_{t} - 1.5741\,a_{t-1} + 0.7995\,a_{t-2} \\ &{}- 0.1297\,a_{t-3} \end{aligned}

Case 5, Lake Shasta: D_{m,t} for the position in the 12-month cycle (\hat\sigma_a^2 = 0.04543).

\begin{aligned} \ln \mathrm{Inflow}_{t} ={}& {-0.1811}\,D_{1,t} + 0.2193\,D_{2,t} - 0.1056\,D_{3,t} + 0.03862\,D_{4,t} - 0.003642\,D_{5,t} - 0.05818\,D_{6,t} \\ &{}+ 0.3601\,D_{7,t} + 0.05300\,D_{8,t} - 0.06760\,D_{9,t} + 0.01282\,D_{10,t} + 0.01876\,D_{11,t} - 0.007467\,D_{12,t} \\ &{}+ 2.4179\,\ln \mathrm{Inflow}_{t-1} - 1.9407\,\ln \mathrm{Inflow}_{t-2} + 0.5173\,\ln \mathrm{Inflow}_{t-3} \\ &{}+ 0.05414\,\sqrt{\mathrm{Precip}_{t}} - 0.1065\,\sqrt{\mathrm{Precip}_{t-1}} + 0.06445\,\sqrt{\mathrm{Precip}_{t-2}} - 0.01130\,\sqrt{\mathrm{Precip}_{t-3}} \\ &{}+ 0.03999\,\mathrm{Temp}_{t} - 0.1128\,\mathrm{Temp}_{t-1} + 0.1044\,\mathrm{Temp}_{t-2} - 0.03175\,\mathrm{Temp}_{t-3} \\ &{}+ a_{t} - 2.2015\,a_{t-1} + 1.6143\,a_{t-2} - 0.3942\,a_{t-3} \end{aligned}

References

• Box, G. E. P., Jenkins, G. M., Reinsel, G. C. and Ljung, G. M. (2015). Time Series Analysis: Forecasting and Control, 5th ed. Wiley. (Series J and Series M.)
• Liu, L.-M. and Hanssens, D. M. (1982). Identification of multiple-input transfer function models. Communications in Statistics – Theory and Methods, 11(3), 297–314.
• Pankratz, A. (1991). Forecasting with Dynamic Regression Models. Wiley.
• Shumway, R. H., Azari, A. S. and Pawitan, Y. (1988). Modeling mortality fluctuations in Los Angeles as functions of pollution and weather effects. Environmental Research, 45, 224–241.
• Shumway, R. H. and Stoffer, D. S. (2017). Time Series Analysis and Its Applications: With R Examples, 4th ed. Springer. (SOI and Recruitment, LA pollution, Lake Shasta.)
• Wei, W. W. S. (2006). Time Series Analysis: Univariate and Multivariate Methods, 2nd ed. Pearson.