Module 6 Lab: Forecasting Fundamentals

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 say, before the future arrives, how wrong you are about to be.

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, and in this module the habit becomes the subject: a forecast is a prediction written down before the future runs.

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.

Where you are on the mixing board: every slider is already on. Module 3 lit the AR and MA blocks, Module 4 chose their orders, Module 5 added the seasonal band. Nothing new comes on today; you use the board, and the question changes from “what is this process?” to “what does it say about next quarter, and how wrong will that be?”

Three questions. What is the best guess? How uncertain is it? And how do you keep score when you are wrong every single time?

1. Setup

What do you need on the table? Three packages, the course helpers, and the series from class: real GDP per person, which you fit on 1990 to 2006 and later let meet the five years that followed. Nothing after 2006 Q4 gets looked at until section 7.

library(ggplot2)
library(forecast)
library(urca)

source("R/forecast_utils.R")     # ts_split() and the course loss functions
source("R/model_selection.R")    # hqic(), the fourth criterion from class
source("R/arma_simulator.R")     # the Module 3 simulator, for the explorations

theme_set(theme_bw(base_size = 13))

gdp_raw <- read.csv("data/A939RX0Q048SBEA.csv", stringsAsFactors = FALSE)
gdp_all <- ts(gdp_raw$A939RX0Q048SBEA, start = c(1947, 1), frequency = 4)

# The worked example works in natural logs, on 1990 Q1 to 2006 Q4 only.
# Everything after 2006 Q4 stays unread until section 7.
y  <- log(window(gdp_all, start = c(1990, 1), end = c(2006, 4)))
dy <- diff(y)                              # quarterly growth, in log points
c(T_levels = length(y), T_growth = length(dy))
T_levels T_growth 
      68       67 
# Two small plotting helpers used all the way down: a sample ACF (or PACF) as a
# data frame, and a correlogram with the +/- 1.96/sqrt(T) reference band.
acf_df <- function(x, lag_max = 24, label = "", partial = FALSE) {
  r <- if (partial) stats::pacf(as.numeric(x), lag.max = lag_max, plot = FALSE)$acf
       else stats::acf(as.numeric(x), lag.max = lag_max, plot = FALSE)$acf[-1]
  data.frame(lag = seq_len(lag_max), rho = as.numeric(r), fit = label)
}
acf_plot <- function(d, T_n, title = NULL, subtitle = NULL, ncol = 2) {
  band <- 1.96 / sqrt(T_n)
  d$fit <- factor(d$fit, levels = unique(d$fit))
  g <- ggplot(d, aes(lag, rho)) +
    geom_hline(yintercept = c(-band, band), linetype = "dashed", colour = "grey55") +
    geom_hline(yintercept = 0, colour = "grey40", linewidth = 0.3) +
    geom_segment(aes(xend = lag, y = 0, yend = rho), colour = "steelblue4",
                 linewidth = 1) +
    labs(title = title, subtitle = subtitle, x = "lag k", y = expression(hat(rho)[k]))
  if (length(levels(d$fit)) > 1) g <- g + facet_wrap(~ fit, ncol = ncol)
  g
}

So what: sixty-eight quarters of one series, ending the year before the Great Recession, and a model you are about to build from them is going to be asked what comes next.

2. The example from class, in three brisk moves

Why GDP per person, and why logs? Because constant percentage growth is a straight line in logs, and the first difference of the log is the quarterly growth rate, the number everyone argues about. Modules 2 through 5 built the whole chain: is it stationary, what does the ACF say, which orders, do the residuals pass. You have done each step before, so this section moves fast: one picture, one table, one prediction per move.

2a. Trend-stationary, or a unit root?

Look at the series before any test touches it. Does the log level hug a straight line, or wander away from one?

level_df <- data.frame(date = as.numeric(time(y)), y = as.numeric(y), t = seq_along(y))
trend_fit <- lm(y ~ t, data = level_df)
level_df$trend <- fitted(trend_fit)

ggplot(level_df, aes(date)) +
  geom_line(aes(y = trend), colour = "grey45", linetype = "dashed") +
  geom_line(aes(y = y), colour = "steelblue4", linewidth = 0.8) +
  labs(title = "Log real GDP per capita, 1990 Q1 to 2006 Q4, with a fitted linear trend",
       x = NULL, y = "log of chained 2017 dollars per person")

round(c(trend_growth_pct_per_year = 400 * unname(coef(trend_fit)["t"])), 2)
trend_growth_pct_per_year 
                     2.17 
acf_plot(acf_df(y, 16), length(y), title = "ACF of the log level: persistence you can see")

The trend grows 2.17% a year and the series rarely strays far from it. The ACF decays slowly, which Module 1 taught you proves nothing on its own. So: test.

Predict first

Sixty-eight quarters, a trend specification, two lags. Name the ADF statistic you expect, to one decimal. Will it clear the 5% critical value of −3.45-3.45? And will KPSS, whose null is the opposite, reject?

adf_level  <- ur.df(y, type = "trend", lags = 2)              # k = 2, chosen by BIC
kpss_short <- ur.kpss(y, type = "tau", lags = "short")
kpss_long  <- ur.kpss(y, type = "tau", lags = "long")
dfgls      <- ur.ers(y, type = "DF-GLS", model = "trend", lag.max = 2)
u_hat      <- residuals(trend_fit)                            # deviations from the trend

round(rbind(
  "ADF, trend, k = 2 (H0: unit root)"      = c(statistic = adf_level@teststat[1],  cv_5pct = adf_level@cval[1, "5pct"]),
  "DF-GLS, trend, k = 2 (H0: unit root)"   = c(statistic = dfgls@teststat[1],      cv_5pct = dfgls@cval[1, "5pct"]),
  "KPSS, trend, short lag (H0: stationary)" = c(statistic = kpss_short@teststat[1], cv_5pct = kpss_short@cval[1, "5pct"]),
  "KPSS, trend, long lag (H0: stationary)"  = c(statistic = kpss_long@teststat[1],  cv_5pct = kpss_long@cval[1, "5pct"])
), 3)
                                        statistic cv_5pct
ADF, trend, k = 2 (H0: unit root)          -3.113  -3.450
DF-GLS, trend, k = 2 (H0: unit root)       -1.965  -3.030
KPSS, trend, short lag (H0: stationary)     0.150   0.146
KPSS, trend, long lag (H0: stationary)      0.085   0.146
round(c(rule_of_thumb_raw = sd(dy) / sd(y), rule_of_thumb_on_deviations = sd(dy) / sd(u_hat)), 3)
          rule_of_thumb_raw rule_of_thumb_on_deviations 
                      0.050                       0.359

Neither null is clearly rejected, and you should sit with that for a second. ADF gives τ=−3.11\tau = -3.11 against −3.45-3.45: not below, so a unit root is not rejected. KPSS sits on its own critical value: over the line in the third decimal with one lag truncation, nowhere near it with the other. The rule of thumb on the raw series is useless here (the trend inflates σy\sigma_y); on the deviations from the trend it is 0.36, still below 0.5, still flagged as persistent. Can 68 quarters tell a trend from a unit root? No. Module 2’s power lesson applies: near the boundary these tests rarely reject, so “never rejects” is weak evidence, not a finding.

A choice, not a finding. Under doubt the course’s rule is to difference, because leaving a unit root in the levels is the more dangerous error and over-differencing leaves a fingerprint you can check: an MA coefficient pinned near −1-1 once you fit. So treat yty_t as I(1)I(1), model the growth rate, and look for that fingerprint.

Modify and re-run

Replace lags = 2 with lags = 8, selectlags = "BIC" and run the ADF line again. Before you do: BIC will still pick two lags, so will τ\tau come out the same, closer to −3.45-3.45, or further away? Why would the same two lags give a different number? (Hint: how many observations does each regression use, and what happened in the first quarters of the window?)

So what: on a sample this short the stationarity verdict is a decision you own, and a forecaster cares because a trend-stationary model’s forecast returns to the line while a unit-root model’s never does.

2b. Was one difference enough?

Here is the series you will actually model: the quarterly growth rate. Before the tests, look.

growth_df <- data.frame(date = as.numeric(time(dy)), g = 100 * as.numeric(dy))
ggplot(growth_df, aes(date, g)) +
  geom_hline(yintercept = 100 * mean(dy), linetype = "dashed", colour = "grey45") +
  geom_line(colour = "steelblue4", linewidth = 0.7) +
  labs(title = "Quarterly growth of real GDP per capita, 1990 Q2 to 2006 Q4",
       x = NULL, y = "percent per quarter")

round(c(mean_pct_per_year = 400 * mean(dy), sd_pct_per_quarter = 100 * sd(dy)), 2)
 mean_pct_per_year sd_pct_per_quarter 
              1.87               0.54

Growth averages 1.87% a year with a standard deviation of 0.54% a quarter, and it crosses its mean constantly. That is what a stationary series looks like. Now the same battery, one level down: drift specification (the trend is gone), KPSS around a constant.

adf_growth  <- ur.df(dy, type = "drift", lags = 1)            # k = 1, chosen by BIC
kpss_growth <- ur.kpss(dy, type = "mu", lags = "short")
rho1 <- acf(dy, plot = FALSE)$acf[2]

round(rbind(
  "ADF, drift, k = 1 (H0: unit root)"    = c(statistic = adf_growth@teststat[1],  cv_5pct = adf_growth@cval[1, "5pct"]),
  "KPSS, level, short (H0: stationary)"  = c(statistic = kpss_growth@teststat[1], cv_5pct = kpss_growth@cval[1, "5pct"])
), 3)
                                    statistic cv_5pct
ADF, drift, k = 1 (H0: unit root)      -3.833  -2.890
KPSS, level, short (H0: stationary)     0.214   0.463
round(c(rule_of_thumb_one_level_down = sd(diff(dy)) / sd(dy), rho_1 = rho1, band = 1.96 / sqrt(length(dy))), 3)
rule_of_thumb_one_level_down                        rho_1 
                       1.186                        0.302 
                        band 
                       0.239

Reject a unit root (τ=−3.83\tau = -3.83 against −2.89-2.89), fail to reject stationarity: the “strong evidence of stationarity” row of Module 2’s table. The rule of thumb one level down is above 1, so a second difference would add variance. And the first look at over-differencing: Module 3 told you that differencing a series that did not need it leaves ρ̂1≈−0.5\hat\rho_1 \approx -0.5; here ρ̂1=+0.3\hat\rho_1 = +0.3. Keep Δyt\Delta y_t.

Predict first

You are about to see the ACF and PACF of growth. With T=67T = 67 the band is ±0.24\pm 0.24. Sketch both: how many lags clear the band in each, and does either one cut off?

both <- rbind(acf_df(dy, 16, "ACF of growth"), acf_df(dy, 16, "PACF of growth", partial = TRUE))
acf_plot(both, length(dy), title = "The series we model")

Both tail off. Two bars stand up in each panel and then everything fades: in the ACF both clear the band, in the PACF the first does and the second (0.22) falls just short of it. The ACF decays no faster than the PACF. Module 3’s table: ACF tails off and PACF cuts off means AR; both tail off means ARMA. Nothing here cuts off cleanly, so AR or ARMA, not MA. Module 5’s rule: start one size bigger and let the criteria and the residuals trim it. Starting model: ARIMA(2,1,2) with drift. The drift is the constant in the growth equation; without it the model would forecast zero growth forever.

So what: one difference turned a borderline series into a textbook stationary one, and the pictures name a class, not a pair of orders.

2c. Which model?

Which (p,q)(p, q)? Fit the three you would write down from the pictures plus the two auto.arima() found, and let four criteria score them. One habit from Module 4: pick your criterion before you look.

Predict first

Five candidates, T=67T = 67. Which criterion do you expect to pick the smallest model, and will all four agree? Name the model you think wins AICc.

fit_drift <- function(p, q) Arima(y, order = c(p, 1, q), include.drift = TRUE)
cands <- list("ARIMA(2,1,2)" = fit_drift(2, 2), "ARIMA(2,1,1)" = fit_drift(2, 1),
              "ARIMA(1,1,1)" = fit_drift(1, 1), "ARIMA(2,1,0)" = fit_drift(2, 0),
              "ARIMA(1,1,0)" = fit_drift(1, 0))

ic_table <- do.call(rbind, lapply(names(cands), function(nm) {
  f  <- cands[[nm]]
  e  <- residuals(f)[-1]           # drop 1990 Q1: a d = 1 fit has no growth to explain there
  pq <- sum(f$arma[1:2])           # fitdf = p + q (Module 5)
  data.frame(model = nm, k = length(coef(f)) + 1, logLik = round(f$loglik, 2),
             AIC = round(f$aic, 2), AICc = round(f$aicc, 2), BIC = round(f$bic, 2),
             HQIC = round(hqic(f), 2),
             LB10_p = round(Box.test(e, lag = 10, type = "Ljung-Box", fitdf = pq)$p.value, 2),
             LB24_p = round(Box.test(e, lag = 24, type = "Ljung-Box", fitdf = pq)$p.value, 2))
}))
ic_table
         model k logLik     AIC    AICc     BIC    HQIC LB10_p LB24_p
1 ARIMA(2,1,2) 6 261.62 -511.25 -509.85 -498.02 -506.01   0.22   0.55
2 ARIMA(2,1,1) 5 261.28 -512.56 -511.58 -501.54 -508.20   0.64   0.60
3 ARIMA(1,1,1) 4 260.28 -512.56 -511.92 -503.74 -509.07   0.56   0.66
4 ARIMA(2,1,0) 4 260.82 -513.63 -512.99 -504.81 -510.14   0.66   0.61
5 ARIMA(1,1,0) 3 259.00 -512.00 -511.61 -505.38 -509.38   0.40   0.66
sapply(ic_table[, c("AIC", "AICc", "BIC", "HQIC")], function(col) ic_table$model[which.min(col)])
           AIC           AICc            BIC           HQIC 
"ARIMA(2,1,0)" "ARIMA(2,1,0)" "ARIMA(1,1,0)" "ARIMA(2,1,0)"

AIC, AICc and HQIC pick ARIMA(2,1,0); BIC, the harshest tariff at this TT, goes one size smaller to ARIMA(1,1,0). Every candidate passes Ljung-Box at both horizons, so the residuals do not choose between them; parsimony does, and the criteria disagree only about the second AR lag. The model you started from is last under every criterion. Look at why:

round(coef(cands[["ARIMA(2,1,2)"]]), 3)
   ar1    ar2    ma1    ma2  drift 
 1.830 -0.861 -1.652  0.652  0.005 
round(Mod(polyroot(c(1, coef(cands[["ARIMA(2,1,2)"]])[c("ma1", "ma2")]))), 3)   # MA roots
[1] 1.000 1.533

An MA root at modulus 1: the over-differencing warning light from Module 3, the model trying to undo the difference. That is the same ambiguity as section 2a in another costume. It is not a verdict for trend-stationarity, but it is a reason not to carry ARIMA(2,1,2) forward. auto.arima() screens such fits out before scoring, which is why it could never have been its answer:

auto.arima(y, d = 1, seasonal = FALSE, stepwise = FALSE, approximation = FALSE)
Series: y 
ARIMA(2,1,0) with drift 

Coefficients:
         ar1     ar2   drift
      0.2317  0.2330  0.0046
s.e.  0.1178  0.1202  0.0011

sigma^2 = 2.716e-05:  log likelihood = 260.82
AIC=-513.63   AICc=-512.99   BIC=-504.81

A very good intern: it agrees with AICc, and you should still read its work. Keep both finalists. The IC table scored them on 1990 to 2006; whether the extra lag earns its keep on 2007 to 2011 is a different question, and the test set in section 7 answers it.

fit_210 <- cands[["ARIMA(2,1,0)"]]
fit_110 <- cands[["ARIMA(1,1,0)"]]
round(rbind("ARIMA(2,1,0)" = c(coef(fit_210), sigma_pct = 100 * sqrt(fit_210$sigma2)),
            "ARIMA(1,1,0)" = c(coef(fit_110)[1], ar2 = NA, coef(fit_110)[2], sigma_pct = 100 * sqrt(fit_110$sigma2))), 4)
                ar1   ar2  drift sigma_pct
ARIMA(2,1,0) 0.2317 0.233 0.0046    0.5211
ARIMA(1,1,0) 0.3008    NA 0.0047    0.5308

So what: the whole Modules 2 to 5 chain ran once on real data and ended where it usually ends, with two defensible models and a scoreboard that cannot separate them. From here on, the question is what they say about the future.

3. What a forecast is

What is the one thing you can be certain of about any forecast you ever make? That it will be wrong. The unemployment rate next month will not be exactly 4.1%. Growth next quarter will not be exactly 0.47%. From Fiedler’s forecasting rules (Across the Board, June 1977): “The moment you forecast you know you’re going to be wrong, you just don’t know when and in which direction.”

So the question is not whether, but how wrong, and whether your wrongness is useful. Which of these two is more likely to be right, and which is more useful: “growth next quarter will be between −3-3% and +3+3%”, or “0.4%, give or take half a point”? The first is almost certainly right and tells nobody anything. The second is almost certainly wrong in the second decimal and is what a decision gets made on.

Three words, then two objects. The origin tt is where you stand. The horizon hh is how far ahead you look. The information set ℱt={yt,yt−1,…}\mathcal{F}_t = \{y_t, y_{t-1}, \ldots\} is everything you know at tt.

  • The point forecast ŷt+h|t=𝔼[yt+h∣ℱt]\hat y_{t+h|t} = \mathbb{E}[y_{t+h} \mid \mathcal{F}_t]: the best guess under squared-error loss, because the conditional mean minimizes expected squared error. Hold on to that qualifier; section 9 is about what happens when you change it.
  • The forecast error et+h|t=yt+h−ŷt+h|te_{t+h|t} = y_{t+h} - \hat y_{t+h|t}: unknowable at tt, which is the point, and observable after the fact, which is why it is an ee and not an ϵ\epsilon.
  • The prediction interval ŷt+h|t±ch\hat y_{t+h|t} \pm c_h with a stated coverage, 80% or 95%: the fine print, where the honesty lives.

Behind both sits a third object, the density forecast: the whole conditional distribution of yt+hy_{t+h}, from which any quantile can be read; the fan chart is one picture of it. The crystal ball is the conditional mean, the error bars are the prediction interval, and a fortune teller who gives you both is an econometrician.

So what: a forecast is a distribution, and the headline number is just its mean under one particular way of scoring mistakes.

4. The recursion

How does a model say anything about t+1t + 1? Take the simplest one, an AR(1):

yt+1=α+ϕyt+ϵt+1. y_{t+1} = \alpha + \phi\, y_t + \epsilon_{t+1}.

At time tt you know yty_t. The only unknown is ϵt+1\epsilon_{t+1}, and its conditional expectation is zero. So ŷt+1|t=α+ϕyt\hat y_{t+1|t} = \alpha + \phi y_t. Two steps ahead, yt+1y_{t+1} is unknown too, so replace it with its forecast: ŷt+2|t=α+ϕŷt+1|t\hat y_{t+2|t} = \alpha + \phi\, \hat y_{t+1|t}. Two rules: replace future innovations with zero; replace future yy’s with their forecasts. Each step feeds the last one forward, which is why it is called recursive.

Build one. This is Module 1’s loop, ϕ=0.7\phi = 0.7, sixty observations you can see and twenty more that stay sealed until section 6.

set.seed(2001)
phi_true <- 0.7; n_obs <- 60; n_future <- 20
eps <- rnorm(100 + n_obs + n_future)
z   <- numeric(length(eps))
for (t in 2:length(z)) z[t] <- phi_true * z[t - 1] + eps[t]      # y_t = 0.7 y_{t-1} + eps_t
y_ar1  <- ts(z[101:(100 + n_obs)])                               # t = 1..60: what you can see
future <- z[(100 + n_obs + 1):length(z)]                         # t = 61..80: sealed for now

fit_ar1 <- Arima(y_ar1, order = c(1, 0, 0))
round(coef(fit_ar1), 3)
      ar1 intercept 
    0.723     0.225

Now the trap, disarmed exactly once. R calls that second coefficient intercept. It is not α̂\hat\alpha. It is the estimated mean μ̂\hat\mu, and α=μ(1−ϕ)\alpha = \mu(1 - \phi). Plug μ̂\hat\mu in where α̂\hat\alpha belongs and your hand forecast disagrees with forecast(); it bites everyone exactly once, usually in a problem set.

Predict first

yT=0.764y_T = 0.764 and μ̂=0.225\hat\mu = 0.225. Write down ŷT+1|T\hat y_{T+1|T} and ŷT+2|T\hat y_{T+2|T} by hand before the chunk does. Then the real question: as h→∞h \to \infty, where does the forecast go? Name the number.

mu_hat    <- unname(coef(fit_ar1)["intercept"])      # the mean, despite the name
phi_hat   <- unname(coef(fit_ar1)["ar1"])
alpha_hat <- mu_hat * (1 - phi_hat)                  # recover alpha-hat
y_T       <- as.numeric(y_ar1)[n_obs]

f1 <- alpha_hat + phi_hat * y_T
f2 <- alpha_hat + phi_hat * f1
f3 <- alpha_hat + phi_hat * f2
trap <- mu_hat + phi_hat * y_T                       # "intercept" used as alpha

fc_ar1 <- forecast(fit_ar1, h = n_future, level = c(80, 95))
data.frame(h = 1:3, by_hand = round(c(f1, f2, f3), 4),
           forecast_pkg = round(as.numeric(fc_ar1$mean)[1:3], 4),
           intercept_as_alpha = round(c(trap, NA, NA), 4))
  h by_hand forecast_pkg intercept_as_alpha
1 1  0.6153       0.6153             0.7784
2 2  0.5075       0.5075                 NA
3 3  0.4295       0.4295                 NA
closed_form <- mu_hat + phi_hat^(1:n_future) * (y_T - mu_hat)   # the recursion, solved
max(abs(closed_form - as.numeric(fc_ar1$mean)))
[1] 5.551115e-17

Same numbers, to machine precision, and the wrong one is off by 0.163 at h=1h = 1. The closed form says where it goes: subtract μ\mu from both sides and the deviation shrinks by ϕ\phi each step,

ŷt+h|t−μ=ϕh(yt−μ), \hat y_{t+h|t} - \mu = \phi^h\,(y_t - \mu),

so the forecast forgets its origin at the rate ϕh\phi^h and lands on μ̂\hat\mu, not on zero. Stationarity is mean reversion; ϕ\phi is how slowly it forgets.

path_df <- data.frame(t = 1:(n_obs + n_future),
                      value = c(as.numeric(y_ar1), as.numeric(fc_ar1$mean)),
                      what = rep(c("observed", "forecast"), c(n_obs, n_future)))
ggplot(path_df, aes(t, value, colour = what)) +
  geom_hline(yintercept = mu_hat, linetype = "dotted", colour = "grey40") +
  geom_vline(xintercept = n_obs + 0.5, colour = "grey70") +
  geom_line() + geom_point(size = 1.2) +
  scale_colour_manual(values = c(observed = "steelblue4", forecast = "firebrick"), name = NULL) +
  labs(title = "An AR(1) forecast is a geometric slide to the mean", x = "t", y = NULL)

The same two rules, bigger models

Does the recipe survive a second lag and an MA term? For an AR(2) the second lag is still observed at h=2h = 2 and becomes a forecast at h=3h = 3. For an ARMA(1,1) the last residual eTe_T enters at h=1h = 1 (it is past, so you use it) and then its explicit term is gone: ϵT+1\epsilon_{T+1} is future, so zero. Verify both against forecast(), to 10−1010^{-10}.

set.seed(1985)
y_ar2   <- ts(arima.sim(model = list(ar = c(0.6, -0.2)), n = 200))
fit_ar2 <- Arima(y_ar2, order = c(2, 0, 0))
ph  <- coef(fit_ar2)[c("ar1", "ar2")]
a2  <- unname(coef(fit_ar2)["intercept"]) * (1 - sum(ph))      # alpha-hat again
yT  <- y_ar2[200]; yT1 <- y_ar2[199]
g1  <- a2 + ph[1] * yT + ph[2] * yT1          # both lags observed
g2  <- a2 + ph[1] * g1 + ph[2] * yT           # first lag is a forecast, second still observed
g3  <- a2 + ph[1] * g2 + ph[2] * g1           # both lags are forecasts
rbind(by_hand = unname(c(g1, g2, g3)), forecast_pkg = as.numeric(forecast(fit_ar2, h = 3)$mean))
                  [,1]       [,2]        [,3]
by_hand      0.4796966 0.03357862 -0.05316904
forecast_pkg 0.4796966 0.03357862 -0.05316904
set.seed(2015)
y_arma   <- ts(arima.sim(model = list(ar = 0.6, ma = 0.4), n = 200))
fit_arma <- Arima(y_arma, order = c(1, 0, 1))
phA <- unname(coef(fit_arma)["ar1"]); thA <- unname(coef(fit_arma)["ma1"])
aA  <- unname(coef(fit_arma)["intercept"]) * (1 - phA)
e_T <- residuals(fit_arma)[200]                              # the last in-sample residual
m1  <- aA + phA * y_arma[200] + thA * e_T     # h = 1: the residual enters, once
m2  <- aA + phA * m1                          # h = 2: epsilon_{T+1} is future, so zero
m3  <- aA + phA * m2
rbind(by_hand = unname(c(m1, m2, m3)), forecast_pkg = as.numeric(forecast(fit_arma, h = 3)$mean))
                 [,1]      [,2]       [,3]
by_hand      -1.86925 -1.354876 -0.9987811
forecast_pkg -1.86925 -1.354876 -0.9987811

The MA term is not gone after h=1h = 1, though. Its contribution θ̂eT\hat\theta e_T is baked into ŷT+1|T\hat y_{T+1|T} and rides forward through the AR recursion, shrinking by ϕ̂\hat\phi each step. A pure MA(qq) forecast goes flat after qq steps; an ARMA fades.

Modify and re-run

In the ar1-draw chunk, change phi_true <- 0.7 to 0.95, then re-run down to the path plot. Before you do: will the forecast reach μ̂\hat\mu sooner or later, and roughly at what hh will it be within 0.1 of it? Why does Module 1’s “ϕ\phi is memory” answer both questions at once?

So what: forecast() is doing nothing you could not do on the board, and every point forecast from a stationary model is a slide toward the mean whose speed is set by ϕ\phi.

5. How wrong: the interval

Why does your weather app give a tighter range for tomorrow than for next Thursday? Because between now and Thursday more unknown things will happen, and each day adds another layer of unresolved randomness. A time-series forecast is the same: every future period adds an ϵ\epsilon you have not seen.

The hh-step error is the sum of those unseen shocks, each weighted by how long the model remembers it:

et+h|t=∑i=0h−1ψiϵt+h−i,Var(et+h|t)=σ2∑i=0h−1ψi2. e_{t+h|t} = \sum_{i=0}^{h-1} \psi_i\, \epsilon_{t+h-i}, \qquad \operatorname{Var}(e_{t+h|t}) = \sigma^2 \sum_{i=0}^{h-1} \psi_i^2 .

The ψi\psi_i are Module 3’s MA(∞\infty) weights, Ψ(L)=Θ(L)/Φ(L)\Psi(L) = \Theta(L)/\Phi(L), and ARMAtoMA() computes them for any model. For an AR(1), ψi=ϕi\psi_i = \phi^i, the same geometry as its ACF, so the sum has a closed form:

Var(et+h|t)=σ21−ϕ2h1−ϕ2→σ21−ϕ2=γ0. \operatorname{Var}(e_{t+h|t}) = \sigma^2\,\frac{1 - \phi^{2h}}{1 - \phi^2} \;\longrightarrow\; \frac{\sigma^2}{1 - \phi^2} = \gamma_0 .

At h=1h = 1 it is σ2\sigma^2, the innovation variance you have tracked since Module 1. As h→∞h \to \infty it is γ0\gamma_0, the unconditional variance: as uncertain as knowing nothing at all.

Predict first

ϕ̂=0.723\hat\phi = 0.723 and σ̂2=0.797\hat\sigma^2 = 0.797. What is γ̂0\hat\gamma_0? And at what hh has the forecast-error variance reached 95% of it? Name both numbers.

psi       <- c(1, ARMAtoMA(ar = phi_hat, ma = numeric(0), lag.max = n_future - 1))   # psi_0 = 1
var_h     <- fit_ar1$sigma2 * cumsum(psi^2)
closed    <- fit_ar1$sigma2 * (1 - phi_hat^(2 * (1:n_future))) / (1 - phi_hat^2)
gamma0_hat <- fit_ar1$sigma2 / (1 - phi_hat^2)

se_pkg <- (as.numeric(fc_ar1$upper[, "95%"]) - as.numeric(fc_ar1$mean)) / qnorm(0.975)
c(psi_sum_vs_closed_form = max(abs(var_h - closed)), vs_forecast_pkg = max(abs(se_pkg^2 - var_h)))
psi_sum_vs_closed_form        vs_forecast_pkg 
          2.220446e-16           6.661338e-16 
round(c(var_h1 = var_h[1], var_h5 = var_h[5], var_h20 = var_h[20], gamma0_hat = gamma0_hat), 3)
    var_h1     var_h5    var_h20 gamma0_hat 
     0.797      1.607      1.673      1.673 
which(var_h >= 0.95 * gamma0_hat)[1]
[1] 5

Three ways, one number: your ψ\psi-sum, your closed form and forecast()’s own band agree. Under Gaussian innovations the error is Gaussian too, so the interval is the point plus or minus zz standard deviations: 1.961.96 for 95%, 1.281.28 for 80%. That is the fan chart you have seen on every central bank’s website.

fan_df <- data.frame(t = n_obs + 1:n_future, mean = as.numeric(fc_ar1$mean),
                     lo80 = as.numeric(fc_ar1$lower[, "80%"]), hi80 = as.numeric(fc_ar1$upper[, "80%"]),
                     lo95 = as.numeric(fc_ar1$lower[, "95%"]), hi95 = as.numeric(fc_ar1$upper[, "95%"]))
obs_df <- data.frame(t = 1:n_obs, y = as.numeric(y_ar1))
ggplot() +
  geom_ribbon(data = fan_df, aes(t, ymin = lo95, ymax = hi95), fill = "steelblue4", alpha = 0.18) +
  geom_ribbon(data = fan_df, aes(t, ymin = lo80, ymax = hi80), fill = "steelblue4", alpha = 0.30) +
  geom_hline(yintercept = mu_hat + c(-1, 1) * 1.96 * sqrt(gamma0_hat), linetype = "dashed", colour = "grey45") +
  geom_line(data = obs_df, aes(t, y), colour = "steelblue4") +
  geom_line(data = fan_df, aes(t, mean), colour = "firebrick", linewidth = 0.9) +
  labs(title = "The fan: 80% and 95% bands; dashed lines are the unconditional 95% limits",
       x = "t", y = NULL)

The fan opens fast and then stops, pinned against the dashed lines. That stop is the whole story of a stationary forecast: the share of the variance the model still explains at horizon hh is ϕ2h\phi^{2h}, and the error carries the rest.

Predict first

At ϕ=0.7\phi = 0.7, at what horizon does the model explain less than 10% of the variance of yt+hy_{t+h}? Name a number of steps.

shares <- data.frame(h = c(1, 2, 4, 8), explained_pct = round(100 * 0.7^(2 * c(1, 2, 4, 8)), 1))
shares$left_over_pct <- 100 - shares$explained_pct
shares
  h explained_pct left_over_pct
1 1          49.0          51.0
2 2          24.0          76.0
3 4           5.8          94.2
4 8           0.3          99.7
c(first_h_below_10pct_at_0.7 = which(0.7^(2 * (1:20)) < 0.10)[1],
  at_the_fitted_phi          = which(phi_hat^(2 * (1:20)) < 0.10)[1])
first_h_below_10pct_at_0.7          at_the_fitted_phi 
                         4                          4

Half the variance is still explained at h=ln0.5/(2lnϕ)h = \ln 0.5 / (2 \ln \phi): about one step at ϕ=0.7\phi = 0.7, about seven at ϕ=0.95\phi = 0.95. Persistence is the dial.

So what: the interval is honest only if its ingredients are, and the three it rests on are Gaussian innovations, the right model, and a known σ2\sigma^2; section 8 breaks the first one on purpose.

6. Stationary forecasts are short-run forecasts

The twenty sealed observations are about to arrive. Before they do: the single path from T=60T = 60 is what you get if you forecast once and never update. A forecaster who re-forecasts every period sees something else.

Predict first

Two questions. How many of the twenty future values will land inside the 95% band you just drew: 20, 19, 18, fewer? And will the one-step forecasts, each made from its own origin, track the series closely, loosely, or not at all?

y_ar1_all <- ts(c(as.numeric(y_ar1), future))          # t = 1..80, the seal broken

# Rolling forecasts with the coefficients FROZEN at the t = 60 fit. Arima(y, model = fit)
# re-applies the fitted coefficients to a longer series without re-estimating; fitted()
# is then the one-step forecast at every t from information through t - 1.
one_step <- as.numeric(window(fitted(Arima(y_ar1_all, model = fit_ar1)), start = n_obs + 1))
two_step <- sapply((n_obs + 1):(n_obs + n_future), function(t)
  as.numeric(forecast(Arima(window(y_ar1_all, end = t - 2), model = fit_ar1), h = 2)$mean)[2])
single   <- as.numeric(fc_ar1$mean)

inside95 <- future >= as.numeric(fc_ar1$lower[, "95%"]) & future <= as.numeric(fc_ar1$upper[, "95%"])
inside80 <- future >= as.numeric(fc_ar1$lower[, "80%"]) & future <= as.numeric(fc_ar1$upper[, "80%"])
c(inside_95 = sum(inside95), inside_80 = sum(inside80), of = n_future)
inside_95 inside_80        of 
       19        15        20 
roll_df <- rbind(data.frame(t = 1:n_obs, value = as.numeric(y_ar1), what = "observed"),
                 data.frame(t = n_obs + 1:n_future, value = future, what = "realised future"),
                 data.frame(t = n_obs + 1:n_future, value = single,   what = "single path from T"),
                 data.frame(t = n_obs + 1:n_future, value = one_step, what = "rolling one-step"),
                 data.frame(t = n_obs + 1:n_future, value = two_step, what = "rolling two-step"))
ggplot(roll_df, aes(t, value, colour = what)) +
  geom_ribbon(data = fan_df, aes(t, ymin = lo95, ymax = hi95), inherit.aes = FALSE,
              fill = "steelblue4", alpha = 0.15) +
  geom_line() +
  scale_colour_manual(values = c("observed" = "steelblue4", "realised future" = "black",
                                 "single path from T" = "grey55", "rolling one-step" = "firebrick",
                                 "rolling two-step" = "darkorange"), name = NULL) +
  coord_cartesian(xlim = c(40, 80)) +
  labs(title = "Forecast once, or re-forecast every step", x = "t", y = NULL) +
  theme(legend.position = "bottom")

round(c(one_step = rmse(future, one_step), two_step = rmse(future, two_step),
        single_path = rmse(future, single), mean_only = rmse(future, rep(mu_hat, n_future))), 3)
   one_step    two_step single_path   mean_only 
      0.983       1.296       1.259       1.264

19 of 20 inside the 95% band. The one-step forecasts track the series one step late, as a one-step AR(1) forecast must; the two-step ones track it less; the single path stops tracking after a few steps and is no better than forecasting μ̂\hat\mu for all twenty. The ranking one-step < two-step is the theory (population RMSEs σ\sigma and σ1+ϕ2\sigma\sqrt{1+\phi^2}); the rest is one draw of twenty.

Two ways to keep score. The single path is a fixed-origin evaluation: fit once, forecast the whole test horizon. The rolling lines are a rolling-origin evaluation: advance the origin one step at a time and re-forecast. If each new observation is added to the training sample, the window is expanding; if the oldest is dropped as the newest arrives, it is a fixed-width rolling window. Both can re-estimate the coefficients at every origin, or freeze them, as this chunk did. Section 7 is an expanding window with frozen coefficients; Module 7 re-estimates.

So what: a stationary model forecasts the next few steps; after that it forecasts the mean, and ϕ\phi sets how many steps “a few” is.

7. The example meets its test set

Your sample ended in 2006 Q4. What came next? Ben Bernanke, testimony before the Joint Economic Committee, 28 March 2007: “At this juncture, however, the impact on the broader economy and financial markets of the problems in the subprime market seems likely to be contained.” You are about to make the forecast the Fed could not, from the same vantage point, with the two models from section 2, and watch 2007 to 2011 arrive one quarter at a time.

Predict first

Write three things in the log. (i) At h=1h = 1, under squared error, which finalist wins: ARIMA(2,1,0) or ARIMA(1,1,0)? (ii) Does either beat a random walk with drift, the forecast that says “growth will be its 1990–2006 average, every quarter”? (iii) How many of the twenty realised quarters land inside the ARIMA(2,1,0)’s 95% band: name a number. The band was drawn by a model that has never seen a recession.

The coefficients stay frozen at their 2006 Q4 values. At every origin the model is applied to the data through that origin (the model = fit trick from section 6) and forecasts one and two quarters ahead. The benchmark is the random walk with drift: ŷt+h|t=yt+hδ‾\hat y_{t+h|t} = y_t + h\bar\delta, with Gaussian bands from the training variance of growth.

y_full <- log(window(gdp_all, start = c(1990, 1), end = c(2011, 4)))     # now the test quarters are read
test_targets <- seq(2007, 2011.75, by = 0.25)                            # 2007 Q1 .. 2011 Q4
q_label <- function(tt) sprintf("%d Q%d", floor(tt + 1e-8), round((tt %% 1) * 4) + 1)
at <- function(x, tt) as.numeric(window(x, start = tt, end = tt))

realised_log    <- sapply(test_targets, function(tt) at(y_full, tt))
realised_growth <- 100 * (realised_log - sapply(test_targets - 0.25, function(tt) at(y_full, tt)))
realised_level  <- exp(realised_log)

# h-step forecasts of every target with the coefficients fixed at the training fit.
rolling_fixed <- function(fit, targets, h, y_series) {
  rows <- t(sapply(targets, function(tt) {
    y_to_origin <- window(y_series, end = tt - h / 4)
    fc <- forecast(Arima(y_to_origin, model = fit), h = h, level = c(80, 95))
    c(point = as.numeric(fc$mean)[h], point1 = as.numeric(fc$mean)[1],
      lo80 = unname(fc$lower[h, 1]), hi80 = unname(fc$upper[h, 1]),
      lo95 = unname(fc$lower[h, 2]), hi95 = unname(fc$upper[h, 2]),
      y_origin = as.numeric(y_to_origin)[length(y_to_origin)])
  }))
  d <- data.frame(target = q_label(targets), rows)
  d$growth_hat <- if (h == 1) 100 * (d$point - d$y_origin) else 100 * (d$point - d$point1)
  d
}
rw_drift_fixed <- function(targets, h, y_series, drift_bar, sd_growth) {
  y_origin <- sapply(targets - h / 4, function(tt) at(y_series, tt))
  point <- y_origin + h * drift_bar; se <- sqrt(h) * sd_growth
  data.frame(target = q_label(targets), point = point, point1 = y_origin + drift_bar,
             lo80 = point - qnorm(0.9) * se, hi80 = point + qnorm(0.9) * se,
             lo95 = point - qnorm(0.975) * se, hi95 = point + qnorm(0.975) * se,
             y_origin = y_origin, growth_hat = 100 * drift_bar)
}
score_track <- function(d, realised_log, realised_growth, targets, worst = targets >= 2008.5 & targets <= 2009.25) {
  in95 <- realised_log >= d$lo95 & realised_log <= d$hi95
  in80 <- realised_log >= d$lo80 & realised_log <= d$hi80
  e    <- realised_growth - d$growth_hat
  c(RMSE_growth = sqrt(mean(e^2)), inside_95 = sum(in95), inside_80 = sum(in80),
    RMSE_2008Q3_2009Q2 = sqrt(mean(e[worst]^2)), inside_95_in_those_4 = sum(in95[worst]))
}

tracks <- list()
for (h in 1:2) {
  tracks[[paste0("ARIMA(2,1,0) h=", h)]] <- rolling_fixed(fit_210, test_targets, h, y_full)
  tracks[[paste0("ARIMA(1,1,0) h=", h)]] <- rolling_fixed(fit_110, test_targets, h, y_full)
  tracks[[paste0("RW + drift h=", h)]]   <- rw_drift_fixed(test_targets, h, y_full, mean(dy), sd(dy))
}
scoreboard <- t(sapply(tracks, score_track, realised_log, realised_growth, test_targets))
round(scoreboard, 3)
                 RMSE_growth inside_95 inside_80 RMSE_2008Q3_2009Q2
ARIMA(2,1,0) h=1       0.767        18        15              1.453
ARIMA(1,1,0) h=1       0.792        17        14              1.491
RW + drift h=1         0.932        16        13              1.874
ARIMA(2,1,0) h=2       0.839        18        15              1.648
ARIMA(1,1,0) h=2       0.897        18        15              1.800
RW + drift h=2         0.932        17        12              1.874
                 inside_95_in_those_4
ARIMA(2,1,0) h=1                    3
ARIMA(1,1,0) h=1                    2
RW + drift h=1                      1
ARIMA(2,1,0) h=2                    2
ARIMA(1,1,0) h=2                    2
RW + drift h=2                      1

Read the first column against your log. At h=1h = 1 the ARIMA(2,1,0) wins under squared error, 0.767 against 0.792, and both beat the random walk’s 0.932: the margin between the two ARIMAs is small next to the gap to the benchmark. The criteria disagreed in section 2c; out of sample the two models are near twins. The 95% band caught 18 of 20, which sounds like a pass until you ask which quarters it missed.

t210 <- tracks[["ARIMA(2,1,0) h=1"]]
readout <- data.frame(target = t210$target, realised = round(realised_growth, 2),
                      forecast = round(t210$growth_hat, 2),
                      lo95 = round(100 * (t210$lo95 - t210$y_origin), 2),
                      hi95 = round(100 * (t210$hi95 - t210$y_origin), 2))
readout$inside <- realised_growth >= 100 * (t210$lo95 - t210$y_origin) & realised_growth <= 100 * (t210$hi95 - t210$y_origin)
readout[5:10, ]
    target realised forecast  lo95 hi95 inside
5  2008 Q1    -0.65     0.40 -0.62 1.43  FALSE
6  2008 Q2     0.37     0.18 -0.84 1.20   TRUE
7  2008 Q3    -0.77     0.18 -0.84 1.20   TRUE
8  2008 Q4    -2.45     0.15 -0.87 1.17  FALSE
9  2009 Q1    -1.34    -0.50 -1.52 0.52   TRUE
10 2009 Q2    -0.38    -0.64 -1.66 0.38   TRUE

2008 Q4: the model said 0.15% with a band from -0.87 to 1.17, and the quarter came in at -2.45%. An error of -2.60 points, outside every model’s band at both horizons. That is the quarter the March 2007 testimony was about, and no model frozen in 2006 could see it. Over the four worst quarters, 2008 Q3 to 2009 Q2, the ARIMA(2,1,0)’s RMSE is 1.45 and only 3 of the four are inside the band. Here is the same thing in dollars.

tail_times <- seq(2003, 2006.75, by = 0.25)
lev_df <- rbind(
  data.frame(date = tail_times, level = sapply(tail_times, function(tt) at(gdp_all, tt)), what = "realised"),
  data.frame(date = test_targets, level = realised_level, what = "realised"))
band_df <- data.frame(date = test_targets, point = exp(t210$point), lo95 = exp(t210$lo95), hi95 = exp(t210$hi95))
ggplot() +
  geom_ribbon(data = band_df, aes(date, ymin = lo95, ymax = hi95), fill = "steelblue4", alpha = 0.2) +
  geom_line(data = band_df, aes(date, point), colour = "steelblue4", linewidth = 0.9) +
  geom_line(data = lev_df, aes(date, level), colour = "black", linewidth = 0.7) +
  geom_point(data = lev_df[lev_df$date >= 2007, ], aes(date, level), colour = "black", size = 1.4) +
  geom_vline(xintercept = 2006.875, colour = "grey60") +
  labs(title = "Real GDP per capita: one-step forecasts (blue) and their 95% band against what happened",
       x = NULL, y = "chained 2017 dollars per person")

The band swallows the recession late: each one-step forecast starts from the latest observed level, so the model is dragged down after each fall, never before it. A one-step forecaster is never far wrong for long and never early.

Modify and re-run

Change inside_95 to the 80% band by swapping lo95/hi95 for lo80/hi80 in score_track(), and re-run the scoreboard. Before you do: an 80% band should miss about four of twenty. Will the random walk’s count fall by more or less than the ARIMAs’ counts, and why would a band that never moves miss differently from one that follows the series?

So what: the model that won the IC table on 1990–2006 also won the test set on 2007–2011, by a little, and the band it drew missed the one quarter that mattered most, because the information set at 2006 Q4 contained no such quarter.

8. The assumption-break: coverage under fat tails

Would you trust a 95% band you had never seen fail? Section 5 said the 95% rests on three things: the right model, a known σ2\sigma^2, and Gaussian innovations. Keep the first two honest and break the third, on data you generated. Simulate an AR(1), fit the correct model, build one-step 80%, 95% and 99% intervals, draw the true yT+1y_{T+1}, and count how often it lands inside. First with Gaussian shocks, then with Student-t3t_3 shocks rescaled to the same variance, so the two runs differ only in the shape of the tails.

coverage_mc <- function(n_reps, n_obs, innov, phi = 0.7, levels = c(80, 95, 99)) {
  hits <- matrix(NA, n_reps, length(levels), dimnames = list(NULL, paste0(levels, "% band")))
  for (r in seq_len(n_reps)) {
    eps <- innov(n_obs + 1); z <- numeric(n_obs + 1)
    for (t in 2:(n_obs + 1)) z[t] <- phi * z[t - 1] + eps[t]
    fit <- Arima(ts(z[1:n_obs]), order = c(1, 0, 0), include.mean = FALSE)
    fc  <- forecast(fit, h = 1, level = levels)
    hits[r, ] <- z[n_obs + 1] >= as.numeric(fc$lower) & z[n_obs + 1] <= as.numeric(fc$upper)
  }
  colMeans(hits)
}
set.seed(2015)
cov_gauss <- coverage_mc(n_reps = 500, n_obs = 200, innov = rnorm)
round(cov_gauss, 3)
80% band 95% band 99% band 
   0.820    0.954    0.990

Close to 80, 95 and 99: when the assumption holds, the band keeps its promise (the small shortfalls are the price you pay for estimating ϕ̂\hat\phi and σ̂2\hat\sigma^2 instead of knowing them). Now break it.

Predict first

Same variance, fatter tails. Does the 80% band’s coverage go up or down under t3t_3? And the 99% band’s? Commit to a direction for each, and if you think they move the same way, say why.

set.seed(2015)
cov_t3 <- coverage_mc(n_reps = 500, n_obs = 200, innov = function(n) rt(n, df = 3) / sqrt(3))
round(rbind(gaussian = cov_gauss, t3_rescaled = cov_t3), 3)
            80% band 95% band 99% band
gaussian        0.82    0.954    0.990
t3_rescaled     0.88    0.940    0.966

They move in opposite directions. The 80% band over-covers (88% here) and the 99% band under-covers (96.6%), with the 95% band sitting near the crossing. No simulation was needed to see it coming: with the parameters known, the true coverage of a Gaussian-width band under rescaled t3t_3 shocks is 88.7%, 95.7%, 97.9% at the three levels. A standardized t3t_3 puts more mass near the center and in the far tails than a Gaussian with the same variance, with less in the shoulders. The 80% band reaches past that fat center, so it catches too much; the 99% band cannot reach the far tails, so it catches too little. The miss rate the 99% band promises is one in a hundred; the exact calculation says it delivers one in about 48, and this run of 500 draws did a little worse than that.

That crossing is a fact about this distribution and these levels, not a law for every fat tail. The general lesson: the shape of your prediction interval is a property of your distributional assumption relative to the data, Module 5’s rule in forecasting clothes, and where it fails is the far tail, which is exactly where risk management lives.

Was 2008 Q4 a fat-tail event?

Section 7’s band was drawn from a model whose innovations were assumed Gaussian. Were the training residuals? And how far out was the miss?

e_train  <- residuals(fit_210)[-1]
kurtosis <- mean(e_train^4) / mean(e_train^2)^2                 # 3 for a Gaussian
z_2008q4 <- (realised_growth[8] - t210$growth_hat[8]) / (100 * sqrt(fit_210$sigma2))
round(c(training_residual_kurtosis = kurtosis, z_score_2008Q4 = z_2008q4,
        gaussian_prob_of_a_miss_that_big = 2 * pnorm(-abs(z_2008q4))), 7)
      training_residual_kurtosis                   z_score_2008Q4 
                       2.9670664                       -4.9917215 
gaussian_prob_of_a_miss_that_big 
                       0.0000006

The training residuals have kurtosis 2.97: not fat-tailed at all, as Gaussian as 67 quarters can look. And 2008 Q4 was a 5-sigma miss, an event the band priced at about one in 1.7 million quarters. The tails the band could not see were not in the 1990–2006 residuals. They were in a regime those years never contained. A fat-tail fix would not have saved this forecast; nothing frozen in 2006 would have.

Modify and re-run

Change df = 3 to df = 5 in the coverage-fat chunk, rescaling by sqrt(5/3) instead of sqrt(3) so the variance stays 1. Before you do: will the 99% band’s shortfall shrink or grow, and will the 80% band’s over-coverage move the same way? Then try df = 30 and say what you expect before you look.

So what: a 95% band is a conditional promise, and the condition it breaks on first is the one you cannot check from inside the training sample.

9. How do we judge ourselves if we are always wrong?

You now have three forecast tracks for the same twenty quarters and a scoreboard that used squared error without asking whether squared error was the right question. Was it? A loss function is the cost of being wrong: g(e)g(e), a function of the forecast error, zero when the forecast is perfect and larger as the miss grows. (It is gg, not LL; LL is the lag operator in this course.) The loss is chosen to reflect your cost function, and chosen before you see the results, because picking the loss after the scores are in is picking the winner.

Think about the Texas grid. You forecast electricity demand and a cold snap hits. Under-forecast, and there is not enough supply: people freeze. Over-forecast, and you ran extra generators: wasted fuel. Both are errors of 500 megawatts; they are nothing like the same error, and a symmetric loss says they are. So ask what each additional unit of miss costs you, and let that marginal cost pick the shape:

What one more unit of error costs Loss Target it optimizes
More than the last unit did (one catastrophe dominates) MSE, 1N∑(yt−ŷt)2\frac{1}{N}\sum (y_t - \hat y_t)^2 conditional mean
The same as the last unit (errors accumulate proportionally) MAE, 1N∑|yt−ŷt|\frac{1}{N}\sum \lvert y_t - \hat y_t \rvert conditional median
Rising for small errors, flat beyond a threshold δ\delta Huber, gδ(e)=12e2g_\delta(e) = \tfrac12 e^2 if |e|≤δ\lvert e\rvert \le \delta, else δ(|e|−12δ)\delta(\lvert e\rvert - \tfrac12\delta) between the two
More in one direction than the other Asymmetric linear, α|e|\alpha\lvert e\rvert if e>0e > 0 (under-prediction), (1−α)|e|(1-\alpha)\lvert e\rvert otherwise the α\alpha-quantile

NN counts test-set observations, never the sample size TT. Draw the four shapes and, underneath them, what the next unit of error costs you.

e <- seq(-4, 4, by = 0.01); delta <- 2; alpha <- 0.75
loss_df <- rbind(
  data.frame(e, value = e^2,                                      loss = "MSE (quadratic)",            panel = "loss g(e)"),
  data.frame(e, value = abs(e),                                   loss = "MAE (absolute)",             panel = "loss g(e)"),
  data.frame(e, value = huber(e, 0, delta),                       loss = "Huber, delta = 2",           panel = "loss g(e)"),
  data.frame(e, value = ifelse(e > 0, alpha, 1 - alpha) * abs(e), loss = "Asymmetric, alpha = 0.75",   panel = "loss g(e)"),
  data.frame(e, value = 2 * abs(e),                               loss = "MSE (quadratic)",            panel = "marginal cost of one more unit"),
  data.frame(e, value = 1,                                        loss = "MAE (absolute)",             panel = "marginal cost of one more unit"),
  data.frame(e, value = pmin(abs(e), delta),                      loss = "Huber, delta = 2",           panel = "marginal cost of one more unit"),
  data.frame(e, value = ifelse(e > 0, alpha, 1 - alpha),          loss = "Asymmetric, alpha = 0.75",   panel = "marginal cost of one more unit"))
ggplot(loss_df, aes(e, value, colour = loss)) +
  geom_line(linewidth = 0.9) +
  facet_wrap(~ panel, ncol = 1, scales = "free_y") +
  labs(title = "Four ways to be wrong, and what the next unit of error costs under each",
       x = "forecast error e = y - y-hat (e > 0 is under-prediction)", y = NULL, colour = NULL) +
  theme(legend.position = "bottom")

The quadratic explodes: off by 4 is sixteen times as bad as off by 1. The absolute grows steadily. Huber rises like the quadratic and then its marginal cost stops at δ\delta, so no single outlier can dominate the evaluation. The asymmetric loss charges 0.750.75 per unit when the actual comes in above the forecast and 0.250.25 when it comes in below. The one loss not on the picture is MAPE, 100N∑|(yt−ŷt)/yt|\frac{100}{N}\sum \lvert (y_t - \hat y_t)/y_t \rvert: scale-free, popular in business, and undefined at yt=0y_t = 0. On a growth rate it is worse than undefined:

round(c(MAPE_of_ARIMA210_on_growth = mape(realised_growth, t210$growth_hat),
        smallest_realised_growth_pct = min(abs(realised_growth))), 2)
  MAPE_of_ARIMA210_on_growth smallest_realised_growth_pct 
                       99.30                         0.08

A 99% error, driven by quarters where growth was a few hundredths of a point and any forecast at all is “off by hundreds of percent”. Said once: never MAPE on a series that passes near zero.

The loss picks the target

Section 3 called the conditional mean the best guess under squared-error loss. Here is why the qualifier was there. Under MAE the best point forecast is the conditional median; under asymmetric loss with weight α\alpha it is the α\alpha-quantile, which is why the Texas operator at α=0.9\alpha = 0.9 forecasts the 90th percentile of demand, deliberately high. For symmetric shocks the mean and the median coincide and you would never notice. Skew the shocks and the loss decides the contest.

Predict first

Centered-exponential shocks: most are small and negative, a few are large and positive. Two one-step forecast tracks on the same fitted AR(1): the model’s conditional mean, and the same track shifted down by the median of the training residuals. Which wins under MSE, and which under MAE? Say both before you look.

set.seed(2020)
n_total <- 600; n_train <- 350
eps_sk <- rexp(n_total) - 1                              # mean 0, skewed right
y_sk <- numeric(n_total); for (t in 2:n_total) y_sk[t] <- 0.6 * y_sk[t - 1] + eps_sk[t]
y_sk <- ts(y_sk); y_sk_train <- window(y_sk, end = n_train); y_sk_test <- window(y_sk, start = n_train + 1)

fit_sk    <- Arima(y_sk_train, order = c(1, 0, 0))
fc_mean   <- window(fitted(Arima(y_sk, model = fit_sk)), start = n_train + 1)   # frozen, one-step
med_shift <- median(residuals(fit_sk))
fc_median <- fc_mean + med_shift
actual_sk <- as.numeric(y_sk_test)

round(c(phi_hat = unname(coef(fit_sk)["ar1"]), median_shift = med_shift,
        share_of_actuals_below_mean_track = mean(actual_sk < as.numeric(fc_mean))), 3)
                          phi_hat                      median_shift 
                            0.578                            -0.391 
share_of_actuals_below_mean_track 
                            0.612 
round(rbind(MSE = c(mean_track = mse(actual_sk, as.numeric(fc_mean)), median_track = mse(actual_sk, as.numeric(fc_median))),
            MAE = c(mean_track = mae(actual_sk, as.numeric(fc_mean)), median_track = mae(actual_sk, as.numeric(fc_median)))), 3)
    mean_track median_track
MSE      1.063        1.207
MAE      0.753        0.727

Each wins under the loss it was built for. Not a paradox: the loss doing its job. MSE asked which track is closest in squared distance, where the occasional big positive shock is expensive, and the mean, which leans toward those shocks, wins. MAE asked which is closest one absolute unit at a time, and the median, which sits where 61% of the actuals actually land (below the mean track), wins. Different questions, different answers.

Within a loss, never between. Two models are compared under one loss, on one series, over the same NN. “A wins on MSE, B wins on MAE” is two questions with two answers, and it is the rule for the next section. Nothing is sacred about these four, either: the notes write a loss for a call center where under-staffing costs more than over-staffing and escalates, and exploration 3 asks you to write yours.

So what: the point forecast you report is the answer to a question the loss function asked, so decide what being wrong costs you before you decide what to forecast.

10. The same forecasts, scored under each loss

Section 7 ranked the three tracks under one loss. Does the ranking survive the others? You have the errors already; this is pure arithmetic.

Predict first

Which loss, if any, makes the random walk with drift the winner at h=1h = 1? Think about what the random walk forecasts (the training average, every quarter) and which direction it is never wrong in. Then name the loss, and if it is the asymmetric one, name the α\alpha at which the flip happens.

g <- realised_growth
score_under <- function(f) c(MSE = mse(g, f), MAE = mae(g, f), "Huber (delta 2)" = mean_huber(g, f, 2),
                              "Asym 0.5" = asym_loss(g, f, 0.5), "Asym 0.7" = asym_loss(g, f, 0.7),
                              "Asym 0.8" = asym_loss(g, f, 0.8))
for (h in 1:2) {
  cat("\nh =", h, "\n")
  tr <- list("ARIMA(2,1,0)" = tracks[[paste0("ARIMA(2,1,0) h=", h)]]$growth_hat,
             "ARIMA(1,1,0)" = tracks[[paste0("ARIMA(1,1,0) h=", h)]]$growth_hat,
             "RW + drift"   = tracks[[paste0("RW + drift h=", h)]]$growth_hat)
  sc <- sapply(tr, score_under)
  print(round(sc, 3))
  print(apply(sc, 1, function(r) names(r)[which.min(r)]))
}

h = 1 
                ARIMA(2,1,0) ARIMA(1,1,0) RW + drift
MSE                    0.589        0.627      0.869
MAE                    0.510        0.543      0.613
Huber (delta 2)        0.285        0.306      0.413
Asym 0.5               0.255        0.272      0.306
Asym 0.7               0.205        0.204      0.209
Asym 0.8               0.180        0.171      0.160
            MSE             MAE Huber (delta 2)        Asym 0.5        Asym 0.7 
 "ARIMA(2,1,0)"  "ARIMA(2,1,0)"  "ARIMA(2,1,0)"  "ARIMA(2,1,0)"  "ARIMA(1,1,0)" 
       Asym 0.8 
   "RW + drift" 

h = 2 
                ARIMA(2,1,0) ARIMA(1,1,0) RW + drift
MSE                    0.704        0.804      0.869
MAE                    0.514        0.574      0.613
Huber (delta 2)        0.335        0.381      0.413
Asym 0.5               0.257        0.287      0.306
Asym 0.7               0.195        0.199      0.209
Asym 0.8               0.164        0.154      0.160
            MSE             MAE Huber (delta 2)        Asym 0.5        Asym 0.7 
 "ARIMA(2,1,0)"  "ARIMA(2,1,0)"  "ARIMA(2,1,0)"  "ARIMA(2,1,0)"  "ARIMA(2,1,0)" 
       Asym 0.8 
 "ARIMA(1,1,0)" 
under_pred <- sapply(list("ARIMA(2,1,0)" = tracks[["ARIMA(2,1,0) h=1"]]$growth_hat,
                          "ARIMA(1,1,0)" = tracks[["ARIMA(1,1,0) h=1"]]$growth_hat,
                          "RW + drift"   = tracks[["RW + drift h=1"]]$growth_hat),
                     function(f) sum(g - f > 0))
under_pred    # quarters each track under-predicted (actual above forecast), h = 1
ARIMA(2,1,0) ARIMA(1,1,0)   RW + drift 
           8            6            5

Under every symmetric loss the ranking does not move: ARIMA(2,1,0), then ARIMA(1,1,0), then the benchmark, at both horizons. Asymmetric loss is a different question. At α=0.5\alpha = 0.5 it is half the MAE and nothing changes. At α=0.7\alpha = 0.7 the ARIMA(1,1,0) edges ahead at h=1h = 1, a near tie. At α=0.8\alpha = 0.8 the random walk with drift wins at h=1h = 1. Why? Think about what you asked each track to do. The random walk forecasts the training average every quarter and never forecasts below it, so when growth rebounds in 2009 and 2010 it is the track that under-predicts least (5 quarters against 8 for the ARIMA(2,1,0), which follows each fall down and then under-shoots the recovery). When under-prediction costs four times what over-prediction costs, the forecast that never goes low is the cheap one. That is “the loss reflects your cost function” on real data: the Texas grid’s row of the table, in GDP.

So what: the winner of a forecasting contest is a property of the loss as much as of the models, so a scoreboard without a stated loss is not a result.

11. The same workflow on the Module 5 handoff

Module 5 ended on a diagnostic-clean ARIMA(1,1,2)(1,0,1)12(1,1,2)(1,0,1)_{12} on the unemployment rate and one question: what does it say about next month? The workflow you have just run on GDP, in four moves: split, fit on training, forecast and overlay, score against a baseline. The baseline is the naive forecast, last value carried forward.

unrate_raw <- read.csv("data/UNRATE.csv", stringsAsFactors = FALSE)
unrate_raw <- unrate_raw[unrate_raw$observation_date >= "1960-01-01" &
                         unrate_raw$observation_date <= "2019-12-01", ]
unrate <- ts(unrate_raw$UNRATE, start = c(1960, 1), frequency = 12)
sp <- ts_split(unrate, test_periods = 24)                   # hold out the last two years
c(train_end = paste(end(sp$train), collapse = "-"), test = paste(c(start(sp$test), end(sp$test)), collapse = "-"))
       train_end             test 
       "2017-12" "2018-1-2019-12"
Predict first

The training set ends in December 2017 at 4.1%. A model with every slider on, against “it stays at 4.1”. Which has the lower MSE over 2018–2019, and by a lot or a little? Write a number for the ratio.

fit_u <- Arima(sp$train, order = c(1, 1, 2), seasonal = list(order = c(1, 0, 1), period = 12))
e_u   <- residuals(fit_u)
round(c(LB10_p = Box.test(e_u, lag = 10, type = "Ljung-Box", fitdf = 5)$p.value,
        LB24_p = Box.test(e_u, lag = 24, type = "Ljung-Box", fitdf = 5)$p.value), 3)
LB10_p LB24_p 
 0.531  0.282

Still clean on the training set alone; if it were not, you would be back in Module 5. Now forecast two years and lay the truth over it.

fc_u    <- forecast(fit_u, h = 24, level = c(80, 95))
naive_u <- naive(sp$train, h = 24, level = c(80, 95))
u_df <- data.frame(date = as.numeric(time(sp$test)), actual = as.numeric(sp$test),
                   sarima = as.numeric(fc_u$mean), naive = as.numeric(naive_u$mean),
                   lo95 = as.numeric(fc_u$lower[, "95%"]), hi95 = as.numeric(fc_u$upper[, "95%"]),
                   lo80 = as.numeric(fc_u$lower[, "80%"]), hi80 = as.numeric(fc_u$upper[, "80%"]))
train_tail <- window(sp$train, start = c(2015, 1))
ggplot() +
  geom_ribbon(data = u_df, aes(date, ymin = lo95, ymax = hi95), fill = "steelblue4", alpha = 0.15) +
  geom_ribbon(data = u_df, aes(date, ymin = lo80, ymax = hi80), fill = "steelblue4", alpha = 0.25) +
  geom_line(data = data.frame(date = as.numeric(time(train_tail)), y = as.numeric(train_tail)), aes(date, y), colour = "grey30") +
  geom_line(data = u_df, aes(date, sarima, colour = "SARIMA point forecast"), linewidth = 0.9) +
  geom_line(data = u_df, aes(date, naive,  colour = "naive: last value"), linetype = "dashed") +
  geom_line(data = u_df, aes(date, actual, colour = "actual"), linewidth = 0.8) +
  scale_colour_manual(values = c("SARIMA point forecast" = "steelblue4", "naive: last value" = "grey45", "actual" = "firebrick"), name = NULL) +
  labs(title = "UNRATE: 24-month forecast from the Module 5 model, with the fan and the truth", x = NULL, y = "percent") +
  theme(legend.position = "bottom")

actual_u <- as.numeric(sp$test)
round(rbind(SARIMA = c(MSE = mse(actual_u, u_df$sarima), MAE = mae(actual_u, u_df$sarima), RMSE = rmse(actual_u, u_df$sarima)),
            naive  = c(MSE = mse(actual_u, u_df$naive),  MAE = mae(actual_u, u_df$naive),  RMSE = rmse(actual_u, u_df$naive))), 4)
          MSE    MAE   RMSE
SARIMA 0.1642 0.3247 0.4052
naive  0.1275 0.3167 0.3571
c(inside_95 = sum(actual_u >= u_df$lo95 & actual_u <= u_df$hi95), of = 24)
inside_95        of 
       24        24

The naive forecast wins on every metric. The model that survived Module 5’s entire gauntlet loses to “carry the last number forward”, because 2018–2019 was a calm glide from 4.1 to 3.5 with no turning point for the model’s machinery to earn its keep on. That is why the random walk is the universal benchmark: on a persistent series in a calm stretch it is not a strawman. All 24 actuals sat inside the 95% band, which is comfort on one path, not proof the band is calibrated. And the margin is 0.008 on MAE, under a hundredth of a point, from one 24-month window: whether that gap means anything is a question for a test, and Module 7 has one.

So what: a model can describe the in-sample dependence adequately and still add nothing over a simpler rule on a stretch with nothing to exploit, and losing to the naive once is information, not an indictment.

12. Explorations

Optional, and outside the assessable surface

Nothing in this section is required, and none of it appears on problem sets or exams. Assessments draw only on the Core tier of the module notes. These are here because the interesting questions usually start where the syllabus stops.

If one of these grabs you, chase it. If a different question grabs you, chase that one instead — that is a better outcome than any of the six below. Your AI tutor’s job in this section is to help you turn a hunch into an experiment you can actually run, then connect what you find back to a Core concept.

1. Can you beat the random walk on unemployment? (guided)

Section 11 was one fixed-origin forecast. Give the model a fairer fight: re-forecast every month of 2019 from an expanding window, coefficients frozen, one step ahead, against the naive and against the seasonal naive (the same month last year).

origins <- seq(2018 + 11/12, 2019 + 10/12, by = 1/12)            # Dec 2018 .. Nov 2019
one_u <- sapply(origins, function(o) as.numeric(forecast(Arima(window(unrate, end = o), model = fit_u), h = 1)$mean))
targets_u <- as.numeric(window(unrate, start = 2019, end = c(2019, 12)))
naive_1   <- as.numeric(window(unrate, start = origins[1], end = origins[12]))
snaive_1  <- as.numeric(window(unrate, start = 2018, end = c(2018, 12)))
round(c(sarima_one_step = rmse(targets_u, one_u), naive = rmse(targets_u, naive_1), seasonal_naive = rmse(targets_u, snaive_1)), 3)
sarima_one_step           naive  seasonal_naive 
          0.097           0.096           0.238

Does one-step-ahead change the verdict? Why would the seasonal naive do so badly on a seasonally adjusted series, and what does that tell you about what Module 5’s seasonal band is absorbing?

2. Was 2007–2011 a fair test, or the cruelest one? (guided)

Move the window. rolling_fixed() takes any set of targets and any fit; give it 2012 Q1 to 2016 Q4, then 2017 Q1 to 2019 Q4, with the 2006 coefficients still frozen, and then with a fresh fit on everything up to the quarter before the window. Does the ARIMA(2,1,0) still beat the ARIMA(1,1,0)? Does either still beat the random walk? Does the band still miss?

y_long <- log(window(gdp_all, start = c(1990, 1), end = c(2019, 4)))
calm <- seq(2012, 2016.75, by = 0.25)
calm_log <- sapply(calm, function(tt) at(y_long, tt))
calm_growth <- 100 * (calm_log - sapply(calm - 0.25, function(tt) at(y_long, tt)))
dy_pre <- diff(window(y_long, end = calm[1] - 0.25))
round(rbind(
  "ARIMA(2,1,0), 2006 coefficients" = score_track(rolling_fixed(fit_210, calm, 1, y_long), calm_log, calm_growth, calm, worst = rep(TRUE, length(calm))),
  "ARIMA(1,1,0), 2006 coefficients" = score_track(rolling_fixed(fit_110, calm, 1, y_long), calm_log, calm_growth, calm, worst = rep(TRUE, length(calm))),
  "RW + drift, 2006 drift"          = score_track(rw_drift_fixed(calm, 1, y_long, mean(dy), sd(dy)), calm_log, calm_growth, calm, worst = rep(TRUE, length(calm)))
)[, 1:3], 3)
                                RMSE_growth inside_95 inside_80
ARIMA(2,1,0), 2006 coefficients       0.445        19        18
ARIMA(1,1,0), 2006 coefficients       0.431        19        18
RW + drift, 2006 drift                0.400        20        19

3. What does your cost function look like? (open)

Pick a decision you actually face: a budget you have to set in advance, a staffing level, how long an exam will take you. Write the cost of under- and over-shooting it as a function of ee, make sure it is zero at e=0e = 0 and grows away from it, and score section 7’s three GDP tracks under it. The notes’ “You Can Write Your Own” is the template. Which track wins under your loss, and would a textbook metric have told you?

4. How fat are the tails you forecast with? (open)

Section 8 used t3t_3 shocks because you could control them. Replace them with the GDP model’s own residuals, resampled with replacement, and run the coverage experiment again. How much coverage does the Gaussian band lose when the innovations are the ones the model actually saw? Then ask the harder question: what would resampling 1990–2006 residuals have told you about 2008?

5. When does re-estimating matter? (guided; Deeper Dive · not assessable)

Section 7 froze the coefficients at 2006 Q4. Re-estimate at every origin instead: replace Arima(y_to_origin, model = fit) in rolling_fixed() with a fresh Arima(y_to_origin, order = c(2, 1, 0), include.drift = TRUE). Does the recession move ϕ̂1\hat\phi_1, ϕ̂2\hat\phi_2 and the drift, and does updating them help or hurt the forecasts during and after it? This is Module 7’s design.

6. Which series do you actually care about? (open, off-syllabus)

Any FRED series you have an opinion about: housing starts, the price of eggs, your state’s unemployment rate. Download it by hand using the GDP file as the template (no key needed), or pull it live with the optional fetch_fred(). Then the whole chain: split, fit, forecast, score, and a baseline you have to beat. Report the band, not just the point.

# Optional. Requires the fredr package and a FRED_KEY in your environment.
source("R/fetch_fred.R")
houst <- fetch_fred("HOUST", start_date = "1990-01-01", frequency = 12)

13. Checks

R/checks.R is the referee: seven assertions about things this module claims are true, in base R so you can read every line. 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")
Module 6 checks
---------------
PASS GDP cache is the 2026-09-30 vintage; 68-quarter window has the quoted growth mean and sd; UNRATE intact  [rows = 318, last = 71212; T = 68, mean growth = 1.871 %/yr, sd = 0.5357 %/qtr; UNRATE window = 720]
PASS Hand-built ADF: trend, k = 2 on 65 rows gives tau = -3.113; growth, drift, k = 1 gives -3.833  [level tau = -3.1128 (n = 65); growth tau = -3.8325 (n = 65)]
PASS AR(1) recursion equals phi^h; variance loop equals the closed form; shares 0.49/0.2401/0.0576/0.0033  [max point gap = 2.8e-17, max variance gap = 4.4e-16, shares = 0.49/0.2401/0.05765/0.003323]
PASS stats::arima on growth: AR(2)+mean logLik = 260.816, AR(1)+mean logLik = 258.998; phi-hats 0.232/0.233 and 0.301  [logLik = 260.8171 / 258.9983; AR(2) phi = 0.2317, 0.2330, mean = 0.00456; AR(1) phi = 0.3008]
PASS Test set 2007-2011, coefficients frozen: AR(2) one-step RMSE = 0.767, RW+drift = 0.932, AR(2) band catches 18 of 20  [RMSE = 0.7674 / 0.9322; inside the 95% band = 18 of 20 (sigma-hat = 0.493 %/qtr)]
PASS huber() is quadratic-then-linear; asym_loss(alpha = 0.5) is MAE/2; MSE picks the mean, MAE the median  [huber gap = 0.0e+00; MSE argmin 0.0157 vs mean 0.0156; MAE argmin -0.3089 vs median -0.3090]
PASS Known AR(1), Gaussian shocks: the 95% one-step band covers within 2 points of 95% over 2,000 draws  [coverage = 0.946]
---------------
7/7 checks passed.

Where you are now

Every slider was already on, and today you used the board. You found that a forecast is a distribution whose headline number is a slide toward the mean at the rate ϕh\phi^h; that its band is honest only on the assumptions it was built from, and the one it breaks on first is the one the training sample cannot show you; that the loss function picks the target and can flip the winner; and that on a calm stretch the naive forecast is a benchmark that earns its reputation.

You now hold one honest forecast object: a point, a band, a loss, and a baseline. The margins you saw today were small — under a hundredth of a point on unemployment, a few hundredths on growth — and from one window each. The open question you carry into Module 7 is the one a statistician asks: does a margin that size mean anything, and when two forecasts disagree, should you pick one or average them? And a question for later still: everything today was one series forecasting itself. What does a forecast from a system of series look like?

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.