Sunspots — Bayesian AR, with a posterior over the cycle period

Companion to sunspots_python.ipynb. A Bayesian AR(p) on the level (no differencing — the series is stationary). The Bayesian payoff: the complex AR roots imply a cycle period, so a posterior over the AR coefficients gives a posterior distribution of the ~11-year cycle length — the uncertainty in the cycle, not just a point value. Forecasts come as posterior-predictive draws (and damp, like any stochastic cycle).

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('sunspots.csv'); y = d['sunspots'].to_numpy(float); yr = d['year'].to_numpy()
H = 40; ytr = y[:-H]; n = len(ytr); p = 9
L = np.column_stack([ytr[p-j: n-j] for j in range(1, p+1)])         # lag design: row t -> [y_{t-1..t-p}]
target = ytr[p:]
with pm.Model() as m:
    mu = pm.Normal('mu', ytr.mean(), 50); phi = pm.Normal('phi', 0, 0.5, shape=p); sigma = pm.HalfNormal('sigma', 50)
    mean = mu*(1 - pt.sum(phi)) + pt.dot(L, phi)
    pm.Normal('obs', mean, sigma, observed=target)
    idata = pm.sample(1500, tune=1500, chains=4, target_accept=0.9, random_seed=1, progressbar=False)
PHI = idata.posterior['phi'].values.reshape(-1, p); MU = idata.posterior['mu'].values.ravel(); SIG = idata.posterior['sigma'].values.ravel()
print('max r_hat %.3f' % float(az.summary(idata, var_names=['phi','mu','sigma'])['r_hat'].max()))

def dom_period(phi):                                                # dominant complex eigenvalue of the companion matrix
    C = np.zeros((p, p)); C[0, :] = phi; C[1:, :-1] = np.eye(p-1)
    ev = np.linalg.eigvals(C); c = ev[np.abs(ev.imag) > 1e-6]
    if len(c) == 0: return np.nan
    dom = c[np.argmax(np.abs(c))]; return 2*np.pi/np.abs(np.angle(dom))
periods = np.array([dom_period(PHI[i]) for i in range(0, len(PHI), 4)]); periods = periods[np.isfinite(periods)]
print('posterior cycle PERIOD: median %.1f yr, 95%% CI [%.1f, %.1f]' %
      (np.median(periods), np.percentile(periods,2.5), np.percentile(periods,97.5)))

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [mu, phi, sigma]

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

There were 23 divergences after tuning. Increase `target_accept` or reparameterize.

max r_hat 1.000
posterior cycle PERIOD: median 10.5 yr, 95% CI [10.1, 10.9]

# posterior-predictive out-of-sample forecast (recursive; the cycle damps)
rng = np.random.default_rng(1); nd = len(PHI)
hist = np.tile(ytr[-p:], (nd, 1))                                   # (nd, p) last p obs, newest last
sims = np.zeros((H, nd))
for h in range(H):
    lag = hist[:, ::-1]                                             # newest first -> y_{t-1..t-p}
    mean = MU*(1 - PHI.sum(1)) + np.einsum('dj,dj->d', lag, PHI)
    nxt = mean + SIG*rng.standard_normal(nd)
    sims[h] = nxt; hist = np.column_stack([hist[:, 1:], nxt])
fm = sims.mean(1); lo = np.percentile(sims,2.5,1); hi = np.percentile(sims,97.5,1)
ytest = y[-H:]; rmse = np.sqrt(np.mean((fm-ytest)**2)); print('OUT-OF-SAMPLE RMSE = %.1f' % rmse)

fig, ax = plt.subplots(1, 2, figsize=(13, 4.4))
ax[0].hist(periods, bins=30, color='steelblue', edgecolor='k', lw=.3, density=True)
ax[0].axvline(np.median(periods), color='firebrick', lw=2, label='median %.1f yr'%np.median(periods))
ax[0].set_xlabel('implied cycle period (years)'); ax[0].set_yticks([]); ax[0].set_title('Posterior of the solar-cycle period'); ax[0].legend()
ax[1].plot(yr[:-H], ytr, color='steelblue', lw=.8, label='train')
ax[1].plot(yr[-H:], ytest, 'o-', color='black', ms=3, lw=.8, label='actual (held out)')
ax[1].plot(yr[-H:], fm, '--', color='firebrick', lw=1.6, label='forecast')
ax[1].fill_between(yr[-H:], lo, hi, color='firebrick', alpha=.2, label='95% predictive')
ax[1].axhline(ytr.mean(), color='gray', ls=':', lw=1); ax[1].set_xlabel('year'); ax[1].set_ylabel('sunspots')
ax[1].set_title('Bayesian AR(%d) forecast (RMSE %.0f) — damps to mean'%(p,rmse)); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('sunspots_pymc.png', dpi=120, bbox_inches='tight'); plt.show()

OUT-OF-SAMPLE RMSE = 35.3

Results

• The Bayesian AR(9) agrees with statsmodels and yields a posterior cycle period centred near ~10–11 years — a genuine credible interval on the solar-cycle length, extracted from the complex roots of each posterior draw.
• The posterior-predictive forecast damps toward the mean as the horizon grows, with the band widening to the unconditional spread — exactly the stochastic-cycle behaviour seen in the frequentist fit, now with honest predictive uncertainty.
• The sampler flagged 23 divergences (of 6,000 draws; r̂ = 1.00) from the awkward AR-on-levels geometry — a small fraction, but enough that the very tight cycle-period interval [10.1, 10.9] yr is best read as slightly optimistic.