Lecture 12: Forecast Evaluation
02 October 2026
By now we can produce forecasts from trends, seasonals, autoregressions, GARCH models, VARs, and from no model at all. Each application ended with a plot of forecast against realization and a number or two. That is not evaluation.
Two questions make it evaluation. Absolute: does a forecast have the properties an optimal forecast must have, so that no cheap fix would improve it? Relative: is one forecast more accurate than another, and is the difference more than luck?
Today: the unforecastability principle and the tests it implies; accuracy measures, benchmarks and forecastability; the Diebold-Mariano test; two applications, a shipping company’s forecasts and quarterly GDP growth; and how to evaluate interval and density forecasts, with the volatility models of a few weeks ago as the test case.
Diebold, Forecasting in Economics, Business, Finance and Beyond, the chapters on evaluating point forecasts and on evaluating interval and density forecasts. The shipping data come from the book’s companion files; the professional forecasts from the Survey of Professional Forecasters.
For a covariance stationary series with Wold representation y_t = \mu + \varepsilon_t + b_1 \varepsilon_{t-1} + b_2 \varepsilon_{t-2} + \cdots, the h-step linear least squares forecast and its error are
y_{t+h,t} = \mu + b_h \varepsilon_t + b_{h+1}\varepsilon_{t-1} + \cdots, \qquad e_{t+h,t} = \varepsilon_{t+h} + b_1 \varepsilon_{t+h-1} + \cdots + b_{h-1}\varepsilon_{t+1}, \qquad \sigma_h^2 = \sigma^2 \Big(1 + \sum_{i=1}^{h-1} b_i^2\Big).
The error is made entirely of shocks that arrive after the forecast is made. Hence the one property from which all others follow, the unforecastability principle:
Optimal forecast errors are unforecastable from information available when the forecast was made.
It holds in great generality: for conditional-mean or linear-projection optimality, for any loss function, and for nonstationary series too. If you could forecast the error, you could add that forecast to the forecast and do better; so an optimal forecast leaves nothing forecastable behind.
Each is a testable restriction on a series we observe: the forecast errors. The tests that follow are all applications of the same idea, with one wrinkle: for h > 1 the errors are serially correlated even when the forecast is perfect, so standard errors must allow for it.
Property (c) is what made the multi-step errors in the ARMA forecasting lecture an MA(h-1).
Regress the forecast errors on a constant and test whether the mean is zero.
In R, arima(e, order = c(0, 0, h - 1)) estimates the mean with MA(h-1) disturbances and reports its standard error; sandwich::NeweyWest() gives the HAC alternative. These slides use a helper, hac_t(e, lag = h - 1), that returns the HAC t statistic of the mean.
“Bias” is a population statement, so even a clearly nonzero sample mean error needs a standard error before it becomes a verdict.
Serial correlation in one-step errors is the most common symptom of a fixable forecast: the error was forecastable from its own past, so a model of the errors would have improved the forecast.
The tests so far use only the errors’ own past. The principle says more: errors should be orthogonal to anything in the information set. Regress them on candidate variables and test that all slopes are zero:
e_{t+h,t} = \alpha_0 + \sum_i \alpha_i x_{it} + u_t, \qquad H_0: \alpha_0 = \alpha_1 = \cdots = 0 .
The most important candidate is the forecast itself, which is certainly in the information set. Equivalently, regress the realization on the forecast, the Mincer-Zarnowitz regression:
y_{t+h} = \beta_0 + \beta_1\, y_{t+h,t} + u_t, \qquad H_0: (\beta_0, \beta_1) = (0, 1) .
Subtracting y_{t+h,t} from both sides shows the two regressions are the same test: (\beta_0, \beta_1) = (0, 1) is (\alpha_0, \alpha_1) = (0, 0).
\beta_0 \neq 0 is bias; \beta_1 < 1 says the forecast over-reacts to its own information and should be shrunk toward the mean. Rejection is constructive: the fitted regression is itself an improved forecast. The Wald test of the joint hypothesis uses HAC standard errors when h > 1.
Mincer and Zarnowitz (1969). mz_test(y, f, lag = h - 1) in these slides reports the coefficients and the Wald test.
Ranking forecasts needs a loss function L(e) (recall the loss-function lecture) and a horizon. With errors e_{t+h,t} = y_{t+h} - y_{t+h,t} over t = 1, \ldots, T:
| Measure | Sample formula | What it measures |
|---|---|---|
| Mean error (bias) | \hat\mu_e = \frac{1}{T}\sum_t e_{t+h,t} | location |
| Error variance | \hat\sigma_e^2 = \frac{1}{T}\sum_t (e_{t+h,t} - \hat\mu_e)^2 | dispersion |
| Mean squared error | \widehat{MSE} = \frac{1}{T}\sum_t e_{t+h,t}^2 = \hat\sigma_e^2 + \hat\mu_e^2 | both, under quadratic loss |
| Root MSE | \sqrt{\widehat{MSE}} | the same, in the units of y |
| Mean absolute error | \widehat{MAE} = \frac{1}{T}\sum_t |e_{t+h,t}| | both, under absolute loss |
| Mean absolute percent error | \frac{1}{T}\sum_t |e_{t+h,t} / y_{t+h}| | scale-free, but fails near y = 0 |
MSE = \text{variance} + \text{bias}^2 is the bias-variance tradeoff: under quadratic loss a biased forecast can beat an unbiased one with more dispersed errors. A small bias for a large variance reduction is a good trade.
MSE squares the units; RMSE and MAE keep them. MAE is less sensitive to a few huge errors.
Accuracy numbers mean little in isolation. Compare them with an allegedly naive benchmark:
R^2_{\text{pred}} = 1 - \frac{\sum_t e_{t,t-1}^2}{\sum_t (y_t - \bar y)^2} \qquad\qquad U = 1 - \frac{\sum_t e_{t,t-1}^2}{\sum_t (y_t - y_{t-1})^2} .
Naive benchmarks are not always naive. Many economic series are nearly random walks, and then no forecaster beats “no change” by much, through no fault of their own: the yen of the unit-root lecture. A skill score near zero may be the honest ceiling.
Theil (1966). Both benchmarks must be computed out of sample too: an expanding-window mean, and the last observation.
Even an optimal forecast makes errors; how large they are is a property of the series. Under quadratic loss, patterned after R^2,
G = 1 - \frac{\operatorname{var}(e_{t+j,t})}{\operatorname{var}(y_{t+j})}, \qquad\text{and more generally}\qquad P(L, \Omega, j, k) = 1 - \frac{E\,L(e_{t+j,t})}{E\,L(e_{t+k,t})}, \quad j \ll k :
the expected loss of an optimal short-horizon forecast relative to that of a long-horizon one. G is the special case of quadratic loss and k = \infty.
Diebold and Kilian (2001). In practice P is estimated from a fitted model’s forecast error variances.
The AR(2) of the cycles lectures, \hat\phi_1 = 1.44, \hat\phi_2 = -0.48, implies \sigma_h^2 = \sigma^2 \sum_{i<h} \psi_i^2 and \operatorname{var}(y) = \sigma^2 \sum_i \psi_i^2, hence P(h) = G at each horizon:
## h=1 h=2 h=4 h=8 h=12 h=20
## 0.96 0.87 0.66 0.34 0.17 0.04Employment is 96 percent forecastable one quarter ahead, two thirds a year ahead, and 4 percent five years ahead. The same calculation for quarterly GDP growth gives about 0.1 at h = 1: GDP growth is nearly unforecastable from its own past, which sets the bar for the application below.
Two forecasts of US inflation: the Survey of Professional Forecasters (S) and the breakeven rate from inflation-indexed bonds (B). Over some sample, \widehat{MSE}_S = 1.80 and \widehat{MSE}_B = 1.92. “S wins.”
Does it? Even if the two population MSEs were identical, one of the two sample MSEs would be smaller. The question in any sample is whether S is truly superior or merely lucky. That is a hypothesis test, with the hypothesis of equal predictive accuracy
H_0: E\big[L(e_{a,t+h,t})\big] = E\big[L(e_{b,t+h,t})\big] \quad\Longleftrightarrow\quad E(d_t) = 0, \qquad d_t = L(e_{a,t+h,t}) - L(e_{b,t+h,t}) .
The population statement is about expected loss; we test it with the sample average of the loss differential d_t, a constructed but observed series.
The forecasts may come from models, surveys, markets or judgment: the test treats the forecast errors as primitives and never asks where they came from.
Assume only that the loss differential is covariance stationary:
\text{Assumption DM:}\qquad E(d_t) = \mu, \qquad \operatorname{cov}(d_t, d_{t-\tau}) = \gamma(\tau), \qquad 0 < \operatorname{var}(d_t) < \infty .
Then, under H_0: \mu = 0,
DM = \frac{\bar d}{\hat\sigma_{\bar d}} \;\overset{d}{\longrightarrow}\; N(0, 1), \qquad \bar d = \frac{1}{T}\sum_{t=1}^{T} d_t ,
where \hat\sigma_{\bar d} is a consistent estimate of the standard deviation of the sample mean, which must allow for serial correlation in d_t: forecast errors overlap when h > 1, and suboptimal forecasts have serially correlated errors at any horizon.
That is all: one assumption, one statistic, one limiting distribution. DM is a t test that the mean of d_t is zero, with the standard error computed robustly. Regress d_t on a constant with HAC standard errors, or with MA(h-1) or AR(p) disturbances chosen by AIC.
Diebold and Mariano (1995). Covariance stationarity is sufficient, not necessary; mixing conditions would do.
Extensions. The regression form invites conditioning: add variables that may explain the loss differential, such as a recession indicator, to ask when one forecast is better (Giacomini and White, 2006). The loss need not be quadratic: |e|, direction-of-change, anything.
Diebold (2015) revisits the test twenty years on. Comparing nested models estimated on the same data changes the limit distribution (West, 1996; Clark and McCracken, 2001); treating forecasts as primitives sidesteps this.
OverSea Services ships cargo across the Atlantic and forecasts weekly volume on each trade lane twice: a quantitative forecast from a model and a judgmental forecast from its sales representatives. Here: the North America to Europe lane, realized volume and the two two-week-ahead forecasts, 499 weeks.
Both forecasts track the realization closely: two weeks ahead is not far. The differences are in the details.
data/shipping_volume.csv, January 1988 to July 1997: Diebold’s textbook application, data from its companion files.
## mean sd RMSE MAE skewness kurtosis
## quantitative -0.027 1.263 1.262 1.003 -0.200 3.181
## judgmental 1.024 1.064 1.476 1.225 -0.106 3.073Quantitative errors are centered on zero; judgmental errors sit about one unit above it (realized volume ran higher than the sales force predicted, a pessimistic bias) but are less dispersed. Both look normal.
In both cases the autocorrelations cut off after displacement 1 (0.52 and 0.50) and the partial autocorrelations oscillate and decay: the MA(1) structure that optimal two-step-ahead errors should have. Property (c) passes.
Two-step errors are MA(1), so regress each error series on a constant with MA(1) disturbances and read off the t statistic of the mean:
bias <- function(e) {
m <- arima(e, order = c(0, 0, 1)) # mean with MA(1) disturbances
mu <- coef(m)["intercept"]; se <- sqrt(m$var.coef["intercept", "intercept"])
round(c(mean = unname(mu), s.e. = se, t = unname(mu / se), theta = unname(coef(m)["ma1"])), 3)
}
rbind(quantitative = bias(sh$eq), judgmental = bias(sh$ej))## mean s.e. t theta
## quantitative -0.027 0.080 -0.333 0.934
## judgmental 1.023 0.067 15.261 0.960## HAC_t_quantitative HAC_t_judgmental
## -0.3819108 17.6025359No bias in the quantitative forecast. A bias of one unit in the judgmental forecast, with a t statistic of 15: sizable and overwhelmingly significant. The HAC version agrees.
Ignoring the MA(1) term would badly understate the standard error of the mean; Diebold reports the same estimates.
Regress realized volume on each forecast, with HAC standard errors for the two-step overlap, and test (\beta_0, \beta_1) = (0, 1):
## beta0 se0 beta1 se1 Wald p
## quantitative 2.162 0.385 0.89 0.019 34.854 0
## judgmental 2.339 0.332 0.93 0.018 352.819 0Both fail, decisively. We expected the judgmental forecast to fail, because it is biased. But until now the quantitative forecast had passed every check; the MZ regression finds a slope well below one: the forecast over-reacts to its own information, and shrinking it toward the mean, 2.2 + 0.89\, \hat y, would have been better. Rejection is constructive.
Diebold runs the same regressions with MA(1) disturbances instead of HAC standard errors; same verdict. The orthogonality tests use only information available at forecast time, so the improvement was feasible.
Under quadratic loss the quantitative forecast wins: RMSE 1.26 against 1.48 (MAE 1.00 against 1.23). The judgmental errors are less dispersed, so the whole gap is the bias. Is the gap real? Form the loss differential:
## mean sd acf1 acf2 acf3
## -0.585 3.416 0.357 -0.069 -0.050
The mean, -0.58, is small relative to the standard deviation, 3.4, but the 499 observations and the MA(1) serial correlation (\hat\rho(1) = 0.36, then nothing) are what the test must account for.
## DM_MA1 DM_HAC
## -2.865698 -3.286535## DM
## -2.915281All three put DM near -3, with a p-value below 0.01. The quantitative forecast is more accurate under quadratic loss, and the verdict is not luck. Under absolute loss (power = 1) the statistic is -3.5.
And yet the loser has the smaller error variance. A bias-corrected judgmental forecast would be a serious competitor, and the two together might beat either. That is next week’s subject.
dm.test() uses a rectangular-window variance estimate with h - 1 autocovariances; the Newey-West version uses Bartlett weights; the MA(1) version models the dependence parametrically. They agree because the loss differential really is MA(1).
Quarterly real GDP growth (annualized percent) and the term spread from the second homework. Recursive one-quarter-ahead forecasts: at each origin T from 1984Q4, estimate on 1959 to T and forecast T + 1; 140 forecasts, 1985Q1 to 2019Q4.
g <- read.csv("data/hw2_gdp_spread.csv") |> mutate(Quarter = yearquarter(quarter))
y <- g$gdp_growth; s <- g$spread
origins <- which(g$Quarter == yearquarter("1984 Q4")):(which(g$Quarter == yearquarter("2019 Q4")) - 1)
fc <- t(sapply(origins, function(T) { # one row of forecasts per origin T
yi <- y[1:T]; si <- s[1:T]
d <- data.frame(y = yi[-(1:2)], l1 = yi[2:(T - 1)], l2 = yi[1:(T - 2)], sp = si[2:(T - 1)])
nd <- data.frame(l1 = yi[T], l2 = yi[T - 1], sp = si[T]) # the regressors for T + 1
f <- function(m) unname(predict(lm(m, d), nd))
c(mean = mean(yi), nochange = yi[T], ar1 = f(y ~ l1), ar2 = f(y ~ l1 + l2),
ar2_spread = f(y ~ l1 + l2 + sp))
}))
ev <- tibble(Quarter = g$Quarter[origins + 1], actual = y[origins + 1]) |> bind_cols(as_tibble(fc))Five model forecasts: the expanding-window mean, no change, AR(1), AR(2), and the AR(2) with the lagged spread, significant in the homework. A sixth is model-free: the SPF median.
data/hw2_gdp_spread.csv. The evaluation stops in 2019: the 2020 quarters (-30 and +35 percent) would account for most of every squared-error measure.
spf <- read.csv("data/spf_median_rgdp_growth.csv") |>
transmute(Quarter = yearquarter(paste(YEAR, QUARTER, sep = " Q")) + 1, # target: survey quarter + 1
spf = DRGDP3)
ev <- ev |> left_join(spf, by = "Quarter")
models <- c("mean", "nochange", "ar1", "ar2", "ar2_spread", "spf")
err <- as_tibble(sapply(models, function(m) ev$actual - ev[[m]])) # errors: actual minus forecast
acc <- function(e) c(ME = mean(e), RMSE = sqrt(mean(e^2)), MAE = mean(abs(e)))
round(t(sapply(err, acc)), 3)## ME RMSE MAE
## mean -0.697 2.330 1.630
## nochange -0.004 2.524 2.036
## ar1 -0.457 2.147 1.557
## ar2 -0.360 2.082 1.525
## ar2_spread -0.732 2.264 1.643
## spf 0.092 2.066 1.571The AR(2) is the most accurate model by RMSE and MAE; adding the spread makes the forecasts worse. The professionals, with every indicator at hand, match the AR(2) and are the only unbiased forecast in the table.
data/spf_median_rgdp_growth.csv: Philadelphia Fed, median forecast DRGDP3. The survey is taken mid-quarter, so its “next quarter” forecast made in quarter t targets growth in t+1. Realizations are the latest vintage; the SPF forecast the first release.
## R2_ar2 R2_ar2_spread R2_spf U_ar2 U_mean
## 0.202 0.056 0.214 0.320 0.148## beta0 se0 beta1 se1 Wald p
## ar2 -0.871 0.830 1.170 0.252 4.344 0.114
## ar2_spread 0.362 0.600 0.676 0.153 30.233 0.000
## spf 0.229 0.729 0.946 0.251 0.277 0.871AR(2) errors: mean -0.36 (HAC t = -2.1), no serial correlation (Q(8): p = 0.85), and t = 0.9 on the lagged spread.
Iterate the recursive AR(2) forecasts to horizons two to four and look at the errors:
hfc <- t(sapply(origins, function(T) { # AR(2) chain-rule forecasts, h = 1, ..., 4
yi <- y[1:T]; b <- coef(lm(yi[-(1:2)] ~ yi[2:(T - 1)] + yi[1:(T - 2)]))
path <- yi[(T - 1):T]
for (h in 1:4) path <- c(path, b[1] + b[2] * path[length(path)] + b[3] * path[length(path) - 1])
path[-(1:2)]
}))
t(sapply(1:4, function(h) { # errors at each horizon, within the sample
ok <- origins + h <= max(origins) + 1; e <- y[origins[ok] + h] - hfc[ok, h]
round(c(h = h, RMSE = sqrt(mean(e^2)), acf = acf(e, 3, plot = FALSE)$acf[-1]), 3)
}))## h RMSE acf1 acf2 acf3
## [1,] 1 2.082 0.077 0.071 -0.012
## [2,] 2 2.183 0.292 0.083 0.001
## [3,] 3 2.298 0.331 0.278 0.018
## [4,] 4 2.315 0.357 0.299 0.109The RMSE rises with the horizon and flattens quickly near the standard deviation of growth (2.25), property (d); the error autocorrelations grow with h as an MA(h-1) allows, property (c). A series this unpredictable reaches its unconditional variance within a year.
dm <- function(a, b) hac_t(err[[a]]^2 - err[[b]]^2, lag = 0) # negative favors the first forecast
round(c(`AR(2) vs mean` = dm("ar2", "mean"), `AR(2) vs no change` = dm("ar2", "nochange"),
`AR(2) vs AR(2)+spread` = dm("ar2", "ar2_spread"), `AR(2) vs AR(1)` = dm("ar2", "ar1"),
`SPF vs AR(2)` = dm("spf", "ar2"), `SPF vs mean` = dm("spf", "mean")), 2)## AR(2) vs mean AR(2) vs no change AR(2) vs AR(2)+spread
## -2.41 -2.65 -2.80
## AR(2) vs AR(1) SPF vs AR(2) SPF vs mean
## -2.22 -0.18 -1.70One-step errors, so lag = 0; a Newey-West lag of 4 changes nothing. Under absolute loss the AR(2)-versus-spread statistic is -1.8: the spread model’s extra errors are a few large ones.
A (1-\alpha) interval forecast is unconditionally calibrated if it brackets the realization (1-\alpha) of the time in the long run. It is conditionally calibrated if it does so at every date, given what was known then.
The volatility lecture made the distinction concrete. Constant-width intervals in a GARCH world can be right on average and wrong at almost every date: too wide in calm periods, too narrow in turbulent ones, with misses clustered in the turbulent ones. The 70-day comparison there was suggestive; a proper evaluation needs a long sequence and a test.
Define the hit sequence of a one-step interval forecast,
I_t = \mathbf{1}\{ y_t \text{ falls inside the interval forecast made at } t-1 \}.
Under correct conditional calibration, I_t \sim \text{iid Bernoulli}(1 - \alpha): the right mean and independence. The hit sequence is the interval forecast’s “error”, and like every one-step error it should be unforecastable.
Christoffersen (1998). For h-step intervals the hits need not be independent but should have (h-1)-dependent structure.
Two likelihood-ratio tests separate the two parts of the hypothesis. With n_1 hits and n_0 misses, and \hat\pi = n_1 / (n_0 + n_1),
LR_{uc} = -2 \ln \frac{(1-p)^{n_0} p^{n_1}}{(1-\hat\pi)^{n_0} \hat\pi^{n_1}} \;\sim\; \chi^2_1 \quad\text{under } H_0: \pi = p = 1 - \alpha ,
tests unconditional coverage. For independence, fit a first-order Markov chain to the hit sequence and test whether the probability of a hit depends on whether the last period was a hit:
LR_{ind} = -2 \ln \frac{(1-\hat\pi)^{n_{00}+n_{10}} \hat\pi^{n_{01}+n_{11}}}{(1-\hat\pi_{01})^{n_{00}} \hat\pi_{01}^{n_{01}} (1-\hat\pi_{11})^{n_{10}} \hat\pi_{11}^{n_{11}}} \;\sim\; \chi^2_1 ,
with n_{ij} the number of transitions from state i to state j and \hat\pi_{ij} = n_{ij} / \sum_j n_{ij}. LR_{cc} = LR_{uc} + LR_{ind} \sim \chi^2_2 tests both at once.
Christoffersen (1998). Under H_0 the hit sequence has no memory, so \hat\pi_{01} = \hat\pi_{11}; a cluster of misses after a miss shows up as \hat\pi_{01} < \hat\pi_{11}.
Both tests are a few lines once the hit sequence is in hand:
lr_uc <- function(hit, p) { # unconditional coverage: is P(hit) = p?
n1 <- sum(hit); n0 <- sum(!hit); pi_hat <- n1 / (n0 + n1)
LR <- -2 * (n0 * log(1 - p) + n1 * log(p) - n0 * log(1 - pi_hat) - n1 * log(pi_hat))
c(LR_uc = LR, p_uc = pchisq(LR, 1, lower.tail = FALSE))
}
lr_ind <- function(hit) { # independence: first-order Markov chain
n <- table(factor(head(hit, -1), c(FALSE, TRUE)), factor(tail(hit, -1), c(FALSE, TRUE)))
p01 <- n[1, 2] / sum(n[1, ]); p11 <- n[2, 2] / sum(n[2, ]); p1 <- sum(n[, 2]) / sum(n)
l0 <- sum(n[, 1]) * log(1 - p1) + sum(n[, 2]) * log(p1) # log likelihood if hits are independent
l1 <- n[1, 1] * log(1 - p01) + n[1, 2] * log(p01) + n[2, 1] * log(1 - p11) + n[2, 2] * log(p11)
c(LR_ind = -2 * (l0 - l1), p_ind = pchisq(-2 * (l0 - l1), 1, lower.tail = FALSE))
}A density forecast p_t(y) is right if it equals the true conditional density f_t(y), which is never observed, not even after the fact. Yet it can be tested. Evaluate the forecast’s cdf at the realization:
z_t = P_t(y_t) = \int_{-\infty}^{y_t} p_t(u)\, du .
If p_t = f_t for every t, the sequence z_t is iid U(0,1) (Rosenblatt, 1952; Diebold, Gunther and Tay, 1998): every quantile of a correct density is hit with the right frequency and without clustering.
Diagnostics, constructive rather than omnibus:
The Kolmogorov-Smirnov test of uniformity is valid but says nothing about why a density forecast fails.
Three forecasters of tomorrow’s return, re-estimated every 250 days and rolled one day at a time over the last 1500 days: a constant-variance normal, a GARCH(1,1) with normal shocks, and one with Student’s t shocks.
r <- 100 * scan("data/nyse_returns.csv", skip = 1, quiet = TRUE)
roll <- function(dist) { # rolling one-step GARCH(1,1) density forecasts
spec <- ugarchspec(variance.model = list(garchOrder = c(1, 1)), mean.model = list(armaOrder = c(0, 0)),
distribution.model = dist)
as.data.frame(ugarchroll(spec, r, n.ahead = 1, forecast.length = 1500, refit.every = 250,
refit.window = "recursive", solver = "hybrid")) # Mu, Sigma, Shape, Realized
}
gn <- roll("norm"); gt <- roll("std")
idx <- (length(r) - 1499):length(r); yr <- r[idx]
cm <- sapply(idx, function(i) mean(r[1:(i - 1)])); cs <- sapply(idx, function(i) sd(r[1:(i - 1)]))
hits <- list(constant = abs(yr - cm) < 1.96 * cs, garch_normal = abs(yr - gn$Mu) < 1.96 * gn$Sigma,
garch_t = abs(yr - gt$Mu) < qdist("std", 0.975, shape = gt$Shape) * gt$Sigma)
round(t(sapply(hits, function(h) c(coverage = mean(h), lr_uc(h, 0.95), lr_ind(h)))), 3)## coverage LR_uc p_uc LR_ind p_ind
## constant 0.881 111.205 0.000 2.472 0.116
## garch_normal 0.934 7.378 0.007 4.331 0.037
## garch_t 0.941 2.602 0.107 3.884 0.049ugarchroll() keeps the parameters fixed between refits but updates the conditional variance every day. The sample is 1988 to 2001; the 1500 evaluation days start in early 1996 and include the 1997 and 1998 bursts.
The independence tests are weaker than the eye: clustering of misses is visible for the constant interval but LR_{ind} is a first-order test on a short series of mostly ones. The histogram and correlograms of the PITs are more informative.
z <- list(constant = pnorm(yr, cm, cs), garch_normal = pnorm(yr, gn$Mu, gn$Sigma),
garch_t = pdist("std", (yr - gt$Mu) / gt$Sigma, shape = gt$Shape))
pit_plot(z$constant, sub = "constant-variance normal") |
pit_plot(z$garch_normal, sub = "GARCH(1,1), normal") |
pit_plot(z$garch_t, sub = "GARCH(1,1), Student's t")
Dashed lines: the band an iid uniform sample of this size respects bin by bin. Constant variance shows the butterfly, too many realizations in both tails: densities too thin-tailed. The GARCH densities are close to uniform; the t version leaves the outer deciles a little full.
Under constant variance the squared PITs are strongly autocorrelated: neglected volatility dynamics show up as forecastable “errors”, as serial correlation does for point forecasts. GARCH halves the autocorrelations but does not remove them: refits only every 250 days leave some variance dynamics uncaptured.
Diebold runs the same sequence on S&P 500 returns; the PITs themselves also keep a small first-order autocorrelation.
library(fpp3); library(sandwich); library(lmtest)
# ---- absolute standards -----------------------------------------------------------------
arima(e, order = c(0, 0, h - 1)) # mean of h-step errors with MA(h-1) disturbances
coeftest(lm(e ~ 1), vcov = NeweyWest(lm(e ~ 1), lag = h - 1, prewhite = FALSE)) # HAC t test of zero mean
hac_t(e, lag = h - 1) # the same, one number (defined in these slides)
correlogram(e, m = 12); Box.test(e, 12, "Ljung-Box") # white noise / MA(h-1) checks
mz_test(y, f, lag = h - 1) # Mincer-Zarnowitz: y on f, Wald test of (0, 1)
# ---- relative standards -----------------------------------------------------------------
c(ME = mean(e), RMSE = sqrt(mean(e^2)), MAE = mean(abs(e))) # accuracy measures
1 - sum(e^2) / sum(e_mean^2); 1 - sum(e^2) / sum(e_nochange^2) # predictive R2; Theil's U
hac_t(e1^2 - e2^2, lag = h - 1) # Diebold-Mariano statistic by regression
forecast::dm.test(e1, e2, h = h, power = 2) # the packaged version; power = 1 for absolute loss
accuracy(fc, data) # fable: RMSE, MAE, MAPE, MASE for a fable of forecasts
# ---- interval and density forecasts -----------------------------------------------------
hit <- abs(y - f) < 1.96 * sigma; mean(hit); lr_uc(hit, 0.95); lr_ind(hit) # coverage and Christoffersen tests
z <- pnorm(y, mu, sigma); pit_plot(z); acf_plot((z - mean(z))^2) # PITs, histogram, dependence
ugarchroll(spec, r, n.ahead = 1, forecast.length = 1500, refit.every = 250) # rolling GARCH density forecastsfable.data/shipping_volume.csv (OverSea, 499 weeks); data/hw2_gdp_spread.csv (GDP growth and term spread, 1959Q1 to 2024Q2); data/spf_median_rgdp_growth.csv (Survey of Professional Forecasters, Federal Reserve Bank of Philadelphia); data/nyse_returns.csv.