US GNP growth — Bayesian AR(1), persistence & forecast¶

Companion to gnp_python.ipynb. A Bayesian AR(1) on quarterly growth: $$y_t = \mu + \phi\,(y_{t-1}-\mu) + \varepsilon_t,\qquad \varepsilon_t\sim N(0,\sigma^2).$$ The persistence $\phi$ has a direct business-cycle reading: the half-life of a growth shock is $\ln 0.5/\ln\phi$ quarters. The Bayesian fit gives a posterior for $\phi$ (and hence the half-life), plus a predictive forecast that — like the frequentist one — reverts quickly to the mean.

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('hamilton_gnp.csv', parse_dates=['date']).dropna(subset=['gnp_growth']).reset_index(drop=True)
y = d['gnp_growth'].to_numpy(float); dates = d['date']; rec = d['recession'].values
H = 16; ytr = y[:-H]; ytest = y[-H:]
with pm.Model() as m:
    mu = pm.Normal('mu', 0.7, 2); phi = pm.Uniform('phi', -0.99, 0.99); sigma = pm.HalfNormal('sigma', 2)
    pm.Normal('y0', mu, sigma/pt.sqrt(1-phi**2), observed=ytr[0])
    pm.Normal('yt', mu + phi*(ytr[:-1]-mu), sigma, observed=ytr[1:])
    idata = pm.sample(1500, tune=1500, chains=4, target_accept=0.9, random_seed=1, progressbar=False)
PHI = idata.posterior['phi'].values.ravel(); MU = idata.posterior['mu'].values.ravel(); SIG = idata.posterior['sigma'].values.ravel()
hl = np.log(0.5)/np.log(np.clip(PHI,1e-3,0.999))
print('phi  = %.2f  [%.2f, %.2f]' % (PHI.mean(), np.percentile(PHI,2.5), np.percentile(PHI,97.5)))
print('mean = %.2f%%/qtr   shock half-life = %.2f quarters (median)  max r_hat %.3f' %
      (MU.mean(), np.median(hl), float(az.summary(idata, var_names=['mu','phi','sigma'])['r_hat'].max())))

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 4 seconds.

phi  = 0.32  [0.14, 0.50]
mean = 0.77%/qtr   shock half-life = 0.61 quarters (median)  max r_hat 1.000
In [2]:
rng = np.random.default_rng(1); nd = len(PHI); prev = np.full(nd, ytr[-1]); sims = np.zeros((H, nd))
for h in range(H):
    prev = MU + PHI*(prev-MU) + SIG*rng.standard_normal(nd); sims[h] = prev
fm = sims.mean(1); lo = np.percentile(sims,2.5,1); hi = np.percentile(sims,97.5,1)
rmse = np.sqrt(np.mean((fm-ytest)**2)); print('OUT-OF-SAMPLE RMSE = %.2f (naive-mean %.2f)' % (rmse, np.sqrt(np.mean((ytest-ytr.mean())**2))))

fig, ax = plt.subplots(1, 2, figsize=(13, 4.3))
ax[0].hist(PHI, bins=40, color='steelblue', edgecolor='k', lw=.3, density=True)
ax[0].axvline(PHI.mean(), color='firebrick', lw=2, label='φ=%.2f (half-life %.1f qtr)'%(PHI.mean(), np.median(hl)))
ax[0].axvline(0, color='gray', ls=':'); ax[0].set_xlabel('AR(1) persistence φ'); ax[0].set_yticks([]); ax[0].set_title('Posterior persistence'); ax[0].legend(fontsize=8)
ax[1].plot(dates[:-H], ytr, color='steelblue', lw=.9, label='train')
ax[1].plot(dates[len(ytr):], ytest, 'o-', color='black', ms=3, lw=.8, label='actual (held out)')
ax[1].plot(dates[len(ytr):], fm, '--', color='firebrick', lw=1.8, label='forecast')
ax[1].fill_between(dates[len(ytr):], 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('growth %')
ax[1].set_title('Bayesian AR(1) forecast (RMSE %.2f) — reverts to mean'%rmse); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('gnp_pymc.png', dpi=120, bbox_inches='tight'); plt.show()

OUT-OF-SAMPLE RMSE = 1.19 (naive-mean 1.20)
No description has been provided for this image

Results¶

  • The posterior puts φ ≈ 0.3 (clearly positive but far from 1) — growth shocks are mildly persistent with a half-life under one quarter. There is genuine business-cycle momentum, but it dissipates fast.
  • The predictive forecast reverts to the mean (~0.7%/qtr) almost immediately and barely beats the naive-mean benchmark — the Bayesian version of the same "GNP growth is hard to predict" conclusion, now with full predictive uncertainty.