1"""H5: HAR forecasts of realized factor variance (5-minute sampled, log) vs GARCH(1,1) on daily FRED factor changes and a 22-day rolling2variance. Horizons 1, 5, 22 days (cumulative variance). OOS 2016-01 → 2026-11, models re-estimated monthly on an expanding window.3Losses: QLIKE and MSE against the realized proxy (sum of 5-min RV over the horizon); Diebold–Mariano with HAC."""4import sys, json, warnings5from pathlib import Path6import numpy as np, pandas as pd7warnings.filterwarnings("ignore")8from arch import arch_model9ROOT = Path(__file__).resolve().parents[1]; CACHE = ROOT / "data/cache"; OUT = ROOT / "results/metrics"10daily = pd.read_parquet(CACHE / "factors_daily_realized.parquet"); daily.index = pd.to_datetime(daily.index)11daily = daily[daily.n_minutes >= 600]12dF = pd.read_csv(OUT / "daily_factors_cmt.csv", index_col=0, parse_dates=True).diff().dropna() * 100 # bp13res = {}14def har_X(lrv):15 return pd.DataFrame({"c": 1.0, "d": lrv, "w": lrv.rolling(5).mean(), "m": lrv.rolling(22).mean()})16def qlike(f, y): return np.log(f) + y / f17def dm_test(l1, l2, h):18 d = np.asarray(l1) - np.asarray(l2); d = d[~np.isnan(d)]; T = len(d)19 lag = max(h - 1, int(np.floor(1.3 * T ** (1 / 2)))) if h > 1 else int(np.floor(1.3 * T ** (1 / 2)))20 mean = d.mean(); var = d.var(ddof=0)21 for k in range(1, lag + 1):22 w = 1 - k / (lag + 1); var += 2 * w * np.mean((d[k:] - mean) * (d[:-k] - mean))23 se = np.sqrt(var / T); return float(mean / se) if se > 0 else np.nan24for fac, k in (("level", "L"), ("slope", "S"), ("curvature", "C")):25 rv = daily[f"RV_{k}_5m"].replace(0, np.nan).dropna() * 1e4 # bp^226 lrv = np.log(rv)27 X = har_X(lrv)28 x = dF[fac].reindex(rv.index).dropna()29 out = {}30 for h in (1, 5, 22):31 target = rv.rolling(h).sum().shift(-h) # realized variance over t+1..t+h32 ltarget = np.log(target)33 dates = rv.index[(rv.index >= "2016-01-01")]34 preds = {"har": [], "garch": [], "roll22": [], "y": [], "date": []}35 months = pd.period_range("2016-01", dates[-1].to_period("M"), freq="M")36 garch_params = None37 for m in months:38 t0 = m.start_time39 tr_idx = rv.index[rv.index < t0]40 te_idx = rv.index[(rv.index >= t0) & (rv.index <= m.end_time)]41 # HAR estimated on data available at t0 (target observed only when t+h < t0)42 Xtr = X.loc[tr_idx].dropna(); ytr = ltarget.loc[Xtr.index]; ok = ytr.notna() & (Xtr.index + pd.Timedelta(days=int(h * 1.6)) < t0)43 b = np.linalg.lstsq(Xtr[ok].values, ytr[ok].values, rcond=None)[0]44 resid_var = np.var(ytr[ok].values - Xtr[ok].values @ b)45 # GARCH(1,1) on daily factor changes available before t046 xs = x.loc[x.index < t0]47 try:48 am = arch_model(xs.values, mean="Constant", vol="GARCH", p=1, q=1, dist="normal", rescale=False)49 fit = am.fit(disp="off", show_warning=False)50 omega, alpha, beta = fit.params["omega"], fit.params["alpha[1]"], fit.params["beta[1]"]51 mu = fit.params["mu"]; cv = fit.conditional_volatility ** 252 last_h, last_e2 = float(cv[-1]), float((xs.values[-1] - mu) ** 2)53 garch_ok = True54 except Exception:55 garch_ok = False56 for d in te_idx:57 if d not in X.index or pd.isna(target.get(d, np.nan)): continue58 xr = X.loc[d]59 if xr.isna().any(): continue60 f_har = float(np.exp(xr.values @ b + 0.5 * resid_var)) # log-normal correction, one-step; for h>1 target is cumulative already61 # rolling 22-day variance of daily factor changes (bp^2), scaled to h days62 xw = x.loc[:d].tail(22); f_roll = float(xw.var(ddof=0) * h) if len(xw) == 22 else np.nan63 # GARCH forecast: iterate from the last available day <= d64 if garch_ok:65 xs2 = x.loc[x.index <= d]66 # recompute conditional variance recursion with fixed params up to d67 e2 = (xs2.values - mu) ** 2; hcur = float(np.var(xs2.values))68 for i in range(1, len(e2)): hcur = omega + alpha * e2[i - 1] + beta * hcur69 fc, tot = hcur, 0.070 for j in range(h):71 fc = omega + (alpha + beta) * fc if j > 0 else omega + alpha * e2[-1] + beta * hcur72 tot += fc73 f_garch = float(tot)74 else: f_garch = np.nan75 preds["har"].append(f_har); preds["garch"].append(f_garch); preds["roll22"].append(f_roll); preds["y"].append(float(target[d])); preds["date"].append(d)76 Pdf = pd.DataFrame(preds).dropna()77 Pdf.to_parquet(CACHE / f"har_preds_{k}_h{h}.parquet")78 y = Pdf.y.values79 L = {m: {"qlike": float(np.mean(qlike(Pdf[m].values, y))), "mse": float(np.mean((Pdf[m].values - y) ** 2)),80 "r2_log": float(1 - np.sum((np.log(Pdf[m].values) - np.log(y)) ** 2) / np.sum((np.log(y) - np.log(y).mean()) ** 2))} for m in ("har", "garch", "roll22")}81 L["dm_qlike_har_vs_garch"] = dm_test(qlike(Pdf.har.values, y), qlike(Pdf.garch.values, y), h)82 L["dm_qlike_har_vs_roll22"] = dm_test(qlike(Pdf.har.values, y), qlike(Pdf.roll22.values, y), h)83 L["dm_mse_har_vs_garch"] = dm_test((Pdf.har.values - y) ** 2, (Pdf.garch.values - y) ** 2, h)84 L["n"] = int(len(Pdf)); L["period"] = [str(Pdf.date.min().date()), str(Pdf.date.max().date())]85 out[f"h{h}"] = L86 print(fac, h, {m: round(L[m]["qlike"], 4) for m in ("har", "garch", "roll22")}, "DM har-garch", round(L["dm_qlike_har_vs_garch"], 2), "DM har-roll", round(L["dm_qlike_har_vs_roll22"], 2), flush=True)87 # full-sample HAR coefficients (h=1) for the table88 Xf = X.dropna(); yf = lrv.shift(-1).loc[Xf.index].dropna(); Xf = Xf.loc[yf.index]89 bf = np.linalg.lstsq(Xf.values, yf.values, rcond=None)[0]90 out["har_coef_h1_full"] = dict(zip(["c", "d", "w", "m"], bf.round(4))); out["har_r2_h1_full"] = float(1 - np.var(yf.values - Xf.values @ bf) / np.var(yf.values))91 res[fac] = out92json.dump(res, open(OUT / "har_forecasts.json", "w"), indent=1)93