"""H5: HAR forecasts of realized factor variance (5-minute sampled, log) vs GARCH(1,1) on daily FRED factor changes and a 22-day rolling
variance. Horizons 1, 5, 22 days (cumulative variance). OOS 2016-01 → 2026-11, models re-estimated monthly on an expanding window.
Losses: QLIKE and MSE against the realized proxy (sum of 5-min RV over the horizon); Diebold–Mariano with HAC."""
import sys, json, warnings
from pathlib import Path
import numpy as np, pandas as pd
warnings.filterwarnings("ignore")
from arch import arch_model
ROOT = Path(__file__).resolve().parents[1]; CACHE = ROOT / "data/cache"; OUT = ROOT / "results/metrics"
daily = pd.read_parquet(CACHE / "factors_daily_realized.parquet"); daily.index = pd.to_datetime(daily.index)
daily = daily[daily.n_minutes >= 600]
dF = pd.read_csv(OUT / "daily_factors_cmt.csv", index_col=0, parse_dates=True).diff().dropna() * 100  # bp
res = {}
def har_X(lrv):
    return pd.DataFrame({"c": 1.0, "d": lrv, "w": lrv.rolling(5).mean(), "m": lrv.rolling(22).mean()})
def qlike(f, y): return np.log(f) + y / f
def dm_test(l1, l2, h):
    d = np.asarray(l1) - np.asarray(l2); d = d[~np.isnan(d)]; T = len(d)
    lag = max(h - 1, int(np.floor(1.3 * T ** (1 / 2)))) if h > 1 else int(np.floor(1.3 * T ** (1 / 2)))
    mean = d.mean(); var = d.var(ddof=0)
    for k in range(1, lag + 1):
        w = 1 - k / (lag + 1); var += 2 * w * np.mean((d[k:] - mean) * (d[:-k] - mean))
    se = np.sqrt(var / T); return float(mean / se) if se > 0 else np.nan
for fac, k in (("level", "L"), ("slope", "S"), ("curvature", "C")):
    rv = daily[f"RV_{k}_5m"].replace(0, np.nan).dropna() * 1e4  # bp^2
    lrv = np.log(rv)
    X = har_X(lrv)
    x = dF[fac].reindex(rv.index).dropna()
    out = {}
    for h in (1, 5, 22):
        target = rv.rolling(h).sum().shift(-h)  # realized variance over t+1..t+h
        ltarget = np.log(target)
        dates = rv.index[(rv.index >= "2016-01-01")]
        preds = {"har": [], "garch": [], "roll22": [], "y": [], "date": []}
        months = pd.period_range("2016-01", dates[-1].to_period("M"), freq="M")
        garch_params = None
        for m in months:
            t0 = m.start_time
            tr_idx = rv.index[rv.index < t0]
            te_idx = rv.index[(rv.index >= t0) & (rv.index <= m.end_time)]
            # HAR estimated on data available at t0 (target observed only when t+h < t0)
            Xtr = X.loc[tr_idx].dropna(); ytr = ltarget.loc[Xtr.index]; ok = ytr.notna() & (Xtr.index + pd.Timedelta(days=int(h * 1.6)) < t0)
            b = np.linalg.lstsq(Xtr[ok].values, ytr[ok].values, rcond=None)[0]
            resid_var = np.var(ytr[ok].values - Xtr[ok].values @ b)
            # GARCH(1,1) on daily factor changes available before t0
            xs = x.loc[x.index < t0]
            try:
                am = arch_model(xs.values, mean="Constant", vol="GARCH", p=1, q=1, dist="normal", rescale=False)
                fit = am.fit(disp="off", show_warning=False)
                omega, alpha, beta = fit.params["omega"], fit.params["alpha[1]"], fit.params["beta[1]"]
                mu = fit.params["mu"]; cv = fit.conditional_volatility ** 2
                last_h, last_e2 = float(cv[-1]), float((xs.values[-1] - mu) ** 2)
                garch_ok = True
            except Exception:
                garch_ok = False
            for d in te_idx:
                if d not in X.index or pd.isna(target.get(d, np.nan)): continue
                xr = X.loc[d]
                if xr.isna().any(): continue
                f_har = float(np.exp(xr.values @ b + 0.5 * resid_var))  # log-normal correction, one-step; for h>1 target is cumulative already
                # rolling 22-day variance of daily factor changes (bp^2), scaled to h days
                xw = x.loc[:d].tail(22); f_roll = float(xw.var(ddof=0) * h) if len(xw) == 22 else np.nan
                # GARCH forecast: iterate from the last available day <= d
                if garch_ok:
                    xs2 = x.loc[x.index <= d]
                    # recompute conditional variance recursion with fixed params up to d
                    e2 = (xs2.values - mu) ** 2; hcur = float(np.var(xs2.values))
                    for i in range(1, len(e2)): hcur = omega + alpha * e2[i - 1] + beta * hcur
                    fc, tot = hcur, 0.0
                    for j in range(h):
                        fc = omega + (alpha + beta) * fc if j > 0 else omega + alpha * e2[-1] + beta * hcur
                        tot += fc
                    f_garch = float(tot)
                else: f_garch = np.nan
                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)
        Pdf = pd.DataFrame(preds).dropna()
        Pdf.to_parquet(CACHE / f"har_preds_{k}_h{h}.parquet")
        y = Pdf.y.values
        L = {m: {"qlike": float(np.mean(qlike(Pdf[m].values, y))), "mse": float(np.mean((Pdf[m].values - y) ** 2)),
                 "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")}
        L["dm_qlike_har_vs_garch"] = dm_test(qlike(Pdf.har.values, y), qlike(Pdf.garch.values, y), h)
        L["dm_qlike_har_vs_roll22"] = dm_test(qlike(Pdf.har.values, y), qlike(Pdf.roll22.values, y), h)
        L["dm_mse_har_vs_garch"] = dm_test((Pdf.har.values - y) ** 2, (Pdf.garch.values - y) ** 2, h)
        L["n"] = int(len(Pdf)); L["period"] = [str(Pdf.date.min().date()), str(Pdf.date.max().date())]
        out[f"h{h}"] = L
        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)
    # full-sample HAR coefficients (h=1) for the table
    Xf = X.dropna(); yf = lrv.shift(-1).loc[Xf.index].dropna(); Xf = Xf.loc[yf.index]
    bf = np.linalg.lstsq(Xf.values, yf.values, rcond=None)[0]
    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))
    res[fac] = out
json.dump(res, open(OUT / "har_forecasts.json", "w"), indent=1)