Module 7 Lab: Forecast Evaluation and Combinations

Econ 6376 — Applied Time Series Econometrics · AI Learning Pack

How to use this lab

Everything here runs as shipped. No API key, no internet, no blanks to fill in. The point is not to get the code working. The point is to find out whether you can tell, before R tells you, which of two forecasts is really better, and whether you should be choosing between them at all.

The one rule from Module 1 still holds:

Write your prediction down before you run the chunk.

Not in your head. In LEARNING_LOG.md, in a comment, on paper. Being wrong on paper is the single most useful thing that can happen to you today.

If you are working with an AI tutor, it will ask for that prediction once, wait one turn, and run the code either way. It is making room for the “huh, I didn’t expect that” moment, because that is the one you remember in November.

Start with a question from class. Ask a room to guess a number, the height of a building or the beans in a jar, and the average guess usually beats most of the people who made one. Why would that be, and should it work for forecasts? Hold your answer; section 7 tests it.

Where you are on the mixing board: every slider of the univariate mean equation is already on, and nothing new switches on today. Module 6 pushed one fitted master equation forward in time. Module 7 puts several forecasts on the table at once (three fitted ARIMA models and three naive rules) and referees between them. Three questions: how do I score several forecasts, is the gap between two of them real, and should I pick one at all?

Notation for this module. From Module 7 on, a forecast carries a tilde: forecaster $m$’s one-step forecast of $y_t$ is $\tilde y^{(m)}_{t \mid t-1}$. A hat means a fitted value or an estimate ($\hat\phi$, and the estimated combination weights $\hat w_0$, $\hat w_m$). Module 6 wrote the same forecast object as $\hat y_{t+h \mid t}$; it is the same thing in new clothes. The forecast error is $e_{m,t} = y_t - \tilde y^{(m)}_{t \mid t-1}$: an out-of-sample $e$, observable once $y_t$ arrives, never an $\epsilon$.

---

1. Setup

Three packages and two course helpers. quadprog is new: it solves one constrained least-squares problem in section 8. It came with forecast, so you probably already have it.

library(ggplot2)
library(forecast)
library(quadprog)

source("R/forecast_utils.R")     # rmse(), mae(), mape(), theil_u2(), mase()
source("R/sarima_simulator.R")   # the Module 5 seasonal simulator

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

# The cached UNRATE pull: published (seasonally adjusted) unemployment, 1960-2019,
# the same window as Modules 4-6. The April 2020 spike would dominate everything.
unrate_raw <- read.csv("data/UNRATE.csv", stringsAsFactors = FALSE)
unrate <- ts(unrate_raw$UNRATE, start = c(1948, 1), frequency = 12)
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 

So what: 720 months of one series, and the last 24 of them are about to be split into two jobs: twelve months to fit combination weights and twelve months to score everything.

---

2. The PS3 design: many origins, one horizon

Module 6 stood at December 2017 and forecast 1, 2, …, 24 months ahead: one origin, many horizons. This module does the opposite: stand at each month in turn and forecast the next one. Many origins, one horizon ($h = 1$ throughout). It is the design the midterm specified and Problem Set 3 uses.

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

• Parameters are estimated once, on $y_1, \ldots, y_T$, and never again.
• The information set rolls. The forecast of $y_t$, $\tilde y_{t \mid t-1}$, uses actuals up to $y_{t-1}$. When $y_t$ arrives it becomes information for the next forecast, but the coefficients do not move.
• The weights window has $T_w = 12$ dates; the holdout has $N = 12$.

2a. One subsetting function, and one trap

# 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)
c(window_keeps = frequency(train_un), brackets_keep = frequency(unrate[1:T_origin]))

 window_keeps brackets_keep 
           12             1 

unrate[1:T] silently drops the monthly frequency. Refiltering with Arima(..., model = fit) survives that, because the fit stores its seasonal period. A fresh fit on such a subset does not: auto.arima() never even searches seasonal models on it.

2b. Three ARIMA tracks, estimated once

fit1 <- Arima(train_un, order = c(1, 1, 2),
              seasonal = list(order = c(1, 0, 1), period = 12))   # Module 5/6 model
fit2 <- Arima(train_un, order = c(1, 1, 0))                        # simpler
fit3 <- auto.arima(train_un)                                        # automatic pick
fits_un <- list(arima_1 = fit1, arima_2 = fit2, arima_3 = fit3)
as.character(fit3)

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

Modify and re-run

Fit auto.arima(unrate[1:T_origin]) instead, on the bracket subset. Before you run it: will it pick (a) the same model, (b) a model with no seasonal terms at all, or (c) an error? Then check frequency() of what you passed it.

2c. Six tracks

Three tracks are ARIMA models. Three are $k$-period average forecasts, $\tilde y_{t \mid t-1} = \frac{1}{k}\sum_{j=1}^{k} y_{t-j}$: $k = 1$ is locf (“last observation carried forward”, the naive random walk), $k = 3$ is avg_k, and $k = t-1$ is meanf, the expanding mean of everything so far. Never call these “moving averages”: in this course a moving average is the MA($q$) block, $\theta_l\epsilon_{t-l}$.

For the ARIMA tracks, Arima(series, model = fit) re-applies the training coefficients to a longer series, and fitted() then returns the one-step forecast at every date from information through $t-1$. The function is called fitted(), but out of sample these are forecasts, so they carry a tilde.

# k-period average forecast of y_t for each t in `dates`;
# 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.
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

(upto(x, max(dates)) is upto(unrate, T + 12) for the weights window and upto(unrate, T + 24) for the holdout. Extending the series never changes an earlier one-step forecast, because a forecast of $y_t$ uses only data before $t$.)

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

c(holdout = isTRUE(all.equal(F_h[, "locf"], as.numeric(unrate)[h_dates - 1])),
  weights = isTRUE(all.equal(F_w[, "locf"], as.numeric(unrate)[w_dates - 1])))

holdout weights 
   TRUE    TRUE 

Predict first

The next plot draws the 24 actuals and the one-step tracks. Before you look: how far apart will the five tracks other than meanf be? (a) Visibly different paths, tenths of a point apart; (b) nearly on top of one another, hundredths apart; (c) one ARIMA track clearly ahead of the rest. And where do you expect meanf, the sixty-year expanding mean, to sit?

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(legend.position = "bottom")

Five tracks hug the actuals and differ by hundredths. meanf is not on the chart: on an I(1) series the expanding mean is the average of a level that has wandered for sixty years, and it sits near 6.0% while unemployment runs below 4%. Remember that track. It decides section 10.

So what: you now hold a $T_w \times M$ and an $N \times M$ matrix of honest one-step forecasts, every column built from the same information at every date, which is what makes the columns comparable at all.

---

3. The course simulation: a world where you know the truth

Every topic here also runs on one simulated DGP, 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 mean is $\mu = 5$. The six simulated tracks mirror the UNRATE pool: arima_1 is ARIMA$(1,0,0)(1,0,0)_{12}$ with a mean (the true specification), arima_2 is ARIMA$(1,0,0)$ with a mean (it misses the seasonal block), arima_3 is ARIMA$(0,1,1)$ (a needless difference), and locf, avg_k ($k = 3$), meanf. Each ARIMA is estimated once on 240 training months.

The DGP is built so that one track is genuinely better: the correctly specified model, with its coefficients still estimated, so it is not an oracle. Section 9 explains why that design choice matters.

# 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]
isTRUE(all.equal(Fs_h[, "locf"], as.numeric(y_sim)[n_train + 12:23]))   # alignment

[1] TRUE

So what: two PS3-shaped datasets now sit side by side, one where the truth is known and one where it never will be. Every tool from here on runs on both.

---

4. Three ways to run an evaluation

As the forecast origin moves forward, do you re-estimate the model, and on which data? Three answers, three schemes:

Scheme	Parameters	Training window at origin $t$	Fits per model
Fixed-parameter one-step (PS3, the midterm, all of this module)	Estimated once on $1{:}T$	Fixed at $1{:}T$; the information set still rolls to $t-1$	1
Expanding window, re-estimated	Re-estimated at every origin	$1{:}(t-1)$	one per origin
Fixed-width rolling window	Re-estimated at every origin	$(t-W){:}(t-1)$	one per origin

Two sentences to keep. 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 the world changed, costly if it did not.

The fixed-parameter scheme is the only one this module builds: one fit per model and the same coefficients at every origin, so differences between tracks are differences between models, not between estimation samples.

Four-step note: this section is a definitional bridge, so it has no Simulate or Real Data step of its own. That is by design.

So what: “out-of-sample” is not one thing; when you read a forecast comparison, the first question is which of these three it ran.

---

5. Relative metrics: Theil’s $U_2$ and MASE

RMSE and MAE are in the units of $y$; MAPE is a percentage. None of them says whether 0.1 points is good. For that you need a relative score, and the two in common use answer different questions:

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

Over the holdout $t = T+13, \ldots, T+24$, with $\tilde y^{(\texttt{locf})}_{t \mid t-1} = y_{t-1}$:

$$U_2^{(m)} = \frac{\text{RMSE}^{(m)}}{\text{RMSE}^{(\texttt{locf})}}, \qquad \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 - y_{t-1}\big|}.$$

In $U_2$ both RMSEs cover the same holdout dates. In MASE the denominator is locf on the training sample only; the naive’s holdout errors never enter. The seasonal variant scales by $\frac{1}{T-s}\sum_{t=s+1}^{T}|y_t - y_{t-s}|$. One notation flag: the helper’s argument m in mase(..., m = 12) is that seasonal lag, which the math writes $s$, because $m$ is the forecaster index.

Compute both by hand once and make 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))))
round(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))), 4)

            U2   MASE
by_hand 1.0106 0.6424
helper  1.0106 0.6424

5a. The scorecard on the course simulation

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

Three passes. The correctly specified arima_1 has the smallest RMSE and MAE, and $U_2 =$ 0.86: on these dates it beat locf. arima_2 shows the metrics disagreeing: $U_2 =$ 0.93 says it beat the naive here, MASE $=$ 1.00 says its holdout error was about the size of the naive’s typical training miss. Both are true; they are different questions. And RMSE and MAE agree on the ranking here only because the gaps are large.

5b. The scorecard on UNRATE

Predict first

2019 was a calm glide in unemployment. For the five credible UNRATE tracks (everything but meanf), which do you expect? (a) MASE below 1 and $U_2$ below 1: the models beat the naive; (b) MASE below 1 but $U_2$ at or near 1; (c) MASE above 1 and $U_2$ above 1. Commit before you run it.

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

Every credible track has MASE between 0.58 and 0.71, comfortably below 1, while its $U_2$ sits between 0.999 and 1.124. 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). MASE below 1 here is a statement that 2019 was easy, not that the models are good. And meanf: RMSE 2.32, $U_2$ 24.2. Module 2’s lesson in new clothes: the series’ order of integration decides which baselines are even credible.

5c. accuracy() does not compute the course’s numbers

# 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)
round(c(acc["Test set", c("RMSE", "MASE", "Theil's U")],
        course_MASE_m1  = mase(test24, fc24$mean, train_un, m = 1),
        course_MASE_m12 = mase(test24, fc24$mean, train_un, m = 12)), 3)

           RMSE            MASE       Theil's U  course_MASE_m1 course_MASE_m12 
          0.405           0.416           3.972           2.510           0.416 

For a monthly ts, accuracy() scales MASE by the seasonal-naive training error, which matches m = 12, not the non-seasonal scale the course uses. On seasonally adjusted UNRATE a 12-month difference measures a year of drift, not seasonality. Its “Theil’s U” works with relative changes; it is a cousin of $U_2$, not the same number. Use mase() and theil_u2().

Modify and re-run

Change the $U_2$ benchmark in metrics_table() from locf to avg_k. Before you run it: which column changes, which stays put, and does any track’s $U_2$ cross 1?

So what: never report MASE alone and call it a win. $U_2$ is the race, MASE the yardstick, and on a calm stretch everything clears the yardstick.

---

6. The Diebold-Mariano test

A metrics table says which number is smaller. It does not say whether the gap is larger than the noise in twelve forecast errors would produce by chance. That was Module 6’s open question.

The taste test. One sip from each of two coffee shops is not a verdict. Twelve paired sips, one each morning, with a running tally of how much better shop A tasted each day: now you can ask whether the average difference is distinguishable from zero. That is a $t$-test on paired differences, and that is all DM is.

With out-of-sample errors $e_{1,t}$, $e_{2,t}$ over $N$ dates and a loss $g(\cdot)$, the loss differential is $d_t = g(e_{1,t}) - g(e_{2,t})$, and

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

where $\hat\sigma_{\bar d}$ is a standard error of $\bar d$ that allows for serial correlation in $d_t$. What h = 1 assumes: with h = 1, dm.test() uses the lag-0 variance only, so it assumes $d_t$ is serially uncorrelated. That is a condition on the loss differential, not merely on the errors: optimal one-step errors make it plausible, but $d_t$ is a nonlinear function of two error series and can be autocorrelated anyway.

Rules of use: out-of-sample errors only (never residuals(fit)); same dates, same loss; $\text{DM} < 0$ means forecast 1 had the smaller average loss. dm.test() also applies the Harvey-Leybourne-Newbold small-sample correction and uses a $t_{N-1}$ reference.

6a. DM really is a $t$-test

Predict first

On the simulation, arima_1 (the correctly specified model) won every metric in 5a. DM on its twelve holdout errors against locf: does it (a) reject at 5%, (b) fail to reject, or (c) depend on which loss you pick? And will the one-sample t.test() on $d_t$ give the same statistic as dm.test()?

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
c(dm.test = unname(dm_sim$statistic), t.test = unname(t.test(d_sim)$statistic))

   dm.test     t.test 
-0.8955949 -0.8955949 

DM $=$ -0.90 (negative: arima_1 had the smaller squared errors), $p =$ 0.390, and the one-sample $t$-test returns the same statistic. Even though arima_1 is the correctly specified model, at $N = 12$ the test does not reject at 5%. Bad luck, or typical?

6b. The assumption-break: make the sample tiny

The correctly specified model has a real expected-loss advantage over locf and over arima_2. How often does DM detect it with $N$ dates? 200 replications: draw a series, estimate the three ARIMAs once on 240 months, build one-step tracks over the next 120, and run DM on the first $N \in \{12, 24, 60, 120\}$. (About half a minute.)

Predict first

At $N = 12$, in what share of replications does DM reject in favour of arima_1 against locf? (a) Above 80%: it is the true model; (b) about half; (c) under 20%. And at $N = 120$?

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 favours 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
data.frame(pair = rownames(power_tab),
           matrix(sprintf("%.1f%%", 100 * power_tab), nrow = 2,
                  dimnames = list(NULL, colnames(power_tab))),
           check.names = FALSE)

                pair  N=12  N=24  N=60 N=120
1    arima_1 vs locf 12.5% 40.5% 76.0% 95.5%
2 arima_1 vs arima_2 10.0% 29.0% 60.5% 86.0%

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\nfavouring arima_1",
       colour = NULL,
       title = "DM rejection rates: a real gap, but N = 12 rarely shows it") +
  theme(legend.position = "bottom")

The table counts only rejections in the right direction (two-sided $p < .05$ and DM $< 0$), so under a calibrated null the benchmark is 2.5%, the dashed line, not 5%. At $N = 12$ DM rejects in favour of arima_1 in 12.5% of replications against locf and 10.0% against arima_2; at $N = 120$, 95.5% and 86.0%.

The lesson is power, with a calibration caveat. DM rejects when the gap is large relative to the variability of $d_t$; at $N = 12$ a real but moderate gap is usually invisible. The caveat: these are rejection rates of the h = 1 test, whose 5% level holds only if $d_t$ is serially uncorrelated. So “DM failed to reject” never means “the forecasts are equally good”; it means “these dates cannot separate them.”

6c. Closing Module 6’s question

arima_1 (the Module 5 model) against locf 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

DM $=$ 0.13, $p =$ 0.90. Asked in the PS3 one-step design, on 2019, the answer is no: the test does not reject on these twelve months, resting on the h = 1 assumption that these $d_t$ are serially uncorrelated. Be careful what this does not claim. It does not test Module 6’s 24-step path: errors from one origin at horizons $1, \ldots, 24$ share the same early shocks and grow in spread with the horizon, so they are not the sequence DM assumes. And it does not say the two forecasts are equally good; 6b just showed how rarely twelve dates separate anything.

Modify and re-run

Run the same test with power = 1 (absolute loss) on UNRATE. Before you run it: does the sign of DM flip, and does the verdict move? Recall from 5b which track had the smaller MAE.

So what: a scorecard is a ranking, not a verdict; DM is the verdict, and with twelve dates it is usually “can’t tell.”

---

7. Why combine?

On UNRATE the five credible tracks are descriptively close, and the one pair we tested does not reject. Picking “the winner” of a twelve-month race is mostly picking noise. A better question: if you cannot pick a winner, why pick at all?

The crowd and the portfolio. Ask a room to guess a quantity and the average of the guesses usually beats most individual guessers, because their errors partly cancel. Forecasts are assets; forecast errors are their risks; a portfolio of imperfectly correlated assets is less risky than its parts. 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,$$

with $\tilde{\mathbf y}_t$ the $M \times 1$ vector of forecasts, $\mathbf w$ the weights, and $w_0$ an intercept that is zero for every method except one. Weights are always $w$, never $\alpha$ or $\beta$. Population-optimal weights are starred, $w^*$; weights estimated from data carry hats, $\hat w_m$.

7a. Two forecasts: the portfolio algebra

With weights $w$ and $1-w$ the combined error is $e_{c,t} = w\,e_{1,t} + (1-w)\,e_{2,t}$ (the subscript $c$ is a label), so

$$\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}),$$

$$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})}, \qquad \text{zero covariance: } w^* = \frac{\text{Var}(e_{1,t})^{-1}}{\text{Var}(e_{1,t})^{-1} + \text{Var}(e_{2,t})^{-1}}.$$

With zero covariance each forecast’s weight is proportional to the inverse of its error variance: the weaker forecast still gets weight, just less. Write the error correlation as $\rho_{1,2}$, with the comma, because $\rho_{12}$ would be the lag-12 autocorrelation.

Predict first

$\text{Var}(e_{1,t}) = 1$, $\text{Var}(e_{2,t}) = 1.5$. As $\rho_{1,2}$ rises from $-0.5$ to $0.9$, what happens to the best mix? (a) It always beats forecast 1 alone, by less as $\rho_{1,2}$ rises; (b) at high enough correlation no mix between 0 and 1 beats forecast 1 alone; (c) the optimal $w$ never leaves $[0, 1]$.

combo_var <- function(w, v1, v2, rho) {
  cv <- rho * sqrt(v1 * v2)
  w^2 * v1 + (1 - w)^2 * v2 + 2 * w * (1 - w) * cv
}
pf <- expand.grid(w = seq(0, 1, by = 0.01), 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: combined error variance as a function of w",
       subtitle = "Var(e1) = 1, Var(e2) = 1.5; dashed line = forecast 1 alone") +
  theme(legend.position = "bottom")

Correlation eats the gain. The gain over the better forecast vanishes once $\rho_{1,2} \ge \sqrt{\text{Var}(e_{1,t})/\text{Var}(e_{2,t})}$, about 0.82 here. Beyond it the unconstrained $w^* > 1$: the formula wants to short the worse forecast (the $\rho_{1,2} = 0.9$ curve is still falling at $w = 1$). Hold that thought for Granger-Ramanathan.

A simulation check that the formula is right:

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)
round(c(w_star = w_star,
        formula = combo_var(w_star, v1, v2, rho),
        simulated = var(w_star * e1 + (1 - w_star) * e2)), 4)

   w_star   formula simulated 
   0.6416    0.7733    0.7681 

Four-step note: this section’s assumption-break is section 9 and its real data are section 10; the two-forecast algebra is the theory both of them test.

So what: averaging is not a hedge against ignorance, it is arithmetic: imperfectly correlated errors partly cancel, and the formula says by how much.

---

8. Four combination methods

All four use only the weights window: the $T_w = 12$ actuals and the $T_w \times M$ matrix of forecasts of those dates. All four are then applied, unchanged, to the holdout as $\hat w_0 + \hat{\mathbf w}'\tilde{\mathbf y}_t$. Each is a small function with one contract: in, the weights-window actuals y and forecast matrix F (F is only an R argument name; the course never uses it as a math symbol); out, list(intercept, weights). The canonical copies live in R/combination_utils.R; here you build them, as in class.

Method 1: simple average

$w_0 = 0$, $w_m = 1/M$. Nothing is estimated, so no hat: 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

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

The $M$-forecast version of the zero-covariance weight. Each weight comes from its own column only. MSE rather than variance, so bias is penalized 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)

$$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},$$

with an intercept and no constraints; OLS delivers $\hat w_0, \hat w_1, \ldots, \hat w_M$. $v_t$ is the regression’s disturbance: not an innovation, not a forecast error.

Predict first

Twelve observations, an intercept and six forecasts. Will the estimated weights look like (a) shares between 0 and 1, (b) a mix of large positive and negative numbers, or (c) R refusing to estimate one of them? (More than one can be true.)

round(coef(lm(ys_w ~ Fs_w)), 3)

(Intercept) Fs_warima_1 Fs_warima_2 Fs_warima_3    Fs_wlocf   Fs_wavg_k 
    261.700       0.779       1.877      -2.043          NA      -0.833 
  Fs_wmeanf 
    -50.827 

The numbers are wild, and one is NA. That is not a bug. 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. With the regression’s own intercept, the two are perfectly collinear and R drops one. Treat the dropped weight as zero (the same as removing the duplicate track); leaving the NA in makes every holdout forecast 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. Module 4’s callback: parameters that fit the sample are not parameters that forecast.

Method 4: constrained least squares

$$\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{s.t.}\quad w_m \ge 0,\;\; \sum_{m=1}^M w_m = 1, \qquad \hat w_0 = 0.$$

A quadratic program. quadprog::solve.QP() takes the quadratic part (Dmat = crossprod(F)), the linear part (dvec = crossprod(F, y)), and the constraints as columns of Amat against bvec, the first meq holding as equalities. Those are R argument names only; they get no math symbols. Look at the constraint matrix:

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

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 

Several weights land exactly on zero: the inequality constraints are doing model selection. CLS cannot short a forecaster, but it has no intercept, so it can reduce a bias by reweighting yet cannot remove a same-direction bias that every track shares. (Without an intercept, arima_2 and locf are not exactly collinear, so CLS runs here.)

The dial

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

Method	$w_0$	$\sum w_m = 1$	$w_m \ge 0$	Estimated from the window	Estimation risk
Simple average	0	yes	yes	nothing	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

Modify and re-run

Add a duplicate track, Fs_dup <- cbind(Fs_w, locf2 = Fs_w[, "locf"]), and run all four methods on it. Before you do: which methods carry on as if nothing happened, which gives a weight of exactly zero, and which stops? Then try the canonical comb_cls() from R/combination_utils.R on it (section 10 shows how to load the canonical copies without overwriting yours) and read its error message.

So what: a combination method is a choice about how much to estimate, from nothing (the simple average) to everything (Granger-Ramanathan).

---

9. The assumption-break: weights that win the window they were fit on

Theory says Granger-Ramanathan should win: it minimizes the weights-window sum of squared errors over all intercepts and weights, so nothing fits that window better. The assumption buried in every optimal-weight formula is that you know the weights. You don’t; you estimate them from twelve observations. Break it.

The DGP is fair to the estimated weights: equal weights are not optimal here (one track is clearly better), so a method that finds arima_1 should beat the simple average. The question is whether it can find it in twelve months.

Stage 1: one realization, $T_w = 12$

Predict first

Score the four combinations from section 8 on both windows. Which method has the best weights-window RMSE, and where does that same method rank on the holdout? (a) Best on both; (b) best in the window, middling on the holdout; (c) best in the window, worst on the holdout.

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 the best weights-window RMSE, 0.282, against 0.877 for the simple average and 0.706 for arima_1 alone. On the holdout it posts 3.009, the worst number in the table. An intercept of 262 and weights from -50.8 to 1.9 are not a discovery about the forecasters; they are a regression with almost as many parameters as data points, fitting noise.

Stage 2: is it a fluke? 200 replications

Each replication simulates a fresh series, estimates the three ARIMAs once, builds the tracks, fits the four combinations on the weights window, and scores the 12-date holdout. Record each method’s holdout RMSE relative to the simple average (below 1: beat it). The same loop sets up Stage 3 by also running weights windows $T_w \in \{24, 60, 120\}$. (About half a minute.)

set.seed(6378)   # Stages 2-3 have their own seed
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))
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")

Granger-Ramanathan’s long right tail is the puzzle in one picture. Median 1.68 times the simple average’s holdout RMSE, 90th percentile 4.07 times, and it beats the simple average in only 16.0% of replications. Inverse-MSE is the simple average with a slight tilt (median 0.99). Constrained LS has a thin tail and, in this DGP, already beats the simple average at the median (0.95): its constraints shrink the estimation error while still letting it find the dominant track. That is evidence for the estimation-error explanation, not against it.

Stage 3: how long a window do estimated weights need?

Predict first

Lengthen the weights window to 24, 60, 120. Does Granger-Ramanathan’s median ratio (a) stay above 1 throughout, (b) cross below 1 somewhere between 24 and 120, or (c) cross below 1 already at 24?

stage3_med <- apply(rel, c(2, 3), median)
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: estimated weights earn their keep only with a longer window") +
  theme(legend.position = "bottom")

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 never quite catches arima_1 alone (0.90) or CLS (0.90). The crossover is real because the optimal weights here are far from $1/M$; with enough data, estimating them pays. (In this sweep the holdout moves with $T_w$; the DGP is stationary, so that is harmless for the ARIMA tracks, though a later holdout gives meanf a slightly longer history.)

The forecast combination puzzle

Stock and Watson (2004) named it: theory says optimally weighted combinations should dominate; in practice the simple average is very hard to beat. The simulation shows why. A combination method is an estimator: its holdout performance is the population gain from its weights minus the estimation error in them. The puzzle is a comparison: attainable gain over equal weights against the cost of estimating the weights. Near-equal optimal weights are the important special case (then there is almost nothing to gain), but not a necessary condition: optimal weights far from $1/M$ can still buy little when the forecasts are largely redundant. With twelve weights-window observations the estimation cost is high, which describes most applied work.

The index fund. Active managers (estimated weights) should beat the index (equal weights) in theory; in practice the fee eats the edge, and here the fee is estimation error. 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.

So what: with twelve observations, estimate as little as you can, unless the gain on the table is large enough to pay the fee.

---

10. The PS3 workflow on a series that isn’t yours

The complete combination step on UNRATE: $N_{\text{total}} =$ 720, $T =$ 696 (train ends December 2017), weights window 2018, holdout 2019, the six tracks from section 2. The metrics table and the DM test of arima_1 against locf are done (5b, 6c). What remains is the combination step.

Predict first

Section 9 says the simple average is hard to beat. The UNRATE pool contains meanf, sitting near 6.0% while the actuals run below 4%. On the 2019 holdout, which wins? (a) The simple average, as the puzzle says; (b) inverse-MSE and CLS; (c) Granger-Ramanathan, because it has an intercept to absorb the bias.

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

Inverse-MSE gives meanf a weight of 0.001, effectively zero. 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. 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 design matrix is about 41,620; anything in the thousands means the solution is extremely sensitive to the data.

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’s holdout RMSE is 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. Picking the holdout winner and then testing it on the same twelve dates would be 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: DM $=$ -9.1, $p =$ 1.9e-06; CLS: DM $=$ -9.1, $p =$ 1.8e-06. A clear rejection at $N = 12$: section 6’s 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 here, these $p$-values assume each $d_t$ is serially uncorrelated (exploration 2 checks it).

The honest reading

• The simple average did not win, and not because the estimated weights got lucky. 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 points off the actuals every month. Here the attainable gain over equal weights is large: inverse-MSE and CLS won by down-weighting meanf, which is exactly what estimated weights are for.
• 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.
• Caveats. One window, and a calm one with no turning point. A different window could rank the methods differently. One realization is one draw; the Monte Carlo of section 9 is the distribution.

Your inline builds against the canonical copies

The course keeps one canonical copy of each function in R/. Load them into a separate environment, so your inline versions are not overwritten, and confirm they agree on UNRATE:

canon <- new.env()
sys.source("R/track_utils.R", envir = canon)
sys.source("R/combination_utils.R", envir = canon)

same_tracks <- isTRUE(all.equal(canon$tracks_on(unrate, fits_un, h_dates), F_h))
same_weights <- sapply(names(ws_un), function(m) {
  f <- get(c(SA = "comb_simple", IMSE = "comb_inv_mse",
             CLS = "comb_cls", GR = "comb_gr")[[m]], envir = canon)
  isTRUE(all.equal(f(y_w, F_w), ws_un[[m]]))
})
c(tracks = same_tracks, same_weights)

tracks     SA   IMSE    CLS     GR 
  TRUE   TRUE   TRUE   TRUE   TRUE 

The canonical copies differ from the inline builds only in their guard rails: comb_gr() warns when it zeroes an aliased track, comb_cls() stops and names exactly collinear tracks, and kavg() refuses dates without enough past.

So what: equal weights are the benchmark to beat when the pool is sound; look at your tracks before you combine them, because a broken member makes the attainable gain large enough to pay for estimating the weights.

---

11. Explorations

Optional, and outside the assessable surface

None of these is required, and none of them is on a problem set or exam. They run on the objects already in memory, and they are ordered from guided to open. Take one that interests you; the best of them is the one you invent.

1. Can you sideline a broken track without estimating anything? (guided; Deeper Dive)

The cross-forecaster median, and a trimmed mean that drops the highest and lowest forecast each month, estimate nothing at all:

round(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]))), 3)

     simple      median     trimmed simple_five 
      0.416       0.099       0.098       0.093 

Why does trimming throw out meanf every month? Test the median against the simple average with dm.test(). Then run the same three on the simulation’s holdout (Fs_h, ys_h), where no track is broken. Do they still help?

2. Does the h = 1 assumption hold on these $d_t$? (guided; Technical Note)

Every UNRATE DM test above assumed serially uncorrelated loss differentials. Look at the autocorrelations of each $d_t$ at lags 1–3, then re-run each test with a Bartlett (Newey-West) variance that allows for serial correlation:

dmdep_pairs <- list(
  "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"]))
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)
}))
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))
}))
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

Does any verdict at 5% move? With twelve points, what size of autocorrelation could the Ljung-Box test even detect? (Use "bartlett", not the default "acf": the truncated sum can turn negative on twelve points, and dm.test() then warns and falls back to h = 1.)

3. Is arima_2 against arima_1 a fair fight for DM? (guided; Deeper Dive)

arima_2 (an AR(1)) is nested in arima_1 (the seasonal AR), and so is locf (the random walk) inside the UNRATE arima_1. Standard DM can be undersized for nested pairs. The Clark-West adjustment tests whether the mean of $z_t = e_{1,t}^2 - \big[e_{2,t}^2 - (\tilde y^{(1)}_{t \mid t-1} - \tilde y^{(2)}_{t \mid t-1})^2\big]$ is positive, one-sided, with forecast 1 from the small model. Add it to the power loop in 6b as a third test. Does it detect arima_1’s advantage over arima_2 more often?

4. Should the weights move? (open; Looking Ahead)

Everything here estimated the weights once. Re-estimate them on a moving 12-month window as each holdout month arrives, on the simulation with $T_w = 120$ of history. Does adaptivity buy anything on a stationary DGP? Now build a DGP with a break halfway through (the better track swaps), and ask again.

5. Can you build a pool that fools the “near-equal weights” story? (open)

The puzzle bites when the attainable gain is small. One way is near-equal optimal weights; another is redundancy. Build a pool of three tracks where two are near-copies of arima_1 and one is locf. Are the optimal weights far from $1/3$? Is the gain over the simple average large or small? Explain with the two-forecast algebra of section 7 and $\rho_{1,2}$.

6. What would you add to the pool? (open, off-syllabus)

Every forecast today used only the series’ own past. Invent a seventh forecaster for the simulated series that uses something else: the true DGP mean, a forecast from a deliberately wrong seasonal period, a forecast that is arima_1 plus noise. Predict what each combination method does with it, then run Stage 1 and see. When does a new forecaster help a combination even though it is worse on its own?

---

12. Checks

R/checks.R is the referee: seven assertions about things this module claims are true. It does not know what you concluded and it cannot grade you. It just refuses to let an incorrect claim about the code stand, whether that claim came from you, your AI tutor, or a typo. If a check fails, CHECKS.md says what it asserts and what a failure usually means.

source("R/checks.R", local = new.env())   # its own environment: your objects stay put

Module 7 checks
---------------
PASS UNRATE cache intact: one blank (2025-10), 720 months in 1960-2019, T = 696 ends December 2017  [rows = 940, blanks = 1, window = 720, T = 696, Dec 2017 = 4.1, Dec 2019 = 3.6]
PASS upto() keeps frequency 12 (brackets give 1); locf = y_{t-1}, avg_k = 3-month mean, meanf = expanding mean; extending the series leaves earlier one-step forecasts unchanged  [refilter gap = 0.0e+00; meanf on the first holdout date = 6.005]
PASS theil_u2() and mase() match their formulas (m = 1 and m = 12), ignore ts dates, and stop on a zero denominator  [U2 = 0.6675, MASE = 0.4990, seasonal MASE = 0.1003]
PASS dm.test(h = 1, power = 2) equals the one-sample t-test on d_t = e1^2 - e2^2, p-value from t with N - 1 = 11 df  [DM = -1.1301, t = -1.1301, p = 0.2825]
PASS Two-forecast w* minimizes Var(e_c) (5 cases); zero covariance gives inverse-variance weights; equal variances halve; w* = 1 at rho_{1,2} = sqrt(V1/V2)  [max |numeric - formula| = 3.3e-16; threshold rho_{1,2} = 0.816]
PASS comb_*: SA = 1/M; IMSE positive, sums to 1, ranks by inverse MSE; CLS on the simplex and beats SA/IMSE in-window; GR = lm(), best in-window; duplicates: GR weight 0 + warning, CLS stops  [in-window SSE: SA 1.023, IMSE 0.464, CLS 0.260, GR 0.183]
PASS UNRATE 2019: auto.arima() picks ARIMA(1,1,3)(2,0,2)[12]; DM arima_1 vs locf p = 0.902; SA holdout RMSE = 0.416; IMSE meanf weight < 0.01; GR intercept = -141; IMSE, CLS reject vs SA, GR does not  [pick = ARIMA(1,1,3)(2,0,2)[12]; p = 0.902; SA RMSE = 0.416; IMSE meanf = 0.0006; GR w0 = -141.45; p vs SA: IMSE 1.9e-06, CLS 1.8e-06, GR 0.403]
---------------
7/7 checks passed.

---

Where you are now

Nothing new switched on today; you refereed between forecasts the board had already produced. You found that $U_2$ and MASE answer different questions and that on a calm stretch everything clears the MASE yardstick; that DM is a $t$-test on loss differentials, assumes at h = 1 that they are serially uncorrelated, and with twelve dates usually says “can’t tell”; that combining works because errors diversify; and that a combination method is an estimator, so with twelve observations the simple average is hard to beat, unless the pool has a broken member, in which case down-weighting it pays.

The question you carry into Module 8: every forecast here used only the series’ own past, and combinations can only reshuffle the information already in the pool. What if another variable could improve the forecast? Module 8 turns on a new slider, $\beta x_t$ and its lags.

Before you close this file: open LEARNING_LOG.md and write down, in your own words, one prediction you got wrong today and what you now think instead. Your words, not your tutor’s — that is the entire value of the exercise.