HAR-RV — forecasting realized volatility (from scratch)

Corsi's (2009) Heterogeneous Autoregressive model

Phase 1 established two facts about realized volatility: it has long memory (autocorrelation decaying over months) and its logarithm is roughly Gaussian. Corsi's (2009) HAR-RV turns both into a forecasting model of disarming simplicity — an ordinary linear regression of today's log-RV on its own recent averages at three horizons: $$\log RV_t = c + \beta_d\,\overline{\log RV}^{(d)}_{t-1} + \beta_w\,\overline{\log RV}^{(w)}_{t-1} + \beta_m\,\overline{\log RV}^{(m)}_{t-1} + \varepsilon_t,$$ where the regressors are the previous day, the average over the previous week (5 days), and the average over the previous month (22 days). The story — the heterogeneous market hypothesis — is that traders act on different horizons (day traders, weekly rebalancers, monthly allocators) whose overlapping memories sum to the slow decay we saw. No fractional integration, no latent state: three averages and OLS, yet it reproduces long memory and out-forecasts GARCH.

We fit it by OLS and by conjugate Bayesian linear regression, then test it out-of-sample against GARCH and a random walk. Data: spx_rv_ret.csv (S&P 500 5-min RV + daily returns, 2000–2013).

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import garch_scratch as G

df = pd.read_csv("spx_rv_ret.csv", parse_dates=["date"])
RV = df["rv"].values * 1e4              # realized variance in %^2
ret = df["ret"].values; dates = df["date"]; l = np.log(RV)

# --- HAR (log-RV) regressors: previous day, weekly (5d) and monthly (22d) averages ---
idx = np.arange(22, len(l))
d = l[idx - 1]
w = np.array([l[t - 5:t].mean()  for t in idx])
m = np.array([l[t - 22:t].mean() for t in idx])
X = np.column_stack([np.ones(len(idx)), d, w, m]); y = l[idx]
ntr = int(0.6 * len(y)); split_i = idx[ntr]                 # 60% train, rest test (incl. 2008)
Xtr, ytr, Xte, yte = X[:ntr], y[:ntr], X[ntr:], y[ntr:]

# --- OLS + conjugate Bayesian posterior (Normal-Inverse-Gamma, flat prior) ---
XtX_inv = np.linalg.inv(Xtr.T @ Xtr)
beta = XtX_inv @ Xtr.T @ ytr
resid = ytr - Xtr @ beta; dof = len(ytr) - 4; s2 = resid @ resid / dof
se = np.sqrt(np.diag(s2 * XtX_inv)); R2 = 1 - resid.var() / ytr.var()
nm = ["const", "beta_d (daily)", "beta_w (weekly)", "beta_m (monthly)"]
print("HAR log-RV, OLS on train (%d obs), in-sample R^2 %.3f:" % (ntr, R2))
for i in range(4): print("  %-18s %7.3f  (se %.3f)" % (nm[i], beta[i], se[i]))
print("  persistence beta_d+beta_w+beta_m = %.3f" % beta[1:].sum())
rng = np.random.default_rng(0); a_ig, b_ig = dof / 2, resid @ resid / 2
bd = np.array([rng.multivariate_normal(beta, (1 / rng.gamma(a_ig, 1 / b_ig)) * XtX_inv) for _ in range(6000)])

fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
for j, lab, c in zip([1, 2, 3], ["daily $\\beta_d$", "weekly $\\beta_w$", "monthly $\\beta_m$"], ["steelblue", "seagreen", "firebrick"]):
    ax[0].hist(bd[:, j], bins=45, density=True, alpha=.6, color=c, label=lab)
ax[0].axvline(0, color="0.4", lw=1, ls="--"); ax[0].set_yticks([]); ax[0].set_xlabel("coefficient")
ax[0].legend(fontsize=9); ax[0].set_title("HAR coefficients (Bayesian posterior)\ndaily + weekly + monthly all contribute")
ax[1].plot(dates.values[idx], np.sqrt(252 * np.exp(y)), color="0.75", lw=.6, label="actual realized vol")
ax[1].plot(dates.values[idx], np.sqrt(252 * np.exp(X @ beta)), color="firebrick", lw=.8, label="HAR fitted")
ax[1].set_ylabel("annualized vol (%)"); ax[1].legend(fontsize=8)
ax[1].set_title("HAR tracks realized volatility (in + out of sample)\n$R^2=%.2f$ on log RV" % R2)
plt.tight_layout(); plt.show()

HAR log-RV, OLS on train (2062 obs), in-sample R^2 0.673:
  const               -0.027  (se 0.014)
  beta_d (daily)       0.255  (se 0.027)
  beta_w (weekly)      0.491  (se 0.044)
  beta_m (monthly)     0.198  (se 0.036)
  persistence beta_d+beta_w+beta_m = 0.944

Reading the figure

• Left — all three horizons matter. Every coefficient is clearly positive (the posteriors sit well away from zero) and they sum to 0.94 — high persistence, but achieved by spreading memory across horizons rather than one autoregressive term. The weekly component carries the most weight here. This overlapping-horizon structure is what reproduces long memory with a 4-parameter linear model.
• Right — it fits. The HAR line tracks realized volatility across 2000–2013, including the 2008 spike, with $R^2=0.67$ on log RV.

# --- out-of-sample: HAR vs GARCH-t vs random walk (forecasting the RV LEVEL) ---
rv_te = np.exp(yte)                                          # actual RV on the test set
f_har = np.exp(Xte @ beta + s2 / 2)                         # log-HAR mean forecast (lognormal bias correction)
th_g, _, _ = G.garch_mle(ret[:split_i], gjr=False, student=True)      # GARCH-t fit on train returns
r2f, _ = G.prepare(ret); neg = (ret < 0).astype(float)
f_garch = G.garch_var(th_g, r2f, float(np.var(ret[:split_i])), neg, False)[idx[ntr:]]
f_rw = RV[idx[ntr:] - 1]                                     # random walk: yesterday's RV
def qlike(a, f): return np.mean(a / f - np.log(a / f) - 1)   # Patton's robust variance-forecast loss
def rmse(a, f): return np.sqrt(np.mean((a - f) ** 2))
print("Out-of-sample (%d test days, incl. 2008):" % len(yte))
res = {}
for tag, f in [("HAR-RV", f_har), ("GARCH-t", f_garch), ("RW-RV", f_rw)]:
    res[tag] = f; print("  %-9s QLIKE %.4f   RMSE %.2f" % (tag, qlike(rv_te, f), rmse(rv_te, f)))

dte = dates.values[idx[ntr:]]; fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
for t, c in [("GARCH-t", "seagreen"), ("RW-RV", "goldenrod"), ("HAR-RV", "firebrick")]:
    ax[0].plot(dte, np.cumsum(rv_te / res[t] - np.log(rv_te / res[t]) - 1), color=c, lw=1.6, label=t)
ax[0].set_ylabel("cumulative QLIKE loss"); ax[0].legend(fontsize=8)
ax[0].set_title("Out-of-sample forecast loss (lower = better)\nHAR-RV wins")
w08 = (pd.to_datetime(dte) >= "2008-06-01") & (pd.to_datetime(dte) <= "2009-06-01")
ax[1].plot(dte[w08], np.sqrt(252 * rv_te[w08]), color="0.7", lw=1, label="actual RVol")
ax[1].plot(dte[w08], np.sqrt(252 * f_har[w08]), color="firebrick", lw=1.5, label="HAR forecast")
ax[1].plot(dte[w08], np.sqrt(252 * f_garch[w08]), color="seagreen", lw=1.3, label="GARCH-t forecast")
ax[1].set_ylabel("annualized vol (%)"); ax[1].legend(fontsize=8)
ax[1].set_title("2008 crisis: HAR adapts faster than GARCH\n(uses realized info, not just returns)")
for a in ax:
    for lab in a.get_xticklabels(): lab.set_rotation(20)
plt.tight_layout(); plt.show()

Out-of-sample (1375 test days, incl. 2008):
  HAR-RV    QLIKE 0.2231   RMSE 2.67
  GARCH-t   QLIKE 0.2962   RMSE 3.12
  RW-RV     QLIKE 0.3405   RMSE 3.19

Reading the figure — HAR vs GARCH out of sample

• HAR forecasts volatility better than GARCH. On 1,375 out-of-sample days, HAR-RV has the lowest QLIKE loss (0.223 vs GARCH-t's 0.296 and a random walk's 0.341), and its cumulative-loss curve stays below both throughout (left). The reason is information: HAR is built from realized (intraday) variance, while GARCH sees only daily returns.
• HAR adapts faster. In the 2008 crisis (right) HAR rises with the actual volatility and — crucially — comes back down with it in 2009, while the GARCH-t forecast lags on both sides: slow to spike, then stuck high for months. A model fed yesterday's realized variance reacts quicker than one inferring volatility from squared returns.
• Fair-fight caveat. HAR uses more information (intraday) than GARCH (daily returns) — which is the point of realized volatility. A head-to-head that also brings in SV and a common scoring rule is what Phase 3 does; here the message is simply that a trivially simple regression on realized measures beats a fitted GARCH.

Results

• HAR-RV is a 4-parameter OLS that captures long memory and beats GARCH out-of-sample on the S&P 500 — Corsi's (2009) result, reproduced from scratch, with a conjugate Bayesian posterior confirming all three horizons contribute.
• Coefficients are stable and interpretable: a daily/weekly/monthly cascade summing to ~0.94 persistence.

References

• Corsi, F. (2009). A simple approximate long-memory model of realized volatility. J. Financial Econometrics 7, 174–196.
• Andersen, Bollerslev, Diebold & Labys (2003). Modeling and forecasting realized volatility. Econometrica 71, 579–625.
• Patton, A. (2011). Volatility forecast comparison using imperfect volatility proxies. J. Econometrics — the QLIKE loss.

Companion notebooks: har_rv_pymc.ipynb (Bayesian HAR in PyMC) and har_rv_R.ipynb (highfrequency::HARmodel).