CO₂ at Mauna Loa — Bayesian structural model, out-of-sample¶

Companion to co2_python.ipynb. A Bayesian structural regression: a quadratic trend (CO₂ is slightly accelerating) + monthly seasonal effects + AR(1) errors, on the raw (additive) scale. $$y_t = \alpha + \beta_1 t + \beta_2 t^2 + \sum_{m=2}^{12} s_m\mathbb{1}[\text{month}=m] + e_t,\quad e_t=\rho e_{t-1}+\varepsilon_t.$$

In [1]:
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('co2.csv', parse_dates=['month']); y = d['co2'].to_numpy(float); month = d['month'].dt.month.to_numpy()
n = len(y); H = 60; ntr = n - H
t = np.arange(n); tm, ts = t[:ntr].mean(), t[:ntr].std(); z = (t-tm)/ts
def design(idx):
    return np.column_stack([np.ones(len(idx)), z[idx], z[idx]**2] + [(month[idx]==m).astype(float) for m in range(2,13)])
Xtr, Xfu = design(np.arange(ntr)), design(np.arange(ntr, n)); k = Xtr.shape[1]; ytr = y[:ntr]
with pm.Model() as m:
    beta = pm.Normal('beta', 0, 50, shape=k); rho = pm.Uniform('rho', -0.99, 0.99); sigma = pm.HalfNormal('sigma', 5)
    mu = pt.dot(Xtr, beta)
    pm.Normal('y0', mu[0], sigma/pt.sqrt(1-rho**2), observed=ytr[0])
    pm.Normal('yt', mu[1:] + rho*(ytr[:-1]-mu[:-1]), sigma, observed=ytr[1:])
    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 7 seconds.

rho 0.919 sigma 0.311 max r_hat 1.000
In [2]:
B = idata.posterior['beta'].values.reshape(-1,k); R = idata.posterior['rho'].values.ravel(); S = idata.posterior['sigma'].values.ravel()
rng = np.random.default_rng(1); nd = len(R)
mu_tr = Xtr@B.T; mu_fu = Xfu@B.T; e = ytr[-1]-mu_tr[-1]; sims = np.zeros((H, nd))
for h in range(H):
    e = R*e + S*rng.standard_normal(nd); sims[h] = mu_fu[h] + e
fm = sims.mean(1); lo = np.percentile(sims,2.5,1); hi = np.percentile(sims,97.5,1); ytest = y[ntr:]
rmse = np.sqrt(np.mean((fm-ytest)**2)); mape = np.mean(np.abs(fm/ytest-1))*100
print('OUT-OF-SAMPLE (5 yr):  RMSE = %.2f ppm   MAPE = %.3f%%' % (rmse, mape))
fig, ax = plt.subplots(figsize=(11, 4.6)); idx = d['month']
ax.plot(idx[ntr-120:ntr], y[ntr-120:ntr], color='steelblue', lw=1, label='train (zoom)')
ax.plot(idx[ntr:], ytest, 'o-', color='black', ms=2.5, lw=.8, label='actual (held out)')
ax.plot(idx[ntr:], fm, '--', color='firebrick', lw=1.6, label='forecast')
ax.fill_between(idx[ntr:], lo, hi, color='firebrick', alpha=.25, label='95% predictive')
ax.set_xlabel('year'); ax.set_ylabel('CO₂ (ppm)'); ax.set_title('Bayesian structural forecast (RMSE %.2f ppm, MAPE %.3f%%)'%(rmse,mape)); ax.legend(fontsize=8)
plt.tight_layout(); plt.savefig('co2_pymc.png', dpi=120, bbox_inches='tight'); plt.show()

OUT-OF-SAMPLE (5 yr):  RMSE = 1.07 ppm   MAPE = 0.245%
No description has been provided for this image

Results¶

  • The quadratic-trend + seasonal + AR(1) model matches the SARIMA: MAPE ≈ 0.25% out of sample, with a full posterior-predictive band. The quadratic term captures the gentle acceleration of the Keeling curve (a purely linear trend would under-forecast).
  • As with AirPassengers, the seasonality is encoded as fixed monthly effects (deterministic) — appropriate because the CO₂ annual cycle is extremely stable.