Realized Volatility

Python · R · the volatility arc, measured

Overview

GARCH and stochastic volatility treat volatility as latent — they infer it from a single daily return. But with intraday prices, a day's variance can be measured directly: sum its squared high-frequency returns to get the realized variance $RV_t=\sum_i r_{t,i}^2$ . This example runs the realized-volatility arc in three phases — measure it, forecast it, and bridge it back to GARCH — and ends with a forecasting tournament that asks whether the intraday information actually pays.

RV_t = \sum_{i} r_{t,i}^2, \qquad \log RV_t = \beta_0 + \beta_d\,\log RV_{t-1}^{(d)} + \beta_w\,\log RV_{t-1}^{(w)} + \beta_m\,\log RV_{t-1}^{(m)} + \varepsilon_t

Realized variance sums a day's squared intraday returns; Corsi's HAR regresses log RV on its own daily, weekly, and monthly averages.

Three phases

Phase 1 — the measures. Realized variance from 1- and 5-minute S&P 500 prices, and the two stylized facts that write the next model: RV has long memory (autocorrelation still $\approx 0.1$ at a 100-day lag — a short-memory GARCH cannot reproduce this), and while RV is wildly right-skewed, $\log RV$ is close to Gaussian. Long memory + log-normality ⇒ regress $\log RV_t$ on its own past over several horizons.

Phase 2 — HAR-RV. Corsi's (2009) Heterogeneous Autoregressive model turns both facts into a forecasting model of disarming simplicity — an OLS of $\log RV$ on its daily, weekly, and monthly averages, mimicking long memory without fractional integration. A 4-parameter regression that beats a fitted GARCH out-of-sample (QLIKE 0.223 vs 0.296) and adapts faster in the 2008 crisis — because it is fed realized information. Built three ways (from-scratch OLS + conjugate Bayes / PyMC NUTS / R highfrequency::HARmodel), plus two extensions: HAR-CJ (jump/continuous split) and HARQ (realized quarticity). The verdict matches the literature — modeling RV's measurement error pays (HARQ $\approx 15\%$ lower QLIKE), separating jumps does not (jumps are real but unpredictable).

Phase 3 — Realized GARCH & the tournament. Hansen, Huang & Shek's (2012) Realized GARCH is the formal bridge: it keeps GARCH's latent $h_t$ but lets the informative $RV$ drive the recursion and adds a measurement equation linking return shocks to log-RV (carrying the leverage via $\tau_1<0$ ). It tracks realized volatility more closely than standard GARCH (correlation 0.75 vs 0.71). The closing forecast horse race scores all four paradigms on the same QLIKE loss — and finds a clean two-camp split: the models that use realized measures (Realized GARCH 0.222, HAR 0.223) forecast $\approx 25\%$ better than the returns-only ones (SV 0.291, GARCH-t 0.296). The verdict of the whole arc: it matters little how the intraday information enters — using it at all is the single biggest improvement to volatility forecasting, worth more than the model class.

Notebooks

Python Phase 1 · Realized measures — RV from intraday prices, long memory & log-normality (Python) Open → Python Phase 2 · HAR-RV from scratch — Corsi (2009), OLS + conjugate Bayes, beats GARCH out-of-sample (Python) Open → Python Phase 2 · HAR-RV in PyMC — Bayesian posterior & model comparison (Python) Open → R Phase 2 · HAR-RV in R — highfrequency::HARmodel (R) Open → Python Phase 2 · HAR extensions — HAR-CJ (jumps) & HARQ (measurement error / realized quarticity) (Python) Open → Python Phase 3 · Realized GARCH — Hansen–Huang–Shek (2012), the GARCH↔realized-volatility bridge (Python) Open → Python Phase 3 · Forecast tournament — HAR & Realized GARCH vs GARCH & SV, one QLIKE scoring rule (Python) Open →

Downloads

realgarch.py Realized GARCH (Hansen–Huang–Shek) — joint return + RV likelihood (NumPy / scipy) spx_rv_ret.csv Daily S&P 500 realized variance + returns (the HAR / Realized-GARCH series) spx_rv.csv S&P 500 realized-volatility series, 2000–2013 spy_5m.csv SPY 5-minute intraday prices (the realized-measure inputs) spy_rm.csv Realized measures (bipower variation, quarticity) for the HAR extensions

The forecast tournament also reuses garch_scratch.py, sv_scratch.py and svgarch_forecast.py from the GARCH and stochastic-volatility examples (bundled in the example folder).

References

Andersen, T. G., Bollerslev, T., Diebold, F. X. & Labys, P. (2003). Modeling and forecasting realized volatility. Econometrica 71(2), 579–625. — realized variance as a model-free volatility measure
Barndorff-Nielsen, O. E. & Shephard, N. (2002). Econometric analysis of realized volatility and its use in estimating stochastic volatility models. Journal of the Royal Statistical Society: Series B 64(2), 253–280. — the theory of realized variance
Barndorff-Nielsen, O. E. & Shephard, N. (2004). Power and bipower variation with stochastic volatility and jumps. Journal of Financial Econometrics 2(1), 1–37. — bipower variation and the jump component
Corsi, F. (2009). A simple approximate long-memory model of realized volatility. Journal of Financial Econometrics 7(2), 174–196. — the HAR-RV model
Bollerslev, T., Patton, A. J. & Quaedvlieg, R. (2016). Exploiting the errors: a simple approach for improved volatility forecasting. Journal of Econometrics 192(1), 1–18. — HARQ (the realized-quarticity / measurement-error correction)
Hansen, P. R., Huang, Z. & Shek, H. H. (2012). Realized GARCH: a joint model for returns and realized measures of volatility. Journal of Applied Econometrics 27(6), 877–906. — the Realized GARCH bridge
Patton, A. J. (2011). Volatility forecast comparison using imperfect volatility proxies. Journal of Econometrics 160(1), 246–256. — the QLIKE loss robust to a noisy RV proxy

Realized GARCH — Source Code

"""
Realized GARCH (Hansen, Huang & Shek 2012), log-linear specification, from scratch.

Standard GARCH updates the latent variance h_t from the noisy squared return r_{t-1}^2.
Realized GARCH instead lets the informative realized measure RV drive it, and ties the latent h_t
to the observed RV_t through a MEASUREMENT equation with a leverage term:

    return:       r_t = sqrt(h_t) z_t,                          z_t ~ N(0,1)
    GARCH eq:     log h_t = omega + beta log h_{t-1} + gamma log RV_{t-1}
    measurement:  log RV_t = xi + phi log h_t + tau1 z_t + tau2 (z_t^2 - 1) + u_t,   u_t ~ N(0, sigma_u^2)

theta = (omega, beta, gamma, xi, phi, tau1, tau2, log sigma_u).  Effective persistence = beta + phi*gamma.
The likelihood is the JOINT density of returns and realized measures.
"""
import numpy as np
from scipy.optimize import minimize


def _unpack(theta):
    omega, beta, gamma, xi, phi, tau1, tau2, lsu = theta
    return omega, beta, gamma, xi, phi, tau1, tau2, np.exp(lsu)


def filter_logh(theta, r, logRV, logh0):
    omega, beta, gamma, xi, phi, tau1, tau2, su = _unpack(theta)
    T = len(r); logh = np.empty(T); logh[0] = logh0
    for t in range(T - 1):
        logh[t + 1] = omega + beta * logh[t] + gamma * logRV[t]      # RV drives the variance recursion
    return logh


def loglik(theta, r, logRV, logh0):
    omega, beta, gamma, xi, phi, tau1, tau2, su = _unpack(theta)
    if su <= 1e-6 or beta < 0 or beta >= 1 or abs(beta + phi * gamma) >= 1:
        return -np.inf
    logh = filter_logh(theta, r, logRV, logh0)
    h = np.exp(logh); z = r / np.sqrt(h)
    u = logRV - xi - phi * logh - tau1 * z - tau2 * (z ** 2 - 1)
    ll_r = -0.5 * (np.log(2 * np.pi) + logh + r ** 2 / h)           # return density
    ll_x = -0.5 * (np.log(2 * np.pi) + 2 * np.log(su) + u ** 2 / su ** 2)   # measurement density
    ll = ll_r + ll_x
    return np.sum(ll) if np.all(np.isfinite(ll)) else -np.inf


def mle(r, RV):
    """Fit Realized GARCH by joint ML. r = returns, RV = realized variance (same units as r^2)."""
    r = np.asarray(r, float); logRV = np.log(np.asarray(RV, float)); logh0 = float(np.log(np.var(r)))
    x0 = np.array([0.1, 0.55, 0.40, 0.0, 1.0, -0.05, 0.05, np.log(0.4)])
    nll = lambda th: -loglik(th, r, logRV, logh0) if np.isfinite(loglik(th, r, logRV, logh0)) else 1e12
    res = minimize(nll, x0, method="Nelder-Mead", options=dict(xatol=1e-7, fatol=1e-7, maxiter=40000))
    return res.x, logh0


def names():
    return ["omega", "beta", "gamma", "xi", "phi", "tau1", "tau2", "sigma_u"]


def params_readable(theta):
    o, b, g, xi, phi, t1, t2, su = _unpack(theta)
    return dict(omega=o, beta=b, gamma=g, xi=xi, phi=phi, tau1=t1, tau2=t2, sigma_u=su,
                persistence=b + phi * g)