HAR extensions — jumps and measurement error¶
HAR-CJ (jumps, Andersen-Bollerslev-Diebold 2007) and HARQ (realized quarticity, Bollerslev-Patton-Quaedvlieg 2016)¶
Plain HAR treats yesterday's realized variance as a single, perfectly-measured number. Two refinements relax that:
- HAR-CJ splits RV into a smooth continuous part and discontinuous jumps (using bipower variation), on the hypothesis that the two have different dynamics.
- HARQ recognises that RV is itself an estimate with sampling error — larger when the day's realized quarticity ($RQ$) is high — and lets the model downweight a noisily-measured RV: the daily coefficient becomes $\beta_d+\beta_{dQ}\sqrt{RQ_{t-1}}$.
Both need intraday measures beyond RV — bipower variation (jumps) and realized quarticity (measurement error). We use highfrequency::SPYRM: daily SPY realized measures (RV, bipower, RQ), 2014–2019, exported to spy_rm.csv.
In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt
df = pd.read_csv("spy_rm.csv", parse_dates=["date"])
# NB: in SPYRM, RV/BPV are in decimal^2 but RQ is in %^4 -> put everything in % units
RV = df["rv"].values * 1e4; BPV = df["bpv"].values * 1e4; RQ = df["rq"].values # %^2, %^2, %^4
dates = df["date"]; n = len(RV)
# --- significant-jump decomposition (Barndorff-Nielsen-Shephard z-test, using RQ) ---
mu1 = np.sqrt(2 / np.pi); theta = mu1 ** -4 + 2 * mu1 ** -2 - 5; M = 78
RJ = (RV - BPV) / RV
Z = np.sqrt(M) * RJ / np.sqrt(theta * np.maximum(1.0, RQ / BPV ** 2))
Jsig = Z > 3.09 # significance level 0.999
J = np.where(Jsig, np.maximum(RV - BPV, 0.0), 0.0) # jump variation (significant only)
C = RV - J # continuous variation
print("significant jumps on %d of %d days (%.1f%%); jumps ~%.1f%% of total variation"
% (Jsig.sum(), n, 100 * Jsig.mean(), 100 * J.sum() / RV.sum()))
idx = np.arange(22, n)
def casc(s): return s[idx-1], np.array([s[t-5:t].mean() for t in idx]), np.array([s[t-22:t].mean() for t in idx])
RVd, RVw, RVm = casc(RV); Cd, Cw, Cm = casc(C); Jd, Jw, Jm = casc(J)
y = RV[idx]; ntr = int(0.6 * len(y)); one = np.ones(len(idx))
sq = np.sqrt(RQ)[idx - 1]; sqz = sq / sq[:ntr].mean() # sqrt(RQ), standardized (1 = average noise)
def ols(X, yy):
b = np.linalg.solve(X.T @ X, X.T @ yy); r = yy - X @ b
se = np.sqrt(np.diag((r @ r / (len(yy) - X.shape[1])) * np.linalg.inv(X.T @ X)))
return b, se
X_har = np.column_stack([one, RVd, RVw, RVm]) # plain HAR (levels)
X_cj = np.column_stack([one, Cd, Cw, Cm, Jd]) # continuous cascade + daily jump
X_q = np.column_stack([one, RVd, (sqz - 1) * RVd, RVw, RVm]) # HARQ (daily RQ interaction)
b_har, _ = ols(X_har[:ntr], y[:ntr]); b_cj, se_cj = ols(X_cj[:ntr], y[:ntr]); b_q, se_q = ols(X_q[:ntr], y[:ntr])
print("HAR-CJ: continuous C_d %.2f C_w %.2f C_m %.2f | jump J_d %.2f (se %.2f)"
% (b_cj[1], b_cj[2], b_cj[3], b_cj[4], se_cj[4]))
print("HARQ: beta_d %.3f beta_dQ %.3f (se %.3f) [<0 => trust a noisy RV less]" % (b_q[1], b_q[2], se_q[2]))
def acf(z, K): z = z - z.mean(); v = z @ z; return np.array([z[:len(z)-k] @ z[k:] / v for k in range(K+1)])
fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
ax[0].plot(dates, np.sqrt(252 * RV), color="0.7", lw=.5, label="realized vol")
ax[0].scatter(dates[Jsig], np.sqrt(252 * RV)[Jsig], color="firebrick", s=18, zorder=3, label="significant jump day")
ax[0].set_ylabel("annualized vol (%)"); ax[0].legend(fontsize=8)
ax[0].set_title("Significant jumps are rare and spiky\n(SPY 2014-2019)")
lags = np.arange(41)
ax[1].plot(lags, acf(C, 40), "o-", color="steelblue", ms=3, label="continuous variation")
ax[1].plot(lags, acf(J, 40), "s-", color="firebrick", ms=3, label="jump variation")
ax[1].axhline(0, color="0.5", lw=.8); ax[1].set_xlabel("lag (days)"); ax[1].set_ylabel("autocorrelation")
ax[1].legend(fontsize=8); ax[1].set_title("Continuous variation persists; jumps do not\n(so jumps don't help forecast)")
plt.tight_layout(); plt.show()
significant jumps on 14 of 1495 days (0.9%); jumps ~0.8% of total variation HAR-CJ: continuous C_d 0.21 C_w 0.24 C_m 0.20 | jump J_d 0.65 (se 0.94) HARQ: beta_d 1.169 beta_dQ -0.107 (se 0.012) [<0 => trust a noisy RV less]
Reading the figure — jumps don't help forecast¶
- Left — jumps are rare. The significance test flags a jump on only 14 of 1,495 days (~1%): a handful of discrete events (the Aug-2015 flash crash, Feb-2018 "volmageddon").
- Right — and transitory. Continuous variation's autocorrelation decays slowly (persistent, like RV), but the jump autocorrelation collapses to zero after one day — jumps do not cluster. Since forecasting is about persistence, jumps carry almost no information about future volatility; in HAR-CJ the daily jump coefficient (0.65) is statistically indistinguishable from zero (se 0.94).
In [2]:
def qlike(a, f): f = np.maximum(f, 1e-6); return np.mean(a / f - np.log(a / f) - 1)
rv_te = y[ntr:]
q = {"HAR": qlike(rv_te, X_har[ntr:] @ b_har), "HAR-CJ": qlike(rv_te, X_cj[ntr:] @ b_cj),
"HARQ": qlike(rv_te, X_q[ntr:] @ b_q)}
print("Out-of-sample QLIKE (%d days): " % len(rv_te) + " ".join("%s %.4f" % (k, v) for k, v in q.items()))
fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
g = np.linspace(np.percentile(sqz, 2), np.percentile(sqz, 98), 100)
ax[0].plot(g, b_q[1] + b_q[2] * (g - 1), color="slateblue", lw=2, label="HARQ effective daily coef")
ax[0].axhline(b_har[1], color="0.6", ls="--", lw=1.2, label="plain HAR daily coef")
ax[0].set_xlabel("$\\sqrt{RQ_{t-1}}$ (relative; 1 = average noise)"); ax[0].set_ylabel("effective daily coefficient")
ax[0].legend(fontsize=8); ax[0].set_title("HARQ trusts yesterday's RV less when it is\nnoisily measured (high realized quarticity)")
kk = list(q); ax[1].bar(kk, [q[k] for k in kk], color=["steelblue", "seagreen", "firebrick"])
for i, k in enumerate(kk): ax[1].text(i, q[k] + 0.0015, "%.4f" % q[k], ha="center", fontsize=10)
ax[1].set_ylabel("out-of-sample QLIKE (lower = better)"); ax[1].set_ylim(min(q.values()) * 0.92, max(q.values()) * 1.05)
ax[1].set_title("Forecast loss: HARQ helps, jumps don't")
plt.tight_layout(); plt.show()
Out-of-sample QLIKE (590 days): HAR 0.2917 HAR-CJ 0.2916 HARQ 0.2490
Reading the figure — HARQ helps¶
- Left — the mechanism. The fitted $\beta_{dQ}=-0.11$ is strongly negative ($t\approx-9$): the effective weight on yesterday's RV shrinks when realized quarticity is high (RV noisy) and grows when RV is cleanly measured. HARQ also weights recent RV more aggressively overall than plain HAR (whose single coefficient is dragged down by the noisy days it cannot distinguish). Trust a precise measurement; discount a noisy one.
- Right — out of sample. HARQ cuts QLIKE from 0.292 to 0.249 (~15% better), while HAR-CJ is essentially tied with plain HAR (0.2916 vs 0.2917). The verdict matches the literature: modeling RV's measurement error pays; separating jumps does not improve point volatility forecasts — because the jump part is unpredictable anyway.
Results¶
- Jumps: real but useless for forecasting. ~1% of days jump, but jumps don't persist, so a jump-aware HAR forecasts no better than plain HAR.
- Measurement error: worth modeling. HARQ's adaptive downweighting of noisy RV is the single most useful HAR refinement here — ~15% lower out-of-sample QLIKE.
References¶
- Andersen, Bollerslev & Diebold (2007). Roughing it up: including jump components in the measurement, modeling, and forecasting of return volatility. Review of Economics and Statistics.
- Barndorff-Nielsen & Shephard (2006). Econometrics of testing for jumps in financial economics using bipower variation. J. Financial Econometrics.
- Bollerslev, Patton & Quaedvlieg (2016). Exploiting the errors: a simple approach for improved volatility forecasting. J. Econometrics 192, 1–18.
Data: SPY 2014–2019 (spy_rm.csv, from highfrequency::SPYRM), the period with bipower + quarticity available. In R, highfrequency::HARmodel(..., type = "HARQ") fits HARQ directly.