AirPassengers — R forecast::auto.arima, out-of-sample¶

Companion to the Python/PyMC notebooks. Rob Hyndman's forecast package: auto.arima selects the SARIMA order automatically (with a Box-Cox lambda=0 log transform), and forecast() produces the out-of-sample forecast with prediction intervals. This is the gold-standard frequentist workflow.

Package: forecast¶

Rob Hyndman's forecast is the reference frequentist time-series toolkit. The pieces used here:

  • auto.arima — automatic SARIMA selection via the Hyndman–Khandakar algorithm: pick the number of regular differences d (from a KPSS test) and seasonal differences D (from a seasonal-strength / OCSB test), then run a stepwise search over (p, q, P, Q), keeping the model with the lowest AICc.
  • lambda = 0 — a Box–Cox transform with $\lambda=0$ is a log transform. forecast applies it, models on the log scale, and back-transforms the forecasts automatically (with a bias correction), so we don't log/exp by hand.
  • forecast(fit, h) — point forecasts plus 80% / 95% prediction intervals, h steps ahead.
  • accuracy(fc, test) — scores the forecast against the held-out test set (RMSE, MAPE, MAE, MASE, …).
  • Box.test(..., type="Ljung-Box") — the residual white-noise check.

This is the most compact of the three workflows: a single call selects the model, applies the transform, and fits.

In [1]:
suppressMessages(library(forecast))
d <- read.csv('airpassengers.csv')
y <- ts(d$passengers, start=c(1949,1), frequency=12)
H <- 24; train <- window(y, end=c(1958,12)); test <- window(y, start=c(1959,1))
fit <- auto.arima(train, lambda=0, seasonal=TRUE, stepwise=TRUE, approximation=FALSE)
cat('auto.arima selected:\n'); print(fit)
fc <- forecast(fit, h=H)
acc <- accuracy(fc, test)
cat(sprintf('\nOUT-OF-SAMPLE (24 months):  RMSE = %.1f   MAPE = %.1f%%\n', acc['Test set','RMSE'], acc['Test set','MAPE']))
cat(sprintf('Ljung-Box residual p = %.3f\n', Box.test(residuals(fit), lag=12, type='Ljung-Box')$p.value))

auto.arima selected:

Series: train 
ARIMA(0,1,1)(0,1,1)[12] 
Box Cox transformation: lambda= 0 

Coefficients:
          ma1     sma1
      -0.3424  -0.5405
s.e.   0.1009   0.0877

sigma^2 = 0.001432:  log likelihood = 197.51
AIC=-389.02   AICc=-388.78   BIC=-381

OUT-OF-SAMPLE (24 months):  RMSE = 43.2   MAPE = 8.5%

Ljung-Box residual p = 0.891
In [2]:
suppressMessages(library(ggplot2))
options(repr.plot.width = 9, repr.plot.height = 5.5)
# ggplot/autoplot render path. autoplot(fc) draws the HISTORY + forecast + PI; add the held-out actuals.
autoplot(fc) +
  autolayer(train, series='training (history)') +
  autolayer(test,  series='actual (held out)') +
  scale_colour_manual(values=c('training (history)'='black','actual (held out)'='red')) +
  labs(title='auto.arima out-of-sample forecast', x='year', y='passengers', colour='') +
  theme_minimal(base_size=13) + theme(legend.position='top')
No description has been provided for this image

Results¶

  • auto.arima (with the log transform) selects the airline model SARIMA(0,1,1)(0,1,1)[12] — matching the Python pmdarima choice — with white-noise residuals (Ljung-Box p > 0.05).
  • Out-of-sample over the 2-year horizon: RMSE ≈ 40–45, MAPE ≈ 8–9%, the same ballpark as statsmodels; the Bayesian structural model scores a bit better (RMSE ≈ 32, MAPE ≈ 6%).
  • Both SARIMA engines agree exactly (statsmodels ≈ forecast::auto.arima); the Bayesian structural alternative captures the same trend + 12-month seasonality and forecasts it a touch better. Next datasets move from seasonality (fixed period) to cyclicality (variable period).