AirPassengers — seasonal ARIMA (Box–Jenkins), with out-of-sample forecasting

statsmodels (frequentist) — the canonical seasonal series

The classic monthly international airline passengers, 1949–1960 (Box & Jenkins). It has a clear upward trend and multiplicative seasonality (the seasonal swings grow with the level) — the textbook case for a seasonal ARIMA (SARIMA) on the log scale. We run the full Box–Jenkins workflow and then forecast the held-out final 2 years to see how it does out of sample.

Packages used

• statsmodels (tsa) — the core engine. SARIMAX fits a seasonal ARIMA (optionally with exogenous regressors) inside a state-space framework, evaluating the exact likelihood with the Kalman filter; get_forecast runs that filter forward to produce forecasts with standard errors and prediction intervals. It also supplies the diagnostic toolkit: adfuller, kpss, STL, acorr_ljungbox.
• pmdarima — auto_arima, the Python port of R's auto.arima. It automates the Box–Jenkins identification step: it searches over ARIMA orders and keeps the model with the lowest AIC (the stepwise Hyndman–Khandakar search), choosing the differencing orders d/D from unit-root tests.

The SARIMA model

A seasonal ARIMA is written SARIMA(p, d, q)(P, D, Q)[s]:

symbol	part	meaning
p / P	AR / seasonal AR	dependence on the last p lags / last P seasonal lags
d / D	differencing / seasonal differencing	d regular differences remove trend; D seasonal differences remove seasonality
q / Q	MA / seasonal MA	dependence on the last q shocks / last Q seasonal shocks
s	season length	12 for monthly data

In backshift-operator form ($B\,y_t = y_{t-1}$): $$\phi_p(B)\,\Phi_P(B^s)\,(1-B)^d(1-B^s)^D\,y_t \;=\; \theta_q(B)\,\Theta_Q(B^s)\,\varepsilon_t .$$ The "airline model" SARIMA(0,1,1)(0,1,1)[12] — the one Box & Jenkins originally fit to this series — is the special case $$(1-B)(1-B^{12})\,y_t = (1+\theta B)(1+\Theta B^{12})\,\varepsilon_t,$$ i.e. one regular difference (removes the trend), one seasonal difference (removes the 12-month pattern), and an MA(1) × seasonal-MA(1) for the short-run and year-over-year shocks — no AR terms at all.

Box–Jenkins, in three steps

1. Identification — stabilise the variance (here: take logs), make the series stationary (choose d, D from ADF/KPSS), and read tentative orders from the ACF/PACF (automated below by auto_arima). Note ADF and KPSS have opposite nulls — ADF's $H_0$ is a unit root (non-stationary), KPSS's $H_0$ is stationarity — so we want ADF p small and KPSS p large.
2. Estimation — fit by maximum likelihood via the Kalman filter.
3. Diagnostics — residuals should be indistinguishable from white noise (Ljung-Box p > 0.05); if not, revise the orders.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.seasonal import STL
from statsmodels.tsa.stattools import adfuller, kpss
from statsmodels.stats.diagnostic import acorr_ljungbox
import pmdarima as pm

d = pd.read_csv('airpassengers.csv', parse_dates=['month'])
y = d.set_index('month')['passengers'].asfreq('MS')
ly = np.log(y)                                            # log -> turns multiplicative seasonality additive
H = 24                                                    # hold out the last 2 years for out-of-sample test
train, test = ly.iloc[:-H], ly.iloc[-H:]
print('total %d months; train %d (%s..%s), test %d (%s..%s)' %
      (len(y), len(train), train.index[0].date(), train.index[-1].date(), len(test), test.index[0].date(), test.index[-1].date()))
# stationarity (on log, log-diff, log-seasonal-diff)
def stat(s, lab):
    s = s.dropna(); ad = adfuller(s)[1]; kp = kpss(s, nlags='auto')[1]
    print('  %-22s ADF p=%.3f  KPSS p=%.3f' % (lab, ad, kp))
print('stationarity tests (ADF H0=unit root; KPSS H0=stationary):')
stat(ly, 'log'); stat(ly.diff(), 'diff'); stat(ly.diff().diff(12), 'diff + seasonal-diff')

total 144 months; train 120 (1949-01-01..1958-12-01), test 24 (1959-01-01..1960-12-01)
stationarity tests (ADF H0=unit root; KPSS H0=stationary):
  log                    ADF p=0.422  KPSS p=0.010
  diff                   ADF p=0.071  KPSS p=0.100
  diff + seasonal-diff   ADF p=0.000  KPSS p=0.100

C:\Users\user\AppData\Local\Temp\ipykernel_28688\1108749064.py:17: InterpolationWarning: The test statistic is outside of the range of p-values available in the
look-up table. The actual p-value is smaller than the p-value returned.

  s = s.dropna(); ad = adfuller(s)[1]; kp = kpss(s, nlags='auto')[1]
C:\Users\user\AppData\Local\Temp\ipykernel_28688\1108749064.py:17: InterpolationWarning: The test statistic is outside of the range of p-values available in the
look-up table. The actual p-value is greater than the p-value returned.

  s = s.dropna(); ad = adfuller(s)[1]; kp = kpss(s, nlags='auto')[1]
C:\Users\user\AppData\Local\Temp\ipykernel_28688\1108749064.py:17: InterpolationWarning: The test statistic is outside of the range of p-values available in the
look-up table. The actual p-value is greater than the p-value returned.

  s = s.dropna(); ad = adfuller(s)[1]; kp = kpss(s, nlags='auto')[1]

Order selection — with a caveat

We let auto_arima search the AR/MA orders but fix d=1, D=1, the differencing the stationarity tests above call for. Left to choose d on its own, auto_arima here picks d=0 — leaving the trend in, so the residuals stay autocorrelated and the forecast degrades. A good reminder that automated selectors still need the analyst's judgement.

# automatic SARIMA order selection -- force d=1, D=1 (the differencing the stationarity tests support)
am = pm.auto_arima(train, seasonal=True, m=12, d=1, D=1, stepwise=True,
                   suppress_warnings=True, error_action='ignore', trace=False)
print('auto_arima chose SARIMA%s x %s[12]   AIC %.1f' % (am.order, am.seasonal_order, am.aic()))
mod = SARIMAX(train, order=am.order, seasonal_order=am.seasonal_order,
              enforce_stationarity=False, enforce_invertibility=False).fit(disp=False)
lb = acorr_ljungbox(mod.resid[13:], lags=[12], return_df=True)['lb_pvalue'].iloc[0]
print('residual Ljung-Box(12) p = %.3f  (>0.05 => residuals look like white noise)' % lb)

auto_arima chose SARIMA(0, 1, 1) x (0, 1, 1, 12)[12]   AIC -389.0
residual Ljung-Box(12) p = 0.900  (>0.05 => residuals look like white noise)

Forecasting out of sample

get_forecast(H) runs the Kalman filter forward H steps, returning the predicted mean and the forecast variance (which grows with the horizon) — hence the widening prediction band. We exponentiate back to passenger counts and score against the held-out actuals with RMSE and MAPE.

# OUT-OF-SAMPLE forecast of the held-out 2 years (back-transformed to passenger counts)
fc = mod.get_forecast(H); fm = np.exp(fc.predicted_mean); ci = np.exp(fc.conf_int(alpha=0.05))
yact = np.exp(test)
rmse = np.sqrt(np.mean((fm.values - yact.values)**2)); mape = np.mean(np.abs(fm.values/yact.values - 1))*100
print('OUT-OF-SAMPLE (24 months):  RMSE = %.1f passengers   MAPE = %.1f%%' % (rmse, mape))

fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
# (A) STL decomposition of the log series
stl = STL(ly, period=12).fit()
ax[0].plot(ly.index, ly, color='steelblue', lw=.8, label='log passengers')
ax[0].plot(ly.index, stl.trend, color='firebrick', lw=1.6, label='STL trend')
ax[0].set_title('Log series + STL trend (seasonality removed)'); ax[0].legend(fontsize=8)
# (B) out-of-sample forecast on the original scale
ax[1].plot(np.exp(train).index, np.exp(train), color='steelblue', lw=1, label='train')
ax[1].plot(yact.index, yact, 'o-', color='black', ms=3, lw=1, label='actual (held out)')
ax[1].plot(fm.index, fm, '--', color='firebrick', lw=1.6, label='forecast')
ax[1].fill_between(ci.index, ci.iloc[:,0], ci.iloc[:,1], color='firebrick', alpha=.2, label='95% PI')
ax[1].set_title('Out-of-sample forecast  (RMSE %.0f, MAPE %.1f%%)' % (rmse, mape)); ax[1].set_ylabel('passengers'); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('airpassengers_py.png', dpi=120, bbox_inches='tight'); plt.show()

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

Spectral view — seasonality as sharp spectral lines

The frequency-domain complement to the SARIMA fit. The periodogram of the (linearly detrended) log series shows a comb of sharp lines at the seasonal harmonics $f=k/12$ cycles per month — a dominant annual line (12-month) plus 6-, 4-, 3-month overtones. Power piling up at exact frequencies is the fingerprint of deterministic seasonality, exactly what the seasonal SARIMA terms $(P,D,Q)_{12}$ encode. It is the sharp-line counterpart to the sunspots broad bump (a fixed period vs a stochastic cycle), and the same lines motivate the Fourier/harmonic-regression alternative to seasonal differencing.

# spectral view: seasonality = SHARP LINES at k/12 (contrast the sunspots broad cyclic bump)
from scipy.signal import periodogram, welch
xtr = train.values
f, Pxx = periodogram(xtr, fs=1.0, detrend="linear")
fw, Pw = welch(xtr, fs=1.0, nperseg=48, detrend="linear")
fig, ax = plt.subplots(figsize=(8.6, 4.7))
for k in range(1, 7): ax.axvline(k/12, color="firebrick", ls="--", lw=1.1, alpha=.7)
ax.plot([], [], "--", color="firebrick", lw=1.1, label="seasonal harmonics k/12")
ax.semilogy(f, Pxx, color="0.55", lw=.9, label="periodogram (log, detrended)")
ax.semilogy(fw, Pw, color="goldenrod", lw=1.4, label="Welch (smoothed)")
ax.set_xlim(0, 0.5); ax.set_ylim(1e-6, None)
ax.set_xlabel("frequency (cycles / month)"); ax.set_ylabel("power (log scale)")
ax.set_title("AirPassengers spectrum: sharp LINES at k/12 = deterministic annual SEASONALITY\n(contrast the sunspots broad bump = a stochastic cycle)")
ax.legend(fontsize=8, loc="upper right")
sec = ax.secondary_xaxis("top", functions=(lambda x: 1/np.maximum(x, 1e-9), lambda p: 1/np.maximum(p, 1e-9)))
sec.set_xlabel("period (months)"); sec.set_xticks([2, 3, 4, 6, 12])
plt.tight_layout(); plt.savefig("airpassengers_spectrum.png", dpi=120, bbox_inches="tight"); plt.show()

Results

• The log transform is essential: raw seasonal swings grow with the level (multiplicative); on the log scale they're constant, so an additive SARIMA fits. Stationarity tests confirm a regular difference (d=1) and a seasonal difference (D=1) are needed.
• auto_arima selects an airline-type SARIMA(p,1,q)(P,1,Q)[12]; residuals pass the Ljung-Box test (white noise), so the model has captured the trend + seasonal structure.
• Out-of-sample: forecasting the held-out final 2 years (a long, 24-step horizon) gives MAPE ≈ 8–9% — the SARIMA tracks both the rising trend and the summer-peak seasonality, with prediction intervals that widen sensibly into the future. (At shorter horizons the error is much smaller; 24-month-ahead is a deliberately demanding test.)
• Seasonality vs cyclicality: here the period is fixed (12 months) — pure seasonality, the easy case. The sunspots/lynx notebooks next show cyclicality (variable period), which SARIMA handles differently (via AR roots, not seasonal differencing).
• Cross-checks: PyMC (Bayesian trend+seasonal+AR) and R (forecast::auto.arima) follow.