GJR-GARCH from scratch — the leverage effect, five methods¶

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

The symmetric GARCH (garch_python.ipynb) treats good and bad news alike. The GJR-GARCH (Glosten–Jagannathan–Runkle) adds a leverage term $\gamma$ that lets negative returns raise volatility more:

$$\sigma_t^2=\omega+\big(\alpha+\gamma\,\mathbf 1[r_{t-1}<0]\big)r_{t-1}^2+\beta\sigma_{t-1}^2.$$

Good news moves vol by $\alpha$, bad news by $\alpha+\gamma$; $\gamma>0$ is the leverage effect. The same five estimators from garch_scratch.py run on the extended parameter vector $(\omega,\alpha,\gamma,\beta)$ (+ $\nu$ for Student-t) — since Gauss–Hermite returns each model's marginal likelihood, we rank GARCH vs GJR × Normal vs t by Bayes factor.

The full model¶

Data. $r_t,\ t=1,\dots,T$ — daily S&P 500 log-returns (in %), demeaned so $\sum_t r_t = 0$. Let $\mathcal F_{t-1}$ denote the information up to $t-1$.

Observation (return) equation. $$r_t = \sigma_t\,\varepsilon_t,\qquad \varepsilon_t \overset{iid}{\sim} D(0,1),$$ with $\sigma_t^2 = \mathrm{Var}(r_t\mid\mathcal F_{t-1})$ the conditional variance and $D$ a zero-mean, unit-variance innovation.

Conditional variance — GJR-GARCH(1,1). $$\boxed{\ \sigma_t^2 \;=\; \omega \;+\; \big(\alpha + \gamma\,\mathbf 1[r_{t-1}<0]\big)\,r_{t-1}^2 \;+\; \beta\,\sigma_{t-1}^2\ }$$ initialized at $\sigma_1^2 = \widehat{\mathrm{Var}}(r)$ (the sample variance). A good-news shock ($r_{t-1}\ge0$) feeds volatility through $\alpha$; a bad-news shock ($r_{t-1}<0$) through $\alpha+\gamma$. Setting $\gamma=0$ recovers the symmetric GARCH(1,1).

Innovation distribution (fit both):

  • Normal: $\varepsilon_t\sim N(0,1)$, i.e. $r_t\mid\mathcal F_{t-1}\sim N(0,\sigma_t^2)$.
  • Student-t: $\varepsilon_t\sim t_\nu$ standardized to unit variance — the conditional density of $r_t$ is Student-$t_\nu$ with scale $s_t=\sigma_t\sqrt{(\nu-2)/\nu}$, so that $\mathrm{Var}(r_t\mid\mathcal F_{t-1})=\sigma_t^2$ exactly, with tail index $\nu>2$.

Parameters. $\theta=(\omega,\alpha,\gamma,\beta)$ for the Gaussian GJR, and $\theta=(\omega,\alpha,\gamma,\beta,\nu)$ for the Student-t GJR:

symbol role
$\omega>0$ baseline variance level
$\alpha\ge0$ ARCH / good-news impact
$\gamma$ leverage — extra impact of bad news ($\gamma>0$ = leverage effect)
$\beta\ge0$ GARCH / volatility persistence
$\nu>2$ Student-$t$ tail (fatter as $\nu\downarrow$)

Admissible region (the support of the prior): $$\omega>0,\quad \alpha\ge0,\quad \beta\ge0,\quad \alpha+\gamma\ge0,\quad \underbrace{\alpha+\tfrac12\gamma+\beta<1}_{\text{covariance stationarity}},\quad \nu>2 .$$ The $\tfrac12\gamma$ appears because a symmetric shock is negative half the time, so the average news-impact coefficient is $\alpha+\tfrac12\gamma$.

Likelihood. Running the recursion gives $\sigma_1^2,\dots,\sigma_T^2$; with $\ell_t$ the conditional log-density of $r_t$ (Normal or standardized-$t$), $$\log L(\theta)=\sum_{t=1}^{T}\ell_t\!\big(r_t;\ \sigma_t^2(\theta),\ \nu\big).$$

Prior. Flat (uniform) on the admissible region above, so the log-posterior equals the log-likelihood there (up to a constant). Every estimator below — the MLE and the four Bayesian methods (Gauss–Hermite quadrature, importance sampling, Metropolis–Hastings, griddy Gibbs) — targets this same posterior; the Gauss–Hermite marginal likelihood $\int L(\theta)\,d\theta$ over the region gives the Bayes factors used for model comparison.

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

# MLE + log marginal likelihood for the 2x2 family (Gauss-Hermite gives logML for the Bayes factors)
fits = {}
for tag, gj, stu, K in [("GARCH-N", False, False, 12), ("GARCH-t", False, True, 9),
                        ("GJR-N", True, False, 8), ("GJR-t", True, True, 7)]:
    mle, se, cov = G.garch_mle(r, gjr=gj, student=stu)
    logml = G.gh_quadrature(r, mle, cov, K=K, gjr=gj)["logml"]
    fits[tag] = dict(mle=mle, se=se, cov=cov, gjr=gj, student=stu, logml=logml, names=G.names(gj, stu))
    print("%-8s logML %.1f   %s" % (tag, logml, "  ".join("%s %.4f" % (n, v) for n, v in zip(G.names(gj, stu), mle))))
print("\nlog Bayes factors vs GARCH-N:  " +
      "  ".join("%s %+.1f" % (t, fits[t]["logml"] - fits["GARCH-N"]["logml"]) for t in fits))

5523 daily returns
GARCH-N  logML -7558.7   omega 0.0135  alpha 0.0881  beta 0.9046

GARCH-t  logML -7359.5   omega 0.0060  alpha 0.0611  beta 0.9358  nu 6.2174
GJR-N    logML -7482.3   omega 0.0186  alpha 0.0078  beta 0.9096  gamma 0.1332

GJR-t    logML -7315.6   omega 0.0119  alpha 0.0083  beta 0.9244  gamma 0.1147  nu 6.7725

log Bayes factors vs GARCH-N:  GARCH-N +0.0  GARCH-t +199.2  GJR-N +76.4  GJR-t +243.1
In [2]:
# five methods on GJR-N -> the leverage parameter gamma, recovered identically
f = fits["GJR-N"]; mle, se, cov = f["mle"], f["se"], f["cov"]
gh  = G.gh_quadrature(r, mle, cov, K=8, gjr=True)
imp = G.importance_sampling(r, mle, cov, N=25000, gjr=True)
mh  = G.metropolis_hastings(r, mle, cov, N=12000, burn=3000, scale=0.55, gjr=True)
gg  = G.griddy_gibbs(r, mle, se, N=2500, burn=600, G=40, gjr=True)
gi = f["names"].index("gamma")
print("gamma (leverage) by method:")
for lab, m in [("MLE", mle), ("Gauss-Hermite", gh["mean"]), ("importance", imp["mean"]),
               ("Metropolis-H", mh["mean"]), ("griddy Gibbs", gg["mean"])]:
    print("  %-14s %.4f" % (lab, m[gi]))
print("news-impact ratio  (alpha+gamma)/alpha = %.1f x  (bad news vs good news)" % ((mle[1]+mle[gi]) / mle[1]))

gamma (leverage) by method:
  MLE            0.1332
  Gauss-Hermite  0.1349
  importance     0.1350
  Metropolis-H   0.1361
  griddy Gibbs   0.1334
news-impact ratio  (alpha+gamma)/alpha = 18.1 x  (bad news vs good news)
In [3]:
fig, ax = plt.subplots(1, 3, figsize=(15, 4.4))
# (A) gamma posterior: five methods overlap
a = ax[0]
a.hist(mh["theta"][:, gi], bins=55, density=True, color="0.82", label="Metropolis-H")
a.hist(gg["theta"][:, gi], bins=55, density=True, histtype="step", color="seagreen", lw=1.6, label="griddy Gibbs")
a.hist(imp["theta"][:, gi], bins=55, weights=imp["w"], density=True, histtype="step", color="goldenrod", lw=1.6, label="importance")
xs = np.linspace(*a.get_xlim(), 300); a.plot(xs, st.norm.pdf(xs, gh["mean"][gi], np.sqrt(gh["cov"][gi, gi])), color="steelblue", lw=2.2, label="Gauss-Hermite")
a.axvline(mle[gi], color="firebrick", ls="--", lw=1.6, label="MLE"); a.axvline(0, color="k", lw=1)
a.set_xlabel(r"$\gamma$  (leverage)"); a.set_yticks([]); a.set_title(r"Leverage $\gamma=%.3f\gg0$ — five methods agree" % mle[gi]); a.legend(fontsize=7.5)
# (B) news-impact curve: symmetric GARCH-t vs asymmetric GJR-t
mT = fits["GARCH-t"]["mle"]; mGt = fits["GJR-t"]["mle"]; gj = fits["GJR-t"]["names"].index("gamma")
x = np.linspace(-6, 6, 400)
nicG = mT[0] + mT[1] * x**2 + mT[2] * (mT[0] / (1 - mT[1] - mT[2]))
sbar = mGt[0] / (1 - mGt[1] - mGt[gj]/2 - mGt[2])
nicJ = mGt[0] + (mGt[1] + mGt[gj] * (x < 0)) * x**2 + mGt[2] * sbar
ax[1].plot(x, nicG, color="steelblue", lw=2, label="GARCH-t (symmetric)")
ax[1].plot(x, nicJ, color="firebrick", lw=2, label="GJR-t (asymmetric)")
ax[1].axvline(0, color="0.7", lw=.7); ax[1].set_xlabel("return shock $r_{t-1}$ (%)"); ax[1].set_ylabel("next-day variance")
ax[1].set_title("News-impact curve: GJR steeper for bad news"); ax[1].legend(fontsize=8)
# (C) model comparison by log marginal likelihood
tags = list(fits); ml = np.array([fits[t]["logml"] for t in tags]); bi = int(np.argmax(ml))
ax[2].barh(range(len(tags)), ml, color=["firebrick" if i == bi else "steelblue" for i in range(len(tags))])
ax[2].set_yticks(range(len(tags))); ax[2].set_yticklabels(tags); ax[2].invert_yaxis()
ax[2].set_xlim(ml.min() - 30, ml.max() + 15)
for i, v in enumerate(ml): ax[2].text(v, i, " %.0f" % v, va="center", fontsize=8)
ax[2].set_xlabel("log marginal likelihood"); ax[2].set_title("GJR-t wins (leverage + fat tails)")
plt.tight_layout(); plt.savefig("garch_asym_python.png", dpi=120, bbox_inches="tight"); plt.show()
No description has been provided for this image
In [4]:
# cross-check: the arch package GJR-GARCH reproduces the from-scratch leverage estimate
from arch import arch_model
aG = arch_model(r, mean="Zero", vol="GARCH", p=1, o=1, q=1, dist="t").fit(disp="off")   # o=1 -> GJR
mGt = fits["GJR-t"]["mle"]; gj = fits["GJR-t"]["names"].index("gamma")
print("arch GJR-t : omega %.4f  alpha %.4f  gamma %.4f  beta %.4f  nu %.1f"
      % (aG.params["omega"], aG.params["alpha[1]"], aG.params["gamma[1]"], aG.params["beta[1]"], aG.params["nu"]))
print("scratch    : omega %.4f  alpha %.4f  gamma %.4f  beta %.4f  nu %.1f"
      % (mGt[0], mGt[1], mGt[gj], mGt[2], mGt[-1]))
print("-> the package and the from-scratch GJR agree; garch_asym_pymc (PyMC) and garch_asym_R (rugarch) close the loop.")

arch GJR-t : omega 0.0120  alpha 0.0084  gamma 0.1148  beta 0.9243  nu 6.8
scratch    : omega 0.0119  alpha 0.0083  gamma 0.1147  beta 0.9244  nu 6.8
-> the package and the from-scratch GJR agree; garch_asym_pymc (PyMC) and garch_asym_R (rugarch) close the loop.

Results¶

  • Leverage is large and unambiguous: $\gamma\approx0.13$ (GJR-N), entirely above 0 and recovered identically by all five estimators. With $\alpha\approx0.008$, a negative shock moves next-day volatility by $\alpha+\gamma\approx0.14$ versus $\approx0.008$ for a positive shock — bad news raises volatility roughly 18× more than good news of the same size (the notebook computes the exact ratio at 18.1). The news-impact curve shows the kink directly.
  • Model comparison (Bayes factors, from the Gauss–Hermite marginal likelihoods): GARCH-N $-7559$ → GARCH-t $-7360$ (fat tails, $+199$) → GJR-N $-7482$ (leverage, $+77$) → GJR-t $-7316$ (both; best overall). Leverage and fat tails each earn their keep, and the asymmetric fat-tailed model wins.
  • Cross-check: the arch package's GJR fit matches the from-scratch $(\omega,\alpha,\gamma,\beta,\nu)$ to the decimal.

This is the from-scratch corner of the Phase-2 (asymmetric) trio; garch_asym_pymc.ipynb fits the same GJR via PyMC, and garch_asym_R.ipynb via rugarch (gjrGARCH / eGARCH).