Python 75.9%
TeX 24%
1#!/usr/bin/env python32"""3================================================================4Auteur : Simon-Pierre Boucher5Contact : contact@spboucher.ai6Projet : Prévision de volatilité réalisée multi-actifs7 (HAR-RV vs GARCH vs Machine Learning)8Fichier : 04_evaluate.py9Description : Étape 04 — Évaluation out-of-sample : assemblage10 des prévisions, combinaisons, pertes QLIKE/MSE,11 Diebold-Mariano vs HAR, Model Confidence Set,12 Mincer-Zarnowitz et forecast encompassing.13Usage : python scripts/04_evaluate.py14================================================================15"""1617from __future__ import annotations1819import logging20import sys21from pathlib import Path2223import numpy as np24import pandas as pd2526sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))2728from wp12 import config, evaluation as ev # noqa: E40229from wp12.models_econ import horizon_target # noqa: E4023031logger = logging.getLogger("04_evaluate")3233FC_DIR = config.path("processed") / "forecasts"34OUT = config.path("reproduced")35HORIZONS = tuple(config.load_config()["evaluation"]["horizons"])36BENCH = "HAR"373839def load_forecasts() -> pd.DataFrame:40 """Concatenate every stored forecast file into one long table."""41 frames = []42 for f in sorted(FC_DIR.glob("*.parquet")):43 df = pd.read_parquet(f)44 frames.append(df)45 out = pd.concat(frames, ignore_index=True)46 out["date"] = pd.to_datetime(out["date"])47 logger.info("loaded %d forecasts, %d models, %d tickers",48 len(out), out["model"].nunique(), out["ticker"].nunique())49 return out505152def attach_actuals(fc: pd.DataFrame, measure: str = "rv5ss") -> pd.DataFrame:53 """Join the realized horizon-average RV target onto the forecast table."""54 panel = pd.read_parquet(config.path("processed") / "panel.parquet")55 targets = []56 for tk, df in panel.groupby("ticker"):57 rv = df.sort_index()[measure]58 for h in HORIZONS:59 t = horizon_target(rv, h).rename("actual").reset_index()60 t["ticker"] = tk61 t["h"] = h62 targets.append(t)63 tgt = pd.concat(targets, ignore_index=True)64 merged = fc.merge(tgt, on=["date", "ticker", "h"], how="inner")65 merged = merged.dropna(subset=["actual", "forecast"])66 merged = merged[merged["actual"] > 0]67 logger.info("merged: %d forecast-actual pairs", len(merged))68 return merged697071def add_combinations(merged: pd.DataFrame) -> pd.DataFrame:72 """Append simple-mean and inverse-MSE combination forecasts."""73 combos = []74 for (tk, h), sub in merged.groupby(["ticker", "h"], observed=True):75 wide = sub.pivot_table(index="date", columns="model", values="forecast")76 actual = sub.drop_duplicates("date").set_index("date")["actual"]77 c = ev.combine_forecasts(wide, actual.reindex(wide.index))78 c = c.stack().rename("forecast").reset_index()79 c.columns = ["date", "model", "forecast"]80 c["ticker"] = tk81 c["h"] = h82 c = c.merge(actual.rename("actual").reset_index(), on="date")83 combos.append(c.dropna(subset=["forecast", "actual"]))84 out = pd.concat([merged] + combos, ignore_index=True)85 logger.info("with combinations: %d rows", len(out))86 return out878889def loss_tables(merged: pd.DataFrame) -> None:90 """Mean losses overall / by class / by ticker, with ratios to HAR."""91 panel_cls = (92 pd.read_parquet(config.path("processed") / "panel.parquet")93 .groupby("ticker")["cls"].first()94 )95 merged = merged.assign(cls=merged["ticker"].map(panel_cls))96 for loss in ("qlike", "mse"):97 lf = ev.LOSSES[loss]98 merged[loss] = lf(merged["actual"].to_numpy(), merged["forecast"].to_numpy())99100 per_asset = merged.groupby(["ticker", "cls", "h", "model"], observed=True)[101 ["qlike", "mse"]].mean().reset_index()102 per_asset.to_csv(OUT / "losses_by_ticker.csv", index=False)103104 def _with_ratio(df: pd.DataFrame, keys: list[str]) -> pd.DataFrame:105 bench = df[df["model"] == BENCH].set_index(keys)["qlike"]106 df = df.set_index(keys)107 df["qlike_ratio"] = df["qlike"] / bench108 return df.reset_index()109110 overall = per_asset.groupby(["h", "model"], observed=True)[111 ["qlike", "mse"]].mean().reset_index()112 _with_ratio(overall, ["h"]).to_csv(OUT / "losses_overall.csv", index=False)113114 by_class = per_asset.groupby(["cls", "h", "model"], observed=True)[115 ["qlike", "mse"]].mean().reset_index()116 _with_ratio(by_class, ["cls", "h"]).to_csv(OUT / "losses_by_class.csv", index=False)117 logger.info("loss tables written")118119120def dm_tables(merged: pd.DataFrame) -> None:121 """Diebold-Mariano tests of every model against the HAR benchmark."""122 rows = []123 for (tk, h), sub in merged.groupby(["ticker", "h"], observed=True):124 wide_f = sub.pivot_table(index="date", columns="model", values="forecast")125 actual = sub.drop_duplicates("date").set_index("date")["actual"]126 if BENCH not in wide_f.columns:127 continue128 lb = ev.qlike(actual.to_numpy(), wide_f[BENCH].to_numpy())129 for m in wide_f.columns:130 if m == BENCH:131 continue132 both = np.isfinite(wide_f[m].to_numpy()) & np.isfinite(wide_f[BENCH].to_numpy())133 la = ev.qlike(actual.to_numpy()[both], wide_f[m].to_numpy()[both])134 stat, p = ev.diebold_mariano(la, lb[both], h=h)135 rows.append((tk, h, m, stat, p))136 dm = pd.DataFrame(rows, columns=["ticker", "h", "model", "dm_stat", "dm_p"])137 dm.to_csv(OUT / "dm_by_ticker.csv", index=False)138 summary = dm.groupby(["h", "model"], observed=True).apply(139 lambda d: pd.Series({140 "n_assets": len(d),141 "median_dm": d["dm_stat"].median(),142 "pct_better": float((d["dm_stat"] < 0).mean()),143 "pct_sig_better": float(((d["dm_stat"] < 0) & (d["dm_p"] < 0.05)).mean()),144 "pct_sig_worse": float(((d["dm_stat"] > 0) & (d["dm_p"] < 0.05)).mean()),145 }), include_groups=False).reset_index()146 summary.to_csv(OUT / "dm_summary.csv", index=False)147 logger.info("DM tables written")148149150def mcs_tables(merged: pd.DataFrame) -> None:151 """90% Model Confidence Set per (ticker, horizon); inclusion shares."""152 cfg = config.load_config()["evaluation"]153 rows = []154 for (tk, h), sub in merged.groupby(["ticker", "h"], observed=True):155 wide_f = sub.pivot_table(index="date", columns="model", values="forecast")156 actual = sub.drop_duplicates("date").set_index("date")["actual"]157 losses = pd.DataFrame({158 m: ev.qlike(actual.to_numpy(), wide_f[m].to_numpy())159 for m in wide_f.columns160 }, index=wide_f.index).dropna()161 if len(losses) < 200:162 continue163 mcs = ev.model_confidence_set(164 losses, alpha=cfg["mcs_alpha"], n_boot=int(cfg["mcs_bootstrap"]))165 mcs["ticker"] = tk166 mcs["h"] = h167 rows.append(mcs)168 allmcs = pd.concat(rows, ignore_index=True)169 allmcs.to_csv(OUT / "mcs_by_ticker.csv", index=False)170 inc = allmcs.groupby(["h", "model"], observed=True)["in_mcs"].mean().rename(171 "mcs_inclusion").reset_index()172 inc.to_csv(OUT / "mcs_inclusion.csv", index=False)173 logger.info("MCS tables written")174175176def mz_and_encompassing(merged: pd.DataFrame) -> None:177 """Mincer-Zarnowitz per model and HAR-vs-ML encompassing tests."""178 rows = []179 for (tk, h, m), sub in merged.groupby(["ticker", "h", "model"], observed=True):180 res = ev.mincer_zarnowitz(sub["actual"].to_numpy(), sub["forecast"].to_numpy(), h=h)181 rows.append((tk, h, m, res["alpha"], res["beta"], res["r2"], res["p_joint"]))182 mz = pd.DataFrame(rows, columns=["ticker", "h", "model", "alpha", "beta", "r2", "p_joint"])183 mz.to_csv(OUT / "mz_by_ticker.csv", index=False)184 mz_sum = mz.groupby(["h", "model"], observed=True).apply(185 lambda d: pd.Series({186 "med_beta": d["beta"].median(), "med_r2": d["r2"].median(),187 "pct_reject": float((d["p_joint"] < 0.05).mean()),188 }), include_groups=False).reset_index()189 mz_sum.to_csv(OUT / "mz_summary.csv", index=False)190191 rows = []192 rivals = [m for m in merged["model"].unique() if m != BENCH]193 for (tk, h), sub in merged.groupby(["ticker", "h"], observed=True):194 wide_f = sub.pivot_table(index="date", columns="model", values="forecast")195 actual = sub.drop_duplicates("date").set_index("date")["actual"]196 if BENCH not in wide_f.columns:197 continue198 for m in rivals:199 if m not in wide_f.columns:200 continue201 res = ev.forecast_encompassing(202 actual.to_numpy(), wide_f[BENCH].to_numpy(), wide_f[m].to_numpy(), h=h)203 rows.append((tk, h, m, res["b1"], res["b2"], res["p_b2"]))204 enc = pd.DataFrame(rows, columns=["ticker", "h", "model", "b1", "b2", "p_b2"])205 enc.to_csv(OUT / "encompassing_by_ticker.csv", index=False)206 enc_sum = enc.groupby(["h", "model"], observed=True).apply(207 lambda d: pd.Series({208 "med_b2": d["b2"].median(),209 "pct_adds_info": float(((d["b2"] > 0) & (d["p_b2"] < 0.05)).mean()),210 }), include_groups=False).reset_index()211 enc_sum.to_csv(OUT / "encompassing_summary.csv", index=False)212 logger.info("MZ and encompassing tables written")213214215def main() -> None:216 """Run the full evaluation battery and persist intermediate tables."""217 config.setup_logging()218 fc = load_forecasts()219 merged = attach_actuals(fc)220 merged = add_combinations(merged)221 merged.to_parquet(config.path("processed") / "merged_forecasts.parquet", index=False)222 loss_tables(merged)223 dm_tables(merged)224 mcs_tables(merged)225 mz_and_encompassing(merged)226 logger.info("evaluation complete")227228229if __name__ == "__main__":230 main()231