HAR-RV in PyMC — the Bayesian view¶

The PPL corner of the HAR trio (from-scratch → PyMC → R)¶

HAR is an ordinary linear regression, so there is no scan and no latent state — PyMC's NUTS samples it in seconds. What Bayes adds over the from-scratch OLS is twofold: the full joint posterior of the coefficients (which should reproduce the conjugate result exactly), and a clean model comparison — is the daily/weekly/monthly cascade really justified, or would a single daily lag do? We answer that with an information criterion (LOO) on the two nested models.

In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import pymc as pm, arviz as az

df = pd.read_csv("spx_rv_ret.csv", parse_dates=["date"])
l = np.log(df["rv"].values * 1e4)                              # log realized variance (%^2)
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])
y = l[idx]; ntr = int(0.6 * len(y))                            # same train split as the from-scratch notebook
d, w, m, y = d[:ntr], w[:ntr], m[:ntr], y[:ntr]

with pm.Model() as har:                                        # HAR is a plain linear regression
    c = pm.Normal("c", 0, 1)
    bd = pm.Normal("bd", 0, 1); bw = pm.Normal("bw", 0, 1); bm = pm.Normal("bm", 0, 1)
    sig = pm.HalfNormal("sigma", 1)
    pm.Normal("y", mu=c + bd * d + bw * w + bm * m, sigma=sig, observed=y)
    idata_har = pm.sample(1500, tune=1000, chains=2, cores=1, target_accept=0.9,
                          random_seed=1, idata_kwargs={"log_likelihood": True})

print("PyMC HAR posterior (train):")
print(az.summary(idata_har, var_names=["c", "bd", "bw", "bm", "sigma"])[["mean", "sd"]].to_string())
print("cross-check from-scratch OLS: bd 0.255  bw 0.491  bm 0.198")

PyMC HAR posterior (train):
          mean      sd
c      -0.0269  0.0131
bd       0.255  0.0264
bw        0.49   0.044
bm       0.199   0.036
sigma   0.5408  0.0083
cross-check from-scratch OLS: bd 0.255  bw 0.491  bm 0.198

The posterior means land on the OLS estimates (as they must for a linear-Gaussian model with flat priors) — a reassuring end-to-end check. Now the Bayesian value-add: is the cascade justified? Fit the nested daily-only model and compare.

In [2]:
with pm.Model() as ar1:                                        # nested: drop weekly + monthly
    c = pm.Normal("c", 0, 1); bd = pm.Normal("bd", 0, 1); sig = pm.HalfNormal("sigma", 1)
    pm.Normal("y", mu=c + bd * d, sigma=sig, observed=y)
    idata_ar1 = pm.sample(1500, tune=1000, chains=2, cores=1, target_accept=0.9,
                          random_seed=1, idata_kwargs={"log_likelihood": True})

cmp = az.compare({"HAR (d+w+m)": idata_har, "AR1 (daily only)": idata_ar1})
ecol = [c for c in cmp.columns if c.startswith("elpd") and c != "elpd_diff"][0]
print("Information-criterion comparison (higher elpd = better):")
print(cmp[["rank", ecol, "elpd_diff", "dse"]].to_string())

fig, ax = plt.subplots(1, 2, figsize=(13, 4.6))
post = idata_har.posterior
ols = {"bd": 0.255, "bw": 0.491, "bm": 0.198}
cols = {"bd": "steelblue", "bw": "seagreen", "bm": "firebrick"}
labs = {"bd": "daily $\\beta_d$", "bw": "weekly $\\beta_w$", "bm": "monthly $\\beta_m$"}
for k in ["bd", "bw", "bm"]:
    ax[0].hist(post[k].values.ravel(), bins=45, density=True, alpha=.55, color=cols[k], label=labs[k])
    ax[0].axvline(ols[k], color=cols[k], ls="--", lw=1.5)
ax[0].set_yticks([]); ax[0].set_xlabel("coefficient"); ax[0].legend(fontsize=9)
ax[0].set_title("PyMC/NUTS HAR posterior\n(dashed = from-scratch OLS — they agree)")
names = list(cmp.index); yb = np.arange(len(names))[::-1]
ax[1].errorbar(cmp[ecol].values, yb, xerr=cmp["se"].values, fmt="o", color="slateblue", capsize=4, ms=8)
ax[1].set_yticks(yb); ax[1].set_yticklabels(names); ax[1].set_ylim(-0.5, len(names) - 0.5)
ax[1].set_xlabel("ELPD (%s)  — higher is better" % ecol.split("_")[-1])
ax[1].set_title("The weekly+monthly cascade is justified\n(HAR beats a daily-only AR(1))")
plt.tight_layout(); plt.show()

Information-criterion comparison (higher elpd = better):
                  rank    elpd  elpd_diff   dse
HAR (d+w+m)          0 -1660.0        0.0   0.0
AR1 (daily only)     1 -1890.0     -230.0  22.0
No description has been provided for this image

Results¶

  • NUTS reproduces the conjugate posterior. The PyMC coefficient posteriors sit exactly on the from-scratch OLS estimates (dashed lines) — as they must for a linear-Gaussian model with flat priors, but a clean end-to-end confirmation that both engines agree.
  • The cascade is justified. The full HAR (daily + weekly + monthly) has an expected log predictive density about 230 higher than a daily-only AR(1) — roughly ten times its standard error. Dropping the weekly and monthly horizons costs real predictive accuracy: the heterogeneous structure earns its keep; it is not decoration.

The from-scratch notebook showed HAR forecasts better than GARCH out of sample; the Bayesian comparison here shows why the HAR form itself is right — each horizon carries genuine information about tomorrow's volatility.

References¶

  • Corsi, F. (2009). A simple approximate long-memory model of realized volatility. J. Financial Econometrics.
  • Vehtari, A., Gelman, A. & Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing.

Companion: har_rv_R.ipynb (highfrequency::HARmodel).