Nile river flow — Bayesian changepoint (estimating the break year)¶

Companion to nile_python.ipynb. Here the break date is not assumed — it is a parameter. A classic Bayesian changepoint model: $$\tau\sim\text{DiscreteUniform},\quad y_t\sim N(\mu_1,\sigma^2)\ \text{if } t<\tau,\ \ N(\mu_2,\sigma^2)\ \text{if } t\ge\tau.$$ NUTS samples $\mu_1,\mu_2,\sigma$; a Metropolis step samples the discrete $\tau$. The posterior of $\tau$ tells us when the regime changed (and how sure we are), and $\mu_2-\mu_1$ gives the size of the drop with full uncertainty.

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('nile.csv'); y = d['flow'].to_numpy(float); year = d['year'].to_numpy()
H = 20; ntr = len(y)-H; ytr = y[:ntr]; idx = np.arange(ntr)
with pm.Model() as m:
    tau = pm.DiscreteUniform('tau', lower=5, upper=ntr-5)
    mu1 = pm.Normal('mu1', 1000, 300); mu2 = pm.Normal('mu2', 1000, 300); sigma = pm.HalfNormal('sigma', 200)
    mu = pm.math.switch(idx >= tau, mu2, mu1)
    pm.Normal('obs', mu, sigma, observed=ytr)
    idata = pm.sample(2000, tune=2000, chains=4, target_accept=0.9, random_seed=1, progressbar=False)
TAU = idata.posterior['tau'].values.ravel(); yearpost = year[0] + TAU
M1 = idata.posterior['mu1'].values.ravel(); M2 = idata.posterior['mu2'].values.ravel()
from scipy import stats
print('break YEAR posterior: mode %d, 95%% CI [%d, %d]' % (int(stats.mode(yearpost, keepdims=True)[0][0]), np.percentile(yearpost,2.5), np.percentile(yearpost,97.5)))
print('level drop mu2-mu1 = %.0f  [%.0f, %.0f]' % ((M2-M1).mean(), np.percentile(M2-M1,2.5), np.percentile(M2-M1,97.5)))

Multiprocess sampling (4 chains in 4 jobs)

CompoundStep

>Metropolis: [tau]

>NUTS: [mu1, mu2, sigma]

Sampling 4 chains for 2_000 tune and 2_000 draw iterations (8_000 + 8_000 draws total) took 7 seconds.

break YEAR posterior: mode 1899, 95% CI [1897, 1900]
level drop mu2-mu1 = -255  [-317, -194]
In [2]:
fig, ax = plt.subplots(1, 2, figsize=(13, 4.4))
vals, cnts = np.unique(yearpost, return_counts=True)
ax[0].bar(vals, cnts/cnts.sum(), color='steelblue', width=.9)
ax[0].axvline(1899, color='firebrick', ls='--', lw=1.5, label='Aswan dam ~1899')
ax[0].set_xlabel('estimated break year τ'); ax[0].set_ylabel('posterior prob'); ax[0].set_title('Posterior of the break YEAR'); ax[0].legend(fontsize=8)
ax[0].set_xlim(1890, 1910)
# fitted two-level mean (posterior) + held-out actuals
mode_tau = int(stats.mode(TAU, keepdims=True)[0][0]); lvl = np.where(np.arange(len(y)) >= mode_tau, M2.mean(), M1.mean())
ax[1].plot(year, y, color='steelblue', lw=1, marker='.', ms=3, label='flow')
ax[1].plot(year, lvl, color='firebrick', lw=2, label='posterior level (μ₁→μ₂)')
ax[1].axvline(year[ntr], color='gray', ls=':', label='train/test split')
ax[1].set_xlabel('year'); ax[1].set_ylabel('flow'); ax[1].set_title('Two-regime mean (drop %.0f units)'%((M1-M2).mean())); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('nile_pymc.png', dpi=120, bbox_inches='tight'); plt.show()
print('OUT-OF-SAMPLE RMSE (forecast = post-break mean): %.0f' % np.sqrt(np.mean((y[ntr:]-M2.mean())**2)))
No description has been provided for this image 
OUT-OF-SAMPLE RMSE (forecast = post-break mean): 127

Results¶

  • The data find the break. The posterior of τ concentrates tightly around 1899 — recovering the Aswan-dam date without being told it. This is the Bayesian advantage over a fixed step dummy: the changepoint and its uncertainty are estimated.
  • The level drop μ₂−μ₁ ≈ −250 units with a credible interval, matching the frequentist intervention estimate.
  • Forecasting the held-out years at the post-break mean μ₂ gives an RMSE in line with the statsmodels intervention model — and a far more honest model than pretending the series has a stochastic trend.

Connection — a changepoint is "Markov-switching with a single switch"¶

This model infers one permanent latent switch (the break year τ). Let the latent state instead flip back and forth over time — governed by a Markov chain (a transition matrix) — and you get Markov-switching / regime-switching, applied to the same Hamilton GNP series in:

  • markov_switching_python.ipynb / bayesian_msar_python.ipynb

It's the same latent-discrete-state machinery; a changepoint is the special absorbing case (once it switches it never returns). See gnp_python.ipynb for the linear benchmark and the full ARIMA → changepoint → Markov-switching ladder.