US GNP growth — R forecast::auto.arima, out-of-sample¶

Companion to the Python/PyMC notebooks. auto.arima on the stationary growth series (d=0), forecast() for the held-out 4 years, ggplot/autoplot with NBER recessions shaded.

In [1]:
suppressMessages({library(forecast); library(ggplot2)})
d <- read.csv('hamilton_gnp.csv'); d <- d[!is.na(d$gnp_growth),]
y <- ts(d$gnp_growth, start=c(1951,2), frequency=4)
H <- 16; train <- head(y, length(y)-H); test <- tail(y, H)
fit <- auto.arima(train, seasonal=FALSE, d=0, stepwise=TRUE, approximation=FALSE)
print(fit)
fc <- forecast(fit, h=H); acc <- accuracy(fc, test)
cat(sprintf('\nOUT-OF-SAMPLE RMSE = %.2f   naive-mean RMSE = %.2f\n',
            acc['Test set','RMSE'], sqrt(mean((as.numeric(test)-mean(train))^2))))

Series: train 
ARIMA(2,0,2) with non-zero mean 

Coefficients:
         ar1      ar2      ma1     ma2    mean
      0.6283  -0.7102  -0.3820  0.8138  0.7778
s.e.  0.1278   0.1308   0.1018  0.1197  0.1164

sigma^2 = 0.9622:  log likelihood = -164.39
AIC=340.79   AICc=341.54   BIC=357.46

OUT-OF-SAMPLE RMSE = 1.21   naive-mean RMSE = 1.20
In [2]:
# recession bands as date ranges for ggplot
yr <- d$date[!is.na(d$date)]; recs <- d$recession
starts <- which(diff(c(0,recs))==1); ends <- which(diff(c(recs,0))==-1)
tt <- as.numeric(time(y)); bands <- data.frame(x1=tt[starts], x2=tt[ends])
options(repr.plot.width = 9, repr.plot.height = 5)
autoplot(fc) +
  geom_rect(data=bands, aes(xmin=x1, xmax=x2, ymin=-Inf, ymax=Inf), inherit.aes=FALSE, fill='grey60', alpha=.3) +
  autolayer(test, series='actual (held out)') +
  geom_hline(yintercept=mean(train), linetype='dotted', colour='grey40') +
  scale_colour_manual(values=c('actual (held out)'='red')) +
  labs(title='auto.arima forecast of US GNP growth (NBER recessions shaded)', x='year', y='growth %/qtr', colour='') +
  theme_minimal(base_size=13) + theme(legend.position='top')
No description has been provided for this image

Business-cycle asymmetry — why a linear ARMA isn't the whole story¶

Recessions are deeper and sharper, not merely lower: split by the NBER flag, recession growth averages well below expansion, and the overall distribution is left-skewed — rare, abrupt contractions a symmetric ARMA cannot represent. This is the nonlinearity behind the Markov-switching companions.

In [3]:
# business-cycle asymmetry: recessions are deeper (left skew) -- a symmetric ARMA can't capture it
options(repr.plot.width = 11, repr.plot.height = 4.5)
g <- d$gnp_growth; rec <- d$recession == 1
skewness <- function(x) { m <- mean(x); mean((x - m)^3) / sd(x)^3 }
par(mfrow = c(1, 2))
brks <- seq(min(g), max(g), length.out = 22)
hist(g[!rec], breaks = brks, col = rgb(.27, .51, .71, .55), border = "white", freq = FALSE, xlim = range(g),
     xlab = "GNP growth (%/qtr)",
     main = sprintf("Recession vs expansion\n(exp mean %.2f, rec mean %.2f)", mean(g[!rec]), mean(g[rec])))
hist(g[rec], breaks = brks, col = rgb(.70, .13, .13, .55), border = "white", freq = FALSE, add = TRUE)
abline(v = 0, lwd = .7)
legend("topright", c("expansion", "recession"),
       fill = c(rgb(.27, .51, .71, .55), rgb(.70, .13, .13, .55)), bty = "n", cex = .8)
hist(g, breaks = 22, col = rgb(.42, .35, .80, .5), border = "white", freq = FALSE,
     xlab = "GNP growth (%/qtr)", main = sprintf("Left-skew %.2f: rare sharp contractions", skewness(g)))
curve(dnorm(x, mean(g), sd(g)), add = TRUE, lty = 2, lwd = 1.6)
par(mfrow = c(1, 1))
No description has been provided for this image

Spectral view — cyclicality as a broad band, not a line¶

The smoothed GNP-growth spectrum is broad and nearly flat (near white noise, echoing the series' weak autocorrelation) with only a gentle low-frequency tilt over the business-cycle band (6–32 quarters) — and no sharp line, the frequency-domain signature of cyclicality (contrast the seasonal combs of AirPassengers / CO₂).

In [4]:
# spectral view: GNP growth near white noise with a gentle business-cycle tilt (no sharp line)
options(repr.plot.width = 8.5, repr.plot.height = 5)
tg <- ts(as.numeric(train), frequency = 1)                 # frequency in cycles/quarter (0..0.5)
sp <- spectrum(tg, spans = c(3, 5), taper = .1, log = "no", plot = FALSE)
plot(sp$freq, sp$spec, type = "l", col = "goldenrod", lwd = 2,
     xlab = "frequency (cycles/quarter)", ylab = "spectral power",
     main = "GNP-growth spectrum: broad, no sharp line (near white noise)")
rect(1/32, 0, 1/6, max(sp$spec) * 1.05, col = rgb(.85, .65, .13, .15), border = NA)
lines(sp$freq, sp$spec, col = "goldenrod", lwd = 2)
abline(v = c(1/32, 1/6), lty = 3, col = "grey50")
legend("topright", c("smoothed periodogram", "business-cycle band (6-32 qtr)"),
       col = c("goldenrod", "grey50"), lty = c(1, 3), lwd = c(2, 1), bty = "n", cex = .8)
No description has been provided for this image

Results¶

  • auto.arima selects a low-order ARMA (here ARMA(2,2); pmdarima/statsmodels preferred a simpler AR(1)). The order is unstable — unsurprising when there is so little signal — but the forecast is the same.
  • Out of sample it fails to beat the naive mean — essentially tied (RMSE 1.21 vs 1.20), the forecast collapsing to the average growth rate within a couple of quarters — all three engines reach the same conclusion (despite picking different ARMA orders): quarterly GNP growth is close to unforecastable.
  • The arc, complete. AirPassengers (fixed seasonality), sunspots (a damping pseudo-cycle), CO₂ (near-deterministic trend+seasonal), and GNP growth (weak, irregular business cycle) span the full spectrum of ARIMA predictability. The GNP case — symmetric ARMA missing the sharp recessions — is exactly what the Markov-switching and (next) GARCH / stochastic-volatility models are built to capture.