GARCH(1,1) from scratch — five ways to estimate the same model¶

The from-scratch corner of the trio (from-scratch · PyMC · R)¶

MLE · Gauss–Hermite quadrature · importance sampling · Metropolis–Hastings · griddy Gibbs¶

One Gaussian GARCH(1,1) likelihood, five estimation routes — the classical MLE and four Bayesian methods that were the workhorses of pre-HMC Bayesian GARCH (Bauwens–Lubrano, Geweke, Ritter–Tanner). The point: with a flat prior on the admissible region they all recover the same posterior.

$$r_t=\sigma_t\varepsilon_t,\quad \sigma_t^2=\omega+\alpha r_{t-1}^2+\beta\sigma_{t-1}^2,\quad \varepsilon_t\sim N(0,1),\qquad \omega>0,\ \alpha,\beta\ge0,\ \alpha+\beta<1.$$

method idea from garch_scratch.py
MLE maximise the log-likelihood; SEs from the numerical Hessian (observed information) garch_mle
Gauss–Hermite deterministic numerical integration of the posterior on a rotated GH tensor grid (Naylor–Smith adaptive GH) — also returns the log marginal likelihood gh_quadrature
importance sampling draw from a heavy-tailed multivariate-$t$ proposal, reweight by posterior/proposal importance_sampling
Metropolis–Hastings random-walk MH on $(\omega,\alpha,\beta)$ metropolis_hastings
griddy Gibbs sample each full conditional by evaluating it on a grid and inverting the CDF (Ritter–Tanner; Bauwens–Lubrano's GARCH sampler) griddy_gibbs

The variance recursion is written as a one-pole IIR filter so every log-likelihood is $O(T)$ (see the module). Data: daily S&P 500 returns, 1987–2009.

In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt, scipy.stats as st, time
import garch_scratch as G
d = pd.read_csv("sp500ret.csv", parse_dates=["date"]); r = d["ret"].values - d["ret"].values.mean()
print("%d daily returns %s..%s" % (len(r), str(d['date'].values[0])[:10], str(d['date'].values[-1])[:10]))

t = time.time(); mle, se, cov = G.garch_mle(r); print("MLE            %4.1fs" % (time.time()-t))
t = time.time(); gh  = G.gh_quadrature(r, mle, cov, K=12); print("Gauss-Hermite  %4.1fs" % (time.time()-t))
t = time.time(); imp = G.importance_sampling(r, mle, cov, N=30000); print("importance     %4.1fs" % (time.time()-t))
t = time.time(); mh  = G.metropolis_hastings(r, mle, cov, N=15000, burn=4000, scale=0.7); print("Metropolis-H   %4.1fs" % (time.time()-t))
t = time.time(); gg  = G.griddy_gibbs(r, mle, se, N=3000, burn=800, G=40); print("griddy Gibbs   %4.1fs" % (time.time()-t))

5523 daily returns 1987-03-10..2009-01-30
MLE             0.0s
Gauss-Hermite   0.1s

importance      1.2s

Metropolis-H    0.8s

griddy Gibbs   19.8s
In [2]:
# posterior mean +/- sd from each method (MLE row = estimate +/- SE)
def sd(th, w=None): return np.sqrt(np.average((th - np.average(th, 0, w))**2, 0, w))
rows = {
 "MLE (max-lik)"      : (mle,        se),
 "Gauss-Hermite"      : (gh["mean"], np.sqrt(np.diag(gh["cov"]))),
 "importance samp."   : (imp["mean"],sd(imp["theta"], imp["w"])),
 "Metropolis-H"       : (mh["mean"], mh["theta"].std(0)),
 "griddy Gibbs"       : (gg["mean"], gg["theta"].std(0)),
}
print("%-18s %-16s %-16s %-16s" % ("method", "omega", "alpha", "beta"))
for k,(m,s) in rows.items():
    print("%-18s %s" % (k, "  ".join("%.4f (%.4f)" % (m[j], s[j]) for j in range(3))))
print("\npersistence alpha+beta = %.4f    GH log marginal likelihood = %.1f" % (mle[1]+mle[2], gh["logml"]))
print("importance-sampling ESS = %.0f / 30000    Metropolis acceptance = %.2f" % (imp["ess"], mh["accept"]))

method             omega            alpha            beta            
MLE (max-lik)      0.0135 (0.0026)  0.0881 (0.0078)  0.9046 (0.0086)
Gauss-Hermite      0.0143 (0.0027)  0.0903 (0.0079)  0.9019 (0.0087)
importance samp.   0.0143 (0.0027)  0.0903 (0.0079)  0.9019 (0.0088)
Metropolis-H       0.0144 (0.0027)  0.0904 (0.0077)  0.9016 (0.0086)
griddy Gibbs       0.0143 (0.0030)  0.0901 (0.0086)  0.9020 (0.0096)

persistence alpha+beta = 0.9926    GH log marginal likelihood = -7558.7
importance-sampling ESS = 17798 / 30000    Metropolis acceptance = 0.58
In [3]:
# THE picture: the five methods land on the same marginal posterior for each parameter
lab = [r"$\omega$", r"$\alpha$", r"$\beta$"]
fig, ax = plt.subplots(1, 3, figsize=(15, 4.4))
for j in range(3):
    a = ax[j]
    a.hist(mh["theta"][:, j], bins=60, density=True, color="0.82", label="Metropolis-H")
    a.hist(gg["theta"][:, j], bins=60, density=True, histtype="step", color="seagreen", lw=1.6, label="griddy Gibbs")
    a.hist(imp["theta"][:, j], bins=60, weights=imp["w"], density=True, histtype="step", color="goldenrod", lw=1.6, label="importance samp.")
    xs = np.linspace(*a.get_xlim(), 300)
    a.plot(xs, st.norm.pdf(xs, gh["mean"][j], np.sqrt(gh["cov"][j, j])), color="steelblue", lw=2.2, label="Gauss-Hermite")
    a.axvline(mle[j], color="firebrick", ls="--", lw=1.6, label="MLE $\\pm$1.96 SE")
    a.axvspan(mle[j]-1.96*se[j], mle[j]+1.96*se[j], color="firebrick", alpha=.10)
    a.set_title(lab[j]); a.set_yticks([])
ax[0].legend(fontsize=7.5, loc="upper left")
plt.suptitle("Five estimation methods, one posterior — GARCH(1,1) on the S&P 500", y=1.02)
plt.tight_layout(); plt.savefig("garch_scratch_posteriors.png", dpi=120, bbox_inches="tight"); plt.show()
No description has been provided for this image
In [4]:
# sanity tie-back: fitted conditional volatility at the posterior mean (they all give ~the same path)
r2, s0 = G.prepare(r); s2 = G.garch_var(mh["mean"], r2, s0)
fig, ax = plt.subplots(figsize=(12, 3.6))
ax.plot(d["date"], np.sqrt(s2), color="steelblue", lw=.7, label="GARCH(1,1) conditional vol (posterior mean)")
ax.plot(d["date"], np.abs(r), color="0.8", lw=.3, zorder=0, label="|return|")
ax.set_ylabel("daily volatility (%)"); ax.set_ylim(0, 8); ax.legend(fontsize=8)
ax.set_title("Fitted volatility — persistence $\\alpha+\\beta$ = %.3f" % (mle[1]+mle[2]))
plt.tight_layout(); plt.savefig("garch_scratch_vol.png", dpi=120, bbox_inches="tight"); plt.show()
No description has been provided for this image

Results¶

  • Five methods, one posterior. MLE, Gauss–Hermite quadrature, importance sampling, Metropolis–Hastings and griddy Gibbs all land on $\omega\approx0.014$, $\alpha\approx0.090$, $\beta\approx0.902$ (persistence $\approx0.993$) with matching posterior SDs — the three panels above are indistinguishable across methods. The MLE sits a hair off the Bayesian posterior mean (the usual small skew of the GARCH posterior), well within the SE.
  • What each buys you. MLE is instant but gives only a point + asymptotic SE. Gauss–Hermite is deterministic (no Monte-Carlo error) and hands back the log marginal likelihood for free — useful for Bayes factors. Importance sampling is trivially parallel and reports an ESS (here ~60% efficiency with a $t_5$ proposal). Metropolis–Hastings and griddy Gibbs scale to models where quadrature's tensor grid explodes; griddy Gibbs — evaluating each full conditional on a grid and inverting the CDF — was the Bayesian-GARCH sampler before HMC.
  • Cost. Quadrature/IS/MH are sub-second to ~1 s; griddy Gibbs is the slowest (a grid of likelihoods per parameter per sweep) but still ~20 s here.

This is the methodological backbone for the arc: the packaged fits — garch_pymc.ipynb (PyMC / NUTS) and garch_R.ipynb (rugarch ML + bayesGARCH MCMC) — are the fast versions of exactly this posterior, and the arch cross-check below confirms the point estimates.

Below, the same five methods with Student-t innovations (one extra parameter $\nu$) — and, because Gauss–Hermite returns the marginal likelihood, a one-line Bayes factor for Normal vs $t$.

Student-t innovations — the same five methods, plus $\nu$, plus a Bayes factor¶

Fat tails are the other headline stylized fact. Replace the Gaussian $\varepsilon_t$ with a standardized Student-t ($\mathrm{Var}=1$), adding degrees of freedom $\nu$, so $\theta=(\omega,\alpha,\beta,\nu)$. Every routine in garch_scratch.py is dimension-general, so the same five estimators now run over four parameters — the Gauss–Hermite grid becomes $K^4$, griddy Gibbs cycles one more full conditional. And since gh_quadrature returns each model's log marginal likelihood, Normal-vs-$t$ is a one-line Bayes factor.

In [5]:
mle_t, se_t, cov_t = G.garch_mle(r, student=True)
gh_t  = G.gh_quadrature(r, mle_t, cov_t, K=9)
imp_t = G.importance_sampling(r, mle_t, cov_t, N=25000)
mh_t  = G.metropolis_hastings(r, mle_t, cov_t, N=12000, burn=3000, scale=0.6)
gg_t  = G.griddy_gibbs(r, mle_t, se_t, N=2500, burn=600, G=40)
rows_t = {"MLE (max-lik)": mle_t, "Gauss-Hermite": gh_t["mean"], "importance samp.": imp_t["mean"],
          "Metropolis-H": mh_t["mean"], "griddy Gibbs": gg_t["mean"]}
print("%-18s %-9s %-9s %-9s %-8s" % ("method", "omega", "alpha", "beta", "nu"))
for k, m in rows_t.items():
    print("%-18s %8.4f %8.4f %8.4f %7.2f" % (k, m[0], m[1], m[2], m[3]))
print("\nnu ~ %.2f (all methods)   persistence alpha+beta = %.4f" % (mle_t[3], mle_t[1] + mle_t[2]))
print("log marginal likelihood:  Normal %.1f    Student-t %.1f    -> log Bayes factor = %+.1f  (decisive for t)"
      % (gh["logml"], gh_t["logml"], gh_t["logml"] - gh["logml"]))

method             omega     alpha     beta      nu      
MLE (max-lik)        0.0060   0.0611   0.9358    6.22
Gauss-Hermite        0.0069   0.0631   0.9329    6.34
importance samp.     0.0068   0.0629   0.9330    6.35
Metropolis-H         0.0068   0.0626   0.9333    6.38
griddy Gibbs         0.0069   0.0627   0.9331    6.34

nu ~ 6.22 (all methods)   persistence alpha+beta = 0.9969
log marginal likelihood:  Normal -7558.7    Student-t -7359.5    -> log Bayes factor = +199.2  (decisive for t)
In [6]:
# nu posterior: the new parameter, recovered identically by all five methods
a = plt.subplots(figsize=(7.6, 4.4))[1]
a.hist(mh_t["theta"][:, 3], bins=60, density=True, color="0.82", label="Metropolis-H")
a.hist(gg_t["theta"][:, 3], bins=60, density=True, histtype="step", color="seagreen", lw=1.6, label="griddy Gibbs")
a.hist(imp_t["theta"][:, 3], bins=60, weights=imp_t["w"], density=True, histtype="step", color="goldenrod", lw=1.6, label="importance samp.")
xs = np.linspace(*a.get_xlim(), 300)
a.plot(xs, st.norm.pdf(xs, gh_t["mean"][3], np.sqrt(gh_t["cov"][3, 3])), color="steelblue", lw=2.2, label="Gauss-Hermite")
a.axvline(mle_t[3], color="firebrick", ls="--", lw=1.6, label="MLE $\\pm$1.96 SE")
a.axvspan(mle_t[3] - 1.96*se_t[3], mle_t[3] + 1.96*se_t[3], color="firebrick", alpha=.10)
a.set_xlabel(r"$\nu$  (Student-t degrees of freedom)"); a.set_yticks([])
a.set_title(r"Fat tails: $\nu \approx %.1f$ — five methods, one posterior (log Bayes factor vs Normal = $+%.0f$)"
            % (mle_t[3], gh_t["logml"] - gh["logml"]))
a.legend(fontsize=8)
plt.tight_layout(); plt.savefig("garch_scratch_t.png", dpi=120, bbox_inches="tight"); plt.show()
No description has been provided for this image

Student-t results¶

  • $\nu \approx 6.3$, recovered identically by all five methods — the returns are conditionally fat-tailed even after GARCH scaling (the same $\nu\approx6$ the arch package gives in garch_python.ipynb).
  • Decisive Bayes factor. Gauss–Hermite gives log marginal likelihoods $-7559$ (Normal) vs $-7360$ (Student-$t$): a log Bayes factor of $+199$ — overwhelming evidence for $t$, at essentially no extra code.
  • Where the fat tails go. Under $t$, the ARCH term shrinks ($\alpha\ \,0.088\!\to\!0.063$) and $\beta$ rises ($0.905\!\to\!0.933$, persistence $0.993\!\to\!0.996$): with $\nu$ explaining the excess kurtosis, the volatility dynamics no longer have to over-react through $\alpha$. This is the classic Normal-vs-$t$ GARCH contrast, here shown five ways with a Bayes factor to settle it.

Cross-check — the x86_64 package (frequentist ML)¶

The from-scratch estimates should match a standard package. x86_64 (Sheppard) fits GARCH(1,1) by ML in one line; it reproduces the from-scratch MLE for both Normal and Student-t innovations.

In [7]:
from arch import arch_model
aN = arch_model(r, mean="Zero", vol="GARCH", p=1, q=1, dist="normal").fit(disp="off")
aT = arch_model(r, mean="Zero", vol="GARCH", p=1, q=1, dist="t").fit(disp="off")
print("arch GARCH-N : omega %.4f  alpha %.3f  beta %.3f          (from-scratch MLE   %.4f %.3f %.3f)"
      % (aN.params["omega"], aN.params["alpha[1]"], aN.params["beta[1]"], mle[0], mle[1], mle[2]))
print("arch GARCH-t : omega %.4f  alpha %.3f  beta %.3f  nu %.1f (from-scratch t-MLE %.4f %.3f %.3f %.1f)"
      % (aT.params["omega"], aT.params["alpha[1]"], aT.params["beta[1]"], aT.params["nu"], mle_t[0], mle_t[1], mle_t[2], mle_t[3]))
print("-> the package and the from-scratch estimator agree; PyMC (garch_pymc) and R (garch_R) close the 4-engine loop.")

arch GARCH-N : omega 0.0135  alpha 0.088  beta 0.904          (from-scratch MLE   0.0135 0.088 0.905)
arch GARCH-t : omega 0.0061  alpha 0.061  beta 0.936  nu 6.2 (from-scratch t-MLE 0.0060 0.061 0.936 6.2)
-> the package and the from-scratch estimator agree; PyMC (garch_pymc) and R (garch_R) close the 4-engine loop.