Module 6: Forecasting Fundamentals

Econ 6376 — Applied Time Series Econometrics

How to Use These Notes

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

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

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

Core Learning Objectives

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

Distinguish point forecasts, interval forecasts, and density forecasts, and state which question each answers.
Define the $h$-step-ahead forecast $\hat{y}_{t+h|t}$ as the conditional mean $\mathbb{E}[y_{t+h} \mid y_t, y_{t-1}, \ldots]$ and explain why the conditional mean minimizes expected squared forecast error.
Compute the recursive $h$-step-ahead forecast for an AR(1) and an AR(2) by hand — including recovering $\hat{\alpha}$ from the mean that R reports — replacing future innovations with zero and future $y$’s with their forecasts.
Explain why random train/test splits are wrong for time series data, implement a proper temporal split using ts_split(), and describe the rolling-window alternative.
Derive the 1-step and $h$-step forecast-error variance for an AR(1) and state the general result using $\psi$-weights: $\text{Var}(e_{t+h|t}) = \sigma^2 \sum_{i=0}^{h-1} \psi_i^2$.
Construct a Gaussian prediction interval $\hat{y}_{t+h|t} \pm z_{\alpha/2} \cdot \sqrt{\text{Var}(e_{t+h|t})}$ and explain why the interval widens with the forecast horizon.
Run a Monte Carlo that checks empirical coverage of Gaussian prediction intervals at the 80%, 95%, and 99% levels under Gaussian innovations, then repeat with fat-tailed innovations and locate where the Gaussian promise fails — in the far tail, not the middle.
Write down MSE, MAE, MAPE, Huber loss, and asymmetric loss; state which statistic each one optimizes; and give a concrete example where asymmetric loss is the right choice.
Use forecast::forecast() to produce point forecasts, prediction intervals, and fan charts from an Arima object, and read the fan chart correctly.
Execute the complete practical workflow on UNRATE: split $\rightarrow$ fit on training $\rightarrow$ forecast into test $\rightarrow$ overlay actuals $\rightarrow$ compute loss $\rightarrow$ compare against a naive random-walk baseline — and interpret the result honestly, including when the baseline is hard to beat.

6.1 Where We Are on the Mixing Board

Modules 1 through 5 built a complete univariate modeling vocabulary. Module 1 introduced memory and the master equation. Module 2 taught you to test for stationarity. Module 3 gave you the AR and MA model families. Module 4 showed you how to pick a model using information criteria. Module 5 showed you how to check a model using residual diagnostics and added the seasonal band. By the end of Module 5, every knob on the master equation mixing board was on.

Term	L3	L4	L5	L6
$\alpha$, $\delta t$	On	On	On	On
$\phi_j y_{t-j}$ (all $p$)	On	On	On	On
$\theta_l \epsilon_{t-l}$ (all $q$)	On	On	On	On
Seasonal block $(P, D, Q)_s$	Off	Off	On	On
$\epsilon_t$	On	On	On	On

Module 6 is the pivot. Until now the course was about building models — test, identify, fit, diagnose, iterate. Starting today the question changes: you have a model that passes diagnostics. Now what?

The answer is forecasting, and this module develops it as a three-part story:

What is the best guess? The point forecast.
How uncertain is that guess? The prediction interval and the fan chart.
How do we keep score? The loss function.

The L5 Handoff

Module 5 ended with a diagnostic-clean SARIMA on UNRATE — the ARIMA$(1,1,2)(1,0,1)_{12}$ that won the IC table and passed Ljung-Box at both course horizons ($h = 10$ and $h = 24$). The residual ACF was within bands, and the lag-12 spike from the Module 4 non-seasonal fit had vanished. (Recall that the airline model was proposed there too, and rejected for this series: UNRATE is seasonally adjusted at the source, so a seasonal difference is one difference too many, and the fit said so by pinning $\hat{\Theta}_1$ at $-1$.) We said: “you have a model, you have checked it — next the question is what you do with it.” This module does exactly that.

Let us reproduce the L5 fit so we have a working model on screen. To pull the series live, store your FRED key in the environment (never hardcode it in course files) and run:

library(fredr)
fredr_set_key(Sys.getenv("FRED_KEY"))

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

So that these notes render identically for everyone — with or without an API key or an internet connection — the code below loads the same series from the course’s cached copy in data/UNRATE.csv (the same cache Modules 2–5 used) and applies the same 1960–2019 pre-COVID sample window Modules 4 and 5 justified (the April 2020 outlier dominates every sample autocovariance; see Module 4, §4.7.1).

library(forecast); library(ggplot2)
source("../helpers/forecast_utils.R")   # ts_split() and the course loss functions

unrate_csv <- read.csv("../data/UNRATE.csv")  # cached FRED pull; see data/README.md
unrate_csv <- subset(unrate_csv,
                     observation_date >= "1960-01-01" &
                     observation_date <= "2019-12-01")

unrate <- ts(unrate_csv$UNRATE, start = c(1960, 1), frequency = 12)

# The L5 winner — ARIMA(1,1,2)(1,0,1)[12]
fit_l5 <- Arima(unrate, order = c(1, 1, 2),
                seasonal = list(order = c(1, 0, 1), period = 12))

checkresiduals(fit_l5)


    Ljung-Box test

data:  Residuals from ARIMA(1,1,2)(1,0,1)[12]
Q* = 21.176, df = 19, p-value = 0.3272

Model df: 5.   Total lags used: 24

Clean residuals. Ljung-Box passes. Every knob on the mixing board is on. The modeling vocabulary is complete. Today’s question is the most natural one in the world: what does this model say about next month?

6.2 What Is a Forecast?

You Are Always Wrong — The Question Is How Wrong

Here is the first thing to internalize about forecasting: your forecast will be wrong. Every single time. The unemployment rate next month will not be exactly 4.1%. GDP growth next quarter will not be exactly 2.3%. Get used to it.

The question is not whether you are wrong — you are — but how wrong, and whether your wrongness is useful. A forecast that says “unemployment will be between 3% and 12%” is technically correct but useless. A forecast that says “4.1% plus or minus 0.3 percentage points” is almost certainly wrong in the third decimal place, but it is wrong in a way that helps people make decisions. Forecasting is the art of being wrong in a useful, quantifiable, and honest way.

In practice, almost everyone focuses on the point forecast — the single number. When the Federal Reserve publishes the Summary of Economic Projections, the headline is “the Fed expects growth of 2.1%.” The uncertainty bands exist in the footnotes. The full density forecast exists in the research division’s working papers. The point estimate is what lands on desks, what moves markets, and what gets quoted in the Wall Street Journal. That is the reality of applied forecasting, and this course respects it.

But point estimates alone are incomplete. A forecast without a measure of uncertainty is a claim without a confidence interval — it tells your audience what you think will happen without telling them how sure you are. The prediction interval is where the honesty lives. In this course we focus on two objects:

The point forecast $\hat{y}_{t+h|t}$: your best guess for $y_{t+h}$ given what you know at time $t$. This is the headline number — the one that goes in the report, the memo, the briefing.
The prediction interval (PI): a range $[\hat{y}_{t+h|t} - c, \; \hat{y}_{t+h|t} + c]$ that contains $y_{t+h}$ with some stated probability (typically 80% or 95%). This is the fine print — the part that keeps you honest.

There is a third object — the density forecast, the entire conditional distribution $f(y_{t+h} \mid y_t, y_{t-1}, \ldots)$ — from which you can read off any quantile or probability statement. Density forecasts appear in central banking research and financial risk management, and we will see them visually as fan charts. But constructing them by hand is beyond the scope of this course, and in most applied settings the point forecast and the interval are what matter.

The crystal ball analogy. A fortune teller who says “you will be happy” is useless. A fortune teller who says “your unemployment rate next month will be 4.1%, and there is a 95% chance it is between 3.8% and 4.4%” is an econometrician. The crystal ball is the conditional mean; the error bars are the prediction interval. The error bars are where the honesty lives.

The Conditional Mean as the Optimal Point Forecast

Under squared-error loss (MSE), the optimal $h$-step-ahead point forecast is the conditional mean:

\[\hat{y}_{t+h|t} = \mathbb{E}[y_{t+h} \mid y_t, y_{t-1}, \ldots, y_1]\]

The conditional mean is the value that minimizes $\mathbb{E}[(y_{t+h} - c)^2 \mid \mathcal{F}_t]$ over all constants $c$. This is the time-series version of the fact that the population mean minimizes expected squared distance.

This matters because when we compute $\hat{y}_{t+h|t}$ for an ARMA model, we are computing a conditional expectation. That is not a modeling choice — it is the mathematical consequence of choosing MSE as our loss function. If you choose a different loss (MAE, asymmetric), you get a different optimal forecast (the conditional median, a conditional quantile). The loss function picks the target. Section 6.7 develops this — and Section 6.7’s simulation makes the mean-vs-median distinction visible on data.

Technical Note — Why the Conditional Mean Minimizes MSE

For any candidate forecast $c$ (measurable with respect to $\mathcal{F}_t$), add and subtract the conditional mean $m = \mathbb{E}[y_{t+h} \mid \mathcal{F}_t]$:

\[\mathbb{E}[(y_{t+h} - c)^2 \mid \mathcal{F}_t] = \underbrace{\mathbb{E}[(y_{t+h} - m)^2 \mid \mathcal{F}_t]}_{\text{irreducible}} + (m - c)^2\]

The cross term vanishes because $\mathbb{E}[y_{t+h} - m \mid \mathcal{F}_t] = 0$. The first term does not depend on $c$; the second is minimized — set to zero — by choosing $c = m$. The same argument with absolute loss (a case split rather than a square expansion) delivers the conditional median as the MAE-optimal forecast. Hamilton Chapter 4.1 develops the full treatment, including forecasts restricted to linear functions of the data.

Forecasting vs. Inference: An Emphasis Shift

The emphasis in Modules 1–5 was on $\hat{\phi}$, $\hat{\theta}$, and the properties of the estimated model. The emphasis starting now is on $\hat{y}_{t+h|t}$ — the model’s predictions. These are related but different activities. A model can have imprecise parameter estimates yet produce useful forecasts (because the forecasts average over the parameter uncertainty). Conversely, a model with sharp parameter estimates can produce terrible forecasts (because the model is misspecified in a way that matters more at longer horizons). The diagnostic toolkit from Module 5 helps separate the two cases, but forecasting introduces its own quality measures — the loss functions we develop in Section 6.7.

6.3 Train/Test Splits for Time Series

Why Random Splits Are Wrong

In cross-sectional analysis, you can shuffle observations and split randomly — every row is exchangeable. In time series, shuffling breaks the temporal order. If you train on January through December and then test on June, you are using December data to predict June. That is forecasting the past.

The Back to the Future analogy. You cannot test your forecast on data you already saw — that is Biff with the sports almanac. The test set has to be data the model has never seen, and for time series that means data from the future relative to the training cutoff. A random split is like shuffling the pages of the almanac — you are still peeking, just in a less obvious way.

The Correct Approach: Temporal Splits

Split the series at a cutoff date. Everything before the cutoff is training; everything after is testing. The test set lives entirely in the future relative to the training set.

# Hold out the last 24 months for testing
sp <- ts_split(unrate, test_periods = 24)
str(sp)

List of 2
 $ train: Time-Series [1:696] from 1960 to 2018: 5.2 4.8 5.4 5.2 5.1 5.4 5.5 5.6 5.5 6.1 ...
 $ test : Time-Series [1:24] from 2018 to 2020: 4 4.1 4 4 3.8 4 3.8 3.8 3.7 3.8 ...

The ts_split() function from helpers/forecast_utils.R returns a list with $train and $test components, both ts objects that preserve the time index.

Two Evaluation Designs

Fixed-origin evaluation. Fit once on training, forecast the entire test horizon. Simple. One estimate, one set of forecasts. This is what we use in this module.

Rolling-window (expanding-window) evaluation. Re-fit the model as each new test observation arrives, forecast one step ahead, then expand the training set and repeat. More realistic. More expensive. Module 7 develops this.

For today we use the simplest version: one split, one fit, one forecast. Module 7 upgrades to rolling windows and formal comparison tests.

6.4 Recursive Forecast Computation

This section shows how to compute $\hat{y}_{t+h|t}$ by hand for AR and ARMA models. The key insight is simple: replace future innovations with zero, and replace future $y$’s with their forecasts.

The AR(1) Case

Start with the simplest possible model:

\[y_{t+1} = \alpha + \phi y_t + \epsilon_{t+1}\]

The 1-step forecast made at time $t$:

\[\hat{y}_{t+1|t} = \mathbb{E}[y_{t+1} \mid \mathcal{F}_t] = \alpha + \phi y_t + \underbrace{\mathbb{E}[\epsilon_{t+1} \mid \mathcal{F}_t]}_{= 0} = \alpha + \phi y_t\]

The key rule: replace future innovations with zero. $\epsilon_{t+1}$ has not happened yet; its conditional expectation is zero.

The 2-step forecast:

\[\hat{y}_{t+2|t} = \mathbb{E}[\alpha + \phi y_{t+1} + \epsilon_{t+2} \mid \mathcal{F}_t] = \alpha + \phi \hat{y}_{t+1|t}\]

The general $h$-step forecast:

\[\hat{y}_{t+h|t} = \alpha + \phi \hat{y}_{t+h-1|t}\]

Each step feeds the previous forecast forward. This is why it is called recursive. The recursion replaces unknown future $y$’s with their forecasts, and unknown future shocks with zero.

Convergence. As $h \rightarrow \infty$, the recursive forecast converges to the unconditional mean $\mu = \alpha / (1 - \phi)$ for any stationary AR(1) ($|\phi| < 1$). The forecast “forgets” its starting point — stationarity means mean-reversion, and the forecast reverts to the long-run mean.

The AR(2) Case

\[y_{t+1} = \alpha + \phi_1 y_t + \phi_2 y_{t-1} + \epsilon_{t+1}\]

1-step: $\hat{y}_{t+1|t} = \alpha + \phi_1 y_t + \phi_2 y_{t-1}$. Both $y_t$ and $y_{t-1}$ are observed.

2-step: $\hat{y}_{t+2|t} = \alpha + \phi_1 \hat{y}_{t+1|t} + \phi_2 y_t$. The first lag uses the forecast; the second lag is still observed.

3-step: $\hat{y}_{t+3|t} = \alpha + \phi_1 \hat{y}_{t+2|t} + \phi_2 \hat{y}_{t+1|t}$. Now both lags are forecasts.

At horizons beyond $p$, every lagged $y$ in the equation is itself a forecast. The recursion feeds on its own output. That is where the uncertainty accumulates.

The General ARMA Case

For ARMA$(p,q)$:

\[y_{t+h} = \alpha + \sum_{j=1}^{p} \phi_j y_{t+h-j} + \sum_{l=1}^{q} \theta_l \epsilon_{t+h-l} + \epsilon_{t+h}\]

The forecast rule applies two replacements:

\[\hat{y}_{t+h|t} = \alpha + \sum_{j=1}^{p} \phi_j \tilde{y}_{t+h-j} + \sum_{l=1}^{q} \theta_l \tilde{\epsilon}_{t+h-l}\]

where:

\[\tilde{y}_{t+h-j} = \begin{cases} \hat{y}_{t+h-j|t} & \text{if } t+h-j > t \text{ (future: use forecast)} \\ y_{t+h-j} & \text{if } t+h-j \leq t \text{ (past: use observed)} \end{cases}\]

\[\tilde{\epsilon}_{t+h-l} = \begin{cases} 0 & \text{if } t+h-l > t \text{ (future: use zero)} \\ e_{t+h-l} & \text{if } t+h-l \leq t \text{ (past: use residual)} \end{cases}\]

Two replacement rules. Future $y$’s get replaced by forecasts. Future $\epsilon$’s get replaced by zero. Past $y$’s are observed data. Past $\epsilon$’s are the residuals $e_t$ from the fitted model.

The MA complication. For an MA or ARMA model, the first few forecast steps can use observed residuals (the ones computed from the fitted model). At horizon $h > q$, all the MA terms have been replaced by zero, and the forecast collapses to a pure AR recursion in the forecasts. This is why MA models produce “flatter” long-horizon forecasts — the MA terms die out quickly.

Worked Example: ARMA(1,1) Forecast by Hand

The replacement rules are clearest when you see them on a model that has both AR and MA terms. Consider an ARMA(1,1):

\[y_{t+1} = \alpha + \phi y_t + \theta \epsilon_t + \epsilon_{t+1}\]

1-step forecast ($h = 1$):

\[\hat{y}_{t+1|t} = \alpha + \phi y_t + \theta e_t\]

Here $y_t$ is observed (past), and $e_t$ is the residual from the fitted model (past) — so we use it. The innovation $\epsilon_{t+1}$ is future, so it becomes zero. This is the one horizon where the MA term contributes: we have a concrete residual $e_t$ to plug in.

2-step forecast ($h = 2$):

\[\hat{y}_{t+2|t} = \alpha + \phi \hat{y}_{t+1|t} + \theta \underbrace{\tilde{\epsilon}_{t+1}}_{= 0}\]

Now $y_{t+1}$ is future, so we replace it with $\hat{y}_{t+1|t}$. And $\epsilon_{t+1}$ is also future — it gets replaced with zero. The MA term vanishes. From $h = 2$ onward, the ARMA(1,1) forecast is a pure AR(1) recursion in the forecasts:

\[\hat{y}_{t+h|t} = \alpha + \phi \hat{y}_{t+h-1|t}, \qquad h \geq 2\]

This is the general pattern for any ARMA($p, q$): the MA terms contribute observed residuals for the first $q$ steps, then drop out. After that, the forecast is driven entirely by the AR structure. The practical consequence is that the MA term improves short-horizon forecasts (where the recent residual carries information) but has no effect on long-horizon forecasts. If two models differ only in their MA specification, their long-horizon forecasts converge.

R’s `intercept` Is a Mean, Not an Intercept

Before you check the recursion against R, one trap to disarm. For a stationary model, Arima() reports a coefficient named intercept — but it is the estimated process mean $\hat{\mu}$, not the master-equation intercept $\hat{\alpha}$. The two are related through the AR coefficients:

\[\alpha = \mu \left(1 - \sum_{j=1}^{p} \phi_j\right)\]

If you plug $\hat{\mu}$ into the recursion where $\hat{\alpha}$ belongs, your hand forecasts will disagree with forecast() — subtly for a near-zero-mean series, badly for anything with a mean far from zero. Recover $\hat{\alpha} = \hat{\mu}(1 - \sum_j \hat{\phi}_j)$ first, as the code below does. (Equivalently, you can run the whole recursion in deviations from $\hat{\mu}$.) The name is a long-standing misnomer inherited from stats::arima(), and it bites everyone exactly once.

Simulation: Verify the Recursion on AR(2)

Simulate an AR(2) with known parameters. Compute the recursion by hand and compare to forecast():

library(forecast)
set.seed(1985)

# Simulate AR(2): phi1 = 0.6, phi2 = -0.2
y_sim <- arima.sim(model = list(ar = c(0.6, -0.2)), n = 200)
y_sim <- ts(y_sim)

# Fit the correct model (we know the truth)
fit_sim <- Arima(y_sim, order = c(2, 0, 0))
coef(fit_sim)

         ar1          ar2    intercept 
 0.473493444 -0.155019783  0.007768097

# Hand recursion for h = 1, 2
mu_hat    <- coef(fit_sim)["intercept"]          # this is mu-hat, NOT alpha-hat
phi1_hat  <- coef(fit_sim)["ar1"]
phi2_hat  <- coef(fit_sim)["ar2"]
alpha_hat <- mu_hat * (1 - phi1_hat - phi2_hat)  # recover alpha-hat (Technical Note)

T_end <- length(y_sim)
y_T   <- y_sim[T_end]
y_Tm1 <- y_sim[T_end - 1]

# 1-step
f1 <- alpha_hat + phi1_hat * y_T + phi2_hat * y_Tm1
# 2-step
f2 <- alpha_hat + phi1_hat * f1 + phi2_hat * y_T

# Compare to forecast()
fc <- forecast(fit_sim, h = 2)
data.frame(
  horizon   = 1:2,
  hand      = c(f1, f2),
  forecast_ = as.numeric(fc$mean)
)

  horizon       hand  forecast_
1       1 0.47969656 0.47969656
2       2 0.03357862 0.03357862

The numbers match — to machine precision, now that $\hat{\alpha}$ has been recovered correctly. forecast() is doing exactly this recursion under the hood, plus computing the prediction intervals we derive in Section 6.5. No magic — just the recursion we wrote on the board, applied at scale.

Simulation: Verify the ARMA(1,1) Recursion

Now verify the ARMA(1,1) case, where the MA term matters at $h = 1$ and disappears at $h = 2$:

set.seed(2015)

# Simulate ARMA(1,1): phi = 0.6, theta = 0.4
y_arma <- arima.sim(model = list(ar = 0.6, ma = 0.4), n = 200)
y_arma <- ts(y_arma)

# Fit the correct model
fit_arma <- Arima(y_arma, order = c(1, 0, 1))
coef(fit_arma)

       ar1        ma1  intercept 
 0.6922884  0.3138359 -0.1976396

mu_hat    <- coef(fit_arma)["intercept"]
phi_hat   <- coef(fit_arma)["ar1"]
theta_hat <- coef(fit_arma)["ma1"]
alpha_hat <- mu_hat * (1 - phi_hat)   # recover alpha-hat again

T_end <- length(y_arma)
y_T   <- y_arma[T_end]
e_T   <- residuals(fit_arma)[T_end]   # the last in-sample residual

# 1-step: MA term uses the observed residual
f1 <- alpha_hat + phi_hat * y_T + theta_hat * e_T

# 2-step: MA term is zero (epsilon_{T+1} is future)
f2 <- alpha_hat + phi_hat * f1

# 3-step: pure AR recursion from here on
f3 <- alpha_hat + phi_hat * f2

# Compare to forecast()
fc_arma <- forecast(fit_arma, h = 3)
data.frame(
  horizon   = 1:3,
  hand      = c(f1, f2, f3),
  forecast_ = as.numeric(fc_arma$mean)
)

  horizon       hand  forecast_
1       1 -1.8692504 -1.8692504
2       2 -1.3548763 -1.3548763
3       3 -0.9987811 -0.9987811

The numbers match again. Notice the key difference from the AR(2) example: the $h = 1$ forecast uses the residual $e_T$, which is information the AR-only forecast would throw away. At $h = 2$ and beyond, the MA term contributes nothing — the forecast is just $\hat{\alpha} + \hat{\phi} \cdot \hat{y}_{t+h-1|t}$. This is why the MA component improves short-horizon accuracy but has no effect on the long-run forecast.

Real Data: The First UNRATE Forecast

Produce the first forecast from the L5 SARIMA on the training set:

sp <- ts_split(unrate, test_periods = 24)

# Re-fit the L5 winner on the training set only
fit_train <- Arima(sp$train, order = c(1, 1, 2),
                   seasonal = list(order = c(1, 0, 1), period = 12))

# Forecast 24 months ahead
fc_unrate <- forecast(fit_train, h = 24)

# First look: point forecasts and prediction intervals
autoplot(fc_unrate) +
  autolayer(sp$test, series = "Actual") +
  ggtitle("UNRATE: 24-month forecast from the L5 SARIMA") +
  theme_bw()

There it is — your first real-data ARMA forecast. The blue line is the point forecast. The shaded bands are the prediction intervals. We will understand those bands in Section 6.5 and grade this forecast in Section 6.8. For now, notice two things: the point forecast hovers near the last observed level while the actuals drift below it — yet stay inside the bands — and the bands get wider as the horizon grows. The widening is a feature, not a bug; whether “inside the bands but persistently below the point forecast” is good enough is exactly what Section 6.8’s scorekeeping is for.

6.5 Forecast Uncertainty and Prediction Intervals

Why Does Uncertainty Grow?

The weather forecast analogy. Your weather app gives a tighter range for tomorrow than for next Thursday. Why? Because between now and next Thursday, more unknown things will happen — storms, pressure changes, fronts that haven’t formed yet. Each day adds another layer of unresolved randomness. Time-series forecasting is exactly the same: each future period adds another $\epsilon$ we don’t know.

The forecast error at horizon $h$ is:

\[e_{t+h|t} = y_{t+h} - \hat{y}_{t+h|t}\]

This error is the sum of all the innovations we could not predict: $\epsilon_{t+1}, \epsilon_{t+2}, \ldots, \epsilon_{t+h}$, weighted by the impulse response of the model (the $\psi$-weights).

The $\psi$-Weight Representation

By the Wold decomposition (Module 3 callback), any stationary ARMA can be written as:

\[y_t = \mu + \sum_{i=0}^{\infty} \psi_i \epsilon_{t-i}, \qquad \psi_0 = 1\]

The $\psi_i$ are the MA($\infty$) coefficients from the Wold decomposition. They measure the impulse response of $y_{t+i}$ to a unit shock $\epsilon_t$, and they are determined entirely by the AR and MA polynomials of the model.

In general, the $\psi$-weights come from expanding the ratio of the MA and AR lag polynomials as a power series:

\[\Psi(L) = \frac{\Theta(L)}{\Phi(L)} = \sum_{i=0}^{\infty} \psi_i L^i\]

where $\Phi(L) = 1 - \phi_1 L - \cdots - \phi_p L^p$ is the AR polynomial and $\Theta(L) = 1 + \theta_1 L + \cdots + \theta_q L^q$ is the MA polynomial (course sign conventions).

AR(1): $\Phi(L) = 1 - \phi L$, $\Theta(L) = 1$. Then $\Psi(L) = 1 / (1 - \phi L) = 1 + \phi L + \phi^2 L^2 + \cdots$, giving $\psi_i = \phi^i$. Geometric decay — the same geometry as the AR(1) ACF from Module 1.

For anything bigger than an AR(1), the expansion is coefficient-matching algebra (the Deeper Dive below works one case by hand), and R does it for you: ARMAtoMA(ar, ma, lag.max) from the stats package returns the $\psi$-weights directly.

Deeper Dive — $\psi$-Weights for the ARMA(1,1) by Hand

Solve for $\psi_1, \psi_2, \ldots$ by matching coefficients in the identity $\Phi(L) \Psi(L) = \Theta(L)$ — or, for a model this small, expand the ratio directly. With $\Phi(L) = 1 - \phi L$ and $\Theta(L) = 1 + \theta L$:

\[\Psi(L) = \frac{1 + \theta L}{1 - \phi L} = (1 + \theta L)(1 + \phi L + \phi^2 L^2 + \cdots)\]

Multiplying out and collecting powers of $L$:

$\psi_0 = 1$
$\psi_1 = \phi + \theta$
$\psi_i = (\phi + \theta) \phi^{i-1}$ for $i \geq 1$

So the first $\psi$-weight incorporates both the AR persistence and the MA effect, and subsequent weights decay geometrically at rate $\phi$. If $\phi = 0.6$ and $\theta = 0.4$, then $\psi_1 = 1.0$, $\psi_2 = 0.6$, $\psi_3 = 0.36$, and so on. Verify against R:

# Psi-weights for ARMA(1,1) with phi = 0.6, theta = 0.4
ARMAtoMA(ar = 0.6, ma = 0.4, lag.max = 5)

[1] 1.0000 0.6000 0.3600 0.2160 0.1296

# The formula (phi + theta) * phi^(i-1) gives: 1.0, 0.6, 0.36, 0.216, 0.1296

For higher-order models the coefficient-matching gets tedious by hand, but the principle is identical, and ARMAtoMA() scales to any $(p, q)$.

Forecast-Error Variance

The $h$-step forecast error can be written as:

\[e_{t+h|t} = \sum_{i=0}^{h-1} \psi_i \epsilon_{t+h-i}\]

since the innovations $\epsilon_{t+1}, \ldots, \epsilon_{t+h}$ are the ones the forecast could not anticipate. The forecast-error variance is therefore:

\[\text{Var}(e_{t+h|t}) = \sigma^2 \sum_{i=0}^{h-1} \psi_i^2\]

Read that formula carefully:

At $h = 1$: $\text{Var}(e_{t+1|t}) = \sigma^2 \psi_0^2 = \sigma^2$. One-step-ahead uncertainty is just the innovation variance — the $\sigma^2$ we have been tracking since Module 1.
At $h = 2$: $\text{Var}(e_{t+2|t}) = \sigma^2(\psi_0^2 + \psi_1^2) = \sigma^2(1 + \psi_1^2)$. Bigger.
As $h \rightarrow \infty$ for a stationary model: the sum converges and the forecast-error variance converges to the unconditional variance $\gamma_0$. The forecast becomes as uncertain as knowing nothing — you are back to the marginal distribution. (For a non-stationary model like our $d = 1$ SARIMA on UNRATE, the sum does not converge — the bands keep growing without bound, which is exactly what the fan chart in Section 6.8 shows.)

The AR(1) Closed Form

For an AR(1), $\psi_i = \phi^i$, so:

\[\text{Var}(e_{t+h|t}) = \sigma^2 \sum_{i=0}^{h-1} \phi^{2i} = \sigma^2 \cdot \frac{1 - \phi^{2h}}{1 - \phi^2}\]

As $h \rightarrow \infty$: $\text{Var}(e_{t+h|t}) \rightarrow \sigma^2 / (1 - \phi^2) = \gamma_0$. The forecast-error variance converges to the unconditional variance — exactly as claimed.

Gaussian Prediction Intervals

Under the assumption $\epsilon_t \sim N(0, \sigma^2)$, the forecast error $e_{t+h|t}$ is a linear combination of Gaussian innovations, hence Gaussian itself:

\[e_{t+h|t} \sim N\left(0, \; \sigma^2 \sum_{i=0}^{h-1} \psi_i^2\right)\]

A $(1 - \alpha) \times 100\%$ prediction interval is:

\[\hat{y}_{t+h|t} \pm z_{\alpha/2} \cdot \sigma \sqrt{\sum_{i=0}^{h-1} \psi_i^2}\]

For the standard 95% PI: $z_{0.025} = 1.96$.

This is what forecast() computes. The point forecast is the conditional mean. The interval is the point forecast plus or minus 1.96 times the square root of the forecast-error variance. Every fan chart you have ever seen is this formula, applied at each horizon and drawn as a band.

Reading the Fan Chart

In the UNRATE forecast from Section 6.4:

Center line: the point forecast $\hat{y}_{t+h|t}$.
Darkest band (80% PI): $\hat{y}_{t+h|t} \pm 1.28 \sigma \sqrt{\sum \psi_i^2}$.
Lighter band (95% PI): $\hat{y}_{t+h|t} \pm 1.96 \sigma \sqrt{\sum \psi_i^2}$.
The widening: each band gets wider as $h$ increases because $\sum_{i=0}^{h-1} \psi_i^2$ grows.

The fan chart is not decorative. It is the single most important output of a forecasting exercise. If your report includes a point forecast without bands, you are hiding the uncertainty — and the uncertainty is the part your audience actually needs to see.

The Key Assumption

The 95% in a 95% PI rests on: $\epsilon_t \sim N(0, \sigma^2)$. If the innovations are not Gaussian — heavier tails, skewness, time-varying variance — the stated coverage is not what the interval actually delivers. The interval is nominal 95% but may under-cover or over-cover in practice — and, as the next section shows, it can do both at once, at different confidence levels. Section 6.6 tests this directly.

6.6 The Assumption-Break: Coverage Under Fat Tails

This is the module’s central demonstration — the direct analog of Module 2’s DF null-distribution simulation and Module 5’s seasonal misspecification experiment. Before you trust a tool on real data, watch how it behaves on data you generated — including data that violates its assumptions.

The Experiment: Gaussian Case

Simulate 1,000 ARMA(1,1) series with Gaussian innovations. For each, fit the correct model, produce 1-step-ahead prediction intervals at the 80%, 95%, and 99% levels, generate $y_{t+1}$ from the true DGP, and check whether it falls inside each interval:

set.seed(2015)
n_reps <- 1000
n_obs  <- 300
phi_true   <- 0.6
theta_true <- 0.4
levels_pct <- c(80, 95, 99)

covered_gaussian <- matrix(NA, n_reps, 3,
                           dimnames = list(NULL, paste0(levels_pct, "% PI")))

for (r in seq_len(n_reps)) {
  # Generate Gaussian innovations
  eps <- rnorm(n_obs + 1, 0, 1)
  y <- numeric(n_obs + 1)
  for (t in 2:(n_obs + 1)) {
    y[t] <- phi_true * y[t - 1] + theta_true * eps[t - 1] + eps[t]
  }

  y_train <- ts(y[1:n_obs])
  y_true  <- y[n_obs + 1]

  fit <- Arima(y_train, order = c(1, 0, 1), include.mean = FALSE)
  fc  <- forecast(fit, h = 1, level = levels_pct)

  covered_gaussian[r, ] <- (y_true >= as.numeric(fc$lower)) &
                           (y_true <= as.numeric(fc$upper))
}

colMeans(covered_gaussian)

80% PI 95% PI 99% PI 
 0.795  0.954  0.987

Roughly 80%, 95%, and 99% — each interval delivers close to its promise when the Gaussian assumption holds. (The slight shortfalls come from estimating the parameters and $\sigma^2$ rather than knowing them; more on that in the lessons below.) The tool works. Now let’s break the assumption.

The Assumption-Break: Fat-Tailed Innovations

Repeat the experiment with $t_3$ innovations (Student-$t$ with 3 degrees of freedom — finite variance but heavy tails). We rescale the $t_3$ draws to have the same variance as the Gaussian case, so the two experiments differ only in tail shape:

set.seed(2015)
covered_fat <- matrix(NA, n_reps, 3,
                      dimnames = list(NULL, paste0(levels_pct, "% PI")))

for (r in seq_len(n_reps)) {
  # Generate t(3) innovations, rescaled to unit variance.
  # Var(t_nu) = nu/(nu-2), so for nu=3, Var = 3.
  # Rescale by sqrt(3) so that Var = 1.
  # This ensures the two experiments differ ONLY in tail shape:
  # same mean, same variance, different kurtosis.
  eps <- rt(n_obs + 1, df = 3) / sqrt(3)
  y <- numeric(n_obs + 1)
  for (t in 2:(n_obs + 1)) {
    y[t] <- phi_true * y[t - 1] + theta_true * eps[t - 1] + eps[t]
  }

  y_train <- ts(y[1:n_obs])
  y_true  <- y[n_obs + 1]

  fit <- Arima(y_train, order = c(1, 0, 1), include.mean = FALSE)
  fc  <- forecast(fit, h = 1, level = levels_pct)

  covered_fat[r, ] <- (y_true >= as.numeric(fc$lower)) &
                      (y_true <= as.numeric(fc$upper))
}

rbind(gaussian = colMeans(covered_gaussian),
      fat_tail = colMeans(covered_fat))

         80% PI 95% PI 99% PI
gaussian  0.795  0.954  0.987
fat_tail  0.880  0.953  0.977

Reading the Table: The Break Is in the Tail, Not the Middle

If you expected every row of the fat-tail line to collapse, look again — the result is stranger and more instructive than a uniform failure.

The nominal 80% interval over-covers — it captures roughly 88% of outcomes. Too wide for its stated job.
The nominal 95% interval is almost exactly right — coverage near 95%, essentially indistinguishable from the Gaussian case.
The nominal 99% interval under-covers — it captures roughly 97.5–98%. That gap looks small until you translate it: the interval promises a miss rate of 1-in-100, and it delivers a miss rate of roughly 1-in-40. Extreme events happen more than twice as often as the interval claims.

Why this pattern? Both innovation distributions have variance 1 by construction, so the width of the fitted intervals is about the same in the two experiments. But a fat-tailed distribution with the same variance as a Gaussian must rearrange its probability mass: more mass piled near the center, and more mass pushed into the extreme tails — the “shoulders” in between are what gets thinned out. A Gaussian-width 80% interval reaches past the fat center, so it over-covers. A Gaussian-width 99% interval cannot reach the fat extremes, so it under-covers. The 95% interval happens to sit near the crossing point where the two distributions’ quantiles agree — a coincidence of where the curves cross, not a property you can rely on.

The practical punchline: fat tails do not hurt you in the middle of the distribution; they kill you in the far tail — exactly where risk management lives. A bank computing a 99% value-at-risk, a grid operator planning for the 1-in-100 demand spike, a fiscal agency stress-testing the extreme scenario: these users live in the tail the Gaussian formula gets most wrong.

Deeper Dive — The Quantile Crossing, Exactly

You can compute where the Gaussian interval succeeds and fails without any simulation. Under the DGP with rescaled $t_3$ innovations, and pretending parameters were known exactly, the 1-step Gaussian interval is $\pm z_{\alpha/2}$ (the innovation SD is 1), and its true coverage is the probability that a scaled $t_3$ lands inside:

nominal <- c(0.80, 0.95, 0.99)
z <- qnorm(1 - (1 - nominal) / 2)
true_coverage <- 1 - 2 * pt(-z * sqrt(3), df = 3)  # unscale: multiply back by sqrt(3)
data.frame(nominal, z, true_coverage = round(true_coverage, 4))

  nominal        z true_coverage
1    0.80 1.281552        0.8869
2    0.95 1.959964        0.9574
3    0.99 2.575829        0.9790

Theory says 88.7%, 95.7%, and 97.9% — matching the Monte Carlo. The Gaussian and unit-variance-$t_3$ quantile functions cross between the 95% and 99% levels: below the crossing the Gaussian interval is too wide, above it too narrow. The Monte Carlo numbers sit slightly below these idealized values because $\hat{\phi}$, $\hat{\theta}$, and $\hat{\sigma}^2$ are estimated, which adds a second, smaller layer of under-coverage that the formula ignores — visible even in the Gaussian row’s 99% entry.

Three Lessons from the Break

Prediction intervals are conditional on distributional assumptions — and the failure is not uniform. The stated coverage is a conditional statement: “correct if the innovations are Gaussian and the model is correctly specified and $\sigma^2$ is correctly estimated.” When the distributional leg breaks, the damage concentrates at the confidence levels furthest from where the quantiles happen to agree — in practice, the far tail.
Fat tails are the most common real-world breach. Financial returns, macro shocks, unemployment spikes — these series routinely exhibit heavier tails than the Gaussian. Matching the variance is not enough: the Gaussian formula gets the whole shape wrong, and the shape error is largest exactly at the extreme quantiles that risk decisions depend on. This is not a defect of the forecast package; it is an honest consequence of the Gaussian assumption built into the formula.
Bootstrap prediction intervals are the pragmatic fix. forecast(..., bootstrap = TRUE) replaces the Gaussian formula with a resampling approach that lets the data’s actual tail behavior drive the PI width. We will not develop bootstrap PIs in this module, but they are the professional response to the problem this simulation exposed.

Callback to Module 5: “Serial correlation in residuals is a property of your model relative to the data.” Here is the forecasting version: the shape of your prediction interval is a property of your distributional assumption relative to the data. If the assumption is wrong, the interval is wrong — and it is most wrong where you can least afford it. The crystal ball has error bars, but the error bars are only as honest as the assumption behind them.

Deeper Dive — Near-Unit-Root Fan Charts

For a quick visual reinforcing the widening-fan story, forecast an AR(1) with $\phi = 0.99$ (near unit root) 60 steps ahead. The fan opens to enormous width because the $\psi$-weights $\phi^i$ decay very slowly:

set.seed(1985)
y_near_ur <- arima.sim(model = list(ar = 0.99), n = 200)
fit_near  <- Arima(ts(y_near_ur), order = c(1, 0, 0))
fc_near   <- forecast(fit_near, h = 60)
autoplot(fc_near) +
  ggtitle("Fan chart: AR(1) with phi = 0.99 — near unit root") +
  theme_bw()

The fan opens like a trumpet. Near the unit root, the $\psi$-weights barely decay, so the forecast-error variance piles up fast. At 60 steps ahead, the PI is so wide it is almost useless. This is why long-horizon forecasts from persistent series are honest about being uninformative — and why you should be suspicious of anyone who gives a confident 5-year unemployment forecast from a univariate model.

6.7 Loss Functions

The Concept: What Does It Cost You to Be Wrong?

A loss function measures the cost of being wrong with your forecast. That cost is not abstract — it depends on the domain, the stakes, and the shape of the consequences. Choosing a loss function is a decision about what kind of errors hurt most, and you make that decision before you see the results, not after.

We use $g(\cdot)$ for loss functions, not $L(\cdot)$, because $L$ is the lag operator in this course. The helper file forecast_utils.R provides mse(), mae(), huber(), mean_huber(), mape(), and asym_loss() so you can compute any of the standard losses without a package dependency.

The Texas grid example. Suppose you are in charge of forecasting energy demand in Texas. A cold snap hits. The electrical grid is strained. If you under-forecast demand — not enough supply — people freeze to death. If you over-forecast demand — too much supply — you wasted money running extra generators. Both are errors, but they are not the same error. The cost of being wrong in one direction (under-forecast, people die) is catastrophically larger than the cost of being wrong in the other direction (over-forecast, wasted fuel). Your loss function must reflect that asymmetry. If you evaluate your forecast with a symmetric loss like MSE, you are implicitly saying that over-forecasting by 500 MW and under-forecasting by 500 MW are equally bad. They are not.

This is the core idea: the loss function encodes what it costs you to be wrong. Different cost structures demand different loss functions. The choice comes first; the formula follows.

Key principle: Like model selection criteria, you cannot compare between loss functions, only within. Two models must be evaluated under the same loss function, with the same dependent variable and the same number of test-set observations. “Model A wins on MSE but Model B wins on MAE” is not a paradox — it is two different questions with two different answers.

How the Cost Structure Maps to the Loss Function

Think about the per-unit cost of being wrong — how much damage does each additional unit of forecast error cause?

The cost of being wrong is increasing explosively. Each additional unit of error does more damage than the last. A forecast that is off by 10 is not just twice as bad as one that is off by 5 — it is four times as bad. One catastrophic miss can dominate everything else. This is the world of MSE: it squares errors, so large errors are penalized quadratically. MSE is the right choice when tail risk is what kills you.
The cost of being wrong is constant per unit. An error of 10 is exactly twice as bad as an error of 5. Every unit of miss costs the same, whether it is the first or the hundredth. This is the world of MAE: it takes absolute values, weighting all errors linearly. MAE is the right choice when errors accumulate proportionally — inventory management, routine staffing, contexts where there is no cliff.
The cost is explosive for small errors but flattens for large ones. You care a lot about getting the forecast right within a tight band, but beyond a certain threshold the damage is already done and additional error doesn’t make it much worse. This is the world of Huber loss: quadratic near zero (like MSE), linear in the tails (like MAE). Huber is the compromise — strict on small mistakes, forgiving on outliers.
The cost is asymmetric. Being wrong in one direction costs more than being wrong in the other. The Texas grid. A central bank under-forecasting inflation. A hospital under-forecasting ER admissions. This is the world of asymmetric loss, which lets you put different weights on over- vs. under-prediction.

The Formulas

In all formulas below, $N$ denotes the number of test-set observations, not the full sample size $T$.

Mean Squared Error (MSE) — explosive per-unit cost:

\[\text{MSE} = \frac{1}{N} \sum_{t=1}^{N} (y_t - \hat{y}_t)^2\]

Optimal forecast under MSE: the conditional mean $\mathbb{E}[y_{t+h} \mid \mathcal{F}_t]$.
Units: squared units of $y$. Often reported as RMSE $= \sqrt{\text{MSE}}$ for interpretability.
R: mse(actual, forecast) and rmse(actual, forecast) from forecast_utils.R.

Mean Absolute Error (MAE) — constant per-unit cost:

\[\text{MAE} = \frac{1}{N} \sum_{t=1}^{N} |y_t - \hat{y}_t|\]

Optimal forecast under MAE: the conditional median of $y_{t+h} \mid \mathcal{F}_t$.
Units: same units as $y$. More robust to outliers than MSE.
R: mae(actual, forecast) from forecast_utils.R.

Mean Absolute Percentage Error (MAPE) — scale-free:

\[\text{MAPE} = \frac{1}{N} \sum_{t=1}^{N} \left| \frac{y_t - \hat{y}_t}{y_t} \right| \times 100\%\]

Scale-free: useful for comparing across series with different units or magnitudes.
Undefined when $y_t = 0$. Asymmetric even when you don’t want it to be (under-predictions are bounded at 100%; over-predictions are unbounded).
Rarely the right choice for macro series; popular in business forecasting.
R: mape(actual, forecast) from forecast_utils.R.

Huber Loss — explosive near zero, constant in the tails:

\[g_\delta(e) = \begin{cases} \frac{1}{2} e^2 & \text{if } |e| \leq \delta \\ \delta \left(|e| - \frac{1}{2}\delta\right) & \text{if } |e| > \delta \end{cases}\]

where $\delta > 0$ is a threshold parameter that controls where the transition from quadratic to linear happens. Smaller $\delta$ makes Huber behave more like MAE; larger $\delta$ makes it behave more like MSE.

R: huber(actual, forecast, delta) returns element-wise losses; mean_huber(actual, forecast, delta) returns the scalar mean.

Asymmetric Loss — directional cost:

\[g(e) = \begin{cases} \alpha \cdot |e| & \text{if } e > 0 \text{ (under-prediction)} \\ (1 - \alpha) \cdot |e| & \text{if } e \leq 0 \text{ (over-prediction)} \end{cases}\]

where $\alpha \in (0, 1)$. Setting $\alpha > 0.5$ penalizes under-prediction more heavily; $\alpha = 0.5$ recovers MAE.

Optimal forecast under asymmetric linear loss: the $\alpha$-quantile of the conditional distribution. This is why the Texas grid operator, if using $\alpha = 0.9$, would forecast the 90th percentile of demand rather than the mean — deliberately over-forecasting to build in a safety margin.
R: asym_loss(actual, forecast, alpha) from forecast_utils.R.

Looking Ahead — Module 7: MASE, the Scale-Free Benchmark Ratio

Module 7 adds one more metric to this menu, and it is worth previewing because it formalizes something we do qualitatively in Section 6.8. The Mean Absolute Scaled Error:

\[\text{MASE} = \frac{\frac{1}{N} \sum_{t=1}^{N} |y_t - \hat{y}_t|}{\frac{1}{T-1} \sum_{t=2}^{T} |y_t - y_{t-1}|}\]

The numerator is the MAE of your forecast on the test set. The denominator is the MAE of the naive random-walk forecast (last-value-carried-forward) computed on the training set. MASE divides your model’s accuracy by the naive benchmark’s accuracy, producing a scale-free ratio:

MASE $< 1$: your model beats the naive benchmark.
MASE $= 1$: your model is exactly as good as repeating the last value.
MASE $> 1$: the naive forecast is better — your model is not earning its complexity.

For seasonal data, the denominator uses the seasonal naive forecast $|y_t - y_{t-s}|$ instead. This is Hyndman and Koehler’s (2006) recommended scale-free metric. Unlike MAPE, MASE is defined when $y_t = 0$, is symmetric, and has a natural interpretation as “how many times better (or worse) than the naive baseline?” Module 7 develops it alongside Theil’s U and the Diebold-Mariano test.

Graphical Comparison

error <- seq(-10, 10, by = 0.01)
loss_mse   <- error^2
loss_mae   <- abs(error)
delta <- 4
loss_huber <- ifelse(abs(error) <= delta,
                     0.5 * error^2,
                     delta * (abs(error) - 0.5 * delta))

ggplot() +
  geom_line(aes(x = error, y = loss_mse, color = "MSE (quadratic)")) +
  geom_line(aes(x = error, y = loss_mae, color = "MAE (linear)")) +
  geom_line(aes(x = error, y = loss_huber, color = "Huber (hybrid)")) +
  scale_color_manual(name = "",
    values = c("MSE (quadratic)" = "blue",
               "MAE (linear)" = "red",
               "Huber (hybrid)" = "purple")) +
  labs(x = "Forecast error", y = "Loss", title = "Loss functions compared") +
  theme_bw() + theme(legend.position = "bottom")

The quadratic (MSE) explodes for large errors — it is the cost curve where tail events are catastrophic. The linear (MAE) grows steadily — every unit of error costs the same. The Huber tracks MSE near zero and flattens to MAE in the tails — the right choice when you care about precision in the normal range but don’t want a single outlier to dominate the evaluation.

Simulation: The Loss Function Picks the Target

The choice of loss function is not just a philosophical preference — it can change which forecast you declare the winner. Section 6.2 made the claim in theory: under MSE the optimal forecast is the conditional mean; under MAE it is the conditional median. For symmetric innovations the two targets coincide and the distinction is invisible. Make the innovations skewed and the targets separate — and then the loss function decides the contest.

Generate an AR(1) whose innovations are centered exponential draws — mean zero, but skewed hard to the right: most shocks are small and negative, a few are large and positive (think of a series that mostly drifts down quietly and occasionally jumps). Fit the correct AR(1) on a training set, then produce two competing one-step-ahead forecast tracks over the test set:

Forecast A (mean-targeting): the model’s conditional mean, straight from the fit.
Forecast B (median-targeting): the same conditional mean, shifted by the median of the training residuals — an estimate of the innovation median, which for a right-skewed distribution sits below the mean.

set.seed(2020)
n_total <- 600
n_train <- 350

# DGP: AR(1) with centered-exponential innovations (mean 0, skewed right)
eps <- rexp(n_total) - 1
y <- numeric(n_total)
for (t in 2:n_total) {
  y[t] <- 0.6 * y[t - 1] + eps[t]
}
y_ts    <- ts(y)
y_train <- window(y_ts, end = n_train)
y_test  <- window(y_ts, start = n_train + 1)

# Fit the correct model on training data only
fit_skew <- Arima(y_train, order = c(1, 0, 0))

# One-step-ahead forecasts across the test set WITHOUT re-estimating:
# Arima(full_series, model = fit) re-applies the training coefficients,
# and fitted() then gives the one-step forecast at every t using only
# information available at t-1. A handy trick worth memorizing.
refit    <- Arima(y_ts, model = fit_skew)
fc_mean  <- window(fitted(refit), start = n_train + 1)

# Median-targeting forecast: shift by the training-residual median
med_shift <- median(residuals(fit_skew))
fc_median <- fc_mean + med_shift
med_shift   # negative: the median of a right-skewed innovation sits below its mean

[1] -0.3911346

actuals <- as.numeric(y_test)

# Score both forecast tracks under both loss functions
rbind(
  MSE = c(mean_fc = mse(actuals, as.numeric(fc_mean)),
          median_fc = mse(actuals, as.numeric(fc_median))),
  MAE = c(mean_fc = mae(actuals, as.numeric(fc_mean)),
          median_fc = mae(actuals, as.numeric(fc_median)))
)

      mean_fc median_fc
MSE 1.0632934 1.2070380
MAE 0.7531036 0.7274896

Read the table by rows. Under MSE, the mean-targeting forecast wins. Under MAE, the median-targeting forecast wins. Same data, same fitted model, two forecast tracks — and each one wins under the loss it was built for. The loss function changed the verdict.

This is not a paradox. It is the loss function doing its job. MSE asked “which forecast is closest in squared distance, where the occasional big positive shock is very expensive?” — and the mean, which leans toward those expensive shocks, wins. MAE asked “which forecast is closest month in, month out, one absolute unit at a time?” — and the median, which sits where most of the observations actually land, wins. Those are different questions with different answers. The right question depends on what it costs you to be wrong — and that is a decision you make before you see the data.

You Can Write Your Own

The five loss functions above are standard — you will find them in every forecasting textbook. But there is nothing sacred about them. They are common because they are mathematically convenient and cover a range of common cost structures. In practice, the cost of being wrong is almost always context-dependent, and context rarely fits neatly into one of five boxes.

If your domain has a specific cost structure, write a loss function that reflects it. A loss function is just a function that takes an actual value and a forecast and returns a non-negative number. The only requirements are: (1) it equals zero when the forecast is perfect, and (2) it increases as the forecast gets worse. Beyond that, the shape is yours to choose.

For example, suppose you are forecasting staffing needs for a call center. Under-staffing by 1 person costs you $500 in lost revenue from dropped calls. Over-staffing by 1 person costs you $100 in idle wages. Under-staffing by 5 costs $5,000 (non-linear — angry customers leave permanently). You could write:

call_center_loss <- function(actual, forecast) {
  e <- actual - forecast
  cost <- ifelse(e > 0,
                 500 * e + 50 * e^2,   # under-staffed: linear + quadratic escalation
                 100 * abs(e))          # over-staffed: linear only
  return(mean(cost))
}

This loss function is not MSE, MAE, Huber, or asymmetric linear — it is a custom function that encodes your cost structure. You evaluate your competing models with call_center_loss(actuals, fc_A$mean) and call_center_loss(actuals, fc_B$mean), and the winner is the model that minimizes your actual cost, not a generic textbook metric.

The standard losses are starting points. If you have domain knowledge about what errors cost, use it. A custom loss function that reflects the real cost structure will almost always be a better evaluation tool than a generic one that does not.

6.8 The Full Workflow on UNRATE

This section executes the complete forecasting workflow end-to-end: split, fit, forecast, overlay actuals, compute loss, and compare against a naive baseline.

Step 1: Split

# Hold out the last 24 months
sp <- ts_split(unrate, test_periods = 24)
cat("Training:", start(sp$train), "to", end(sp$train), "\n")

Training: 1960 1 to 2017 12

cat("Testing: ", start(sp$test),  "to", end(sp$test),  "\n")

Testing:  2018 1 to 2019 12

Step 2: Fit on Training Only

# The L5 winner, re-fit on the training set
fit_train <- Arima(sp$train, order = c(1, 1, 2),
                   seasonal = list(order = c(1, 0, 1), period = 12))
summary(fit_train)

Series: sp$train 
ARIMA(1,1,2)(1,0,1)[12] 

Coefficients:
         ar1      ma1     ma2    sar1     sma1
      0.9061  -0.9295  0.2095  0.5283  -0.7921
s.e.  0.0272   0.0451  0.0401  0.0735   0.0524

sigma^2 = 0.02531:  log likelihood = 292.43
AIC=-572.86   AICc=-572.73   BIC=-545.59

Training set error measures:
                        ME      RMSE       MAE         MPE    MAPE      MASE
Training set -0.0003949227 0.1584036 0.1215638 -0.03216234 2.06971 0.1555944
                   ACF1
Training set 0.01084174

checkresiduals(fit_train)


    Ljung-Box test

data:  Residuals from ARIMA(1,1,2)(1,0,1)[12]
Q* = 22.046, df = 19, p-value = 0.282

Model df: 5.   Total lags used: 24

We re-fit on training only. The test set has not touched the estimation. The diagnostics still pass — if they did not, we would go back to Module 5 and reconsider the model.

Step 3: Forecast into the Test Period

fc <- forecast(fit_train, h = 24)

# Overlay actuals
autoplot(fc) +
  autolayer(sp$test, series = "Actual", color = "red") +
  ggtitle("UNRATE forecast: ARIMA(1,1,2)(1,0,1)[12], 24-month horizon") +
  ylab("UNRATE (%)") +
  theme_bw()

Walk through the plot:

The blue line (point forecast) stays close to the last observed level, drifting slightly upward.
The red line (actuals) glides slowly downward — 2018–2019 was a long, calm decline in unemployment.
At short horizons (1–3 months), the actuals are close to the point forecast.
At longer horizons, the actuals stray further below the point forecast — but remain comfortably inside the 95% band the whole way.

Step 4: Compute Loss and Compare to Baseline

actuals  <- as.numeric(sp$test)
fc_point <- as.numeric(fc$mean)

# --- Naive random-walk baseline ---
# The random walk says: best guess for all future periods is the last
# training observation. This is the simplest possible benchmark.
naive_fc <- rep(tail(sp$train, 1), length(sp$test))

cat("=== SARIMA (1,1,2)(1,0,1)[12] ===\n")

=== SARIMA (1,1,2)(1,0,1)[12] ===

cat("  MSE: ", mse(actuals, fc_point), "\n")

  MSE:  0.1642191

cat("  MAE: ", mae(actuals, fc_point), "\n")

  MAE:  0.3247277

cat("  RMSE:", rmse(actuals, fc_point), "\n\n")

  RMSE: 0.4052396

cat("=== Naive random walk ===\n")

=== Naive random walk ===

cat("  MSE: ", mse(actuals, naive_fc), "\n")

  MSE:  0.1275

cat("  MAE: ", mae(actuals, naive_fc), "\n")

  MAE:  0.3166667

cat("  RMSE:", rmse(actuals, naive_fc), "\n")

  RMSE: 0.3570714

Look carefully, because real data just handed us a humbling verdict: the naive random walk beats the diagnostic-clean SARIMA on every metric — narrowly on MAE, more visibly on MSE and RMSE. The model that survived Module 5’s entire diagnostic gauntlet loses to “carry forward the last number.”

Three honest readings, all of which belong in your toolkit:

The naive benchmark is genuinely hard to beat on a persistent series in a calm period. UNRATE is one of the most persistent series in macroeconomics, and 2018–2019 was a placid, nearly straight glide from 4.1% down to 3.5%. Last-value-carried-forward is a superb forecast for exactly that environment. This is why the naive random walk is the universal benchmark — not because it is a strawman, but because it usually isn’t.
Diagnostics and forecast accuracy answer different questions. Module 5 certified that the SARIMA describes the in-sample dependence structure adequately. That is not the same claim as “it forecasts this particular 24-month window better than a simpler rule.” A model can be adequate as a description and still add little forecasting value over a stretch with no turning points and no seasonal surprises to exploit — the two situations where its extra machinery would earn its keep.
One window is one draw. The margins here are small — a few hundredths on MAE — and they come from a single 24-month evaluation window. Is that gap statistically meaningful, or is it noise? You cannot tell by staring at two numbers. That question — “is this forecast significantly better than that one?” — is exactly what Module 7’s Diebold-Mariano test answers.

What the SARIMA still delivers that the naive forecast cannot: honest uncertainty. The naive rule emits a number with no bands; the SARIMA’s prediction intervals contained every one of the 24 test observations. If your deliverable is a decision under uncertainty rather than a single number in a spreadsheet, that is not a consolation prize — it is the product.

The Full Fan Chart

fc_fan <- forecast(fit_train, h = 24, fan = TRUE)
autoplot(fc_fan) +
  autolayer(sp$test, series = "Actual") +
  ggtitle("Fan chart: UNRATE") +
  theme_bw()

The darker the band, the more confident the interval. The outermost whisper of shading is the 99% PI. If the actual wanders outside even that, you have a genuine surprise — though remember Section 6.6: for real series with fat-tailed shocks, that outermost band is exactly where the Gaussian formula is least trustworthy.

Common Pitfalls and Misconceptions

“The point forecast is the most likely outcome.” Not exactly. Under MSE, the point forecast is the conditional mean. Under symmetric unimodal distributions, the mean equals the mode (most likely), but this is not guaranteed for skewed distributions. The forecast minimizes expected loss, which is a different statement than “most likely.”
“A wider prediction interval means a worse model.” Not necessarily. A wider PI means more honest uncertainty. A model with a narrow PI that under-covers is worse than a model with a wide PI that covers correctly. Width reflects the forecast-error variance, which depends on the horizon, the persistence of the series, and $\sigma^2$.
“The 95% PI means 95 out of 100 future observations will fall inside.” Close, but the coverage guarantee is conditional on the model being correct and the innovations being Gaussian. If either assumption fails, actual coverage drifts from nominal — and Section 6.6 showed the drift is worst in the far tail: under fat-tailed innovations the 95% interval held up while the 99% interval missed more than twice as often as promised.
“I should always use MSE because it is the most common.” MSE is the default, not the universal answer. If your application penalizes large errors disproportionately, MSE is appropriate. If your application treats all errors proportionally, use MAE. If your application has asymmetric costs, use asymmetric loss. Huber is a useful compromise. The choice is domain-specific.
“Random train/test splits work fine for time series if the sample is large enough.” No. Size does not fix the temporal-ordering problem. A random split allows information from the future to leak into the training set, making out-of-sample evaluation meaningless.
“Long-horizon forecasts from ARMA models are useful.” For stationary models, the forecast converges to the unconditional mean and the PI converges to the unconditional variance. Beyond a certain horizon, the forecast is uninformative — it is just the average. The useful horizon depends on the persistence ($\phi$ near 1 extends it; $\phi$ near 0 shortens it).
Plugging R’s intercept into the recursion as $\hat{\alpha}$. For stationary fits, Arima()’s intercept is the estimated mean $\hat{\mu}$. Recover $\hat{\alpha} = \hat{\mu}(1 - \sum_j \hat{\phi}_j)$ before running the hand recursion, or your numbers will not match forecast().
“My model lost to the naive forecast, so the modeling was pointless.” Losing to the naive over one calm window is information, not indictment — it tells you the window had nothing for your model’s extra machinery to exploit, and it reminds you that the margin needs a formal test (Module 7) before you conclude anything. And the model still produced calibrated uncertainty bands, which the naive rule cannot.

Connection to Enders

Forecasting with ARMA models: Enders Chapter 2, pp. 68–78
Forecast-error variance and prediction intervals: Enders Chapter 2, pp. 73–76
The Wold decomposition and $\psi$-weights: Enders Chapter 2, pp. 51–55

A convention warning when you cross-reference. As in Modules 3–5, Enders writes AR coefficients as $a_i$ where we write $\phi_j$, and works in the characteristic-root-inside-the-unit-circle framing where the course states stationarity as roots of $\Phi(L)$ lying outside the unit circle — same condition, reciprocal roots.

Hyndman & Athanasopoulos (Forecasting: Principles and Practice, 3rd ed.) covers loss functions, train/test splits, and the forecast package in detail. Chapter 5 (evaluating forecast accuracy, including MASE) and Chapter 9 (ARIMA forecasting) are particularly relevant. Note: Hyndman uses a strict stationarity definition (stronger than the course’s weak/covariance version); the forecasting material is compatible regardless.

Hamilton (Time Series Analysis) covers the forecast-error variance derivation in Chapter 4.3 with full matrix algebra. The conditional-mean optimality result is in Chapter 4.1.

Practice Problems

Core Practice

Simulate an AR(1) with $\phi = 0.7$, $\sigma = 1$, $T = 300$. Fit the correct AR(1). Compute $\hat{y}_{T+h|T}$ by hand for $h = 1, 2, 3, 4, 5$ using the recursive formula — remember that Arima()’s intercept is $\hat{\mu}$, so recover $\hat{\alpha} = \hat{\mu}(1 - \hat{\phi})$ first. Compare your answers to forecast(fit, h = 5)$mean. Report both columns. At what horizon does the forecast effectively reach the unconditional mean?
Repeat Problem 1 for an ARMA(1,1) with $\phi = 0.6$, $\theta = 0.4$. Explain why the MA term affects only the $h = 1$ forecast and disappears from $h = 2$ onward.
Run the coverage Monte Carlo from Section 6.6 (1,000 reps, ARMA(1,1), 1-step-ahead PIs at the 80%, 95%, and 99% levels) for Gaussian innovations and for $t_5$ innovations rescaled to unit variance (less fat-tailed than $t_3$; note $\text{Var}(t_5) = 5/3$). Report the six coverage rates in a table. At which levels does the $t_5$ case deviate from nominal, in which direction, and how does the pattern compare to the $t_3$ results in the notes?
For the cached UNRATE series (1960–2019, as in these notes), split with test_periods = 36. Fit the L5 winner ARIMA$(1,1,2)(1,0,1)_{12}$ on training. Produce a 36-month forecast. Compute MSE, MAE, and RMSE on the test set. Also compute the same metrics for a naive random-walk forecast. Which forecast wins on which metrics, and — given Section 6.8 — should the outcome change your opinion of the model? Produce a fan chart with fan = TRUE and write two sentences interpreting the picture.
Using the same UNRATE split from Problem 4, also fit a non-seasonal ARIMA(1,1,1) on training and forecast 36 months. Compare it to the seasonal model on MSE and MAE. Does the seasonal model win on both metrics? Write one paragraph explaining why or why not.
Conceptual. Explain in three to four sentences why a random train/test split is invalid for time series forecasting. Give a concrete example of information leakage that a random split would cause.
Conceptual. A hospital administrator needs to forecast daily ER admissions for staffing. Under-forecasting (too few staff) has a much higher cost than over-forecasting (too many staff). Which loss function should they use, and what does the optimal forecast under that loss represent?
Challenge. Write a function coverage_mc(phi, theta, sigma, n_obs, n_reps, levels, innov_fun) that takes an innovation-generating function (e.g., rnorm or function(n) rt(n, df = 3) / sqrt(3)) and a vector of confidence levels, and returns the empirical coverage of each level’s PI. Use it to produce a table of coverage rates at the 80%, 95%, and 99% levels for Gaussian, $t_5$, and $t_3$ innovations.

Key Takeaways

A forecast is a conditional mean. Under MSE, the optimal point forecast is $\hat{y}_{t+h|t} = \mathbb{E}[y_{t+h} \mid \mathcal{F}_t]$. Computing it is a recursion: replace future innovations with zero, future $y$’s with their forecasts.
Forecast uncertainty grows with the horizon. The forecast-error variance $\sigma^2 \sum_{i=0}^{h-1} \psi_i^2$ increases with $h$ because each step adds an unresolved shock. For stationary models, it converges to the unconditional variance.
Prediction intervals assume Gaussianity — and fat tails break the promise where it matters most. Under fat-tailed innovations with matched variance, nominal 95% coverage roughly survives, the 80% interval over-covers, and the 99% interval misses more than twice as often as promised. The failure lives in the far tail. Bootstrap PIs are the pragmatic fix.
Loss functions are a choice, not a default. MSE optimizes the mean; MAE optimizes the median; Huber compromises; asymmetric loss optimizes a quantile. Under skewed errors, mean-targeting and median-targeting forecasts genuinely trade places depending on the loss. Choose the loss before you compute the forecast.
Train/test splits must respect time order. Random splits leak future information. Use ts_split() or rolling windows.
A forecast needs a baseline — and expect the baseline to be tough. Raw loss numbers are uninformative without a comparison. The naive random walk is the universal benchmark precisely because, on persistent series in calm periods, it often wins. Whether a margin is real is a statistical question — Module 7’s Diebold-Mariano test answers it.
Fan charts are the professional output. Point forecast without uncertainty bands is a headline without the story — and calibrated bands are the part of the product a naive rule can never supply.

Looking Ahead — Module 7

You now have a forecast object: point forecasts, prediction intervals, a fan chart, and loss numbers against a baseline. Section 6.8 ended on an unresolved note that is really a question for a statistician: the naive forecast won by a few hundredths — does that margin mean anything?

Module 7 (Forecast Evaluation and Combinations) picks up exactly there:

The Diebold-Mariano test: a formal answer to “is forecast A significantly better than forecast B?” — the hypothesis test that Section 6.8’s eyeball comparison was begging for.
Evaluation metrics in full: RMSE, MAE, MAPE, Theil’s U, and MASE (previewed in Section 6.7) as a working toolkit.
Rolling-window evaluation: replacing this module’s single fixed-origin split with re-fitting as each observation arrives — more data points for the comparison, and the design the DM test wants.
Forecast combinations: when two forecasts disagree, don’t pick — average. Equal weights, MSE-inverse weights, OLS (Granger-Ramanathan) weights, and the striking empirical regularity that the simple average is almost impossible to beat.

The theme of this module was producing one honest forecast. The theme of the next is refereeing between several.

Term	L3	L4	L5	L6
\(\alpha\), \(\delta t\)	On	On	On	On
\(\phi_j y_{t-j}\) (all \(p\))	On	On	On	On
\(\theta_l \epsilon_{t-l}\) (all \(q\))	On	On	On	On
Seasonal block \((P, D, Q)_s\)	Off	Off	On	On
\(\epsilon_t\)	On	On	On	On