GARCH-in-Mean in PyMC (NUTS)¶

The Bayesian-PPL corner of the GARCH-M trio (from-scratch → PyMC → R)¶

The from-scratch notebook sampled the risk premium $\lambda$ with a hand-written random-walk Metropolis; here PyMC's NUTS samples the same GARCH-M-$t$ model, with the variance recursion written as a pytensor.scan. Because the mean $\mu+\lambda\sqrt{h_t}$ feeds back into $h_t$, the scan carries the conditional variance forward one step at a time. Stationarity is imposed by a reparameterisation (persistence $\sim\text{Beta}$, split into $\alpha,\beta$); $\lambda$ gets a weakly-informative $N(0,0.5)$ prior so the data speak.

In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import pymc as pm, pytensor, pytensor.tensor as pt, arviz as az

R = pd.read_csv("sp500ret.csv")["ret"].values.astype(float); v = float(np.var(R))
def step(r_prev, h_prev, mu, lam, omega, alpha, beta):
    eps = r_prev - mu - lam * pt.sqrt(h_prev)                     # in-mean residual
    return omega + alpha * eps ** 2 + beta * h_prev

with pm.Model() as model:
    mu = pm.Normal("mu", 0.0, 0.2); lam = pm.Normal("lam", 0.0, 0.5)      # lambda = price of risk
    persist = pm.Beta("persist", 20, 1.5); frac = pm.Beta("frac", 1.5, 10)
    alpha = pm.Deterministic("alpha", persist * frac); beta = pm.Deterministic("beta", persist * (1 - frac))
    uncvar = pm.HalfNormal("uncvar", sigma=v); omega = pm.Deterministic("omega", uncvar * (1 - persist))
    nu = pm.Deterministic("nu", 2.0 + pm.Gamma("nu_m2", 2.0, 0.2))
    h_rest, _ = pytensor.scan(step, sequences=[pt.as_tensor(R[:-1])], outputs_info=[uncvar],
                              non_sequences=[mu, lam, omega, alpha, beta])
    h = pt.concatenate([pt.stack([uncvar]), h_rest]); sd = pt.sqrt(h)
    pm.StudentT("obs", nu=nu, mu=mu + lam * sd, sigma=sd * pt.sqrt((nu - 2.0) / nu), observed=R)
    idata = pm.sample(1000, tune=1000, chains=2, target_accept=0.95, random_seed=1)

print(az.summary(idata, var_names=["mu", "lam", "alpha", "beta", "persist", "nu"]).to_string())
lam_d = idata.posterior["lam"].values.ravel()
print("\nNUTS lambda mean %.4f  95%% [%.4f, %.4f]  P(lambda>0) %.2f"
      % (lam_d.mean(), np.quantile(lam_d, .025), np.quantile(lam_d, .975), (lam_d > 0).mean()))

fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
ax[0].hist(lam_d, bins=55, density=True, color="slateblue", alpha=.75, label="PyMC NUTS posterior")
mle_l, mle_s = 0.0218, 0.0362; xx = np.linspace(lam_d.min(), lam_d.max(), 200)
ax[0].plot(xx, np.exp(-0.5*((xx-mle_l)/mle_s)**2)/(mle_s*np.sqrt(2*np.pi)), color="firebrick", lw=2, label="from-scratch MLE $\\pm$se")
ax[0].axvline(0, color="0.4", lw=1.2, ls="--"); ax[0].set_yticks([]); ax[0].set_xlabel("$\\lambda$"); ax[0].legend(fontsize=8)
ax[0].set_title("PyMC/NUTS posterior of the risk premium $\\lambda$\nP($\\lambda>0$) = %.2f" % (lam_d > 0).mean())
eng = [("from-scratch MLE", mle_l, [mle_l-1.96*mle_s, mle_l+1.96*mle_s]),
       ("from-scratch Metropolis", 0.0296, [-0.0379, 0.0951]),
       ("PyMC NUTS", lam_d.mean(), [np.quantile(lam_d, .025), np.quantile(lam_d, .975)])]
y = np.arange(len(eng))[::-1]
for (nm, m, ci), yi in zip(eng, y):
    ax[1].plot(ci, [yi, yi], color="slateblue", lw=4, alpha=.5); ax[1].plot(m, yi, "o", color="navy", ms=8)
ax[1].axvline(0, color="firebrick", lw=1.4, ls="--"); ax[1].set_yticks(y); ax[1].set_yticklabels([e[0] for e in eng])
ax[1].set_ylim(-0.6, len(eng)-0.4); ax[1].set_xlabel("$\\lambda$  (risk premium)")
ax[1].set_title("Three engines, one (weak) premium\n(all agree; all straddle 0)")
plt.tight_layout(); plt.show()

           mean      sd eti89_lb eti89_ub  ess_bulk  ess_tail r_hat mcse_mean  mcse_sd
mu         0.04   0.031  -0.0082    0.087       993       857  1.00     0.001  0.00082
lam       0.025   0.038   -0.034    0.086       995       909  1.00    0.0012  0.00098
alpha    0.0633  0.0067    0.053    0.074      1510      1224  1.01   0.00017  0.00013
beta     0.9315  0.0075     0.92     0.94      1497      1234  1.01   0.00019  0.00014
persist  0.9949  0.0021     0.99        1       955      1103  1.00   6.6e-05  5.6e-05
nu          6.4    0.51      5.7      7.2      1509      1116  1.01     0.013     0.01

NUTS lambda mean 0.0249  95% [-0.0478, 0.1023]  P(lambda>0) 0.75
No description has been provided for this image

Results¶

  • Converged and clean. All $\hat R\le1.01$ with ~900+ effective samples; the scan-based GARCH-M-$t$ posterior is well-behaved under NUTS.
  • Three engines, one answer. The NUTS posterior of $\lambda$ (mean $0.025$, 95% $[-0.048,\,0.102]$, $P(\lambda>0)=0.75$) sits right on top of the from-scratch MLE and the hand-written Metropolis — the right panel stacks all three, and every interval straddles zero. Independent samplers, one weak-but-positive premium.
  • The variance parameters are textbook. Persistence $\alpha+\beta\approx0.995$, $\nu\approx6.4$ — the same near-integrated, fat-tailed volatility the plain-GARCH notebooks found; adding the in-mean term barely moves them, consistent with $\lambda$ being small.

Together with garchm_python (from-scratch MLE + Metropolis) and garchm_R (rugarch), the risk-return tradeoff is now estimated four ways across three engines, all agreeing: positive in sign, weak at daily frequency, and — as the frequency comparison showed — more convincingly positive only at lower frequencies.