Module 7: Forecast Evaluation and Combinations

Econ 6376 — Applied Time Series Econometrics

How to Use These Notes

This chapter has one main reading path and three optional layers.

Core material is the main text. It contains the concepts, notation, interpretations, and applied skills expected of everyone.
Deeper Dive sections explain why a result works or develop it more fully. They are useful, but they can be skipped on a first reading.
Technical Note sections state qualifications or formal details that matter for precise reasoning.
Looking Ahead sections introduce an idea that will be taught formally in a later module.

If you missed lecture, read the main text, run the Core code, and complete the Core Practice problems. Then return to the optional sections that address your questions or interests.

Core Learning Objectives

By the end of this module, you should be able to:

Describe the PS3 one-step design — train $1{:}T$ , weights window $(T{+}1){:}(T{+}12)$ , holdout $(T{+}13){:}(T{+}24)$ with $T = N_{\text{total}} - 24$ — and explain why the information set rolls while the parameters do not. Produce one-step tracks with Arima(y, model = fit) + fitted() and the three $k$ -period average forecasts (locf, avg_k, meanf), with correct alignment.
Name the three evaluation schemes — fixed-parameter one-step (the PS3 design), expanding window re-estimated, fixed-width rolling window — and say what each holds fixed and what each costs.
Compute and interpret RMSE, MAE, MAPE, MASE, and Theil’s $U_2$ , and state precisely what question each relative metric answers ( $U_2$ : did the model beat the naive forecast in this holdout? MASE: how large is the holdout error relative to the naive forecast’s typical training error?).
State the Diebold-Mariano test — loss differential $d_t = g(e_{1,t}) - g(e_{2,t})$ , $H_0\!: \mathbb{E}[d_t] = 0$ , $\text{DM} = \bar d/\hat\sigma_{\bar d}$ with a standard error that allows for serial correlation in $d_t$ (at h = 1, dm.test() instead assumes serially uncorrelated $d_t$ ), $N(0,1)$ under the null — run dm.test(e1, e2, h = 1, power = 2) on out-of-sample errors, read the output, and explain why a short evaluation window often fails to reject.
Write the two-forecast combined error variance, explain diversification with it, and derive that with zero error covariance the optimal weight is inversely proportional to each forecast’s error variance.
Write down and compute, inline, the four combination methods — simple average, inverse-MSE weights, Granger-Ramanathan (unconstrained least squares with an intercept), and constrained least squares (no intercept, $w_m \ge 0$ , $\sum_m w_m = 1$ ) — and state what each estimates and what each restricts.
Run the assumption-break simulation and explain the forecast combination puzzle as a trade-off between the population gain from optimal weights and the estimation error in those weights — including when equal weights are not hard to beat.
Execute the full PS3 workflow on UNRATE and interpret it honestly: six tracks, a metrics table, DM on the Module 6 question, four combinations, holdout RMSE, and DM of each estimated-weight combination against the simple average.
State the forward pointer: every forecast in this module used only the series’ own past; Module 8 asks whether another variable can help.

7.1 Where We Are on the Mixing Board

The Module 6 Handoff

Module 6 ended on an unresolved note. We fit the diagnostic-clean SARIMA from Module 5 — ARIMA $(1,1,2)(1,0,1)_{12}$ on UNRATE — forecast 2018–2019, and lost to the naive random walk, “carry the last value forward,” by a few hundredths. Module 6’s last sentence was a question for a statistician: does that margin mean anything?

This module answers it, and then asks a better question. The arc has three parts:

How do I score several forecasts? Two relative metrics, Theil’s $U_2$ and MASE.
Is the difference between two forecasts real? The Diebold-Mariano test.
Should I even pick one? Forecast combinations — and the surprising reason the simplest combination is so hard to beat.

On the mixing board, nothing new switches on. Module 6 used one fitted master equation to project forward; Module 7 puts several fitted equations (and a few naive rules) on the table at once and referees between them.

Term	L4	L5	L6	L7
$\alpha$ , $\delta t$	On	On	On	On
$\phi_j y_{t-j}$ (all $p$ )	On	On	On	On
$\theta_l \epsilon_{t-l}$ (all $q$ )	On	On	On	On
Seasonal block $(P, D, Q)_s$	Off	On	On	On
$\epsilon_t$	On	On	On	On
What we do with it	Select	Diagnose	Forecast	Evaluate and combine $M$ forecasts

The PS3 Design: Many Origins, One Horizon

The whole module runs on one design — the same one the midterm specified and Problem Set 3 uses. We restate it here rather than assume you have built it.

|<---------- train ---------->|<-- weights -->|<-- holdout -->|
1                            T  T+1       T+12 T+13        T+24
   estimate every model here    fit comb. weights  score everything
                             ^
                    forecast origin T = N_total - 24

Parameters are estimated once, on $y_1, \ldots, y_T$ , and never again.
The information set rolls. Every forecast is one step ahead, $h = 1$ : the forecast of $y_t$ , written $\tilde y_{t \mid t-1}$ , uses actuals up to $y_{t-1}$ . When $y_{t}$ arrives, it becomes information for the forecast of $y_{t+1}$ — but the coefficients do not move.
The weights window ( $T_w = 12$ dates) is where we will estimate combination weights in Part B. The holdout ( $N = 12$ dates) is where everything gets scored.

Module 6’s design was one origin, many horizons: stand at December 2017 and forecast 1, 2, …, 24 months ahead. Module 7’s is many origins, one horizon: stand at each month in turn and forecast the next one. Keep the two straight — they produce different error series, and (as §7.4 explains) different statistical properties.

Six Forecast Tracks

We compare $M = 6$ forecasters, indexed $m = 1, \ldots, M$ . Forecaster $m$ ’s one-step forecast of $y_t$ is $\tilde y^{(m)}_{t \mid t-1}$ , and its forecast error is

$e_{m,t} = y_t - \tilde y^{(m)}_{t \mid t-1}.$

(This $e$ is an out-of-sample forecast error, an observable quantity — not the innovation $\epsilon_t$ , and not an in-sample residual. A word on notation: from this module on, a forecast carries a tilde, $\tilde y$ , and a hat means a fitted value or an estimate, such as $\hat\phi$ ; Module 6’s $\hat y_{t+h \mid t}$ is the same object we now write $\tilde y_{t+h \mid t}$ .) Three tracks come from ARIMA models; three are $k$ -period average forecasts,

$\tilde y_{t \mid t-1} = \frac{1}{k}\sum_{j=1}^{k} y_{t-j},$

with $k = 1$ (locf, “last observation carried forward” — the naive random walk), a chosen $k$ (avg_k), and $k = t - 1$ (meanf, the expanding mean of everything so far). Do not call these “moving averages.” In this course moving average means the MA( $q$ ) block of the master equation, $\theta_l \epsilon_{t-l}$ — a completely different object.

Track	UNRATE specification	Why it is in the pool
`arima_1`	ARIMA $(1,1,2)(1,0,1)_{12}$	The Module 5/6 model
`arima_2`	ARIMA $(1,1,0)$	A simpler credible alternative
`arima_3`	`auto.arima()` on train	The automatic pick
`locf`	$k = 1$	The naive benchmark
`avg_k`	$k = 3$	Smooths month-to-month sampling noise
`meanf`	$k = t-1$	The long-run mean

Why $k = 3$ for UNRATE? Module 2 established that UNRATE behaves like an I(1) series: almost all of the information about next month sits in the latest value. Averaging over a long window would drag the forecast back toward levels the series left long ago; averaging over three months only smooths the household survey’s month-to-month sampling noise. The choice of $k$ is a statement about persistence.

Loading UNRATE and Building the Tracks

To pull the series live, store your FRED key in the environment (never hardcode it) and run:

library(fredr)
fredr_set_key(Sys.getenv("FRED_KEY"))
unrate_raw <- fredr(series_id = "UNRATE",
                    observation_start = as.Date("1960-01-01"),
                    observation_end   = as.Date("2019-12-31"))
unrate <- ts(unrate_raw$value, start = c(1960, 1), frequency = 12)

So that these notes render identically with or without a key, we load the course’s cached copy instead — the same 1960–2019 window as Modules 4–6 (the April 2020 spike would dominate everything; see Module 4, §4.7.1).

library(forecast); library(ggplot2); library(quadprog)
source("../helpers/data_utils.R")      # load_cached_fred()
source("../helpers/simulators.R")      # sarima_simulator()
source("../helpers/forecast_utils.R")  # rmse(), mae(), mape(), mase(), theil_u2()

fmt <- function(x, d = 2) formatC(x, format = "f", digits = d, big.mark = ",")

unrate <- load_cached_fred("../data/UNRATE.csv")
unrate <- window(unrate, start = c(1960, 1), end = c(2019, 12))

N_total  <- length(unrate)
T_origin <- N_total - 24        # last training date: December 2017
c(N_total = N_total, T = T_origin)

N_total       T 
    720     696

One small function does all the subsetting. It uses window(), not square brackets:

# First n observations of a ts, KEEPING its frequency and dates.
upto <- function(x, n) window(x, end = time(x)[n])

train_un <- upto(unrate, T_origin)
frequency(train_un)            # still 12

[1] 12

frequency(unrate[1:T_origin])  # `[` returns a plain vector: frequency 1

[1] 1

That second line is a trap worth one sentence. unrate[1:T] silently drops the monthly frequency. With Arima(..., model = fit) the stored seasonal period survives, so the midterm’s y[1:(T + 12)] idiom works for refiltering — but any fresh fit on such a subset (auto.arima(), a seasonal Arima()) silently loses the season. Use window(), or rebuild with ts(..., frequency = 12).

Now estimate the three ARIMA models once, on train:

fit1 <- Arima(train_un, order = c(1, 1, 2),
              seasonal = list(order = c(1, 0, 1), period = 12))   # Module 5/6
fit2 <- Arima(train_un, order = c(1, 1, 0))                        # simpler
fit3 <- auto.arima(train_un)                                        # automatic
fit3

Series: train_un 
ARIMA(1,1,3)(2,0,2)[12] 

Coefficients:
         ar1      ma1     ma2     ma3    sar1     sar2     sma1    sma2
      0.8902  -0.9112  0.1766  0.0535  0.7527  -0.2728  -0.9904  0.2575
s.e.  0.0324   0.0484  0.0499  0.0408  0.3005   0.1922   0.3008  0.2588

sigma^2 = 0.02525:  log likelihood = 294.76
AIC=-571.52   AICc=-571.26   BIC=-530.62

fits_un <- list(arima_1 = fit1, arima_2 = fit2, arima_3 = fit3)

On the training sample auto.arima() picks ARIMA(1,1,3)(2,0,2)[12] — a richer model than the Module 5 choice.

The track builder. For the ARIMA tracks we use the Module 6 idiom (§6.7): Arima(series, model = fit) re-applies the training coefficients to a longer series without re-estimating, and fitted() then returns, at every date $t$ , the one-step forecast built from information up to $t-1$ (the function is called fitted(), but out of sample these are forecasts, so they carry a tilde, $\tilde y^{(m)}_{t \mid t-1}$ ). For the $k$ -period averages we just average the previous $k$ actuals. (upto(), kavg(), and tracks_on() are built inline here; their canonical copies live in helpers/track_utils.R, which Module 8 and Problem Set 3 can source.)

# k-period average forecast of y_t for each t in `dates`:
#   mean of the k previous actuals; k = Inf gives the expanding mean (meanf).
kavg <- function(x, k, dates) {
  x <- as.numeric(x)
  sapply(dates, function(t) {
    if (is.infinite(k)) mean(x[1:(t - 1)]) else mean(x[(t - k):(t - 1)])
  })
}

# The six one-step tracks for the dates in `dates`.
# ARIMA parameters stay at their training values (model = fit);
# only the information set rolls forward.
tracks_on <- function(x, fits, dates, k = 3) {
  xs <- upto(x, max(dates))
  os <- function(fit) as.numeric(fitted(Arima(xs, model = fit)))[dates]
  cbind(arima_1 = os(fits[[1]]),
        arima_2 = os(fits[[2]]),
        arima_3 = os(fits[[3]]),
        locf    = kavg(x, 1,   dates),
        avg_k   = kavg(x, k,   dates),
        meanf   = kavg(x, Inf, dates))
}

w_dates <- T_origin + 1:12      # weights window: 2018
h_dates <- T_origin + 13:24     # holdout:        2019

F_w <- tracks_on(unrate, fits_un, w_dates)   # T_w x M
F_h <- tracks_on(unrate, fits_un, h_dates)   # N   x M
y_w <- as.numeric(unrate)[w_dates]
y_h <- as.numeric(unrate)[h_dates]
round(head(F_h, 3), 3)

     arima_1 arima_2 arima_3 locf avg_k meanf
[1,]   3.926   3.911   3.922  3.9 3.833 6.005
[2,]   3.988   4.011   3.991  4.0 3.900 6.002
[3,]   3.871   3.779   3.879  3.8 3.900 5.999

Note upto(x, max(dates)): for the weights window that is upto(unrate, T + 12), and for the holdout upto(unrate, T + 24) — exactly the PS3 calls. (PS3 is entitled to see the holdout actuals as they arrive; the midterm, which stopped at the weights window, was not.) Because a one-step forecast at $t$ uses only data before $t$ , extending the series never changes an earlier one-step forecast.

The alignment check. Before computing a single metric, confirm that locf on date $t$ equals $y_{t-1}$ exactly. If it does not, every number downstream is wrong.

all.equal(F_h[, "locf"], as.numeric(unrate)[h_dates - 1])

[1] TRUE

all.equal(F_w[, "locf"], as.numeric(unrate)[w_dates - 1])

[1] TRUE

Here are the six tracks against the actuals:

dates_all <- c(w_dates, h_dates)
tr_df <- data.frame(
  time  = rep(as.numeric(time(unrate))[dates_all], 6),
  value = c(rbind(F_w, F_h)),
  track = rep(colnames(F_w), each = 24))
ggplot() +
  geom_line(data = subset(tr_df, track != "meanf"),
            aes(time, value, colour = track), linewidth = 0.6) +
  geom_line(aes(x = as.numeric(time(unrate))[dates_all],
                y = as.numeric(unrate)[dates_all]),
            colour = "black", linewidth = 1.1) +
  geom_vline(xintercept = mean(as.numeric(time(unrate))[T_origin + 12:13]),
             linetype = "dashed") +
  annotate("text", x = 2018.4, y = 4.25, label = "weights window") +
  annotate("text", x = 2019.4, y = 4.25, label = "holdout") +
  labs(x = NULL, y = "UNRATE (%)", colour = "Track",
       title = "UNRATE 2018-2019: actuals (black) and five one-step tracks",
       subtitle = paste0("meanf is off the chart, near ",
                         fmt(mean(F_h[, "meanf"]), 1), "%")) +
  theme_bw() + theme(legend.position = "bottom")

Five of the six tracks hug the actuals; they differ by hundredths. The sixth, meanf, is not on the chart at all: on an I(1) series, the expanding mean is a sixty-year average of a level that has wandered, and it sits near 6.0% while unemployment runs below 4%. Keep that track in mind — it will matter a great deal in §7.8.

The Course Simulation

Every topic in this module also runs on one simulated DGP, so we can see what the tools do when we know the truth. It is a stationary seasonal AR:

$(1 - 0.7L)(1 - 0.5L^{12})\,y_t = \alpha + \epsilon_t, \qquad \epsilon_t \sim N(0,1), \qquad \alpha = 5(1-0.7)(1-0.5) = 0.75,$

so the process mean is $\mu = 5$ . The six simulated tracks mirror the UNRATE pool: arima_1 = ARIMA $(1,0,0)(1,0,0)_{12}$ with a mean (the true specification), arima_2 = ARIMA $(1,0,0)$ with a mean (it misses the seasonal block), arima_3 = ARIMA $(0,1,1)$ (a needless difference — a plausible student error), and locf, avg_k ( $k = 3$ ), meanf. Each ARIMA is estimated once on a 240-month training sample. (In the simulation arima_3 is a fixed specification rather than auto.arima(), so that every replication of the Monte Carlos below compares the same six forecasters.) Section 7.7 explains why the DGP was built this way.

# The course DGP: (1 - 0.7L)(1 - 0.5L^12) y_t = 0.75 + eps_t,  mean 5
course_dgp <- function(n) {
  sarima_simulator(n = n, alpha = 0.75, phi = 0.7, Phi = 0.5, s = 12, sigma = 1)
}

# Estimate the three simulated ARIMA tracks once, on the training sample
fit_sim <- function(train) list(
  arima_1 = Arima(train, order = c(1, 0, 0),
                  seasonal = list(order = c(1, 0, 0), period = 12),
                  method = "CSS-ML"),                       # true specification
  arima_2 = Arima(train, order = c(1, 0, 0), method = "CSS-ML"),  # no season
  arima_3 = Arima(train, order = c(0, 1, 1), method = "CSS-ML"))  # over-differenced

n_train <- 240
set.seed(6376)
y_sim    <- course_dgp(n_train + 24)          # train 240 | weights 12 | holdout 12
fits_sim <- fit_sim(upto(y_sim, n_train))
Fs_w <- tracks_on(y_sim, fits_sim, n_train + 1:12)
Fs_h <- tracks_on(y_sim, fits_sim, n_train + 13:24)
ys_w <- as.numeric(y_sim)[n_train + 1:12]
ys_h <- as.numeric(y_sim)[n_train + 13:24]
all.equal(Fs_h[, "locf"], as.numeric(y_sim)[n_train + 12:23])   # alignment

[1] TRUE

We now have two complete PS3-style datasets: twelve weights-window dates and twelve holdout dates, six tracks each — one simulated, one real.

7.2 Three Ways to Run an Evaluation

Module 6 (§6.3) named the evaluation designs in passing. Here we pin them down, because “out-of-sample evaluation” means different things in different papers, and the difference matters for what your numbers mean. The question is: as the forecast origin moves forward, do you re-estimate the model, and on which data?

# Dot diagram (after Hyndman's "rolling forecasting origin" figure):
# one row per forecast origin, one dot per time point.
T0 <- 10; n_orig <- 6; W <- 6
roles <- c("Estimation sample", "Information only (not re-estimated)",
           "Forecast target", "Not used")
dots <- do.call(rbind, lapply(c("Fixed-parameter one-step (PS3)",
                                "Expanding window, re-estimated",
                                "Fixed-width rolling window"), function(sch) {
  do.call(rbind, lapply(1:n_orig, function(i) {
    target <- T0 + i
    tt <- 1:(T0 + n_orig)
    role <- ifelse(tt == target, roles[3], roles[4])
    if (sch == "Fixed-parameter one-step (PS3)") {
      role[tt <= T0] <- roles[1]
      role[tt > T0 & tt < target] <- roles[2]
    } else if (sch == "Expanding window, re-estimated") {
      role[tt < target] <- roles[1]
    } else {
      role[tt < target & tt >= target - W] <- roles[1]
    }
    data.frame(scheme = sch, origin = i, time = tt, role = role)
  }))
}))
dots$scheme <- factor(dots$scheme, levels = unique(dots$scheme))
dots$role   <- factor(dots$role, levels = roles)
ggplot(dots, aes(time, origin, colour = role)) +
  geom_point(size = 3) +
  facet_wrap(~ scheme, ncol = 1) +
  scale_y_reverse(breaks = 1:n_orig) +
  scale_colour_manual(values = c("steelblue", "lightblue", "firebrick", "grey85")) +
  labs(x = "Time", y = "Forecast origin", colour = NULL,
       title = "Three evaluation schemes, one-step-ahead forecasts") +
  theme_bw() + theme(legend.position = "bottom")

Scheme	Parameters	Training window at origin $t$	Fits per model	Used in
Fixed-parameter one-step	Estimated once on $1{:}T$	Fixed at $1{:}T$ for estimation; the information set still rolls to $t-1$	1	PS3, the midterm, all of Module 7
Expanding window, re-estimated	Re-estimated at every origin	$1{:}(t-1)$	one per origin	`tsCV()` default (Deeper Dive below)
Fixed-width rolling window	Re-estimated at every origin	$(t-W){:}(t-1)$	one per origin	`tsCV(..., window = W)` (Deeper Dive below)

Two sentences to remember. Re-estimation buys robustness to change at the cost of more computation and noisier parameters. Fixed-width windows forget old regimes on purpose — useful if you believe the world changed, costly if it did not, because each fit uses fewer observations.

The fixed-parameter scheme is the cheapest and the cleanest for teaching: one fit per model, and every forecast uses the same coefficients, so differences between tracks are differences between models, not between estimation samples. It is the only scheme this module builds.

Four-step note: this section is a definitional bridge — Core has no Math, Simulate, or Real Data step here, by design (design principles §2, “When to deviate”). The optional tsCV() Deeper Dive below applies the two re-estimating schemes to UNRATE, but it is not part of Core.

Beyond the Fixed-Parameter Scheme

The two re-estimating schemes are worth seeing once in code, even though this module does not build on them. The optional Deeper Dive below runs both on UNRATE.

Deeper Dive — `tsCV()`: Re-Estimating at Every Origin

The forecast package implements both re-estimating schemes in one function, tsCV(y, forecastfunction, h). It walks the origin forward, calls your forecast function on the data up to that origin, and records the $h$ -step error. We restrict it to the origins of the last ten years (2010–2019) so it runs in well under a minute, and compare arima_1’s specification against the naive forecast. One quirk: with window = W, tsCV() starts at origin W + initial, so the rolling call needs initial = n0 - 119 (not n0) to cover the same origins as the other two.

ff_sarima <- function(x, h) {
  forecast(Arima(x, order = c(1, 1, 2),
                 seasonal = list(order = c(1, 0, 1), period = 12)), h = h)
}
ff_naive <- function(x, h) rwf(x, h = h)

n0 <- length(unrate) - 120
e_expand <- tsCV(unrate, ff_sarima, h = 1, initial = n0)                    # expanding
e_roll   <- tsCV(unrate, ff_sarima, h = 1, window = 120, initial = n0 - 119)  # W = 120
e_naive  <- tsCV(unrate, ff_naive,  h = 1, initial = n0)

# Pair the errors date by date: keep only origins where all three exist
ok <- complete.cases(e_expand, e_roll, e_naive)
n_tscv <- sum(ok)
data.frame(scheme = c("expanding SARIMA", "rolling SARIMA (W = 120)", "naive"),
           n_errors = n_tscv,
           RMSE = c(sqrt(mean(e_expand[ok]^2)),
                    sqrt(mean(e_roll[ok]^2)),
                    sqrt(mean(e_naive[ok]^2))))

                    scheme n_errors      RMSE
1         expanding SARIMA      119 0.1356625
2 rolling SARIMA (W = 120)      119 0.1375083
3                    naive      119 0.1489459

dm_tscv <- dm.test(e_expand[ok], e_naive[ok], h = 1, power = 2)
dm_tscv


    Diebold-Mariano Test

data:  e_expand[ok]e_naive[ok]
DM = -2.3478, Forecast horizon = 1, Loss function power = 2, p-value =
0.02055
alternative hypothesis: two.sided

Three points.

More evaluation points. Re-estimating schemes produce an error at every origin you can afford to fit — here 119, against the PS3 holdout’s 12. That is the cheapest way to buy statistical power for the Diebold-Mariano test of §7.4. Over these ten years the re-estimated SARIMA beats the naive with a DM $p$ -value of 0.021 — compare the verdict on twelve months in §7.4. (More dates is not the only difference: the window and the scheme differ too.)
Robustness versus precision. The rolling window adapts if the dynamics change, but each fit uses only 120 months; the expanding window uses everything and adapts slowly.
Alignment. tsCV() stores each error on the origin date: element $t$ holds $y_{t+h} - \tilde y_{t+h \mid t}$ . Shift by $h$ before lining its errors up against the PS3 tracks, which are stored on the target date.

None of this is assessable, and PS3 does not ask for tsCV().

7.3 Relative Metrics: Theil’s $U_2$ and MASE

Concept: Two Different Questions

RMSE, MAE, and MAPE (Module 6, §6.7) are absolute scores: RMSE and MAE are in the units of $y$ , and MAPE is a percentage. An RMSE of 0.1 percentage points sounds good — but is it? It depends entirely on how hard the series is to forecast. A model is worth building only if it does better than something trivial, so we need relative scores. There are two in common use, and they answer different questions. Conflating them is the single most common error in this part of the module.

Theil’s $U_2$ asks: did the model beat the naive forecast in this holdout? It is a head-to-head race: the model and locf, on the same twelve dates.
MASE (mean absolute scaled error) asks: how big was the holdout error relative to the naive forecast’s typical one-step miss in training? It is a yardstick, not a race. The yardstick was measured on a different stretch of data.

Math

For model $m$ , over the $N = 12$ holdout dates $t = T+13, \ldots, T+24$ , with the naive forecast $\tilde y^{(\texttt{locf})}_{t \mid t-1} = y_{t-1}$ :

$U_2^{(m)} = \frac{\text{RMSE}^{(m)}}{\text{RMSE}^{(\texttt{locf})}} = \sqrt{\frac{\sum_{t=T+13}^{T+24}\big(y_t - \tilde y^{(m)}_{t \mid t-1}\big)^2}{\sum_{t=T+13}^{T+24}\big(y_t - \tilde y^{(\texttt{locf})}_{t \mid t-1}\big)^2}},$

$\text{MASE}^{(m)} = \frac{\frac{1}{N}\sum_{t=T+13}^{T+24} \big|y_t - \tilde y^{(m)}_{t \mid t-1}\big|}{\frac{1}{T-1}\sum_{t=2}^{T} \big|y_t - \tilde y^{(\texttt{locf})}_{t \mid t-1}\big|}.$

The two metrics use the naive differently. In $U_2^{(m)}$ the numerator and denominator cover the same holdout dates. In $\text{MASE}^{(m)}$ the numerator is model $m$ on the holdout, and the denominator is locf on the training sample $t = 2, \ldots, T$ only.

Reading $U_2$ . $U_2 < 1$ : the model beat locf on these dates. $U_2 = 1$ : a tie. $U_2 > 1$ : it lost. Both RMSEs are computed on exactly the same dates, so nothing else differs.

Reading MASE. The denominator is the mean absolute error the naive forecast made, one step ahead, across the training sample (Hyndman and Koehler 2006). MASE $< 1$ says your holdout errors were smaller than the naive’s typical training-sample miss. It does not say you beat the naive in the holdout — the naive’s holdout errors never enter the formula. If the holdout was an unusually calm stretch, everything — including the naive — will have MASE below 1.

Seasonal MASE. For seasonal data the yardstick can instead be the seasonal-naive forecast $\tilde y^{(\text{snaive})}_{t \mid t-1} = y_{t-s}$ , still measured on training only:

$\text{MASE}_s^{(m)} = \frac{\frac{1}{N}\sum_{t=T+13}^{T+24} \big|y_t - \tilde y^{(m)}_{t \mid t-1}\big|}{\frac{1}{T-s}\sum_{t=s+1}^{T} \big|y_t - y_{t-s}\big|}.$

Why not just MAPE? MAPE divides by $y_t$ , so it is undefined when $y_t = 0$ and explodes near it. MASE scales by a training-sample error instead, so it is defined for any series that is not constant. UNRATE never touches zero, so MAPE is fine here; on an inflation rate or an interest-rate spread it would not be.

The helpers in forecast_utils.R are theil_u2(actual, forecast, benchmark) and mase(actual, forecast, train, m = 1). One notation flag: the helper’s argument m is the seasonal lag (m = 12 gives $\text{MASE}_s$ for monthly data) — in the math we write that lag as $s$ , because $m$ is the forecaster index in this module.

Let’s compute both by hand once for arima_1 on UNRATE, and confirm the helpers agree:

f_a1   <- F_h[, "arima_1"]
f_locf <- F_h[, "locf"]

u2_hand   <- sqrt(mean((y_h - f_a1)^2)) / sqrt(mean((y_h - f_locf)^2))
mase_hand <- mean(abs(y_h - f_a1)) / mean(abs(diff(as.numeric(train_un))))

rbind(by_hand = c(U2 = u2_hand, MASE = mase_hand),
      helper  = c(U2 = theil_u2(y_h, f_a1, f_locf),
                  MASE = mase(y_h, f_a1, train_un)))

              U2      MASE
by_hand 1.010649 0.6423775
helper  1.010649 0.6423775

Simulate: The Metrics Table on the Course DGP

One function builds the whole scorecard. Each row is a track; locf is always the $U_2$ benchmark.

Four-step note: no new assumption-break here; the §7.4 break (shrink the evaluation sample) covers Part A.

metrics_table <- function(actual, Fmat, train) {
  bench <- Fmat[, "locf"]
  t(apply(Fmat, 2, function(f) c(
    RMSE = rmse(actual, f), MAE = mae(actual, f), MAPE = mape(actual, f),
    MASE = mase(actual, f, train), U2 = theil_u2(actual, f, bench))))
}
met_sim <- metrics_table(ys_h, Fs_h, upto(y_sim, n_train))
round(met_sim, 3)

         RMSE   MAE   MAPE  MASE    U2
arima_1 1.120 0.887 20.094 0.831 0.861
arima_2 1.208 1.072 23.237 1.004 0.929
arima_3 1.372 1.302 27.260 1.219 1.054
locf    1.301 1.224 25.745 1.146 1.000
avg_k   1.793 1.550 32.269 1.451 1.378
meanf   1.701 1.420 29.363 1.329 1.308

Read it in three passes.

The correctly specified model wins this realization, by every metric. arima_1 (the true specification, with coefficients estimated from 240 months) has the smallest RMSE and MAE, and $U_2 =$ 0.86: on these twelve dates it beat locf.
$U_2$ and MASE can disagree. Look at arima_2: $U_2 =$ 0.93 says it beat the naive in this holdout, while MASE $=$ 1.00 says its holdout error was about the size of the naive’s typical training miss. Both are true. They answer different questions.
Check whether the metrics agree on the ranking. Here RMSE and MAE put the six tracks in the same order, because the gaps are large. When two tracks are close, squared and absolute loss can rank them differently — Module 6’s lesson that the loss function picks the question. (UNRATE, below, is such a case: locf has the smallest MAE, while on RMSE it and arima_3 are essentially tied — see the $U_2$ column there — so neither metric orders the close tracks with any confidence.)

Real Data: The Metrics Table on UNRATE

met_un <- metrics_table(y_h, F_h, train_un)
round(met_un, 3)

         RMSE   MAE   MAPE   MASE     U2
arima_1 0.097 0.083  2.258  0.642  1.011
arima_2 0.097 0.079  2.127  0.607  1.016
arima_3 0.096 0.080  2.168  0.618  0.999
locf    0.096 0.075  2.030  0.580  1.000
avg_k   0.108 0.092  2.478  0.709  1.124
meanf   2.315 2.312 63.103 17.875 24.182

credible <- c("arima_1", "arima_2", "arima_3", "locf", "avg_k")

This is the contrast we warned about, in the wild. Every credible track has MASE between 0.58 and 0.71 — comfortably below 1 — while its $U_2$ sits between 1.00 and 1.12. In words: small errors by the training yardstick, and no better than locf on these dates (the closest, arima_3 at 0.999, is a tie for practical purposes, not a win). 2019 was a calm glide in unemployment, so every reasonable forecast had small errors — including the naive. MASE $< 1$ here is a statement about 2019 being easy, not about the models being good.

And then there is meanf: RMSE 2.32, MASE 17.9, $U_2$ 24.2. Recall Module 2: the series’ order of integration decides which baselines are even credible. An expanding mean is a sensible forecast for a stationary series (it converges to $\mu$ ); for an I(1) series there is no $\mu$ to converge to.

Other Theil Statistics, and `accuracy()` Defaults

forecast::accuracy() reports MASE and a “Theil’s U,” and both differ from the course definitions in ways that bite:

# Module 6's fixed-origin 24-step forecast of 2018-2019, scored by accuracy()
fc24   <- forecast(fit1, h = 24)
test24 <- window(unrate, start = time(unrate)[T_origin + 1])
acc    <- accuracy(fc24, test24)
acc["Test set", c("RMSE", "MASE", "Theil's U")]

     RMSE      MASE Theil's U 
0.4052396 0.4156320 3.9722351

c(course_MASE_m1  = mase(test24, fc24$mean, train_un, m = 1),
  course_MASE_m12 = mase(test24, fc24$mean, train_un, m = 12))

 course_MASE_m1 course_MASE_m12 
       2.510409        0.415632

MASE. For any ts with frequency above 1, accuracy() scales by the seasonal-naive training MAE by default — our m = 12 version, 0.416 here, against 2.510 for the non-seasonal scale the course uses. UNRATE is seasonally adjusted at the source, so a 12-month difference measures a year of drift, not seasonality; the non-seasonal scale is the honest yardstick. Pass d = 1, D = 0 to accuracy() to get it, or use mase().
Theil’s U. accuracy()’s “Theil’s U” works with relative changes, $(\tilde y_{t \mid T} - y_{t-1})/y_{t-1}$ against $(y_t - y_{t-1})/y_{t-1}$ , and drops the first test date. It is a cousin of $U_2$ , not the same number.
accuracy() is also designed around fixed-origin forecast objects, not the one-step tracks we build. The course computes MASE and $U_2$ with the helpers, so the definition is the one above.

Deeper Dive — Theil’s $U_1$

Older literature sometimes reports a different Theil statistic,

$U_1 = \frac{\sqrt{\frac1N\sum_t \big(y_t - \tilde y_{t \mid t-1}\big)^2}}{\sqrt{\frac1N\sum_t y_t^2} + \sqrt{\frac1N\sum_t \tilde y_{t \mid t-1}^2}},$

which is bounded between 0 (perfect) and 1. The bound sounds attractive, but $U_1$ has no benchmark built in and depends on the level of the series, so a model can score near 0 just because $y_t$ is large. The course uses $U_2$ . Not assessable.

7.4 The Diebold-Mariano Test

Concept: Scoring Is Not Testing

A metrics table tells you which number is smaller. It does not tell you whether the gap is larger than what the noise in twelve forecast errors would produce by chance. Module 6’s question was exactly this: two RMSEs that differ in the second decimal place — is that a real difference?

The taste test. Two coffee shops, one sip each, and you prefer shop A. Is shop A better? One sip is not a verdict. Twelve paired sips, one from each shop each morning, with a running tally of how much better A tasted each day — now you can ask whether the average difference is distinguishable from zero. That is a $t$ -test on the paired differences, and it is all the Diebold-Mariano (DM) test is: a $t$ -test on the loss differentials.

Math

Take two forecasts with out-of-sample errors $e_{1,t}$ and $e_{2,t}$ over an evaluation window of $N$ dates, and a loss function $g(\cdot)$ (Module 6: squared error $g(e) = e^2$ , absolute error $g(e) = |e|$ ). Define the loss differential

$d_t = g(e_{1,t}) - g(e_{2,t}).$

If the two forecasts are equally accurate, the loss differential is zero on average:

$H_0\!: \mathbb{E}[d_t] = 0, \qquad \bar d = \frac{1}{N}\sum_t d_t, \qquad \text{DM} = \frac{\bar d}{\hat\sigma_{\bar d}} \;\stackrel{a}{\sim}\; N(0,1) \text{ under } H_0.$

Here $\hat\sigma_{\bar d}$ is the standard error of $\bar d$ — a standard error that allows for serial correlation in $d_t$ . That clause matters: loss differentials are a time series, and if they are autocorrelated the naive standard error $\hat\sigma_d/\sqrt N$ — where $\hat\sigma_d$ is the sample standard deviation of $d_t$ — is the wrong one. (Why they can be autocorrelated, and what the correction looks like, is the Technical Note below.)

What h = 1 assumes. With h = 1, dm.test() uses the lag-0 variance only, so $\hat\sigma_{\bar d}$ becomes the familiar $\hat\sigma_d/\sqrt N$ (after the small-sample correction). That version does not allow for serial correlation: it assumes that $d_t$ is serially uncorrelated. The assumption is about the loss differential, not merely about the forecast errors. Optimal one-step forecasts have serially uncorrelated errors, and that makes the assumption plausible — but it does not guarantee it, because $d_t = g(e_{1,t}) - g(e_{2,t})$ is a nonlinear function of two error series and can be autocorrelated even when neither error is. A misspecified competitor makes it more likely. (§7.8’s Technical Note checks this for the UNRATE tests; optional.)

Three rules of use:

Out-of-sample errors only. Feed dm.test() forecast errors from dates the models did not train on — never residuals(fit). In-sample residuals reward overfitting.
Same dates, same loss. Both error series must cover the same $N$ dates.
The sign tells you who won. $\text{DM} < 0$ means $\bar d < 0$ : forecast 1 had the smaller average loss.

In R:

dm.test(e1, e2, h = 1, power = 2)   # power = 2: squared-error loss; power = 1: absolute

h is the forecast horizon — always 1 in this module’s forecasts — and power picks the loss. dm.test() applies the Harvey-Leybourne-Newbold (1997) small-sample correction and compares the statistic with a $t_{N-1}$ distribution rather than $N(0,1)$ ; name it, don’t derive it (Enders eq. 2.62 and the text after it).

Worked Example: DM Really Is a $t$ -Test

On the course simulation, the correctly specified model arima_1 (its coefficients estimated from the 240 training months) against locf over the twelve holdout dates:

es_h <- ys_h - Fs_h                          # N x M matrix of holdout errors
dm_sim <- dm.test(es_h[, "arima_1"], es_h[, "locf"], h = 1, power = 2)
dm_sim


    Diebold-Mariano Test

data:  es_h[, "arima_1"]es_h[, "locf"]
DM = -0.89559, Forecast horizon = 1, Loss function power = 2, p-value =
0.3896
alternative hypothesis: two.sided

d_sim <- es_h[, "arima_1"]^2 - es_h[, "locf"]^2   # the loss differentials
t.test(d_sim)$statistic                            # an ordinary one-sample t-test

         t 
-0.8955949

Read the dm.test() output line by line: DM is the statistic (-0.90; negative, so arima_1 had the smaller squared errors); Forecast horizon = 1 and Loss function power = 2 echo your choices; the $p$ -value is 0.390; the alternative is two-sided. And the one-sample $t$ -test on $d_t$ returns the same statistic: at $h = 1$ , after the small-sample correction, DM is literally a $t$ -test on the loss differentials.

Even though arima_1 is the correctly specified model and won on all five metrics in §7.3, at $N = 12$ the test does not reject at 5%. Is that bad luck, or is it typical? Simulate.

Simulate — The Assumption-Break: Make the Sample Tiny

The DGP was built so that the correctly specified model — arima_1, with its coefficients estimated from 240 observations — has a real expected-loss advantage over locf and over arima_2. Estimated coefficients carry error, so that advantage is not guaranteed in every replication, but it is there on average. The question is how often DM detects it with $N$ evaluation dates. Monte Carlo: $R = 200$ replications; each draws a fresh series, estimates the three ARIMAs once on 240 training months, builds one-step tracks over the next 120 months, and runs DM on the first $N \in \{12, 24, 60, 120\}$ of them.

set.seed(6377)
R <- 200
N_grid <- c(12, 24, 60, 120)
pairs  <- c("arima_1 vs locf", "arima_1 vs arima_2")
rej <- array(NA, c(R, length(N_grid), 2),
             dimnames = list(NULL, paste0("N=", N_grid), pairs))

for (r in seq_len(R)) {
  y    <- course_dgp(n_train + 120)
  fits <- fit_sim(upto(y, n_train))                 # estimated once
  e    <- as.numeric(y)[n_train + 1:120] -
          tracks_on(y, fits, n_train + 1:120)       # 120 x 6 one-step errors
  for (j in seq_along(N_grid)) {
    ee <- e[1:N_grid[j], ]
    for (k in 1:2) {
      rival <- c("locf", "arima_2")[k]
      tst <- dm.test(ee[, "arima_1"], ee[, rival], h = 1, power = 2)
      # count a rejection only when it favors arima_1 (DM < 0)
      rej[r, j, k] <- tst$p.value < 0.05 && tst$statistic < 0
    }
  }
}
power_tab <- t(apply(rej, c(2, 3), mean))   # share of directional rejections (p < .05, DM < 0)
power_tab

                    N=12  N=24  N=60 N=120
arima_1 vs locf    0.125 0.405 0.760 0.955
arima_1 vs arima_2 0.100 0.290 0.605 0.860

pw <- data.frame(N = rep(N_grid, each = 2), pair = rep(pairs, 4),
                 rate = c(power_tab))
ggplot(pw, aes(N, rate, colour = pair)) +
  geom_line(linewidth = 1) + geom_point(size = 2.5) +
  # two-sided p < .05 AND DM < 0: under a calibrated null this happens 2.5% of the time
  geom_hline(yintercept = 0.025, linetype = "dashed") +
  annotate("text", x = 120, y = 0.025, vjust = -0.6, hjust = 1, size = 3.5,
           label = "2.5%: directional rejections under a calibrated null") +
  scale_x_continuous(breaks = N_grid) + ylim(0, 1) +
  labs(x = "Evaluation dates N",
       y = "Share of directional rejections in favor of arima_1\n(two-sided p < .05, DM < 0)",
       colour = NULL,
       title = "DM rejection rates: a real gap, but N = 12 rarely shows it") +
  theme_bw() + theme(legend.position = "bottom")

(The test is two-sided; the table counts only rejections in the right direction, with arima_1 the winner. Under a calibrated null, half of the 5% two-sided rejections land on each side, so the benchmark for these directional rates is 2.5%, the dashed line, not 5%.) At $N = 12$ — the PS3 holdout — DM rejects in favor of arima_1 in 12.5% of replications for arima_1 vs locf and 10.0% for arima_1 vs arima_2. At $N = 120$ , 95.5% and 86.0%. The correctly specified model has the advantage on average; with twelve dates, the test usually cannot detect it.

The main lesson is power — with a calibration caveat. State it precisely: DM rejects when the loss gap is large relative to the variability of $d_t$ . The caveat: these rates are rejection rates of the h = 1 test, whose 5% level is calibrated only if $d_t$ is serially uncorrelated. locf and arima_2 leave structure in their errors that arima_1 captures, so that is an assumption about these $d_t$ , not a guarantee; the table measures power to the extent the assumption holds. At $N = 12$ a real but moderate gap is usually invisible. A large gap still shows up at $N = 12$ — §7.8 produces exactly that, a rejection on twelve dates (its $p$ -value is printed there). So “DM failed to reject” never means “the forecasts are equally good”; it means “these dates cannot separate them.”

Real Data: Closing Module 6’s Open Question

Now the question Module 6 left open. arima_1 (the Module 5 model) against locf (the naive) on the 2019 holdout:

e_h <- y_h - F_h                       # N x M holdout errors on UNRATE
dm_un <- dm.test(e_h[, "arima_1"], e_h[, "locf"], h = 1, power = 2)
dm_un


    Diebold-Mariano Test

data:  e_h[, "arima_1"]e_h[, "locf"]
DM = 0.12594, Forecast horizon = 1, Loss function power = 2, p-value =
0.902
alternative hypothesis: two.sided

$\text{DM} =$ 0.13, $p =$ 0.90. Module 6 asked whether the naive’s win over the Module 5 model meant anything. Asked in the PS3 one-step design, on 2019, the answer is no — the test does not reject on these twelve months. (The verdict rests on the h = 1 assumption that these $d_t$ are serially uncorrelated; §7.8 checks it and re-runs the test with a serial-correlation-robust variance.)

Be careful about what this does and does not claim. It does not test Module 6’s 24-step path. That was a different design — one origin, many horizons — and DM does not apply to it. Errors made from a single origin at horizons $1, 2, \ldots, 24$ all share the same early shocks, so they are serially correlated, and their spread grows with the horizon. They are not the comparable sequence of loss differentials DM assumes, so a DM $p$ -value computed on them cannot be trusted. The valid design is the one used here: rolling origins, one horizon. (The Technical Note shows that even its correction for serially correlated errors does not cover the single-origin case.) And it does not say the two forecasts are equally good; the power simulation just showed that twelve dates rarely separate even a correctly specified estimated model from a misspecified one. (The tsCV() Deeper Dive in §7.2, with 119 re-estimated origins, does reject — $p =$ 0.021 — with more dates (and a different window and scheme).)

Technical Note — Why $\hat\sigma_{\bar d}$ Must Allow for Serial Correlation

From Module 6 (§6.5), the $h$ -step forecast error from a correctly specified model is

$e_{t+h \mid t} = \sum_{i=0}^{h-1}\psi_i\,\epsilon_{t+h-i}.$

Two consecutive $h$ -step errors, made from origins $t$ and $t+1$ , share $h-1$ of the same innovations. So even optimal $h$ -step errors are serially correlated up to lag $h-1$ — an MA( $h-1$ ) — and the loss differentials built from them inherit that dependence. At $h = 1$ the errors share nothing and, for a well-specified model, are serially uncorrelated; at $h > 1$ they are not. Uncorrelated errors are not the end of the story, though: the test needs $d_t = g(e_{1,t}) - g(e_{2,t})$ to be serially uncorrelated, and a nonlinear function of two error series can be autocorrelated even when each error series is not (for instance, when both errors’ volatility is persistent, or when one forecaster is misspecified).

With autocorrelated $d_t$ , the variance of the mean is not $\gamma_0^{(d)}/N$ but the long-run variance

$\text{Var}(\bar d) \approx \frac{1}{N}\sum_{k=-(N-1)}^{N-1} \gamma_k^{(d)},$

where $\gamma_k^{(d)}$ is the lag- $k$ autocovariance of $d_t$ . Enders eq. (2.62) truncates the sum at a generic $q$ lags; Diebold and Mariano recommend $q = h - 1$ . Truncation has a defect: the estimated sum can be negative. Newey and West (1987) fix this by down-weighting the higher-lag autocovariances with the Bartlett kernel, $\hat J = \sum_{|k| \le S} \big(1 - \tfrac{|k|}{S+1}\big)\hat\gamma_k^{(d)}$ , which is guaranteed non-negative (Mikusheva, MIT 14.384 Lecture 3, develops the kernel estimator and bandwidth choice). In dm.test(), h = h uses the truncated sum with lags up to $h-1$ (varestimator = "acf", the default) or a Bartlett-weighted one (varestimator = "bartlett"), followed by the Harvey-Leybourne-Newbold correction $\sqrt{(N + 1 - 2h + h(h-1)/N)/N}$ and a $t_{N-1}$ reference. At $h = 1$ only the lag-0 variance is used — which is why the $t$ -test check above matched exactly, and why h = 1 assumes rather than allows for serially uncorrelated $d_t$ . So the h argument does double duty: it is the forecast horizon, and it sets how many autocovariance lags enter $\hat\sigma_{\bar d}$ . That gives a practical check of the h = 1 assumption: look at the ACF of $d_t$ , then re-run the test with a variance that does allow for serial correlation (h = 2 or more, varestimator = "bartlett") and see whether the verdict moves. Raising h buys a serial-correlation-robust variance at the price of a different small-sample factor; one caution: with varestimator = "acf" the truncated sum can come out negative, and dm.test() then warns and falls back to h = 1 (the output then reports horizon 1), so prefer "bartlett" for the sensitivity check. The §7.8 Technical Note runs that check on the UNRATE tests.

Enders also notes a shortcut: regress $d_t$ on a constant and use a HAC (Newey-West) standard error on the constant. That gives the same test; we mention it and do not run it, since the sandwich package is not part of the course toolkit.

One caution for anyone tempted to run DM on Module 6’s design: errors from a single origin at horizons $1, \ldots, 24$ are not a stationary series at all — their variance grows with the horizon and every one of them shares the early innovations. Neither the $h = 1$ formula nor the MA( $h-1$ ) correction describes that object. Rolling origins, one horizon, is the design DM was built for. Not assessable.

Deeper Dive — Nested Models: Clark-McCracken and the Clark-West Fix

DM’s normal (or $t$ ) reference distribution can fail when one model is nested in the other (Clark and McCracken 2001; Enders’ example is an AR(1) inside an ARMA(2,1)). If the small model is the truth, the large one can only lose — its extra parameters are estimated noise — so under the null the loss differential is not centered where the test assumes, and DM is undersized and loses power. Clark and West (2007) adjust the larger model’s loss: with forecast 1 from the small model and forecast 2 from the large one,

$z_t = e_{1,t}^2 - \Big[e_{2,t}^2 - \big(\tilde y^{(1)}_{t \mid t-1} - \tilde y^{(2)}_{t \mid t-1}\big)^2\Big],$

and test whether the mean of $z_t$ is positive, one-sided. Note that in the course simulation arima_2 (an AR(1) with mean) is nested in arima_1 (the seasonal AR), so the second row of the power table is a nested comparison. There the large model is the correctly specified one — the alternative, not the null. Clark-West is designed to address the size and power problems of standard DM in nested comparisons like this one; we have not run it in our simulation, so the table says nothing about how much it would change. The same nesting holds on UNRATE: locf is the random walk ARIMA $(0,1,0)$ , which is nested in arima_1’s ARIMA $(1,1,2)(1,0,1)_{12}$ . So the §7.4 “no” for arima_1 against locf comes from a test with low power at $N = 12$ , applied to a nested pair for which standard DM is not the purpose-built tool — weak evidence of equal accuracy. Not assessable; PS3 does not ask for it.

7.5 Why Combine?

Concept: From Selection to Combination

Section 7.4 left us in an awkward spot. On UNRATE the five credible tracks are descriptively close — their RMSEs differ by hundredths — and the one comparison we tested, arima_1 against locf, does not reject; picking “the winner” of a twelve-month race is mostly picking noise. There is a better question: if you cannot pick a winner, why pick at all? Each forecast carries some information, and the loser may carry information the winner lacks.

The portfolio analogy. Forecasts are assets; forecast errors are their risks. A portfolio of imperfectly correlated assets has lower variance than its components, because some losses offset others. A combination of imperfectly correlated forecasts has smaller error variance for the same reason. The idea is old: Bates and Granger (1969) proposed weighting forecasts by their accuracy; Granger and Ramanathan (1984) recast the weights as a regression; Diebold (1989) connected combination to forecast encompassing; Diebold and Pauly (1990) proposed shrinking estimated weights toward the simple average; and Stock and Watson (2004) named the empirical regularity that drives Part B of this module — the forecast combination puzzle.

A combined forecast is

$\tilde y^{(c)}_{t \mid t-1} = w_0 + \sum_{m=1}^{M} w_m\,\tilde y^{(m)}_{t \mid t-1} = w_0 + \mathbf{w}'\tilde{\mathbf{y}}_t,$

where $\tilde{\mathbf y}_t$ is the $M \times 1$ vector of the $M$ forecasts of $y_t$ , $\mathbf w$ the $M \times 1$ weight vector, and $w_0$ an intercept that is zero for every method except one. (Weights are always $w$ in this course — never $\alpha$ or $\beta$ , which are taken. The cardinal DGP-vs-estimation rule applies to them too: optimal population weights are starred, $w^*$ ; weights estimated from data carry hats, $\hat w_0$ , $\hat w_m$ , $\hat{\mathbf w}$ ; unhatted $w$ is a generic or fixed weight.)

Math: Two-Forecast Portfolio Algebra

Combine two forecasts with weights $w$ and $1-w$ : $\tilde y^{(c)}_{t \mid t-1} = w\,\tilde y^{(1)}_{t \mid t-1} + (1-w)\,\tilde y^{(2)}_{t \mid t-1}$ . Because the weights sum to one, $y_t = w\,y_t + (1-w)\,y_t$ , so the combined error $e_{c,t} = y_t - \tilde y^{(c)}_{t \mid t-1}$ (the subscript $c$ is a label for “combination,” not a forecaster index) is the same combination of the individual errors:

$e_{c,t} = w\,e_{1,t} + (1-w)\,e_{2,t},$

$\text{Var}(e_{c,t}) = w^2\,\text{Var}(e_{1,t}) + (1-w)^2\,\text{Var}(e_{2,t}) + 2w(1-w)\,\text{Cov}(e_{1,t}, e_{2,t}).$

This is Enders eqs. (2.71)–(2.72). Minimize over $w$ — set the derivative to zero, $2w\,\text{Var}(e_{1,t}) - 2(1-w)\,\text{Var}(e_{2,t}) + 2(1-2w)\,\text{Cov}(e_{1,t},e_{2,t}) = 0$ — and solve:

$w^* = \frac{\text{Var}(e_{2,t}) - \text{Cov}(e_{1,t}, e_{2,t})}{\text{Var}(e_{1,t}) + \text{Var}(e_{2,t}) - 2\,\text{Cov}(e_{1,t}, e_{2,t})}.$

This is Enders eq. (2.73). With zero covariance (the form Enders states in the text after it, unnumbered), divide top and bottom by $\text{Var}(e_{1,t})\,\text{Var}(e_{2,t})$ :

$w^* = \frac{\text{Var}(e_{1,t})^{-1}}{\text{Var}(e_{1,t})^{-1} + \text{Var}(e_{2,t})^{-1}}.$

Read the results.

Diversification works unless the errors are perfectly correlated (for weights between 0 and 1, with equal error variances). With equal error variances $V$ (a scalar shorthand for $\text{Var}(e_{1,t}) = \text{Var}(e_{2,t})$ , used only in this bullet) and zero covariance, $w^* = 0.5$ and $\text{Var}(e_{c,t}) = 0.25V + 0.25V = 0.5V$ : averaging two equally good, uncorrelated forecasts halves the error variance. (Enders works the same example and states the composite variance as a quarter of either forecast’s; each forecast contributes $0.25V$ , but the total is $0.5V$ .)
With zero covariance, each forecast’s weight is proportional to the inverse of its error variance. The less accurate forecast still gets weight — just less. This is the motivation for Method 2 in §7.6.
Correlation eats the gain. Among weights between 0 and 1, the best combination’s advantage over the better forecast alone shrinks as the error correlation rises, and disappears entirely once $\text{Cov}(e_{1,t}, e_{2,t}) \ge \text{Var}(e_{1,t})$ (forecast 1 the better one). Dividing both sides by $\sqrt{\text{Var}(e_{1,t})\,\text{Var}(e_{2,t})}$ puts the same threshold in terms of the error correlation $\rho_{1,2}$ (defined below the bullets): the gain vanishes once $\rho_{1,2} \ge \sqrt{\text{Var}(e_{1,t})/\text{Var}(e_{2,t})}$ , about 0.82 for the variances in the figure below. Beyond that point the unconstrained formula gives $w^* > 1$ : it shorts the worse forecast (in the figure below, the $\rho_{1,2} = 0.9$ curve is still falling at $w = 1$ ).

(We write Var and Cov, as Enders does, rather than $\sigma_1^2$ and $\sigma_{12}$ : in this course $\sigma^2$ is reserved for the innovation variance.)

A picture makes the last point concrete. Take $\text{Var}(e_{1,t}) = 1$ , $\text{Var}(e_{2,t}) = 1.5$ and vary the correlation between the two forecasts’ errors, $\rho_{1,2} = \text{Cov}(e_{1,t}, e_{2,t})/\sqrt{\text{Var}(e_{1,t})\text{Var}(e_{2,t})}$ (two indices separated by a comma, to keep it apart from the lag- $k$ autocorrelation $\rho_k$ — written without the comma, $\rho_{12}$ would read as the lag-12 autocorrelation):

combo_var <- function(w, v1, v2, rho) {
  cv <- rho * sqrt(v1 * v2)
  w^2 * v1 + (1 - w)^2 * v2 + 2 * w * (1 - w) * cv
}
w_grid <- seq(0, 1, by = 0.01)
pf <- expand.grid(w = w_grid, rho = c(-0.5, 0, 0.5, 0.9))
pf$var <- combo_var(pf$w, 1, 1.5, pf$rho)
ggplot(pf, aes(w, var, colour = factor(rho))) +
  geom_line(linewidth = 1) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  labs(x = "Weight w on forecast 1", y = "Var(combined error)",
       colour = "Error correlation",
       title = "Diversification: the combined error variance as a function of w",
       subtitle = "Var(e1) = 1, Var(e2) = 1.5; dashed line = forecast 1 alone") +
  theme_bw() + theme(legend.position = "bottom")

For every correlation below about 0.8 the curve dips below the dashed line — some mix of the two beats the better forecast alone. At $\rho_{1,2} = 0.9$ the optimum is $w^* > 1$ : with highly correlated errors, the formula wants to short the worse forecast. Hold that thought; it is exactly what unconstrained regression weights do in §7.6.

A quick simulation check that the formula is right — draw correlated errors and compare the empirical variance of the combination with the formula:

set.seed(6379)
v1 <- 1; v2 <- 1.5; rho <- 0.3
e1 <- rnorm(1e5)
e2 <- sqrt(v2) * (rho * e1 + sqrt(1 - rho^2) * rnorm(1e5))
cv <- rho * sqrt(v1 * v2)
w_star <- (v2 - cv) / (v1 + v2 - 2 * cv)
c(w_star = w_star,
  formula = combo_var(w_star, v1, v2, rho),
  simulated = var(w_star * e1 + (1 - w_star) * e2))

   w_star   formula simulated 
0.6416308 0.7733040 0.7681394

Four-step note: this section’s assumption-break is shared with §7.7, and its real-data step with §7.8 — the two-forecast algebra is the theory both of them test.

Convention flag. Enders writes the individual forecasts as $f_{it}$ and the regression weights as $\alpha_0, \ldots, \alpha_n$ , and defines the error as $y - f$ in the derivation but $f - y$ in the spread example. The course uses $\tilde y^{(m)}_{t \mid t-1}$ , $w_0$ and $w_m$ , and $e = y - \tilde y$ throughout. The sign convention does not change any variance or weight.

7.6 Four Combination Methods

All four methods use only the weights window: the $T_w = 12$ actuals $y_t$ there and the $T_w \times M$ matrix of forecasts of those dates. All four are then applied, unchanged, to the holdout forecasts as $\hat w_0 + \hat{\mathbf w}'\tilde{\mathbf y}_t$ (hats because, except for the simple average, the weights are estimated from the window). We build each one inline as a small function with a common contract: input the weights-window actuals y and forecast matrix F (F is just an R argument name; in math the letter is reserved for the Granger $F$ -test of Module 9, so the forecast matrix never gets a symbol of its own here); output list(intercept, weights). The canonical copies live in helpers/combination_utils.R, which Problem Set 3 can source.

We demonstrate on the course simulation’s weights window (Fs_w, ys_w). Section 7.7 scores these exact weights on the holdout.

Four-step note: this section is Concept and Math with the inline builds; its assumption-break is §7.7 and its Real Data step is §7.8.

Method 1 — Simple Average

$w_0 = 0, \qquad w_m = \frac{1}{M}.$

Nothing is estimated — which is why these weights carry no hat. That is its whole point: it is the zero-estimation-error benchmark.

comb_simple <- function(y, F) {
  M <- ncol(F)                       # y unused: kept for a common interface
  list(intercept = 0, weights = setNames(rep(1 / M, M), colnames(F)))
}

Method 2 — Inverse-MSE Weights

The $M$ -forecast version of the zero-covariance weight from §7.5 (Bates and Granger 1969; Enders eq. 2.74):

$\hat w_m = \frac{1/\text{MSE}_m}{\sum_{j=1}^{M} 1/\text{MSE}_j}, \qquad \text{MSE}_m = \frac{1}{T_w}\sum_{t \in \text{weights}} e_{m,t}^2, \qquad \hat w_0 = 0.$

$M$ numbers are estimated, but each from its own column only; the weights are non-negative and sum to one by construction. We use MSE rather than the error variance, so a biased forecaster is penalized for its bias too.

comb_inv_mse <- function(y, F) {
  mse <- colMeans((y - F)^2)         # y recycles down each column
  if (any(mse == 0)) stop("a track with zero MSE is a leak, not a forecast")
  list(intercept = 0, weights = (1 / mse) / sum(1 / mse))
}
round(comb_inv_mse(ys_w, Fs_w)$weights, 3)

arima_1 arima_2 arima_3    locf   avg_k   meanf 
  0.271   0.141   0.109   0.093   0.159   0.229

Method 3 — Granger-Ramanathan (Unconstrained Least Squares)

Regress the actuals on the forecasts, with an intercept and no constraints (Granger and Ramanathan 1984; Enders eq. 2.75):

$y_t = w_0 + \sum_{m=1}^{M} w_m\,\tilde y^{(m)}_{t \mid t-1} + v_t, \qquad t \in \text{weights window},$

where $v_t$ is the regression’s disturbance (not an innovation, and not one of the forecast errors $e_{m,t}$ ). OLS on the twelve weights-window dates delivers the estimates $\hat w_0, \hat w_1, \ldots, \hat w_M$ .

The intercept absorbs any systematic bias; the weights can be negative and need not sum to one. It is one lm() call:

coef(lm(ys_w ~ Fs_w))

(Intercept) Fs_warima_1 Fs_warima_2 Fs_warima_3    Fs_wlocf   Fs_wavg_k 
261.7004758   0.7789967   1.8768134  -2.0428886          NA  -0.8332475 
  Fs_wmeanf 
-50.8267308

Two things to see. First, the numbers are wild — large positive and negative weights and an intercept far from zero, on a series whose mean is 5. Second, one coefficient is NA. That is not a bug, and it is worth understanding. In this DGP, arima_2 is an AR(1) with a mean, so its one-step forecast is exactly $\hat\mu(1-\hat\phi) + \hat\phi\,y_{t-1}$ — an intercept plus a multiple of locf. Once the regression has its own intercept, locf adds nothing that arima_2 does not already supply, the two are perfectly collinear, and R drops one (NA). The fitted values are unaffected; the right move is to treat the dropped weight as zero, which is the same as removing the duplicate track. If we left the NA in, every holdout forecast w0 + F %*% w would be NA.

comb_gr <- function(y, F) {
  b <- coef(lm(y ~ F))
  w <- setNames(b[-1], colnames(F))
  w[is.na(w)] <- 0      # aliased (exactly collinear) track: weight 0
  list(intercept = unname(b[1]), weights = w)
}
gr_s <- comb_gr(ys_w, Fs_w)
round(c(w0 = gr_s$intercept, gr_s$weights), 3)

     w0 arima_1 arima_2 arima_3    locf   avg_k   meanf 
261.700   0.779   1.877  -2.043   0.000  -0.833 -50.827

Say it out loud: with $M = 6$ and $T_w = 12$ , Granger-Ramanathan estimates seven numbers from twelve observations. (In the simulation, the collinearity leaves six identifiable; on UNRATE all seven are estimated.) Module 4’s callback: extra parameters that fit the sample are not the same as parameters that forecast.

Method 4 — Constrained Least Squares

No intercept; weights non-negative and summing to one:

$\hat{\mathbf w} = \arg\min_{\mathbf w}\; \sum_{t \in \text{weights}} \big(y_t - \mathbf w'\tilde{\mathbf y}_t\big)^2 \quad \text{subject to}\quad w_m \ge 0,\;\; \sum_{m=1}^M w_m = 1, \qquad \hat w_0 = 0.$

This is a quadratic program. Expanding the sum of squares,

$\sum_t (y_t - \mathbf w'\tilde{\mathbf y}_t)^2 = \sum_t y_t^2 - 2\,\mathbf w'\Big(\sum_t \tilde{\mathbf y}_t y_t\Big) + \mathbf w'\Big(\sum_t \tilde{\mathbf y}_t\tilde{\mathbf y}_t'\Big)\mathbf w,$

so the objective is fully described by an $M \times M$ quadratic part, $\sum_t \tilde{\mathbf y}_t\tilde{\mathbf y}_t'$ , and an $M \times 1$ linear part, $\sum_t \tilde{\mathbf y}_t y_t$ . quadprog::solve.QP() solves exactly this kind of problem. It takes the quadratic part as Dmat (in code, crossprod(F)), the linear part as dvec (crossprod(F, y)), and the constraints as a matrix Amat and a vector bvec, read column by column as “column $\cdot$ weights $\ge$ entry of bvec”, with the first meq columns holding as equalities. (Those are R argument names only; we give them no math symbols, because $D$ , $d$ , $A$ , and $b$ already mean other things in this course.) The constraint matrix stacks the sum-to-one column and the $M$ non-negativity columns — look at them:

M <- ncol(Fs_w)
Amat <- cbind(rep(1, M), diag(M))   # col 1: sum-to-one; cols 2..M+1: w_m >= 0
bvec <- c(1, rep(0, M))
Amat

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    1    1    0    0    0    0    0
[2,]    1    0    1    0    0    0    0
[3,]    1    0    0    1    0    0    0
[4,]    1    0    0    0    1    0    0
[5,]    1    0    0    0    0    1    0
[6,]    1    0    0    0    0    0    1

bvec

[1] 1 0 0 0 0 0 0

Column 1 of Amat with bvec[1] = 1 and meq = 1 says $\sum_m \hat w_m = 1$ exactly; columns 2 through $M+1$ with zeros in bvec say each $\hat w_m \ge 0$ .

comb_cls <- function(y, F) {
  M    <- ncol(F)
  D    <- crossprod(F)                # M x M:  F'F
  d    <- crossprod(F, y)             # M x 1:  F'y
  Amat <- cbind(rep(1, M), diag(M))
  bvec <- c(1, rep(0, M))
  w <- solve.QP(D, d, Amat, bvec, meq = 1)$solution
  w <- pmax(w, 0); w <- w / sum(w)    # clip round-off like -1e-17 to zero
  list(intercept = 0, weights = setNames(w, colnames(F)))
}
cls_s <- comb_cls(ys_w, Fs_w)
round(cls_s$weights, 3)

arima_1 arima_2 arima_3    locf   avg_k   meanf 
  0.657   0.000   0.000   0.000   0.000   0.343

The weights read as shares: “66% arima_1, 34% meanf.” Several weights land exactly on zero: the inequality constraints are doing model selection. Constrained least squares is stable, interpretable, and cannot short a forecaster — but it has no separate intercept for bias adjustment. Reweighting can reduce bias, but it cannot eliminate a common same-direction bias shared by every forecast.

(solve.QP needs the quadratic part Dmat to be positive definite. Exactly collinear tracks — say, two identical columns — break that, and the helper version stops with a message naming them. Without an intercept, arima_2 and locf are not exactly collinear, so CLS runs here.)

Summary

Method	$w_0$	$\sum w_m = 1$	$w_m \ge 0$	Numbers estimated from the window	Estimation risk
Simple average	0	yes	yes	0	none
Inverse-MSE	0	yes	yes	$M$ , one column at a time	low
Constrained LS	0	yes	yes	$M$ , jointly, boxed in	moderate
Granger-Ramanathan	free	no	no	$M + 1$ , jointly, unrestricted	high

ws_sim <- list(SA = comb_simple(ys_w, Fs_w), IMSE = comb_inv_mse(ys_w, Fs_w),
               CLS = comb_cls(ys_w, Fs_w),   GR = comb_gr(ys_w, Fs_w))
round(sapply(ws_sim, function(w) c(w0 = w$intercept, w$weights)), 3)

           SA  IMSE   CLS      GR
w0      0.000 0.000 0.000 261.700
arima_1 0.167 0.271 0.657   0.779
arima_2 0.167 0.141 0.000   1.877
arima_3 0.167 0.109 0.000  -2.043
locf    0.167 0.093 0.000   0.000
avg_k   0.167 0.159 0.000  -0.833
meanf   0.167 0.229 0.343 -50.827

Read down the table as a dial: the simple average estimates nothing, inverse-MSE estimates a little, CLS estimates more but inside a box, Granger-Ramanathan estimates everything. Section 7.7 asks which end of the dial wins on data the weights have not seen.

Technical Note — Full-Covariance Bates-Granger Is a Restricted Regression

For $M$ forecasts, let $\mathbf e_t = y_t\mathbf 1 - \tilde{\mathbf y}_t$ be the $M$ -vector of one-step errors at date $t$ and $\boldsymbol\Sigma_e = \mathbb E[\mathbf e_t\mathbf e_t']$ their $M \times M$ second-moment matrix. (Collision warning: bold $\mathbf e_t$ is also the VAR residual vector of Module 11, and unsubscripted $\boldsymbol\Sigma$ the VAR innovation covariance of Module 12; the subscript $e$ marks the forecasting objects.) Minimizing $\mathbf w'\boldsymbol\Sigma_e\mathbf w$ subject to $\mathbf 1'\mathbf w = 1$ gives

$\mathbf w^* = \frac{\boldsymbol\Sigma_e^{-1}\mathbf 1}{\mathbf 1'\boldsymbol\Sigma_e^{-1}\mathbf 1}$

(Bates and Granger 1969 for $M = 2$ — the $w^*$ of §7.5 — and Newbold and Granger 1974 for $M > 2$ ; Wang, Hyndman, Li and Kang 2023, eq. 1). The equivalence with a regression takes three lines. When $\mathbf 1'\mathbf w = 1$ , $y_t - \mathbf w'\tilde{\mathbf y}_t = \mathbf w'(y_t\mathbf 1 - \tilde{\mathbf y}_t) = \mathbf w'\mathbf e_t$ . So $\sum_t (y_t - \mathbf w'\tilde{\mathbf y}_t)^2 = \mathbf w'\big(\sum_t \mathbf e_t\mathbf e_t'\big)\mathbf w = T_w\,\mathbf w'\hat{\boldsymbol\Sigma}_e\mathbf w$ . Least squares with no intercept and a sum-to-one constraint therefore minimizes the sample version of the same objective — it is full-covariance Bates-Granger (for unbiased forecasts, whose error second-moment matrix is the error covariance matrix Bates and Granger work with) with $\hat{\boldsymbol\Sigma}_e$ plugged in, $\hat{\mathbf w} = \hat{\boldsymbol\Sigma}_e^{-1}\mathbf 1/(\mathbf 1'\hat{\boldsymbol\Sigma}_e^{-1}\mathbf 1)$ (Granger-Ramanathan’s restricted regression; Wang et al. eq. 2). Drop the off-diagonal elements of $\hat{\boldsymbol\Sigma}_e$ and you get Method 2 (with MSE on the diagonal); add $w_m \ge 0$ and you get Method 4. Not assessable.

Deeper Dive — Three Regressions, Enders’ Fix, and Robust Combinations

The three Granger-Ramanathan regressions (Wang et al. 2023, §2.2, eqs. 2–4): (2) no intercept, $\sum_m w_m = 1$ — full-covariance Bates-Granger; (3) no intercept, unconstrained; (4) intercept, unconstrained — our Method 3. Fit all three on the simulated weights window and score them on the holdout. For this comparison we drop locf from the pool: in the course DGP arima_2’s forecast is a constant plus $\hat\phi$ times locf, so with both in the pool the two columns already span a constant, and regression (3) would be regression (4) in disguise — the same fit, with the intercept hidden inside the arima_2 and locf weights.

pool5 <- setdiff(colnames(Fs_w), "locf")        # five tracks: no hidden constant
F5_w <- Fs_w[, pool5]; F5_h <- Fs_h[, pool5]
M5 <- ncol(F5_w)
# (2): equality constraint only, via solve.QP with meq = 1 and no inequalities
w2 <- solve.QP(crossprod(F5_w), crossprod(F5_w, ys_w),
               matrix(1, M5, 1), 1, meq = 1)$solution
# (3): no intercept, no constraints
w3 <- coef(lm(ys_w ~ F5_w - 1))
# (4): Method 3
w4 <- comb_gr(ys_w, F5_w)
gr3 <- cbind(`(2) sum-to-one` = w2, `(3) no intercept` = w3,
             `(4) GR` = w4$weights)
rownames(gr3) <- pool5
gr3_rmse <- c(rmse(ys_h, F5_h %*% w2), rmse(ys_h, F5_h %*% w3),
              rmse(ys_h, w4$intercept + F5_h %*% w4$weights))
round(rbind(w0 = c(0, 0, w4$intercept), gr3, holdout_RMSE = gr3_rmse), 3)

             (2) sum-to-one (3) no intercept  (4) GR
w0                    0.000            0.000 261.700
arima_1               0.650            0.372   0.779
arima_2               4.072            3.372   1.877
arima_3              -4.029           -3.415  -2.043
avg_k                 1.279            0.144  -0.833
meanf                -0.972            0.706 -50.827
holdout_RMSE          1.557            2.478   3.009

The three now genuinely differ: holdout RMSEs of 1.557, 2.478, and 3.009 for regressions (2), (3), and (4). In this realization each restriction removed let the weights chase the twelve weights-window points harder — and the holdout paid for it.

Enders’ fix. In his spread example (§13, eq. 2.76), Enders fits regression (2), finds negative weights, and successively drops the forecaster with the most negative weight and re-estimates until every weight is non-negative. Regression (2) has no intercept, so here we can start from the full six-track pool:

keep <- colnames(Fs_w)
repeat {
  Fk <- Fs_w[, keep, drop = FALSE]
  wk <- solve.QP(crossprod(Fk), crossprod(Fk, ys_w),
                 matrix(1, length(keep), 1), 1, meq = 1)$solution
  if (all(wk >= -1e-10)) break
  keep <- keep[-which.min(wk)]
}
setNames(round(wk, 3), keep)

arima_1   meanf 
  0.657   0.343

rmse(ys_h, Fs_h[, keep, drop = FALSE] %*% wk)

[1] 1.237391

It often lands close to CLS — both end with a few non-negative weights that sum to one — but it gets there by deleting forecasters one at a time instead of solving one constrained problem.

Robust combinations. The cross-forecaster median, and a trimmed mean that drops the highest and lowest of the six forecasts each date, estimate nothing at all — and, unlike the simple average, they are immune to a single broken track:

c(simple  = rmse(ys_h, rowMeans(Fs_h)),
  median  = rmse(ys_h, apply(Fs_h, 1, median)),
  trimmed = rmse(ys_h, apply(Fs_h, 1, function(v) mean(sort(v)[2:5]))))

  simple   median  trimmed 
1.284367 1.238565 1.278935

Keep the median and trimmed mean in mind for §7.8, where one track is broken.

One-line mentions. Eigenvector weights (Hsiao and Wan 2014), who also give conditions under which the simple average is optimal. Information-criterion weights (Enders §13, who writes SBC — the course’s BIC): $\hat w_m \propto \exp[(\text{BIC}_{\min} - \text{BIC}_m)/2]$ , normalized to sum to one, with $\text{BIC}_{\min}$ the best model’s value. Shrinkage of estimated weights toward $1/M$ (Diebold and Pauly 1990; Stock and Watson 2004) — a formal version of “estimate a little, not a lot.”

None of this is assessable.

7.7 The Assumption-Break: Weights That Win the Window They Were Fit On

Concept: A Combination Method Is an Estimator

Theory says Granger-Ramanathan should win. It minimizes the weights-window sum of squared errors over all intercepts and weights, so no other set of weights — equal, inverse-MSE, or constrained — can fit that window better. But “the fit” is in-sample. Granger-Ramanathan’s weights are estimates, made from twelve observations, and estimates carry error. The assumption we are about to break is the one buried in every optimal-weight formula: that you know the weights. Shrink the weights window and watch what happens.

The DGP was built so the experiment is fair. If equal weights were already optimal, every estimated weight would be pure noise and the simple average would win trivially. So the course DGP has a clearly dominant track — arima_1, the correctly specified model (its coefficients still estimated, so it is not an oracle) — and the optimal weights are far from $1/M$ . A weighting scheme that finds arima_1 should beat the simple average. The question is whether it can find it with twelve observations.

Stage 1 — One Realization, $T_w = 12$

Score the weights from §7.6 on both windows:

apply_w <- function(w, F) as.numeric(w$intercept + F %*% w$weights)
stage1 <- rbind(
  cbind(weights_window = apply(Fs_w, 2, rmse, actual = ys_w),
        holdout        = apply(Fs_h, 2, rmse, actual = ys_h)),
  t(sapply(ws_sim, function(w) c(rmse(ys_w, apply_w(w, Fs_w)),
                                 rmse(ys_h, apply_w(w, Fs_h))))))
round(stage1, 3)

        weights_window holdout
arima_1          0.706   1.120
arima_2          0.979   1.208
arima_3          1.115   1.372
locf             1.208   1.301
avg_k            0.922   1.793
meanf            0.769   1.701
SA               0.877   1.284
IMSE             0.804   1.281
CLS              0.681   1.237
GR               0.282   3.009

Granger-Ramanathan has by far the best weights-window RMSE — 0.282, against 0.877 for the simple average and 0.706 for the correctly specified arima_1 alone. And on the holdout it posts 3.009: the worst number in the table. The weights printed in §7.6 explain why — an intercept of 262 and weights like -50.8 and 1.9 are not a discovery about the forecasters; they are a regression with almost as many parameters as data points, fitting the noise in twelve months. The other three methods sit close to one another, and CLS — which found arima_1 — is the best of them.

One realization could be a fluke. Stage 2 says it is not.

Stage 2 — Monte Carlo, $R = 200$

For each replication: simulate a fresh series, estimate the three ARIMAs once, build the six tracks, fit the four combinations on the weights window, and score the 12-date holdout. We record each method’s holdout RMSE relative to the simple average (below 1: beat the simple average). The same loop also runs Stage 3, so we set it up for weights windows $T_w \in \{12, 24, 60, 120\}$ now: the training window stays fixed at 240 and the ARIMAs are fit once; the weights window extends to $T_w$ ; the 12-date holdout follows it.

set.seed(6378)   # Stages 2-3 have their own seed (Stage 1 is the §7.1 realization, seed 6376)
R <- 200
Tw_grid <- c(12, 24, 60, 120)
methods <- c("IMSE", "CLS", "GR")
rel <- array(NA, c(R, length(Tw_grid), length(methods) + 1),
             dimnames = list(NULL, paste0("Tw=", Tw_grid), c(methods, "arima_1 alone")))

for (r in seq_len(R)) {
  y    <- course_dgp(n_train + 120 + 12)
  fits <- fit_sim(upto(y, n_train))
  Fa   <- tracks_on(y, fits, n_train + 1:132)     # all one-step tracks
  ya   <- as.numeric(y)[n_train + 1:132]
  for (j in seq_along(Tw_grid)) {
    wi <- 1:Tw_grid[j]; hi <- Tw_grid[j] + 1:12   # weights window, then holdout
    ws <- list(SA   = comb_simple(ya[wi], Fa[wi, ]),
               IMSE = comb_inv_mse(ya[wi], Fa[wi, ]),
               CLS  = comb_cls(ya[wi], Fa[wi, ]),
               GR   = comb_gr(ya[wi], Fa[wi, ]))
    r_h <- sapply(ws, function(w) rmse(ya[hi], apply_w(w, Fa[hi, ])))
    rel[r, j, ] <- c(r_h[methods], rmse(ya[hi], Fa[hi, "arima_1"])) / r_h["SA"]
  }
}

summ <- function(x) c(median = median(x), mean = mean(x),
                      p90 = unname(quantile(x, 0.9)), share_beating_SA = mean(x < 1))
stage2 <- t(apply(rel[, "Tw=12", methods], 2, summ))
# display copy: ratios to 2 dp, the share as a percentage to 1 dp (R = 200, so exact)
data.frame(round(stage2[, c("median", "mean", "p90")], 2),
           pct_beating_SA = sprintf("%.1f%%", 100 * stage2[, "share_beating_SA"]))

     median mean  p90 pct_beating_SA
IMSE   0.99 0.98 1.01          74.5%
CLS    0.95 0.95 1.12          64.0%
GR     1.68 2.22 4.07          16.0%

s2 <- data.frame(method = factor(rep(methods, each = R), levels = methods),
                 rel = c(rel[, "Tw=12", methods]))
ggplot(s2, aes(method, rel)) +
  geom_boxplot(fill = "grey90", outlier.alpha = 0.4) +
  geom_hline(yintercept = 1, colour = "firebrick", linetype = "dashed") +
  scale_y_log10() +
  labs(x = NULL, y = "Holdout RMSE / simple-average RMSE (log scale)",
       title = "Stage 2: 200 replications, weights fit on T_w = 12",
       subtitle = "Below the dashed line: beat the simple average") +
  theme_bw()

Granger-Ramanathan’s long right tail is the puzzle in one picture. Its median holdout RMSE is 1.68 times the simple average’s, its 90th percentile 4.07 times, and it beats the simple average in only 16.0% of replications. So the Stage 1 realization was typical: GR lost to the simple average 84.0% of the time.

The other two methods tell the rest of the story. Inverse-MSE is the simple average with a slight tilt (median 0.99): it leans toward arima_1 and away from the weak tracks, estimating so little that it almost never does much worse. Constrained LS sits between the two on risk — a thin right tail — and in this DGP it already beats the simple average at the median (0.95), because its constraints shrink the estimation error while still letting it find the dominant track. That is evidence for the estimation-error explanation of the puzzle, not against it.

Stage 3 — The Sweep: How Long a Window Do Estimated Weights Need?

stage3_med  <- apply(rel, c(2, 3), median)
stage3_mean <- apply(rel, c(2, 3), mean)
round(stage3_med, 2)

       IMSE  CLS   GR arima_1 alone
Tw=12  0.99 0.95 1.68          0.90
Tw=24  0.98 0.93 1.12          0.90
Tw=60  0.99 0.94 0.98          0.92
Tw=120 0.98 0.90 0.93          0.90

s3 <- data.frame(Tw = rep(Tw_grid, 4),
                 method = rep(colnames(stage3_med), each = length(Tw_grid)),
                 med = c(stage3_med))
ggplot(s3, aes(Tw, med, colour = method)) +
  geom_line(linewidth = 1) + geom_point(size = 2.5) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  scale_x_continuous(breaks = Tw_grid) +
  labs(x = "Weights-window length T_w", y = "Median holdout RMSE / simple average",
       colour = NULL,
       title = "Stage 3: Granger-Ramanathan's weights earn their keep only with a longer window") +
  theme_bw() + theme(legend.position = "bottom")

(A note on design: in this sweep the holdout dates move with $T_w$ , a deliberate departure from the PS3 split; the shortcut saves refitting the ARIMAs four times per replication. Because the DGP is stationary, moving the holdout does not change the distribution of $y_t$ . It is not completely innocuous, though: a later holdout also lengthens the history behind meanf, the expanding mean, so its forecasts are a little more precise at $T_w = 120$ than at $T_w = 12$ — even here. The ARIMA tracks, whose coefficients stay fixed, are unaffected.)

As the window grows, Granger-Ramanathan’s median falls from 1.68 at $T_w = 12$ to 1.12 at 24 and 0.98 at 60 — crossing the simple average between $T_w = 24$ and $60$ — and 0.93 at 120. It keeps improving but never quite catches the correctly specified arima_1 alone (0.90) or CLS (0.90). The crossover is real because the DGP was built so that the optimal weights are far from $1/M$ ; with enough data, estimating them pays. Note that the longer window is a Granger-Ramanathan requirement: CLS (0.95) and inverse-MSE (0.99) already sit below the simple average at $T_w = 12$ .

Reproducibility note. Every simulation in these notes runs inline and is seeded, so every render is identical: seed 6376 for the course realization of §7.1 (Stage 1 and the §7.3, §7.4, and §7.6 examples), 6377 for the §7.4 power table, 6378 for Stages 2–3, and 6379 for the §7.5 variance check. Every simulation or data number quoted in the prose is printed by the render from the objects in the code, not typed by hand, so if you change a seed or a setting, those numbers update with the code. Two kinds of statement do not update by themselves: the words around the numbers (a “crossover” or a “rejection” was checked against the render for these seeds, and could need rewording under others), and the forward references in §7.4, which describe §7.8’s results (its rejection on twelve dates, and its check of the h = 1 assumption) in words because they have not yet been computed at that point.

Three Lessons

A combination method is an estimator. Its holdout performance is the population gain from its weights minus the estimation error in them. The weights window is where the estimation error is paid.
The puzzle is a comparison: attainable gain against estimation cost (Smith and Wallis 2009; Claeskens, Magnus, Vasnev and Wang 2016; Wang et al. 2023, §2.6). It bites when the improvement that optimal weights could deliver over equal weights is small relative to the cost of estimating them. Near-equal optimal weights are the important special case — then there is almost nothing to gain — but not the only one: the optimal weights can differ substantially from $1/M$ and still buy little when the forecasts are largely redundant (Elliott 2011, reviewed in Wang et al. 2023, §2.6; on redundant pools, §2.5). With twelve weights-window observations, the estimation cost is high, which describes most applied work — and Problem Set 3.
Restrictions are a dial between the two. $1/M$ estimates nothing, inverse-MSE a little, CLS more inside a box, Granger-Ramanathan everything. Moving down the list pays only when the window is long or the gain is large.

The Forecast Combination Puzzle

Stock and Watson (2004) gave it a name: theory says optimally weighted combinations should dominate; in practice, the simple average is very hard to beat. The simulation shows the explanation (Timmermann 2006; Smith and Wallis 2009): the optimal weights have to be estimated, and with short windows the estimation error swamps the gain. The governing comparison is the attainable improvement over equal weights versus the cost of estimating the weights. The most important case in which the gain is small is optimal weights that are near equal, so that there is little to swamp; but it is a special case, not a necessary condition — substantially different optimal weights can still yield little gain when the forecasts are redundant, carrying much the same information (Elliott 2011, reviewed in Wang et al. 2023, §2.6; on redundant pools, §2.5). Enders’ spread example illustrates the near-equal case: his seven forecasts had nearly identical error variances, and on its single forecast date the simple and variance-weighted averages landed closest — one date, so an illustration rather than evidence.

The index fund analogy. Active managers (estimated weights) should beat the index (equal weights) in theory — they pick the winners. In practice, fees eat the edge, and here the fee is estimation error. The index is very hard to beat. But an index fund that holds a bankrupt stock at full weight is not a good index. Equal weights are robust to estimation error, not to a broken member of the pool. Section 7.8 shows what that means on real data.

Looking Ahead — Rolling-Window Combination Weights

Everything here estimated the weights once, on a fixed window. An alternative re-estimates them on a moving window as each holdout date arrives (Stock and Watson 1998), and a related literature lets the weights drift over time (Diebold and Pauly 1987). That trades estimation error — each re-fit uses only the latest window — for adaptivity to structural change: the combination-weight analogue of the fixed-width rolling scheme in §7.2. It is not built in this course. Not assessable.

7.8 The PS3 Workflow on a Series That Isn’t Yours

This is the complete Problem Set 3 workflow, run end to end on UNRATE: $N_{\text{total}} =$ 720 months (1960–2019), $T =$ 696 (train ends December 2017), weights window = 2018, holdout = 2019, and the six tracks built in §7.1. The metrics table and the DM test of arima_1 against locf are already done (§7.3 and §7.4). What remains is the combination step.

Step 1: Fit the Four Methods on 2018

ws_un <- list(SA = comb_simple(y_w, F_w), IMSE = comb_inv_mse(y_w, F_w),
              CLS = comb_cls(y_w, F_w),   GR = comb_gr(y_w, F_w))
w_tab <- sapply(ws_un, function(w) c(w0 = w$intercept, w$weights))
round(w_tab, 3)

           SA  IMSE  CLS       GR
w0      0.000 0.000 0.00 -141.450
arima_1 0.167 0.235 0.69    1.224
arima_2 0.167 0.158 0.00   -5.825
arima_3 0.167 0.229 0.00    1.088
locf    0.167 0.179 0.00    5.098
avg_k   0.167 0.198 0.31   -2.605
meanf   0.167 0.001 0.00   24.803

kappa_w <- kappa(cbind(1, F_w), exact = TRUE)   # condition number, GR design matrix
kappa_w

[1] 41619.94

Read across the rows.

Inverse-MSE gives meanf a weight of 0.001 — effectively zero — and spreads the rest between 0.16 and 0.24 across the five credible tracks.
CLS puts 0.69 on arima_1 and 0.31 on avg_k, and zero everywhere else.
Granger-Ramanathan returns an intercept of -141 and a weight of 24.8 on meanf, with large offsetting weights elsewhere. Why? Over twelve months, the expanding mean of more than 696 observations barely moves, so the meanf column is almost a constant — almost collinear with the intercept. The condition number of the regression’s design matrix (intercept plus six tracks; kappa(..., exact = TRUE)) is about 41,620; anything in the thousands means the least-squares solution is extremely sensitive to the data. The regression is trading an enormous intercept against an enormous meanf weight to fit twelve points.

Step 2: Score Everything on 2019

comb_h <- sapply(ws_un, function(w) apply_w(w, F_h))   # N x 4 combined forecasts
rmse_un <- c(apply(F_h, 2, rmse, actual = y_h),
             apply(comb_h, 2, rmse, actual = y_h))
round(rmse_un, 3)

arima_1 arima_2 arima_3    locf   avg_k   meanf      SA    IMSE     CLS      GR 
  0.097   0.097   0.096   0.096   0.108   2.315   0.416   0.093   0.094   0.463

The simple average has a holdout RMSE of 0.416 — roughly four times the credible single tracks’ 0.096 to 0.108. Inverse-MSE (0.093) and CLS (0.094) match or edge out the best single tracks. Granger-Ramanathan (0.463) is worse even than the simple average.

Step 3: Test Each Estimated-Weight Combination Against the Simple Average

Three tests, declared in advance and all reported. We do not pick the holdout winner and then test it on the same twelve dates — that is choosing the hypothesis after seeing the data.

ec_h <- y_h - comb_h
dm_comb <- t(sapply(c("IMSE", "CLS", "GR"), function(m) {
  tst <- dm.test(ec_h[, m], ec_h[, "SA"], h = 1, power = 2)
  c(DM = unname(tst$statistic), p_value = unname(tst$p.value))
}))
signif(dm_comb, 3)

        DM  p_value
IMSE -9.10 1.87e-06
CLS  -9.12 1.83e-06
GR    0.87 4.03e-01

Inverse-MSE against the simple average: $\text{DM} =$ -9.1, $p =$ 1.9e-06. CLS: $\text{DM} =$ -9.1, $p =$ 1.8e-06. A clear rejection at $N = 12$ — the §7.4 power lesson from the other side. A gap this large relative to the noise in $d_t$ shows up even in twelve months. Granger-Ramanathan against the simple average: $p =$ 0.40, no detectable difference. Like every h = 1 result in this module, these $p$ -values assume that each $d_t$ is serially uncorrelated; the Technical Note below checks that assumption, and The Honest Reading reports what it found.

Technical Note — Checking the `h = 1` Assumption on $d_t$

The four PS3-design UNRATE tests were run at h = 1: arima_1 against locf (§7.4) and the three combinations against the simple average (above). For each, we look at the autocorrelations of the squared-error loss differential $d_t$ at lags 1–3, with a Ljung-Box test of all three, and then re-run the DM test with a Bartlett (Newey-West) variance that allows for serial correlation in $d_t$ up to lag $h - 1$ , for h = 2, 3, 4. (We use "bartlett" rather than the default truncated "acf" sum because the truncated sum can turn negative on twelve points, in which case dm.test() warns and falls back to h = 1 (the output then reports horizon 1).)

# The h = 1 DM tests assume d_t is serially uncorrelated. Check it.
dmdep_pairs <- list(                          # (error 1, error 2) for each test
  "arima_1 vs locf" = cbind(e_h[, "arima_1"], e_h[, "locf"]),
  "IMSE vs SA"      = cbind(ec_h[, "IMSE"],   ec_h[, "SA"]),
  "CLS vs SA"       = cbind(ec_h[, "CLS"],    ec_h[, "SA"]),
  "GR vs SA"        = cbind(ec_h[, "GR"],     ec_h[, "SA"]))

# (1) ACF of d_t = e1^2 - e2^2 at lags 1-3, and a Ljung-Box test of all three
dmdep_acf <- t(sapply(dmdep_pairs, function(e) {
  d <- e[, 1]^2 - e[, 2]^2
  r <- acf(d, lag.max = 3, plot = FALSE)$acf[2:4]
  c(lag1 = r[1], lag2 = r[2], lag3 = r[3],
    LB_p = Box.test(d, lag = 3, type = "Ljung-Box")$p.value)
}))

# (2) Sensitivity: DM p-values with a Bartlett (HAC) variance, lags up to h - 1
dmdep_sens <- t(sapply(dmdep_pairs, function(e) {
  p <- function(h) unname(dm.test(e[, 1], e[, 2], h = h, power = 2,
                                  varestimator = "bartlett")$p.value)
  c(h1_lag0 = p(1), bartlett_h2 = p(2), bartlett_h3 = p(3), bartlett_h4 = p(4))
}))
stopifnot(all.equal(unname(dmdep_sens[, "h1_lag0"]),        # column 1 reproduces
                    unname(c(dm_un$p.value, dm_comb[, "p_value"]))))  # the tests above
round(dmdep_acf, 3)

                  lag1   lag2   lag3  LB_p
arima_1 vs locf -0.326 -0.113 -0.082 0.580
IMSE vs SA      -0.200 -0.279  0.026 0.587
CLS vs SA       -0.223 -0.253  0.028 0.605
GR vs SA         0.097 -0.047  0.050 0.973

signif(dmdep_sens, 3)

                 h1_lag0 bartlett_h2 bartlett_h3 bartlett_h4
arima_1 vs locf 9.02e-01    8.91e-01    8.85e-01    8.79e-01
IMSE vs SA      1.87e-06    1.54e-06    6.33e-07    6.07e-07
CLS vs SA       1.83e-06    1.30e-06    5.47e-07    5.43e-07
GR vs SA        4.03e-01    4.64e-01    5.07e-01    5.56e-01

Three readings.

No autocorrelation is detectable — but twelve points cannot detect much. The Ljung-Box $p$ -values run from 0.58 to 0.97. The sample autocorrelations are not all small, though: the lag-1 values run from -0.33 to 0.10, and with $N = 12$ an autocorrelation of that size is well inside the noise. The check is consistent with uncorrelated $d_t$ ; it does not establish it.
The verdicts survive the robust variance. arima_1 against locf has $p$ between 0.88 and 0.89 across the three Bartlett variances; inverse-MSE and CLS against the simple average stay below 2e-06; Granger-Ramanathan stays between 0.46 and 0.56. None of the four verdicts at 5% changes.
Direction. Where the lag-1 autocorrelation is negative, the robust $p$ -values here came out smaller than the lag-0 ones, so in these data the lag-0 standard error was, if anything, conservative. That is an observation about these four tests, not a theorem — raising h also changes the Harvey-Leybourne-Newbold small-sample factor.

Not assessable beyond the Core statement in §7.4: at h = 1, dm.test() assumes serially uncorrelated $d_t$ .

The Honest Reading

This is not the textbook story, and the notes would be dishonest to pretend otherwise.

The simple average did not win, and the reason is not that the estimated weights got lucky. Here the attainable gain over equal weights is large enough to pay for estimating the weights. The pool contains meanf, a forecaster with no business on an I(1) series, and $1/M$ gives it a sixth of the vote — enough to drag the average about 0.39 percentage points off the actuals every month. Equal weights are robust to estimation error, not to a broken member. Inverse-MSE and CLS won by down-weighting meanf, which is exactly what estimated weights are for: the gain from sidelining one broken member dwarfs the estimation cost of these restricted weights.
The tests rest on an assumption about $d_t$ , and the check qualifies rather than confirms it. Every DM test here used h = 1, which assumes serially uncorrelated loss differentials. The autocorrelations of the four $d_t$ series are not detectably different from zero (Ljung-Box $p \ge$ 0.58), but twelve dates give that check little power; re-running the tests with a serial-correlation-robust variance, no verdict changes (Technical Note above). So the UNRATE verdicts are robust to that check — not definitive.
Granger-Ramanathan lost for the Stage 1 reason, amplified by near-collinearity: seven numbers from twelve observations, two of which — the intercept and the nearly constant meanf column — are almost the same regressor.
Is estimation the only way to sideline a broken track? That question goes beyond Core; the optional Deeper Dive at the end of this section takes it up.
Caveats. One window. 2018–2019 was a calm glide with no turning point, so every credible track had small errors. A different window could rank the methods differently. One realization is one draw; the Monte Carlo of §7.7 is the distribution.

What this means for Problem Set 3. Look at your six tracks before you combine them. If one of them is broken on your series — meanf on an I(1) series is the usual suspect — say so, explain why from the series’ properties, and expect it to drag the simple average. Your scorecard still reports every track and every combination; the reading is where the judgment goes.

Deeper Dive — Robust Combinations on UNRATE

c(simple       = rmse(y_h, rowMeans(F_h)),
  median       = rmse(y_h, apply(F_h, 1, median)),
  trimmed      = rmse(y_h, apply(F_h, 1, function(v) mean(sort(v)[2:5]))),
  simple_five  = rmse(y_h, rowMeans(F_h[, credible])))

     simple      median     trimmed simple_five 
 0.41601502  0.09869660  0.09751126  0.09250425

The median and the trimmed mean (drop the highest and lowest forecast each month) land with the credible tracks, without estimating a single weight: meanf is always the highest forecast, so trimming throws it out. The last entry — a simple average over the five credible tracks — makes the point directly: once the broken member is out of the pool, equal weights are as good as anything else on this window. Not assessable.

Common Pitfalls and Misconceptions

“MASE $< 1$ means my model beat the naive.” No. MASE compares your holdout error with the naive’s training error. The head-to-head on the same dates is $U_2 < 1$ . On UNRATE 2019, every credible track had MASE at or below 0.71 and $U_2$ between 0.999 and 1.12.
“DM failed to reject, so the models are equally good.” It means the gap is not detectable with these dates. At $N = 12$ , DM rarely detects even a correctly specified estimated model’s advantage over a misspecified one.
“DM always fails with twelve observations.” No. Power depends on the gap relative to the noise in $d_t$ . On UNRATE, inverse-MSE against the simple average rejected with $p =$ 1.9e-06 at $N = 12$ .
“I can feed residuals(fit) to dm.test().” Out-of-sample forecast errors only. In-sample residuals reward the model that overfits.
“Granger-Ramanathan weights are optimal, so use them.” Optimal in population. With seven parameters from twelve observations, mostly noise: they win the window they were fit on and lose the next one.
“A negative weight is a bug.” It is what unconstrained least squares does with correlated forecasts. It is a warning about estimation error, not an insight about the forecaster. (And an NA weight is exact collinearity — treat it as zero.)
“The simple average always wins.” Not when the pool contains a broken member — meanf on an I(1) series. Equal weights are robust to estimation error, not to a bad forecaster.
“avg_k is a moving-average model.” It is a $k$ -period average forecast. In this course a moving average is the MA( $q$ ) block, $\theta_l\epsilon_{t-l}$ .
“Re-fit the ARIMA on y[1:(T+12)] before forecasting the holdout.” Parameters are estimated once, on train; Arima(y, model = fit) reuses them while the information set rolls.
Subsetting a ts with square brackets. y[1:T] drops the frequency; a fresh auto.arima() on it never searches seasonal models. Use window().

Connection to Enders and Other Texts

Diebold-Mariano: Enders Chapter 2, Section 9 (“Properties of Forecasts”), eqs. (2.61)–(2.62): the loss differential, the truncated long-run variance, the Harvey-Leybourne-Newbold version with a $t_{H-1}$ reference, the regress-on-a-constant shortcut, and the Clark-McCracken nested-model caveat with the Clark-West adjustment.
Combining forecasts: Enders Chapter 2, Section 13, eqs. (2.71)–(2.76): the two-forecast variance, the optimal and inverse-variance weights, the Granger-Ramanathan regression, and the spread example, which illustrates the near-equal-weights case of the puzzle (seven nearly identical error variances); for one forecast date, the simple and variance-weighted averages landed closest, ahead of the regression and SBC (= BIC) weights.

Convention flags when you cross-reference Enders.

Enders	Course	Note
$f_{it}$	$\tilde y^{(m)}_{t\mid t-1}$	individual forecasts
$\alpha_0, \alpha_1, \ldots, \alpha_n$	$w_0, w_1, \ldots, w_M$	GR regression weights; $\alpha$ is the master-equation intercept in this course
$n$ forecasts	$M$ forecasts
$H$ (number of forecast errors)	$N$	evaluation-window size
$e_{it} = y_t - f_{it}$ , then $f_{it} - y_t$ in the spread example	$e_{m,t} = y_t - \tilde y^{(m)}_{t\mid t-1}$ throughout	sign flips do not change variances or weights
$a_i$ , characteristic roots inside	$\phi_j$ , roots of $\Phi(L)$ outside	as in Modules 3–6
Harvey, Leybourne & Newbold cited as (1998)	Harvey, Leybourne & Newbold (1997)	Enders’ DM discussion cites a different HLN paper (1998); the small-sample correction used by `dm.test()` is from HLN (1997), IJF
SBC	BIC	same criterion, different name

Hyndman (2014 Forecasting: Principles and Practice workshop notes): §2.6 (“Evaluating forecast accuracy”) defines MASE and §5.1 (“Cross-validation”) the rolling forecasting origin used in the §7.2 figure. Hyndman and Athanasopoulos’s textbook covers the same ground in its chapter on evaluating forecast accuracy; note that the forecast package’s accuracy() defaults to the seasonal MASE scale for seasonal ts objects (§7.3).

Wang, Hyndman, Li and Kang (2023), “Forecast combinations: an over 50-year review” (International Journal of Forecasting): §2.2 for the optimal weights and the three Granger-Ramanathan regressions (eqs. 1–4), §2.5 on which forecasts to combine (diversity and redundant pools), §2.6 for the combination puzzle. They write $N$ for the number of forecasts — the course’s $M$ .

Elliott, G. (2011), “Averaging and the optimal combination of forecasts,” working paper, University of California, San Diego: bounds on the gain of optimal weights over equal weights, and why that gain is often too small to cover estimation error (reviewed in Wang et al. 2023, §2.6).

On shrinking weights toward $1/M$ : Diebold and Pauly (1990), “The use of prior information in forecast combination,” International Journal of Forecasting 6(4), 503–508; and Stock and Watson (2004). Diebold (1989) is the reference for combination and forecast encompassing.

Hamilton (Time Series Analysis), §10.5 (“The Sample Mean of a Vector Process”), develops the long-run variance of a sample mean and the Newey-West estimator (eq. 10.5.15) behind the Technical Note in §7.4.

Practice Problems

Core Practice

On the course DGP with a different seed, build the six tracks and report RMSE, MAE, MASE, and $U_2$ on the holdout. Find a track with MASE $< 1$ and $U_2 > 1$ (or explain why none exists), and say in one sentence what each of the two numbers means for that track.
On the same realization, run dm.test() for arima_1 against each of the other five tracks. Which differences are significant at 5%? Relate your answer to the size of each RMSE gap.
Re-run the §7.4 power Monte Carlo with power = 1 (absolute loss). Does DM gain or lose power at each $N$ ? Why might it?
Two forecasts have error variances 1 and 4 and error correlation $\rho_{1,2} = 0.25$ . Compute the optimal two-forecast weight and the combined error variance; compare with the simple average and with the zero-covariance (inverse-variance) weight.
Compute all four combinations on the course DGP with $T_w = 12$ and with $T_w = 60$ . Report the weights and the holdout RMSEs. In which case do the Granger-Ramanathan weights look sensible?
On UNRATE, repeat §7.8 with avg_k at $k = 6$ and $k = 12$ (use the k argument of tracks_on()). Do the combination rankings change? Justify your choice of $k$ from the series’ persistence.
Conceptual. A colleague’s Granger-Ramanathan combination has an intercept of -141 and a weight of 25 on one forecast. Explain what produced it, and what you would report instead.
Conceptual. In three to four sentences, explain why the simple average is hard to beat, and name one situation in which it is not. Use the phrase “estimation error” once.

Beyond Core (optional)

Re-run the tsCV() Deeper Dive (§7.2) and compare the DM verdict on the ten years of re-estimated origins with the 12-date verdict of §7.4. Why can they differ?
Implement Enders’ drop-the-most-negative rule on UNRATE and compare it with CLS.
Add the median and trimmed mean to the UNRATE holdout table and test each against the simple average with DM.

Key Takeaways

Many origins, one horizon. The PS3 design estimates each model once and rolls the information set: every forecast is one step ahead, and the three $k$ -period average forecasts (locf, avg_k, meanf) sit beside three ARIMA tracks.
$U_2$ asks “did you beat the naive here?”; MASE asks “how big was your miss relative to the naive’s typical miss in training?” Report both, and never read MASE $< 1$ as a win.
A metrics table is a ranking, not a verdict. The Diebold-Mariano test is a $t$ -test on the loss differentials $d_t$ , with a standard error that allows for serial correlation; at h = 1, dm.test() uses the lag-0 version, which assumes $d_t$ is serially uncorrelated. Whether it rejects depends on the gap relative to the noise in $d_t$ — and twelve dates buy little power. On UNRATE 2019 the five credible tracks are descriptively close (RMSE), and the one track-vs-track comparison tested (arima_1 vs locf) does not reject.
Combine because errors diversify. With zero error covariance, the optimal weight is proportional to the inverse error variance.
A combination method is an estimator. With twelve weights-window observations, estimate as little as you can: Granger-Ramanathan’s seven numbers win the window and lose the holdout.
Equal weights are the benchmark to beat — when the pool is sound. The forecast combination puzzle is the cost of estimating the weights outweighing the attainable gain over equal weights — most often because the optimal weights are near equal or the forecasts are redundant. A broken pool member (like meanf on an I(1) series) makes the attainable gain large, so estimated weights (or, see the Deeper Dive, robust ones) can pay for themselves.
Every forecast here used only the series’ own past. Module 8 asks whether another variable can help: the ADL model.

Looking Ahead — Module 8

Look back at what every forecast in this module used: the series’ own past, and nothing else. Three ARIMA models, three averages of past values, and four ways of combining them — and on UNRATE in 2019 the five credible tracks were descriptively close in RMSE, the one track-vs-track comparison tested (arima_1 vs locf, carrying the last value forward) did not reject, and the best combinations only matched the best single tracks’ RMSE. Combinations can only reshuffle the information already in the pool. The combination puzzle hands us a sharper question: what if an exogenous regressor could improve the forecast? Could another variable — initial unemployment-insurance claims, say — carry information that no univariate track has?

Module 8 (ADL Models and Dynamic Multipliers) adds that information directly. It turns on a new slider on the mixing board — $\beta x_t$ and its lags — and asks how a change in $x$ today works its way into $y$ over time. The theme of this module was refereeing between univariate forecasts. The theme of the next is letting a second variable into the room.

How to Use These Notes

Core Learning Objectives

7.1 Where We Are on the Mixing Board

The Module 6 Handoff

The PS3 Design: Many Origins, One Horizon

Six Forecast Tracks

Loading UNRATE and Building the Tracks

The Course Simulation

7.2 Three Ways to Run an Evaluation

Beyond the Fixed-Parameter Scheme

Deeper Dive — tsCV(): Re-Estimating at Every Origin

7.3 Relative Metrics: Theil’s U2U_2 and MASE

Concept: Two Different Questions

Math

Simulate: The Metrics Table on the Course DGP

Real Data: The Metrics Table on UNRATE

Other Theil Statistics, and accuracy() Defaults

Deeper Dive — Theil’s U1U_1

7.4 The Diebold-Mariano Test

Concept: Scoring Is Not Testing

Math

Worked Example: DM Really Is a tt-Test

Simulate — The Assumption-Break: Make the Sample Tiny

Real Data: Closing Module 6’s Open Question

Technical Note — Why σ̂d‾\hat\sigma_{\bar d} Must Allow for Serial Correlation

Deeper Dive — Nested Models: Clark-McCracken and the Clark-West Fix

7.5 Why Combine?

Concept: From Selection to Combination

Math: Two-Forecast Portfolio Algebra

7.6 Four Combination Methods

Method 1 — Simple Average

Method 2 — Inverse-MSE Weights

Method 3 — Granger-Ramanathan (Unconstrained Least Squares)

Method 4 — Constrained Least Squares

Summary

Technical Note — Full-Covariance Bates-Granger Is a Restricted Regression

Deeper Dive — Three Regressions, Enders’ Fix, and Robust Combinations

7.7 The Assumption-Break: Weights That Win the Window They Were Fit On

Concept: A Combination Method Is an Estimator

Stage 1 — One Realization, Tw=12T_w = 12

Stage 2 — Monte Carlo, R=200R = 200

Stage 3 — The Sweep: How Long a Window Do Estimated Weights Need?

Three Lessons

The Forecast Combination Puzzle

Looking Ahead — Rolling-Window Combination Weights

7.8 The PS3 Workflow on a Series That Isn’t Yours

Step 1: Fit the Four Methods on 2018

Step 2: Score Everything on 2019

Step 3: Test Each Estimated-Weight Combination Against the Simple Average

Technical Note — Checking the h = 1 Assumption on dtd_t

The Honest Reading

Deeper Dive — Robust Combinations on UNRATE

Common Pitfalls and Misconceptions

Connection to Enders and Other Texts

Practice Problems

Core Practice

Beyond Core (optional)

Key Takeaways

Looking Ahead — Module 8

Deeper Dive — `tsCV()`: Re-Estimating at Every Origin

7.3 Relative Metrics: Theil’s $U_2$ and MASE

Other Theil Statistics, and `accuracy()` Defaults

Deeper Dive — Theil’s $U_1$

Worked Example: DM Really Is a $t$ -Test

Technical Note — Why $\hat\sigma_{\bar d}$ Must Allow for Serial Correlation

Stage 1 — One Realization, $T_w = 12$

Stage 2 — Monte Carlo, $R = 200$

Technical Note — Checking the `h = 1` Assumption on $d_t$