GARCH-in-Mean — the risk-return tradeoff (from scratch)¶
Does higher volatility earn a higher expected return?¶
Finance theory says risk should be paid for: a more volatile asset ought to offer a higher expected return. GARCH-in-Mean (Engle, Lilien & Robins 1987) tests this directly by letting the conditional volatility enter the mean equation: $$r_t = \mu + \lambda\,\sqrt{h_t} + \varepsilon_t,\qquad \varepsilon_t=\sqrt{h_t}\,z_t,\qquad h_t=\omega+(\alpha+\gamma\,\mathbf 1[\varepsilon_{t-1}<0])\varepsilon_{t-1}^2+\beta h_{t-1},$$ with $z_t\sim N(0,1)$ or a standardized $t_\nu$. The new parameter $\lambda$ is the price of risk: $\lambda>0$ is the risk-return tradeoff — when the model expects a turbulent day ($\sqrt{h_t}$ high) it also expects a higher return.
Unlike a plain GARCH, the mean now depends on $h_t$, which depends on past residuals, which depend on the mean — so the likelihood must be evaluated by a genuine forward recursion (it cannot be vectorised like an ordinary GARCH filter). We estimate by maximum likelihood (numerical-Hessian standard errors) and by a random-walk Metropolis sampler for the posterior of $\lambda$. Everything is in garchm_scratch.py; the data is the S&P 500 daily series used throughout the arc.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import garchm_scratch as M
from scipy.stats import chi2
# --- can the estimator recover a KNOWN price of risk? ---
def simulate(T, mu, lam, omega, alpha, beta, seed=0):
rng = np.random.default_rng(seed); h = np.empty(T); r = np.empty(T); eps = np.empty(T)
h[0] = omega / (1 - alpha - beta)
for t in range(T):
eps[t] = np.sqrt(h[t]) * rng.standard_normal()
r[t] = mu + lam * np.sqrt(h[t]) + eps[t] # in-mean term lambda * sqrt(h_t)
if t + 1 < T:
h[t + 1] = omega + alpha * eps[t] ** 2 + beta * h[t]
return r
rsim = simulate(6000, 0.02, 0.10, 0.02, 0.08, 0.90, seed=1)
th, se = M.mle(rsim, inmean="sd")
print("simulation recovery (T=6000):")
for nm, v, s, tr in zip(M.names(), th, se, [0.02, 0.10, 0.02, 0.08, 0.90]):
print(" %-7s %8.4f (se %.4f) true %.3f" % (nm, v, s, tr))
simulation recovery (T=6000): mu 0.0177 (se 0.0486) true 0.020 lambda 0.0964 (se 0.0548) true 0.100 omega 0.0222 (se 0.0044) true 0.020 alpha 0.0860 (se 0.0084) true 0.080 beta 0.8921 (se 0.0108) true 0.900
Estimating the price of risk on the S&P 500¶
Three specifications: Normal errors, Student-$t$ errors, and GJR (with leverage) plus Student-$t$. In each, $\lambda$ is the in-mean coefficient; a likelihood-ratio test compares the model to its $\lambda=0$ restriction (plain GARCH).
R = pd.read_csv("sp500ret.csv")["ret"].values.astype(float)
specs = [("GARCH-M-N", False, False), ("GARCH-M-t", False, True), ("GJR-M-t", True, True)]
fits = {}
print("Risk premium lambda (in-mean on conditional SD), S&P 500 daily:")
for tag, gjr, st in specs:
th, se = M.mle(R, gjr=gjr, student=st, inmean="sd")
ll = M.loglik(th, R, gjr, "sd"); fits[tag] = dict(th=th, se=se)
print(" %-11s lambda %6.4f se %.4f t %5.2f AIC %.0f" % (tag, th[1], se[1], th[1] / se[1], -2 * ll + 2 * len(th)))
th_t = fits["GARCH-M-t"]["th"]
llr = M.loglik_restricted(R, gjr=False, student=True, inmean="sd") # re-fit the lambda=0 model
lr = 2 * (M.loglik(th_t, R, False, "sd") - llr)
print(" LR test (in-mean vs re-fit lambda=0, GARCH-M-t): LR %.2f p %.3f" % (lr, 1 - chi2.cdf(lr, 1)))
Risk premium lambda (in-mean on conditional SD), S&P 500 daily: GARCH-M-N lambda 0.0582 se 0.0399 t 1.46 AIC 15087 GARCH-M-t lambda 0.0218 se 0.0362 t 0.60 AIC 14684 GJR-M-t lambda 0.0359 se 0.0360 t 1.00 AIC 14602 LR test (in-mean vs re-fit lambda=0, GARCH-M-t): LR 0.36 p 0.546
The posterior of the risk premium¶
Maximum likelihood gives a point estimate and a standard error; a Bayesian random-walk Metropolis sampler gives the whole posterior of $\lambda$ — and directly the probability $P(\lambda>0)$ that the risk-return tradeoff is positive.
post = M.metropolis(R, th_t, fits["GARCH-M-t"]["se"], gjr=False, student=True,
inmean="sd", N=10000, burn=3000, scale=0.35, seed=1)
lam_draws = post["theta"][:, 1]; pgt0 = (lam_draws > 0).mean()
print("Bayesian GARCH-M-t (accept %.2f): lambda mean %.4f 95%% [%.4f, %.4f] P(lambda>0) %.2f"
% (post["accept"], lam_draws.mean(), np.quantile(lam_draws, .025), np.quantile(lam_draws, .975), pgt0))
fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
ax[0].hist(lam_draws, bins=50, density=True, color="steelblue", alpha=.75, label="Bayesian posterior")
mle_l, mle_s = th_t[1], fits["GARCH-M-t"]["se"][1]
xx = np.linspace(lam_draws.min(), lam_draws.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="MLE $\\pm$ se")
ax[0].axvline(0, color="0.4", lw=1.2, ls="--")
ax[0].set_title("Posterior of the risk premium $\\lambda$ (GARCH-M-t)\nP($\\lambda>0$) = %.2f — positive but uncertain" % pgt0)
ax[0].set_xlabel("$\\lambda$ (price of risk, per unit daily vol)"); ax[0].set_yticks([]); ax[0].legend(fontsize=8)
labels = [s[0] for s in specs]; lams = [fits[l]["th"][1] for l in labels]; ses = [fits[l]["se"][1] for l in labels]
y = np.arange(len(labels))[::-1]
ax[1].errorbar(lams, y, xerr=1.96 * np.array(ses), fmt="o", color="seagreen", capsize=4, ms=7, label="MLE $\\pm$1.96 se")
ax[1].plot([np.quantile(lam_draws, .025), np.quantile(lam_draws, .975)], [y[1], y[1]], color="steelblue", lw=6, alpha=.35, solid_capstyle="round", label="Bayes 95% (GARCH-M-t)")
ax[1].axvline(0, color="firebrick", lw=1.4, ls="--")
ax[1].set_yticks(y); ax[1].set_yticklabels(labels); ax[1].set_ylim(-0.6, len(labels) - 0.4)
ax[1].set_xlabel("$\\lambda$ (risk premium)"); ax[1].legend(fontsize=8, loc="lower right")
ax[1].set_title("Risk premium across specifications\n(all intervals straddle 0; shrinks under fat tails)")
plt.tight_layout(); plt.show()
Bayesian GARCH-M-t (accept 0.25): lambda mean 0.0296 95% [-0.0379, 0.0951] P(lambda>0) 0.79
Does frequency matter? The tradeoff at monthly frequency¶
The premium $\lambda\sqrt{h_t}$ is a drift term: over a horizon $\tau$ it accumulates like $\tau$, while the return's noise grows only like $\sqrt{\tau}$. So the signal-to-noise ratio for detecting it improves at lower frequencies — the mean effect is swamped by noise in daily data but has more room to show up monthly. This is the classic finding of French, Schwert & Stambaugh (1987) and Ghysels, Santa-Clara & Valkanov (2005). We re-fit the (correctly tail-specified) GARCH-M-$t$ to a long monthly S&P 500 series (spx_monthly.csv, from Yahoo Finance).
rm = pd.read_csv("spx_monthly.csv")["ret"].values.astype(float) # long monthly S&P 500 (yfinance ^GSPC)
th_m, se_m = M.mle(rm, gjr=False, student=True, inmean="sd")
post_m = M.metropolis(rm, th_m, se_m, gjr=False, student=True, inmean="sd", N=10000, burn=3000, scale=0.35, seed=1)
pgt0_m = (post_m["theta"][:, 1] > 0).mean() # monthly posterior P(lambda>0)
print("GARCH-M-t risk premium by frequency (posterior probability it is positive):")
print(" daily (n=%d) lambda %.4f P(lambda>0) %.2f" % (len(R), fits["GARCH-M-t"]["th"][1], pgt0))
print(" monthly (n=%d) lambda %.4f P(lambda>0) %.2f" % (len(rm), th_m[1], pgt0_m))
lam = [fits["GARCH-M-t"]["th"][1], th_m[1]]; pg = [pgt0, pgt0_m]
fig, ax = plt.subplots(1, 2, figsize=(12, 4.5))
labels = ["daily\n(1987-2009)", "monthly\n(1985-2026)"]; x = np.arange(2); cols = ["steelblue", "firebrick"]
ax[0].bar(x, lam, 0.55, color=cols)
for xi, v in zip(x, lam): ax[0].text(xi, v + 0.004, "%.3f" % v, ha="center", fontsize=10)
ax[0].set_xticks(x); ax[0].set_xticklabels(labels); ax[0].set_ylabel("estimated risk premium $\\lambda$")
ax[0].set_title("The premium is ~10x larger at monthly frequency\n(GARCH-M-t)")
ax[1].bar(x, pg, 0.55, color=cols)
for xi, v in zip(x, pg): ax[1].text(xi, v + 0.012, "%.2f" % v, ha="center", fontsize=11)
ax[1].axhline(0.5, color="0.3", ls="--", lw=1.3, label="0.5 (no information)")
ax[1].set_xticks(x); ax[1].set_xticklabels(labels); ax[1].set_ylabel("posterior P($\\lambda>0$)"); ax[1].set_ylim(0, 1.0)
ax[1].set_title("...so we are more confident the premium is positive\n(but not decisively — even monthly)")
ax[1].legend(fontsize=8, loc="lower right")
plt.tight_layout(); plt.show()
GARCH-M-t risk premium by frequency (posterior probability it is positive): daily (n=5523) lambda 0.0218 P(lambda>0) 0.79 monthly (n=498) lambda 0.2114 P(lambda>0) 0.94
The premium is larger and clearer at monthly frequency — $\hat\lambda$ rises roughly tenfold ($0.02\to0.21$; though part of that raw jump is mechanical, since $\lambda$ multiplies $\sqrt{h_t}$ and monthly volatility is on a larger scale), and — the scale-free evidence — our posterior confidence that it is positive climbs from $P(\lambda>0)=0.79$ (daily) to $0.94$ (monthly). But even $0.94$ is not decisive — well short of the $\gtrsim0.975$ one would want before declaring the tradeoff established. Four decades of monthly data move the needle without settling the question; this is why the literature reaches for very long samples (Lundblad 2007 uses nearly a century).
(Aside: the Normal model reports a larger $\lambda$ at daily frequency — but that is spurious, the same fat tails being misread as a premium; once the $t$ tails are modeled it collapses, which is why this uses GARCH-M-$t$.)
Results¶
- The estimator is sound — on simulated data it recovers a known $\lambda=0.10$ (and the variance parameters) from 6,000 observations.
- The risk premium has the right sign but is statistically weak. On the S&P 500, $\hat\lambda>0$ in every specification (posterior $P(\lambda>0)=0.79$), but no interval excludes zero and the likelihood-ratio test for the in-mean term is insignificant (LR 0.36, $p=0.55$). This is the well-known risk-return tradeoff "puzzle": the positive relation theory predicts is real, but a tiny mean effect is swamped by daily-return noise. (Lower frequency helps — see the monthly comparison above.)
- Fat tails absorb the premium. Under Normal errors $\hat\lambda=0.058$ ($t=1.46$); allowing Student-$t$ innovations shrinks it to $0.022$ ($t=0.60$). Part of what a Gaussian model reads as a risk premium is really just heavy tails — once the tails are modeled, the premium fades.
- Leverage is robust; the premium is not. In GJR-M-t the leverage term stays firmly positive ($\gamma\approx0.11$, as quantified in
garch_asym_python), while $\lambda$ stays small and insignificant. The asymmetry of volatility is a firm feature of the data; the mean risk premium is not.
The economic message: GARCH-M is the right tool to test the risk-return tradeoff, and it confirms the direction — but daily equity data cannot resolve the magnitude. Lower-frequency (monthly) data, longer samples, and models that separate the premium from fat tails and leverage do better.
Reference¶
- Engle, R. F., Lilien, D. M. & Robins, R. P. (1987). Estimating time varying risk premia in the term structure: the ARCH-M model. Econometrica 55, 391–407.
Companion notebooks: garchm_pymc.ipynb (PyMC/NUTS posterior of $\lambda$) and garchm_R.ipynb (rugarch with archm=TRUE).