AirPassengers — Bayesian time series (trend + seasonal + AR errors), out-of-sample

Companion to airpassengers_python.ipynb. A SARIMA is one way to handle trend + seasonality; the Bayesian route here is an interpretable structural regression: on the log scale, $$\log y_t = \alpha + \beta\,t + \sum_{m} s_m\,\mathbb{1}[\text{month}=m] + e_t,\qquad e_t = \rho\,e_{t-1} + \varepsilon_t,$$ a linear trend + monthly seasonal effects + AR(1) errors (the autocorrelation a plain regression would miss). NUTS gives full posterior forecast distributions for the held-out 2 years.

The model in detail

This is a structural / regression-with-ARMA-errors model — an interpretable alternative to SARIMA. On the log scale, $$\log y_t = \underbrace{\alpha + \beta\,t}_{\text{linear trend}} \;+\; \underbrace{\textstyle\sum_{m=2}^{12} s_m\,\mathbb{1}[\text{month}=m]}_{\text{11 monthly seasonal effects}} \;+\; e_t,\qquad e_t=\rho\,e_{t-1}+\varepsilon_t,\;\;\varepsilon_t\sim N(0,\sigma^2).$$

• Trend — a deterministic linear slope on (standardised) time.
• Seasonality — 11 monthly dummy coefficients (January is the baseline); these are fixed (deterministic) seasonal effects.
• AR(1) errors — the regression residuals are serially correlated; a plain OLS fit would have autocorrelated errors and therefore wrong standard errors and intervals. The AR(1) term absorbs that correlation, and $\rho$ measures the month-to-month error persistence.

Priors. beta ~ N(0, 5) (weakly informative on the log scale), rho ~ Uniform(-0.99, 0.99) (enforces a stationary AR(1)), sigma ~ HalfNormal(1).

Likelihood (conditional AR(1)). For $t\ge 1$ we use $y_t\mid y_{t-1}\sim N\!\big(\mu_t+\rho(y_{t-1}-\mu_{t-1}),\,\sigma^2\big)$, where $\mu_t$ is the trend+seasonal mean; the first observation gets the stationary marginal $N\!\big(\mu_0,\,\sigma^2/(1-\rho^2)\big)$. Plugging in the observed lag makes this fully vectorised (no slow scan loop).

Package

PyMC builds the model symbolically and draws the posterior with NUTS (the No-U-Turn Hamiltonian Monte Carlo sampler); arviz checks convergence (r̂ ≈ 1.00 means the four chains agree).

Deterministic vs stochastic seasonality

SARIMA removes seasonality by differencing ($1-B^{12}$), which lets the seasonal pattern drift over time — stochastic seasonality. Here the seasonal effects are fixed dummies — deterministic seasonality, simpler and appropriate when the seasonal shape is stable (as it is for this series). Both are legitimate; they forecast similarly here.

import numpy as np, pandas as pd, pymc as pm, pytensor.tensor as pt, arviz as az, matplotlib.pyplot as plt
d = pd.read_csv('airpassengers.csv', parse_dates=['month'])
y = d['passengers'].to_numpy(float); ly = np.log(y); month = d['month'].dt.month.to_numpy()
n = len(y); H = 24; ntr = n - H
t = np.arange(n); tstd = (t - t[:ntr].mean()) / t[:ntr].std()
def design(idx):                                          # [1, t, 11 month dummies]
    D = np.column_stack([np.ones(len(idx)), tstd[idx]] + [(month[idx]==m).astype(float) for m in range(2,13)])
    return D
Xtr = design(np.arange(ntr)); Xfu = design(np.arange(ntr, n)); k = Xtr.shape[1]
ytr = ly[:ntr]

with pm.Model() as m:
    beta = pm.Normal('beta', 0, 5, shape=k); rho = pm.Uniform('rho', -0.99, 0.99); sigma = pm.HalfNormal('sigma', 1.0)
    mu = pt.dot(Xtr, beta)
    pm.Normal('y0', mu[0], sigma/pt.sqrt(1-rho**2), observed=ytr[0])                  # stationary start
    pm.Normal('yt', mu[1:] + rho*(ytr[:-1] - mu[:-1]), sigma, observed=ytr[1:])       # AR(1)-errors conditional
    idata = pm.sample(1500, tune=1500, chains=4, target_accept=0.95, random_seed=1, progressbar=False)
print('rho = %.3f   sigma = %.3f   max r_hat %.3f' %
      (float(idata.posterior['rho'].mean()), float(idata.posterior['sigma'].mean()),
       float(az.summary(idata, var_names=['beta','rho','sigma'])['r_hat'].max())))

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [beta, rho, sigma]

Sampling 4 chains for 1_500 tune and 1_500 draw iterations (6_000 + 6_000 draws total) took 6 seconds.

rho = 0.798   sigma = 0.036   max r_hat 1.000

Out-of-sample forecast (simulated)

There is no closed-form predictive distribution, so we simulate it. For each posterior draw we extrapolate the trend + seasonal mean over the next 24 months and propagate the AR(1) error forward, $e_{T+h}=\rho\,e_{T+h-1}+\varepsilon$, then exponentiate back to counts. Averaging over draws gives the posterior-mean forecast; the 2.5/97.5 percentiles give a genuine posterior-predictive interval that carries both parameter and innovation uncertainty (not a Gaussian asymptotic approximation).

# out-of-sample posterior-predictive forecast: extrapolate trend+seasonal, propagate AR(1) errors
post = idata.posterior; B = post['beta'].values.reshape(-1,k); R = post['rho'].values.ravel(); S = post['sigma'].values.ravel()
rng = np.random.default_rng(1); ndraw = len(R)
mu_tr_last = (Xtr @ B.T)                                  # (ntr, ndraw) fitted means on train
mu_fu = (Xfu @ B.T)                                       # (H, ndraw) future means
e_prev = ytr[-1] - mu_tr_last[-1]                         # last residual per draw
sims = np.zeros((H, ndraw))
for h in range(H):
    e_prev = R*e_prev + S*rng.standard_normal(ndraw)
    sims[h] = mu_fu[h] + e_prev
ypred = np.exp(sims)                                      # back to passenger counts (full predictive)
fm = ypred.mean(1); lo = np.percentile(ypred,2.5,1); hi = np.percentile(ypred,97.5,1)
yact = y[ntr:]
rmse = np.sqrt(np.mean((fm-yact)**2)); mape = np.mean(np.abs(fm/yact-1))*100
print('OUT-OF-SAMPLE (24 months):  RMSE = %.1f   MAPE = %.1f%%' % (rmse, mape))

fig, ax = plt.subplots(figsize=(11, 4.6)); idx = d['month']
ax.plot(idx[:ntr], y[:ntr], color='steelblue', lw=1, label='train')
ax.plot(idx[ntr:], yact, 'o-', color='black', ms=3, lw=1, label='actual (held out)')
ax.plot(idx[ntr:], fm, '--', color='firebrick', lw=1.6, label='posterior-mean forecast')
ax.fill_between(idx[ntr:], lo, hi, color='firebrick', alpha=.2, label='95% predictive')
ax.set_ylabel('passengers'); ax.set_title('Bayesian trend+seasonal+AR(1) forecast  (RMSE %.0f, MAPE %.1f%%)'%(rmse,mape)); ax.legend(fontsize=8)
plt.tight_layout(); plt.savefig('airpassengers_pymc.png', dpi=120, bbox_inches='tight'); plt.show()

OUT-OF-SAMPLE (24 months):  RMSE = 31.5   MAPE = 5.8%

Results

• The Bayesian structural model (linear trend + monthly seasonal effects + AR(1) errors, ρ ≈ 0.6–0.8) forecasts the held-out 2 years with accuracy comparable to the SARIMA, and gives a full predictive distribution — the 95% band comes from propagating both parameter uncertainty and the AR innovations.
• It encodes seasonality as fixed monthly effects (deterministic), whereas SARIMA uses seasonal differencing (stochastic seasonality) — two valid choices; here the seasonal pattern is stable, so both work well.
• Bayesian advantage: the forecast intervals are genuine posterior-predictive intervals (no Gaussian/asymptotic approximation), and every component (trend slope, each month, ρ) has an interpretable posterior.