Realized volatility — measuring what GARCH and SV only infer

Phase 1: realized measures and the stylized facts of high-frequency volatility

GARCH and SV treat volatility as latent: they infer it from a single daily return. But if we have intraday prices, we can measure a day's variance directly — sum its squared high-frequency returns. That is the realized variance, $$RV_t=\sum_{i=1}^{M} r_{t,i}^2,\qquad r_{t,i}=\text{the $i$-th intraday return on day $t$},$$ which, as the sampling grid gets finer, converges to the day's integrated variance — the quantity GARCH/SV can only estimate. This is the insight of Andersen & Bollerslev (1998), "Answering the Skeptics": judged against a daily squared return (a hopelessly noisy proxy) GARCH looked weak; judged against realized variance it forecasts well.

This notebook builds RV from real SPY intraday data and establishes the stylized facts that the next phase (HAR-RV) will model.

Data: spy_5m.csv / spy_1m.csv — SPY 5-min and 1-min bars pulled from Yahoo Finance (a ~60-day / ~7-day snapshot, for constructing RV); spx_rv.csv — a long series of S&P 500 5-min realized variance, 2000–2013 (from the midasr package, originally the Oxford-Man Realized Library), for the stylized facts and the modeling to come.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy.stats import norm, skew

# --- realized variance from 5-minute SPY prices (a 60-day snapshot) ---
df5 = pd.read_csv("spy_5m.csv", parse_dates=["dt"]); df5["day"] = df5["dt"].dt.date
rows = []; prev = None
for day, g in df5.groupby("day"):
    p = g["close"].values; r = np.diff(np.log(p)) * 100                # intraday 5-min log returns (%)
    rv = np.sum(r ** 2)                                                # realized variance
    bv = (np.pi / 2) * np.sum(np.abs(r[1:]) * np.abs(r[:-1]))          # bipower variation (jump-robust)
    rd = np.log(p[-1] / prev) * 100 if prev is not None else np.nan    # daily close-to-close return
    rows.append((day, rv, bv, rd)); prev = p[-1]
D = pd.DataFrame(rows, columns=["day", "rv", "bv", "rd"]).dropna().reset_index(drop=True)
D["rvol"] = np.sqrt(252 * D["rv"]); D["jump"] = np.maximum(D["rv"] - D["bv"], 0)
print("5-min realized variance, %d trading days:  mean annualized RVol %.1f%%" % (len(D), D["rvol"].mean()))
print("  a single day's r^2 is unbiased but ~%.0f%%-noisy as a vol proxy; 5-min RV (78 obs/day) is only ~%.0f%%"
      % (100 * np.sqrt(2), 100 * np.sqrt(2 / 78)))

fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
hiday = D.loc[D["rvol"].idxmax(), "day"]; g = df5[df5["day"] == hiday]
rr = np.diff(np.log(g["close"].values)) * 100
ax[0].bar(range(len(rr)), rr ** 2, color="steelblue", width=1.0)
axt = ax[0].twinx(); axt.plot(np.cumsum(rr ** 2), color="firebrick", lw=2)
ax[0].set_xlabel("5-min bar within the day"); ax[0].set_ylabel("squared 5-min return", color="steelblue")
axt.set_ylabel("cumulative RV  $\\to$", color="firebrick")
ax[0].set_title("One day (%s): RV = sum of intraday squared returns" % hiday)
ax[1].plot(D["day"], np.abs(D["rd"]) * np.sqrt(252), color="0.7", lw=.9, label="|daily return| (annualized)")
ax[1].plot(D["day"], D["rvol"], color="firebrick", lw=1.8, label="realized vol  $\\sqrt{RV}$ (annualized)")
ax[1].set_ylabel("annualized volatility (%)"); ax[1].legend(fontsize=8)
ax[1].set_title("Realized vol is a PRECISE daily estimate;\na single |return| is a noisy one (Andersen-Bollerslev 1998)")
for lab in ax[1].get_xticklabels(): lab.set_rotation(25)
plt.tight_layout(); plt.show()

5-min realized variance, 59 trading days:  mean annualized RVol 9.4%
  a single day's r^2 is unbiased but ~141%-noisy as a vol proxy; 5-min RV (78 obs/day) is only ~16%

Reading the figure

• Left — how RV is built. On a single day the squared 5-min returns (bars) accumulate into the day's realized variance (red line). Instead of one number per day ($r_t^2$), we now have ~78 intraday observations of variance, and their sum estimates the day's true variance far more precisely.
• Right — why it matters (Andersen–Bollerslev 1998). Both $\sqrt{RV_t}$ (red) and $|r_t|$ (grey) are unbiased estimates of the day's volatility, but $r_t^2$ is a single draw with enormous sampling error (relative error ~$\sqrt2\approx141\%$), while 5-min RV averages 78 observations (relative error ~$\sqrt{2/78}\approx16\%$). The red line is a smooth, trustworthy volatility series; the grey one is the noisy proxy that made people doubt GARCH. This is the whole reason realized volatility exists.

(A note on sampling frequency: sampling *too finely lets bid–ask bounce and other microstructure noise inflate RV — the "volatility signature plot" problem — so 5 minutes is the standard compromise. For an asset as liquid as SPY the noise is already negligible at the minute level, so we simply sample at 5-min.)*

# intraday pattern (diurnal U-shape) and the jump component
df5["r5"] = df5.groupby("day")["close"].transform(lambda s: np.log(s).diff()) * 100
df5["tod"] = df5["dt"].dt.strftime("%H:%M")
uu = df5.dropna(subset=["r5"]).groupby("tod")["r5"].apply(lambda z: np.sqrt(np.mean(z.values ** 2)))
njump = int((D["jump"] / D["rv"] > 0.25).sum())
print("intraday U-shape RMS 5-min return:  open %.3f%%   midday %.3f%%   close %.3f%%"
      % (uu.iloc[0], uu.iloc[len(uu) // 2], uu.iloc[-1]))
print("jumps: %d of %d days have RV-BV > 25%% of RV;  jumps are ~%.1f%% of total variation on average"
      % (njump, len(D), 100 * (D["jump"] / D["rv"]).mean()))

tod = list(uu.index); fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
ax[0].plot(range(len(uu)), uu.values, color="steelblue", lw=2)
ax[0].fill_between(range(len(uu)), uu.values, color="steelblue", alpha=.15)
ticks = [i for i, t in enumerate(tod) if t.endswith(":00")]
ax[0].set_xticks(ticks); ax[0].set_xticklabels([tod[i] for i in ticks])
ax[0].set_xlabel("time of day (ET)"); ax[0].set_ylabel("RMS 5-min return (%)")
ax[0].set_title("Intraday volatility is U-shaped\n(high at the open & close, quiet midday)")
ax[1].fill_between(D["day"], D["bv"], D["rv"], where=(D["rv"] > D["bv"]), color="firebrick", alpha=.35, label="jump variation (RV$-$BV)")
ax[1].plot(D["day"], D["rv"], color="firebrick", lw=1.3, label="RV (total)")
ax[1].plot(D["day"], D["bv"], color="steelblue", lw=1.3, label="bipower variation (continuous)")
ax[1].set_ylabel("daily variation (%$^2$)"); ax[1].legend(fontsize=8)
ax[1].set_title("Bipower variation isolates the continuous part;\nRV $-$ BV flags jumps")
for lab in ax[1].get_xticklabels(): lab.set_rotation(25)
plt.tight_layout(); plt.show()

intraday U-shape RMS 5-min return:  open 0.093%   midday 0.046%   close 0.110%
jumps: 4 of 59 days have RV-BV > 25% of RV;  jumps are ~8.6% of total variation on average

Reading the figure

• Left — intraday structure. Volatility is not uniform through the day: it is high right after the open (overnight news being priced) and into the close, and quiet at midday — the classic U-shape. (The midday spikes are the 2 pm ET data/FOMC releases.) A model of daily RV can ignore this, but it is what the intraday data reveals.
• Right — jumps. Bipower variation (Barndorff-Nielsen–Shephard) multiplies adjacent absolute returns, so a lone large move contributes little — it estimates only the continuous part of variation. The gap $RV-BV$ (shaded) is therefore the jump contribution: discontinuous moves from news. Here ~9% of total variation is jumps, concentrated on a handful of days — the raw material for jump-aware models (HAR-J) later.

# stylized facts on a long series: S&P 500 realized variance, 2000-2013
L = pd.read_csv("spx_rv.csv", parse_dates=["date"]); x = L["rv"].values * 1e4    # daily RV in %^2
L["rvol"] = np.sqrt(252 * L["rv"]) * 100                                          # annualized realized vol
def acf(z, K):
    z = z - z.mean(); v = np.sum(z * z)
    return np.array([np.sum(z[:len(z) - k] * z[k:]) / v for k in range(K + 1)])
a_rv = acf(x, 120); phi = a_rv[1]; logrv = np.log(x)
print("long memory:  ACF(RV) lag1 %.2f  lag22 %.2f  lag100 %.2f   (an AR(1) with phi=%.2f is %.3f at lag100)"
      % (a_rv[1], a_rv[22], a_rv[100], phi, phi ** 100))
print("distribution: skew(RV) %.1f  ->  skew(log RV) %.2f   (log-RV is close to Gaussian)" % (skew(x), skew(logrv)))

fig, ax = plt.subplots(1, 3, figsize=(15, 4.1))
ax[0].plot(L["date"], L["rvol"], color="firebrick", lw=.6)
ax[0].set_ylabel("annualized realized vol (%)"); ax[0].set_title("S&P 500 realized volatility, 2000-2013\n(5-min RV; the 2008 crisis spike)")
lags = np.arange(121)
ax[1].bar(lags, a_rv, color="steelblue", width=1.0, label="ACF of RV")
ax[1].plot(lags, phi ** lags, color="firebrick", lw=2, ls="--", label="AR(1) $\\phi^k$ (exponential)")
ax[1].set_xlabel("lag (days)"); ax[1].set_ylabel("autocorrelation"); ax[1].legend(fontsize=8)
ax[1].set_title("Long memory:\nRV autocorrelation decays slowly")
ax[2].hist(logrv, bins=45, density=True, color="steelblue", alpha=.75)
xx = np.linspace(logrv.min(), logrv.max(), 200)
ax[2].plot(xx, norm.pdf(xx, logrv.mean(), logrv.std()), color="firebrick", lw=2, label="Normal fit")
ax[2].set_xlabel("log RV"); ax[2].set_yticks([]); ax[2].legend(fontsize=8)
ax[2].set_title("log-RV is approximately Gaussian\n(skew %.2f vs %.1f for RV)" % (skew(logrv), skew(x)))
plt.tight_layout(); plt.show()

long memory:  ACF(RV) lag1 0.69  lag22 0.38  lag100 0.11   (an AR(1) with phi=0.69 is 0.000 at lag100)
distribution: skew(RV) 10.5  ->  skew(log RV) 0.51   (log-RV is close to Gaussian)

Reading the figure — and what it means for modeling

• Left — the series. S&P 500 realized volatility over 2000–2013: quiet years, the 2008 spike to ~140% annualized, and a slow return to calm. This is the target series the next phase will forecast.
• Middle — long memory. The autocorrelation of RV decays very slowly — still ~0.1 at a lag of 100 trading days — whereas an AR(1) fit to the lag-1 correlation (red dashed) would have vanished by lag 15. Volatility today is correlated with volatility months ago. A short-memory model (like GARCH's geometric decay) cannot reproduce this; it is the single most important fact about RV.
• Right — log-normality. RV itself is wildly right-skewed (skew ~10), but $\log RV$ is close to Gaussian (skew ~0.5). So volatility models should work with log RV.

These two facts write the next model for us. Long memory + log-normality $\Rightarrow$ regress $\log RV_t$ on its own past over several horizons. That is exactly Corsi's (2009) HAR-RV — a simple linear cascade of daily, weekly, and monthly averaged RV that mimics long memory without fractional integration. Phase 2 builds it three ways (from-scratch / PyMC / R) and pits it against GARCH.

References

• Andersen, T. G. & Bollerslev, T. (1998). Answering the skeptics: yes, standard volatility models do provide accurate forecasts. International Economic Review 39, 885–905.
• Andersen, T. G., Bollerslev, T., Diebold, F. X. & Labys, P. (2003). Modeling and forecasting realized volatility. Econometrica 71, 579–625.
• Barndorff-Nielsen, O. E. & Shephard, N. (2004). Power and bipower variation with stochastic volatility and jumps. J. Financial Econometrics 2, 1–37.
• Corsi, F. (2009). A simple approximate long-memory model of realized volatility. J. Financial Econometrics 7, 174–196.

Data: intraday from Yahoo Finance (spy_5m.csv, spy_1m.csv); long RV series spx_rv.csv (S&P 500 5-min RV 2000–2013, via midasr::rvsp500, originally the Oxford-Man Realized Library).