Core pipeline: API client, RV measures, econ/ML/eval/robustness/backtest modules, tests

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Simon-Pierre Boucher committed 1 mo ago (Aug 10, 2026) parent 28bfb36

16 changed files +3,014 −1

modified config.yaml +1 −1

		`@@ -14,7 +14,7 @@ api:`
14	14	max_retries: 5
15	15	backoff_base_sec: 1.0 # exponential backoff: base * 2**attempt
16	16	timeout_sec: 60
17		− page_limit: 50000 # max rows per request page
	17	+ page_limit: 100000 # max rows per request page (validated live)
18	18	cache_dir: "data/raw" # parquet cache, one file per (asset, ticker, timeframe)
19	19
20	20	# ---------------------------------------------------------------- universe

added scripts/01_download_data.py +48 −0

		`@@ -0,0 +1,48 @@`
	1	+#!/usr/bin/env python3
	2	+"""
	3	+================================================================
	4	+Auteur : Simon-Pierre Boucher
	5	+Contact : contact@spboucher.ai
	6	+Projet : Prévision de volatilité réalisée multi-actifs
	7	+ (HAR-RV vs GARCH vs Machine Learning)
	8	+Fichier : 01_download_data.py
	9	+Description : Étape 01 — Téléchargement des barres 1-min et daily
	10	+ pour tout l'univers (cache parquet data/raw).
	11	+Usage : python scripts/01_download_data.py [--refresh]
	12	+================================================================
	13	+"""
	14	+
	15	+from __future__ import annotations
	16	+
	17	+import argparse
	18	+import logging
	19	+import sys
	20	+from pathlib import Path
	21	+
	22	+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
	23	+
	24	+from wp12 import config, data_pipeline # noqa: E402
	25	+from wp12.api_client import HFMarketDataClient # noqa: E402
	26	+
	27	+logger = logging.getLogger("01_download")
	28	+
	29	+
	30	+def main() -> None:
	31	+ """Download 1-min and daily bars for every instrument in the universe."""
	32	+ ap = argparse.ArgumentParser()
	33	+ ap.add_argument("--refresh", action="store_true", help="force re-download")
	34	+ args = ap.parse_args()
	35	+
	36	+ config.setup_logging()
	37	+ logging.getLogger("urllib3").setLevel(logging.WARNING)
	38	+ client = HFMarketDataClient()
	39	+ todo = config.universe(include_predictors=True)
	40	+ for i, inst in enumerate(todo, 1):
	41	+ logger.info("[%d/%d] %s (%s)", i, len(todo), inst["ticker"], inst["asset"])
	42	+ bars = data_pipeline.download_instrument(client, inst, refresh=args.refresh)
	43	+ for tf, df in bars.items():
	44	+ logger.info(" %s: %d rows", tf, len(df))
	45	+
	46	+
	47	+if __name__ == "__main__":
	48	+ main()

added scripts/02_build_rv.py +43 −0

		`@@ -0,0 +1,43 @@`
	1	+#!/usr/bin/env python3
	2	+"""
	3	+================================================================
	4	+Auteur : Simon-Pierre Boucher
	5	+Contact : contact@spboucher.ai
	6	+Projet : Prévision de volatilité réalisée multi-actifs
	7	+ (HAR-RV vs GARCH vs Machine Learning)
	8	+Fichier : 02_build_rv.py
	9	+Description : Étape 02 — Construction des mesures de volatilité
	10	+ réalisée journalières et du panel multi-actifs.
	11	+Usage : python scripts/02_build_rv.py [--refresh]
	12	+================================================================
	13	+"""
	14	+
	15	+from __future__ import annotations
	16	+
	17	+import argparse
	18	+import logging
	19	+import sys
	20	+from pathlib import Path
	21	+
	22	+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
	23	+
	24	+from wp12 import config, data_pipeline # noqa: E402
	25	+
	26	+logger = logging.getLogger("02_build_rv")
	27	+
	28	+
	29	+def main() -> None:
	30	+ """Build daily realized measures for the universe and save the panel."""
	31	+ ap = argparse.ArgumentParser()
	32	+ ap.add_argument("--refresh", action="store_true", help="recompute RV files")
	33	+ args = ap.parse_args()
	34	+
	35	+ config.setup_logging()
	36	+ logging.getLogger("urllib3").setLevel(logging.WARNING)
	37	+ panel = data_pipeline.build_panel(refresh=args.refresh)
	38	+ logger.info("panel summary:\n%s",
	39	+ panel.groupby("ticker")["rv5ss"].agg(["count", "mean"]))
	40	+
	41	+
	42	+if __name__ == "__main__":
	43	+ main()

added scripts/03_run_models.py +153 −0

		`@@ -0,0 +1,153 @@`
	1	+#!/usr/bin/env python3
	2	+"""
	3	+================================================================
	4	+Auteur : Simon-Pierre Boucher
	5	+Contact : contact@spboucher.ai
	6	+Projet : Prévision de volatilité réalisée multi-actifs
	7	+ (HAR-RV vs GARCH vs Machine Learning)
	8	+Fichier : 03_run_models.py
	9	+Description : Étape 03 — Prévisions out-of-sample roulantes pour
	10	+ toutes les familles de modèles (HAR, GARCH, ML,
	11	+ deep, pooled), parallélisées par actif.
	12	+Usage : python scripts/03_run_models.py --family econ\|ml\|deep\|pooled\|all
	13	+ [--tickers SPY,BTC] [--workers 8]
	14	+================================================================
	15	+"""
	16	+
	17	+from __future__ import annotations
	18	+
	19	+import argparse
	20	+import logging
	21	+import sys
	22	+from concurrent.futures import ProcessPoolExecutor, as_completed
	23	+from pathlib import Path
	24	+
	25	+import pandas as pd
	26	+
	27	+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
	28	+
	29	+from wp12 import config # noqa: E402
	30	+
	31	+logger = logging.getLogger("03_run_models")
	32	+
	33	+FC_DIR = config.path("processed") / "forecasts"
	34	+
	35	+HORIZONS = tuple(config.load_config()["evaluation"]["horizons"])
	36	+WINDOW = int(config.load_config()["evaluation"]["estimation_window_days"])
	37	+REFIT = int(config.load_config()["evaluation"]["reestimation_freq_days"])
	38	+
	39	+
	40	+def _load_rv(ticker: str) -> pd.DataFrame:
	41	+ return pd.read_parquet(config.path("processed") / f"rv_{ticker}.parquet")
	42	+
	43	+
	44	+def run_econ(ticker: str) -> str:
	45	+ """All HAR-family and GARCH-family forecasts for one ticker."""
	46	+ from wp12 import models_econ as me
	47	+
	48	+ df = _load_rv(ticker)
	49	+ parts = []
	50	+ for spec in ["HAR", "HAR-J", "HAR-CJ", "SHAR", "HARQ", "LogHAR"]:
	51	+ parts.append(me.rolling_har_forecasts(
	52	+ df, spec, horizons=HORIZONS, window=WINDOW, refit_every=REFIT))
	53	+ for kind in ["GARCH", "GJR", "EGARCH"]:
	54	+ parts.append(me.rolling_garch_forecasts(
	55	+ df, kind, horizons=HORIZONS, window=WINDOW, refit_every=REFIT))
	56	+ parts.append(me.rolling_realgarch_forecasts(
	57	+ df, horizons=HORIZONS, window=WINDOW, refit_every=REFIT))
	58	+ out = pd.concat(parts, ignore_index=True)
	59	+ out["ticker"] = ticker
	60	+ f = FC_DIR / f"econ_{ticker}.parquet"
	61	+ out.to_parquet(f, index=False)
	62	+ return f"{ticker}: {len(out)} econ forecasts"
	63	+
	64	+
	65	+def run_ml(ticker: str) -> str:
	66	+ """Tabular ML forecasts (2 feature sets x 6 estimators) for one ticker."""
	67	+ from wp12 import models_ml as ml
	68	+
	69	+ panel = pd.read_parquet(config.path("processed") / "panel.parquet")
	70	+ frames = ml.ml_features(panel)
	71	+ frame = frames[ticker]
	72	+ parts = []
	73	+ for name in ["LASSO", "Ridge", "ElasticNet", "RF", "XGBoost", "LightGBM"]:
	74	+ for fs in ["har", "extended"]:
	75	+ parts.append(ml.rolling_ml_forecasts(
	76	+ frame, name, feature_set=fs, horizons=HORIZONS,
	77	+ window=WINDOW, refit_every=REFIT))
	78	+ out = pd.concat(parts, ignore_index=True)
	79	+ out["ticker"] = ticker
	80	+ f = FC_DIR / f"ml_{ticker}.parquet"
	81	+ out.to_parquet(f, index=False)
	82	+ return f"{ticker}: {len(out)} ML forecasts"
	83	+
	84	+
	85	+def run_deep(ticker: str) -> str:
	86	+ """LSTM and Transformer forecasts for one ticker."""
	87	+ from wp12 import models_ml as ml
	88	+
	89	+ panel = pd.read_parquet(config.path("processed") / "panel.parquet")
	90	+ frames = ml.ml_features(panel)
	91	+ frame = frames[ticker]
	92	+ parts = [
	93	+ ml.rolling_deep_forecasts(frame, arch, horizons=HORIZONS, window=WINDOW)
	94	+ for arch in ["LSTM", "Transformer"]
	95	+ ]
	96	+ out = pd.concat(parts, ignore_index=True)
	97	+ out["ticker"] = ticker
	98	+ f = FC_DIR / f"deep_{ticker}.parquet"
	99	+ out.to_parquet(f, index=False)
	100	+ return f"{ticker}: {len(out)} deep forecasts"
	101	+
	102	+
	103	+def run_pooled() -> str:
	104	+ """Pooled multi-asset LightGBM (single job, all tickers jointly)."""
	105	+ from wp12 import models_ml as ml
	106	+
	107	+ panel = pd.read_parquet(config.path("processed") / "panel.parquet")
	108	+ frames = ml.ml_features(panel)
	109	+ out = ml.rolling_pooled_lgbm(
	110	+ frames, horizons=HORIZONS, window=WINDOW, refit_every=REFIT)
	111	+ f = FC_DIR / "pooled.parquet"
	112	+ out.to_parquet(f, index=False)
	113	+ return f"pooled: {len(out)} forecasts"
	114	+
	115	+
	116	+def main() -> None:
	117	+ """Dispatch model runs across tickers with a process pool."""
	118	+ ap = argparse.ArgumentParser()
	119	+ ap.add_argument("--family", default="all",
	120	+ choices=["econ", "ml", "deep", "pooled", "all"])
	121	+ ap.add_argument("--tickers", default=None, help="comma-separated subset")
	122	+ ap.add_argument("--workers", type=int, default=8)
	123	+ args = ap.parse_args()
	124	+
	125	+ config.setup_logging()
	126	+ FC_DIR.mkdir(parents=True, exist_ok=True)
	127	+ tickers = (
	128	+ args.tickers.split(",") if args.tickers
	129	+ else [i["ticker"] for i in config.universe()]
	130	+ )
	131	+ jobs: list[tuple] = []
	132	+ if args.family in ("econ", "all"):
	133	+ jobs += [(run_econ, tk) for tk in tickers]
	134	+ if args.family in ("ml", "all"):
	135	+ jobs += [(run_ml, tk) for tk in tickers]
	136	+ if args.family in ("deep", "all"):
	137	+ jobs += [(run_deep, tk) for tk in tickers]
	138	+
	139	+ if jobs:
	140	+ with ProcessPoolExecutor(max_workers=args.workers) as pool:
	141	+ futs = {pool.submit(fn, tk): (fn.__name__, tk) for fn, tk in jobs}
	142	+ for fut in as_completed(futs):
	143	+ name, tk = futs[fut]
	144	+ try:
	145	+ logger.info("done %s %s -> %s", name, tk, fut.result())
	146	+ except Exception as exc: # noqa: BLE001
	147	+ logger.error("FAILED %s %s: %s", name, tk, exc)
	148	+ if args.family in ("pooled", "all"):
	149	+ logger.info(run_pooled())
	150	+
	151	+
	152	+if __name__ == "__main__":
	153	+ main()

added src/wp12/__init__.py +15 −0

		`@@ -0,0 +1,15 @@`
	1	+"""
	2	+================================================================
	3	+Auteur : Simon-Pierre Boucher
	4	+Contact : contact@spboucher.ai
	5	+Projet : Prévision de volatilité réalisée multi-actifs
	6	+ (HAR-RV vs GARCH vs Machine Learning)
	7	+Fichier : __init__.py
	8	+Description : Package WP12 — pipeline de prévision de volatilité
	9	+ réalisée multi-actifs (UQO Working Paper No. 12).
	10	+================================================================
	11	+"""
	12	+
	13	+__version__ = "0.1.0"
	14	+__author__ = "Simon-Pierre Boucher"
	15	+__email__ = "contact@spboucher.ai"

added src/wp12/api_client.py +256 −0

		`@@ -0,0 +1,256 @@`
	1	+"""
	2	+================================================================
	3	+Auteur : Simon-Pierre Boucher
	4	+Contact : contact@spboucher.ai
	5	+Projet : Prévision de volatilité réalisée multi-actifs
	6	+ (HAR-RV vs GARCH vs Machine Learning)
	7	+Fichier : api_client.py
	8	+Description : Client réutilisable pour l'API HF Market Data —
	9	+ rate limiting, retry avec backoff exponentiel,
	10	+ pagination et cache local parquet.
	11	+================================================================
	12	+"""
	13	+
	14	+from __future__ import annotations
	15	+
	16	+import io
	17	+import logging
	18	+import time
	19	+from pathlib import Path
	20	+from typing import Any
	21	+
	22	+import pandas as pd
	23	+import requests
	24	+
	25	+from . import config
	26	+
	27	+logger = logging.getLogger(__name__)
	28	+
	29	+_RETRYABLE = {429, 500, 502, 503, 504}
	30	+
	31	+
	32	+class HFMarketDataClient:
	33	+ """Client for https://www.hfmarketdata.io with caching and retries.
	34	+
	35	+ Parameters
	36	+ ----------
	37	+ base_url : str, optional
	38	+ API root; defaults to the value in ``config.yaml``.
	39	+ cache_dir : pathlib.Path, optional
	40	+ Directory for the local parquet cache; defaults to ``data/raw``.
	41	+
	42	+ Notes
	43	+ -----
	44	+ Bars are fetched in ascending order and paginated by advancing the
	45	+ ``start`` parameter one minute past the last received timestamp until
	46	+ an empty page is returned — this is robust to any server-side row cap.
	47	+ """
	48	+
	49	+ def __init__(self, base_url: str \| None = None, cache_dir: Path \| None = None) -> None:
	50	+ cfg = config.load_config()["api"]
	51	+ self.base_url = (base_url or cfg["base_url"]).rstrip("/")
	52	+ self.cache_dir = Path(cache_dir) if cache_dir else config.path("raw")
	53	+ self.min_interval = 1.0 / float(cfg.get("rate_limit_per_sec", 4))
	54	+ self.max_retries = int(cfg.get("max_retries", 5))
	55	+ self.backoff_base = float(cfg.get("backoff_base_sec", 1.0))
	56	+ self.timeout = float(cfg.get("timeout_sec", 60))
	57	+ self.page_limit = int(cfg.get("page_limit", 50_000))
	58	+ self._last_request_ts = 0.0
	59	+ self._session = requests.Session()
	60	+
	61	+ # ------------------------------------------------------------- low level
	62	+ def _throttle(self) -> None:
	63	+ """Sleep as needed to respect the configured request rate."""
	64	+ wait = self._last_request_ts + self.min_interval - time.monotonic()
	65	+ if wait > 0:
	66	+ time.sleep(wait)
	67	+ self._last_request_ts = time.monotonic()
	68	+
	69	+ def _get(self, path: str, params: dict[str, Any] \| None = None) -> requests.Response:
	70	+ """Issue a GET with throttling and exponential-backoff retries.
	71	+
	72	+ Parameters
	73	+ ----------
	74	+ path : str
	75	+ Path relative to the API root (e.g. ``/v1/status``).
	76	+ params : dict, optional
	77	+ Query-string parameters.
	78	+
	79	+ Returns
	80	+ -------
	81	+ requests.Response
	82	+ The successful response (status code < 400).
	83	+
	84	+ Raises
	85	+ ------
	86	+ requests.HTTPError
	87	+ If the request keeps failing after ``max_retries`` attempts.
	88	+ """
	89	+ url = f"{self.base_url}{path}"
	90	+ last_exc: Exception \| None = None
	91	+ for attempt in range(self.max_retries + 1):
	92	+ self._throttle()
	93	+ try:
	94	+ resp = self._session.get(url, params=params, timeout=self.timeout)
	95	+ if resp.status_code in _RETRYABLE:
	96	+ raise requests.HTTPError(f"HTTP {resp.status_code}", response=resp)
	97	+ resp.raise_for_status()
	98	+ return resp
	99	+ except (requests.ConnectionError, requests.Timeout, requests.HTTPError) as exc:
	100	+ resp_obj = getattr(exc, "response", None)
	101	+ if resp_obj is not None and resp_obj.status_code not in _RETRYABLE:
	102	+ raise
	103	+ last_exc = exc
	104	+ delay = self.backoff_base * (2**attempt)
	105	+ logger.warning(
	106	+ "GET %s failed (%s), retry %d/%d in %.1fs",
	107	+ path, exc, attempt + 1, self.max_retries, delay,
	108	+ )
	109	+ time.sleep(delay)
	110	+ raise requests.HTTPError(f"GET {url} failed after {self.max_retries} retries") from last_exc
	111	+
	112	+ def get_json(self, path: str, params: dict[str, Any] \| None = None) -> Any:
	113	+ """GET a JSON endpoint and return the decoded payload."""
	114	+ return self._get(path, params).json()
	115	+
	116	+ # ------------------------------------------------------------- endpoints
	117	+ def status(self) -> dict[str, Any]:
	118	+ """Return the ``/v1/status`` dataset inventory."""
	119	+ return self.get_json("/v1/status")
	120	+
	121	+ def tickers(self, asset: str, search: str \| None = None, limit: int = 100) -> list[str]:
	122	+ """List tickers available for an asset class.
	123	+
	124	+ Parameters
	125	+ ----------
	126	+ asset : str
	127	+ API asset class (``stock``, ``etf``, ``futures``, ``crypto``,
	128	+ ``index``, ``fx``).
	129	+ search : str, optional
	130	+ Substring filter applied server-side.
	131	+ limit : int
	132	+ Maximum number of tickers returned.
	133	+ """
	134	+ params: dict[str, Any] = {"limit": limit}
	135	+ if search:
	136	+ params["search"] = search
	137	+ payload = self.get_json(f"/v1/{asset}/tickers", params)
	138	+ return payload.get("tickers", payload)
	139	+
	140	+ def fetch_bars(
	141	+ self,
	142	+ asset: str,
	143	+ ticker: str,
	144	+ timeframe: str = "1min",
	145	+ adjustment: str \| None = None,
	146	+ start: str \| None = None,
	147	+ end: str \| None = None,
	148	+ ) -> pd.DataFrame:
	149	+ """Download OHLCV bars, transparently handling pagination.
	150	+
	151	+ Parameters
	152	+ ----------
	153	+ asset, ticker : str
	154	+ Instrument identification.
	155	+ timeframe : str
	156	+ ``1min``, ``5min``, ``30min``, ``1hour`` or ``1day``.
	157	+ adjustment : str, optional
	158	+ Price-adjustment scheme; defaults to the class setting in
	159	+ ``config.yaml``.
	160	+ start, end : str, optional
	161	+ Inclusive date bounds (``YYYY-MM-DD`` or full datetime).
	162	+
	163	+ Returns
	164	+ -------
	165	+ pandas.DataFrame
	166	+ Columns ``datetime`` (naive exchange-local timestamps),
	167	+ ``open``, ``high``, ``low``, ``close``, ``volume``; sorted,
	168	+ duplicate timestamps dropped.
	169	+ """
	170	+ adjustment = adjustment or config.adjustment_for(asset)
	171	+ frames: list[pd.DataFrame] = []
	172	+ cursor = start
	173	+ n_pages = 0
	174	+ while True:
	175	+ params: dict[str, Any] = {
	176	+ "timeframe": timeframe,
	177	+ "adjustment": adjustment,
	178	+ "order": "asc",
	179	+ "limit": self.page_limit,
	180	+ "format": "csv",
	181	+ }
	182	+ if cursor:
	183	+ params["start"] = cursor
	184	+ if end:
	185	+ params["end"] = end
	186	+ resp = self._get(f"/v1/bars/{asset}/{ticker}", params)
	187	+ page = pd.read_csv(io.StringIO(resp.text)) if resp.text.strip() else pd.DataFrame()
	188	+ if page.empty:
	189	+ break
	190	+ page["datetime"] = pd.to_datetime(page["datetime"])
	191	+ frames.append(page)
	192	+ n_pages += 1
	193	+ last_dt = page["datetime"].iloc[-1]
	194	+ logger.debug("%s/%s %s page %d: %d rows, up to %s",
	195	+ asset, ticker, timeframe, n_pages, len(page), last_dt)
	196	+ if len(page) < self.page_limit:
	197	+ break
	198	+ cursor = (last_dt + pd.Timedelta(minutes=1)).strftime("%Y-%m-%d %H:%M:%S")
	199	+ if not frames:
	200	+ logger.warning("%s/%s %s: no data returned", asset, ticker, timeframe)
	201	+ return pd.DataFrame(columns=["datetime", "open", "high", "low", "close", "volume"])
	202	+ df = pd.concat(frames, ignore_index=True)
	203	+ df = (
	204	+ df.drop(columns=["ticker"], errors="ignore")
	205	+ .drop_duplicates(subset="datetime")
	206	+ .sort_values("datetime")
	207	+ .reset_index(drop=True)
	208	+ )
	209	+ if "volume" not in df.columns: # some index series carry no volume
	210	+ df["volume"] = float("nan")
	211	+ logger.info("%s/%s %s: %d bars (%s -> %s)", asset, ticker, timeframe,
	212	+ len(df), df["datetime"].iloc[0], df["datetime"].iloc[-1])
	213	+ return df
	214	+
	215	+ # ------------------------------------------------------------- cache
	216	+ def _cache_file(self, asset: str, ticker: str, timeframe: str) -> Path:
	217	+ safe = ticker.replace("/", "-")
	218	+ return self.cache_dir / f"{asset}_{safe}_{timeframe}.parquet"
	219	+
	220	+ def get_bars(
	221	+ self,
	222	+ asset: str,
	223	+ ticker: str,
	224	+ timeframe: str = "1min",
	225	+ start: str \| None = None,
	226	+ end: str \| None = None,
	227	+ refresh: bool = False,
	228	+ ) -> pd.DataFrame:
	229	+ """Return bars from the local parquet cache, downloading if absent.
	230	+
	231	+ Parameters
	232	+ ----------
	233	+ asset, ticker, timeframe : str
	234	+ Instrument identification and bar frequency.
	235	+ start, end : str, optional
	236	+ Date bounds used when the series must be downloaded; when the
	237	+ cache is hit the full cached range is returned (filter at the
	238	+ call site if needed).
	239	+ refresh : bool
	240	+ Force a re-download even when a cache file exists.
	241	+
	242	+ Returns
	243	+ -------
	244	+ pandas.DataFrame
	245	+ Same layout as :meth:`fetch_bars`.
	246	+ """
	247	+ f = self._cache_file(asset, ticker, timeframe)
	248	+ if f.exists() and not refresh:
	249	+ logger.debug("cache hit: %s", f.name)
	250	+ return pd.read_parquet(f)
	251	+ df = self.fetch_bars(asset, ticker, timeframe, start=start, end=end)
	252	+ if not df.empty:
	253	+ f.parent.mkdir(parents=True, exist_ok=True)
	254	+ df.to_parquet(f, index=False)
	255	+ logger.info("cached %s (%.1f MB)", f.name, f.stat().st_size / 1e6)
	256	+ return df

added src/wp12/backtest.py +210 −0

		`@@ -0,0 +1,210 @@`
	1	+"""
	2	+================================================================
	3	+Auteur : Simon-Pierre Boucher
	4	+Contact : contact@spboucher.ai
	5	+Projet : Prévision de volatilité réalisée multi-actifs
	6	+ (HAR-RV vs GARCH vs Machine Learning)
	7	+Fichier : backtest.py
	8	+Description : Valeur économique des prévisions — stratégie de
	9	+ volatility timing (cible de vol, cap de levier,
	10	+ coûts de transaction), gains d'utilité mean-
	11	+ variance, et backtesting de la VaR (tests de
	12	+ Kupiec et Christoffersen).
	13	+================================================================
	14	+"""
	15	+
	16	+from __future__ import annotations
	17	+
	18	+import logging
	19	+
	20	+import numpy as np
	21	+import pandas as pd
	22	+from scipy import stats
	23	+
	24	+logger = logging.getLogger(__name__)
	25	+
	26	+ANN = 252
	27	+
	28	+
	29	+def volatility_timing(
	30	+ ret: pd.Series,
	31	+ pred_var: pd.Series,
	32	+ vol_target_ann: float = 0.10,
	33	+ leverage_cap: float = 3.0,
	34	+ tc_bps: float = 5.0,
	35	+) -> pd.DataFrame:
	36	+ """Volatility-timing strategy driven by one-day-ahead variance forecasts.
	37	+
	38	+ The weight on the risky asset is
	39	+ ``w_t = min(cap, sigma_target / sigma_hat_{t+1})`` where
	40	+ ``sigma_hat`` is the annualized forecast volatility; the position is
	41	+ financed at zero rate and transaction costs are charged on turnover.
	42	+
	43	+ Parameters
	44	+ ----------
	45	+ ret : pandas.Series
	46	+ Daily close-to-close log returns of the asset.
	47	+ pred_var : pandas.Series
	48	+ Daily variance forecasts indexed by forecast origin (the
	49	+ position for day t+1 uses the forecast dated t).
	50	+ vol_target_ann : float
	51	+ Annualized volatility target.
	52	+ leverage_cap : float
	53	+ Maximum absolute weight.
	54	+ tc_bps : float
	55	+ One-way transaction cost in basis points of traded notional.
	56	+
	57	+ Returns
	58	+ -------
	59	+ pandas.DataFrame
	60	+ Daily strategy ledger: ``weight``, ``gross``, ``net`` returns.
	61	+ """
	62	+ sigma_ann = np.sqrt(pred_var.clip(lower=1e-12) * ANN)
	63	+ w = (vol_target_ann / sigma_ann).clip(upper=leverage_cap)
	64	+ # weight decided at origin t applies to return of t+1
	65	+ w_applied = w.shift(1)
	66	+ aligned = pd.concat([ret.rename("ret"), w_applied.rename("w")], axis=1).dropna()
	67	+ gross = aligned["w"] * aligned["ret"]
	68	+ turnover = aligned["w"].diff().abs().fillna(aligned["w"].abs())
	69	+ net = gross - turnover * tc_bps / 1e4
	70	+ return pd.DataFrame({"weight": aligned["w"], "gross": gross, "net": net})
	71	+
	72	+
	73	+def perf_stats(returns: pd.Series) -> dict[str, float]:
	74	+ """Annualized performance statistics of a daily return series."""
	75	+ r = returns.dropna()
	76	+ if len(r) < 60:
	77	+ return {k: np.nan for k in ("ann_ret", "ann_vol", "sharpe", "max_dd")}
	78	+ mu = r.mean() * ANN
	79	+ sig = r.std() * np.sqrt(ANN)
	80	+ cum = r.cumsum()
	81	+ dd = float((cum - cum.cummax()).min())
	82	+ return {
	83	+ "ann_ret": float(mu),
	84	+ "ann_vol": float(sig),
	85	+ "sharpe": float(mu / sig) if sig > 0 else np.nan,
	86	+ "max_dd": dd,
	87	+ }
	88	+
	89	+
	90	+def utility_gain(
	91	+ net_a: pd.Series, net_b: pd.Series, risk_aversion: float = 5.0
	92	+) -> float:
	93	+ """Annualized management fee an investor would pay to switch B -> A.
	94	+
	95	+ Mean-variance utility (Fleming, Kirby & Ostdiek 2001):
	96	+ ``Delta`` solves ``U(r_A - Delta) = U(r_B)`` with
	97	+ ``U(r) = E[r] - (gamma/2) Var[r]``.
	98	+
	99	+ Parameters
	100	+ ----------
	101	+ net_a, net_b : pandas.Series
	102	+ Daily net strategy returns (A = candidate, B = benchmark).
	103	+ risk_aversion : float
	104	+ Relative risk-aversion coefficient gamma.
	105	+
	106	+ Returns
	107	+ -------
	108	+ float
	109	+ Annualized utility gain (fee), in return units.
	110	+ """
	111	+ a, b = net_a.dropna(), net_b.dropna()
	112	+ common = a.index.intersection(b.index)
	113	+ a, b = a[common], b[common]
	114	+ ua = a.mean() - 0.5 * risk_aversion * a.var()
	115	+ ub = b.mean() - 0.5 * risk_aversion * b.var()
	116	+ return float((ua - ub) * ANN)
	117	+
	118	+
	119	+# ---------------------------------------------------------------- VaR backtest
	120	+def var_forecast(
	121	+ pred_var: pd.Series, level: float = 0.01
	122	+) -> pd.Series:
	123	+ """Gaussian Value-at-Risk from a variance forecast (positive number)."""
	124	+ z = stats.norm.ppf(level)
	125	+ return -z * np.sqrt(pred_var.clip(lower=1e-12))
	126	+
	127	+
	128	+def kupiec_test(violations: np.ndarray, level: float) -> tuple[float, float]:
	129	+ """Kupiec (1995) proportion-of-failures LR test.
	130	+
	131	+ Returns
	132	+ -------
	133	+ tuple
	134	+ ``(LR statistic, p-value)`` under chi2(1).
	135	+ """
	136	+ n = len(violations)
	137	+ x = int(violations.sum())
	138	+ if n == 0:
	139	+ return np.nan, np.nan
	140	+ pi_hat = x / n
	141	+ if x in (0, n):
	142	+ lr = -2 * n * (np.log(1 - level) if x == 0 else np.log(level))
	143	+ else:
	144	+ ll0 = (n - x) * np.log(1 - level) + x * np.log(level)
	145	+ ll1 = (n - x) * np.log(1 - pi_hat) + x * np.log(pi_hat)
	146	+ lr = -2 * (ll0 - ll1)
	147	+ return float(lr), float(stats.chi2.sf(lr, df=1))
	148	+
	149	+
	150	+def christoffersen_test(violations: np.ndarray, level: float) -> tuple[float, float]:
	151	+ """Christoffersen (1998) conditional-coverage LR test (chi2(2)).
	152	+
	153	+ Combines unconditional coverage with first-order independence of the
	154	+ violation indicator sequence.
	155	+ """
	156	+ v = violations.astype(int)
	157	+ n = len(v)
	158	+ if n < 2:
	159	+ return np.nan, np.nan
	160	+ pairs = np.column_stack([v[:-1], v[1:]])
	161	+ n00 = int(((pairs[:, 0] == 0) & (pairs[:, 1] == 0)).sum())
	162	+ n01 = int(((pairs[:, 0] == 0) & (pairs[:, 1] == 1)).sum())
	163	+ n10 = int(((pairs[:, 0] == 1) & (pairs[:, 1] == 0)).sum())
	164	+ n11 = int(((pairs[:, 0] == 1) & (pairs[:, 1] == 1)).sum())
	165	+ pi01 = n01 / max(n00 + n01, 1)
	166	+ pi11 = n11 / max(n10 + n11, 1)
	167	+ pi = (n01 + n11) / max(n00 + n01 + n10 + n11, 1)
	168	+
	169	+ def _safe_log(x: float) -> float:
	170	+ return np.log(x) if x > 0 else 0.0
	171	+
	172	+ ll_ind = (n00 * _safe_log(1 - pi01) + n01 * _safe_log(pi01)
	173	+ + n10 * _safe_log(1 - pi11) + n11 * _safe_log(pi11))
	174	+ ll_null = (n00 + n10) * _safe_log(1 - pi) + (n01 + n11) * _safe_log(pi)
	175	+ lr_ind = -2 * (ll_null - ll_ind)
	176	+ lr_uc, _ = kupiec_test(violations, level)
	177	+ lr_cc = lr_uc + lr_ind
	178	+ return float(lr_cc), float(stats.chi2.sf(lr_cc, df=2))
	179	+
	180	+
	181	+def var_backtest(
	182	+ ret: pd.Series, pred_var: pd.Series, level: float
	183	+) -> dict[str, float]:
	184	+ """Full VaR backtest for one model / asset / level.
	185	+
	186	+ Parameters
	187	+ ----------
	188	+ ret : pandas.Series
	189	+ Daily returns.
	190	+ pred_var : pandas.Series
	191	+ One-day variance forecasts indexed by origin (applied to t+1).
	192	+ level : float
	193	+ VaR tail level (0.01 or 0.05).
	194	+
	195	+ Returns
	196	+ -------
	197	+ dict
	198	+ Violation rate and Kupiec / Christoffersen p-values.
	199	+ """
	200	+ var = var_forecast(pred_var, level).shift(1)
	201	+ aligned = pd.concat([ret.rename("ret"), var.rename("var")], axis=1).dropna()
	202	+ viol = (aligned["ret"] < -aligned["var"]).to_numpy()
	203	+ lr_uc, p_uc = kupiec_test(viol, level)
	204	+ lr_cc, p_cc = christoffersen_test(viol, level)
	205	+ return {
	206	+ "n": float(len(viol)),
	207	+ "viol_rate": float(viol.mean()) if len(viol) else np.nan,
	208	+ "kupiec_p": p_uc,
	209	+ "christoffersen_p": p_cc,
	210	+ }

added src/wp12/config.py +115 −0

		`@@ -0,0 +1,115 @@`
	1	+"""
	2	+================================================================
	3	+Auteur : Simon-Pierre Boucher
	4	+Contact : contact@spboucher.ai
	5	+Projet : Prévision de volatilité réalisée multi-actifs
	6	+ (HAR-RV vs GARCH vs Machine Learning)
	7	+Fichier : config.py
	8	+Description : Chargement de config.yaml, chemins du projet et
	9	+ helpers d'accès à l'univers d'actifs.
	10	+================================================================
	11	+"""
	12	+
	13	+from __future__ import annotations
	14	+
	15	+import logging
	16	+import os
	17	+from functools import lru_cache
	18	+from pathlib import Path
	19	+from typing import Any
	20	+
	21	+import yaml
	22	+
	23	+ROOT = Path(os.environ.get("WP12_ROOT", Path(__file__).resolve().parents[2]))
	24	+CONFIG_PATH = ROOT / "config.yaml"
	25	+
	26	+
	27	+def setup_logging(level: int = logging.INFO) -> None:
	28	+ """Configure root logging once, with a compact console format.
	29	+
	30	+ Parameters
	31	+ ----------
	32	+ level : int
	33	+ Logging level applied to the root logger.
	34	+ """
	35	+ logging.basicConfig(
	36	+ level=level,
	37	+ format="%(asctime)s %(levelname)-7s %(name)s: %(message)s",
	38	+ datefmt="%H:%M:%S",
	39	+ force=False,
	40	+ )
	41	+
	42	+
	43	+@lru_cache(maxsize=1)
	44	+def load_config() -> dict[str, Any]:
	45	+ """Load and cache the central YAML configuration.
	46	+
	47	+ Returns
	48	+ -------
	49	+ dict
	50	+ Parsed contents of ``config.yaml`` at the repository root.
	51	+ """
	52	+ with open(CONFIG_PATH, encoding="utf-8") as fh:
	53	+ return yaml.safe_load(fh)
	54	+
	55	+
	56	+def path(key: str) -> Path:
	57	+ """Resolve a path from the ``paths`` section of the config.
	58	+
	59	+ Parameters
	60	+ ----------
	61	+ key : str
	62	+ Key inside the ``paths`` mapping (e.g. ``"raw"``, ``"figures"``).
	63	+
	64	+ Returns
	65	+ -------
	66	+ pathlib.Path
	67	+ Absolute path anchored at the repository root.
	68	+ """
	69	+ p = ROOT / load_config()["paths"][key]
	70	+ p.mkdir(parents=True, exist_ok=True)
	71	+ return p
	72	+
	73	+
	74	+def universe(include_predictors: bool = False) -> list[dict[str, str]]:
	75	+ """Flatten the asset universe into a list of instrument records.
	76	+
	77	+ Parameters
	78	+ ----------
	79	+ include_predictors : bool
	80	+ If True, append predictor-only instruments (e.g. VIX).
	81	+
	82	+ Returns
	83	+ -------
	84	+ list of dict
	85	+ Records with keys ``ticker``, ``asset`` and ``cls`` (asset class
	86	+ label used throughout the paper: equity, fx, crypto, futures,
	87	+ index).
	88	+ """
	89	+ cfg = load_config()
	90	+ out: list[dict[str, str]] = []
	91	+ for cls, items in cfg["universe"].items():
	92	+ for it in items:
	93	+ out.append({"ticker": it["ticker"], "asset": it["asset"], "cls": cls})
	94	+ if include_predictors:
	95	+ for cls, items in cfg.get("predictors", {}).items():
	96	+ for it in items:
	97	+ out.append({"ticker": it["ticker"], "asset": it["asset"], "cls": cls})
	98	+ return out
	99	+
	100	+
	101	+def adjustment_for(asset: str) -> str:
	102	+ """Return the price-adjustment scheme configured for an asset class.
	103	+
	104	+ Parameters
	105	+ ----------
	106	+ asset : str
	107	+ API asset class (``stock``, ``etf``, ``futures``, ``crypto``,
	108	+ ``index``, ``fx``).
	109	+
	110	+ Returns
	111	+ -------
	112	+ str
	113	+ Adjustment keyword understood by the HF Market Data API.
	114	+ """
	115	+ return load_config()["adjustment"][asset]

added src/wp12/data_pipeline.py +136 −0

		`@@ -0,0 +1,136 @@`
	1	+"""
	2	+================================================================
	3	+Auteur : Simon-Pierre Boucher
	4	+Contact : contact@spboucher.ai
	5	+Projet : Prévision de volatilité réalisée multi-actifs
	6	+ (HAR-RV vs GARCH vs Machine Learning)
	7	+Fichier : data_pipeline.py
	8	+Description : Orchestration données — téléchargement (via
	9	+ api_client), construction des mesures RV
	10	+ journalières par actif et assemblage du panel
	11	+ multi-actifs (data/processed).
	12	+================================================================
	13	+"""
	14	+
	15	+from __future__ import annotations
	16	+
	17	+import logging
	18	+from pathlib import Path
	19	+
	20	+import pandas as pd
	21	+
	22	+from . import config, realized_vol
	23	+from .api_client import HFMarketDataClient
	24	+
	25	+logger = logging.getLogger(__name__)
	26	+
	27	+
	28	+def download_instrument(
	29	+ client: HFMarketDataClient,
	30	+ inst: dict[str, str],
	31	+ timeframes: tuple[str, ...] = ("1min", "1day"),
	32	+ refresh: bool = False,
	33	+) -> dict[str, pd.DataFrame]:
	34	+ """Download (or load from cache) all timeframes for one instrument.
	35	+
	36	+ Parameters
	37	+ ----------
	38	+ client : HFMarketDataClient
	39	+ API client with parquet cache.
	40	+ inst : dict
	41	+ Instrument record with keys ``ticker``, ``asset``, ``cls``.
	42	+ timeframes : tuple of str
	43	+ Bar frequencies to materialize.
	44	+ refresh : bool
	45	+ Force re-download.
	46	+
	47	+ Returns
	48	+ -------
	49	+ dict
	50	+ Mapping timeframe -> bars DataFrame.
	51	+ """
	52	+ cfg = config.load_config()["sample"]
	53	+ out = {}
	54	+ for tf in timeframes:
	55	+ out[tf] = client.get_bars(
	56	+ inst["asset"], inst["ticker"], tf,
	57	+ start=cfg["start"], end=cfg["end"], refresh=refresh,
	58	+ )
	59	+ return out
	60	+
	61	+
	62	+def rv_file(ticker: str) -> Path:
	63	+ """Path of the processed daily-RV parquet for one ticker."""
	64	+ return config.path("processed") / f"rv_{ticker}.parquet"
	65	+
	66	+
	67	+def build_instrument_rv(
	68	+ client: HFMarketDataClient,
	69	+ inst: dict[str, str],
	70	+ refresh: bool = False,
	71	+) -> pd.DataFrame:
	72	+ """Build (and cache) the daily realized-measure panel for one ticker.
	73	+
	74	+ Parameters
	75	+ ----------
	76	+ client : HFMarketDataClient
	77	+ API client.
	78	+ inst : dict
	79	+ Instrument record (``ticker``, ``asset``, ``cls``).
	80	+ refresh : bool
	81	+ Recompute even if the processed file exists.
	82	+
	83	+ Returns
	84	+ -------
	85	+ pandas.DataFrame
	86	+ Daily measures indexed by date, with ``ticker`` and ``cls``
	87	+ columns added.
	88	+ """
	89	+ f = rv_file(inst["ticker"])
	90	+ if f.exists() and not refresh:
	91	+ return pd.read_parquet(f)
	92	+ cfg = config.load_config()
	93	+ bars = client.get_bars(
	94	+ inst["asset"], inst["ticker"], "1min",
	95	+ start=cfg["sample"]["start"], end=cfg["sample"]["end"],
	96	+ )
	97	+ daily = realized_vol.build_daily_rv(
	98	+ bars, cls=inst["cls"],
	99	+ jump_alpha=cfg["realized_vol"]["jump_test_alpha"],
	100	+ )
	101	+ if daily.empty:
	102	+ logger.warning("empty RV panel for %s", inst["ticker"])
	103	+ return daily
	104	+ daily["ticker"] = inst["ticker"]
	105	+ daily["cls"] = inst["cls"]
	106	+ daily.to_parquet(f)
	107	+ logger.info("%s: %d RV days -> %s", inst["ticker"], len(daily), f.name)
	108	+ return daily
	109	+
	110	+
	111	+def build_panel(refresh: bool = False) -> pd.DataFrame:
	112	+ """Assemble the multi-asset daily RV panel for the whole universe.
	113	+
	114	+ Parameters
	115	+ ----------
	116	+ refresh : bool
	117	+ Recompute per-ticker RV files even when cached.
	118	+
	119	+ Returns
	120	+ -------
	121	+ pandas.DataFrame
	122	+ Long panel (date x ticker) of daily realized measures, also
	123	+ written to ``data/processed/panel.parquet``.
	124	+ """
	125	+ client = HFMarketDataClient()
	126	+ frames = []
	127	+ for inst in config.universe(include_predictors=True):
	128	+ df = build_instrument_rv(client, inst, refresh=refresh)
	129	+ if not df.empty:
	130	+ frames.append(df)
	131	+ panel = pd.concat(frames).sort_index()
	132	+ out = config.path("processed") / "panel.parquet"
	133	+ panel.to_parquet(out)
	134	+ logger.info("panel: %d rows, %d tickers -> %s",
	135	+ len(panel), panel["ticker"].nunique(), out.name)
	136	+ return panel

added src/wp12/evaluation.py +292 −0

		`@@ -0,0 +1,292 @@`
	1	+"""
	2	+================================================================
	3	+Auteur : Simon-Pierre Boucher
	4	+Contact : contact@spboucher.ai
	5	+Projet : Prévision de volatilité réalisée multi-actifs
	6	+ (HAR-RV vs GARCH vs Machine Learning)
	7	+Fichier : evaluation.py
	8	+Description : Évaluation des prévisions — pertes QLIKE/MSE,
	9	+ test de Diebold-Mariano (HAC), Model Confidence
	10	+ Set (bootstrap stationnaire), régressions de
	11	+ Mincer-Zarnowitz, tests d'encompassing et
	12	+ combinaisons de prévisions.
	13	+================================================================
	14	+"""
	15	+
	16	+from __future__ import annotations
	17	+
	18	+import logging
	19	+
	20	+import numpy as np
	21	+import pandas as pd
	22	+from scipy import stats
	23	+
	24	+logger = logging.getLogger(__name__)
	25	+
	26	+EPS = 1e-12
	27	+
	28	+
	29	+# ---------------------------------------------------------------- losses
	30	+def qlike(actual: np.ndarray, forecast: np.ndarray) -> np.ndarray:
	31	+ """QLIKE loss (per observation), robust ranking under noisy proxies.
	32	+
	33	+ Parameters
	34	+ ----------
	35	+ actual, forecast : numpy.ndarray
	36	+ Realized and predicted variance (both strictly positive).
	37	+
	38	+ Returns
	39	+ -------
	40	+ numpy.ndarray
	41	+ ``a/f - log(a/f) - 1`` — non-negative, zero iff perfect.
	42	+ """
	43	+ a = np.maximum(actual, EPS)
	44	+ f = np.maximum(forecast, EPS)
	45	+ ratio = a / f
	46	+ return ratio - np.log(ratio) - 1.0
	47	+
	48	+
	49	+def mse(actual: np.ndarray, forecast: np.ndarray) -> np.ndarray:
	50	+ """Squared-error loss (per observation) on variance levels."""
	51	+ return (actual - forecast) ** 2
	52	+
	53	+
	54	+LOSSES = {"qlike": qlike, "mse": mse}
	55	+
	56	+
	57	+# ---------------------------------------------------------------- DM test
	58	+def newey_west_var(x: np.ndarray, lags: int) -> float:
	59	+ """Newey-West long-run variance of the mean of a series."""
	60	+ x = x - x.mean()
	61	+ n = len(x)
	62	+ v = float(x @ x) / n
	63	+ for k in range(1, min(lags, n - 1) + 1):
	64	+ w = 1.0 - k / (lags + 1.0)
	65	+ v += 2.0 * w * float(x[k:] @ x[:-k]) / n
	66	+ return v
	67	+
	68	+
	69	+def diebold_mariano(
	70	+ loss_a: np.ndarray, loss_b: np.ndarray, h: int = 1
	71	+) -> tuple[float, float]:
	72	+ """Diebold-Mariano test of equal predictive accuracy.
	73	+
	74	+ Parameters
	75	+ ----------
	76	+ loss_a, loss_b : numpy.ndarray
	77	+ Per-observation losses of models A and B (aligned).
	78	+ h : int
	79	+ Forecast horizon; the HAC bandwidth is ``h - 1``.
	80	+
	81	+ Returns
	82	+ -------
	83	+ tuple
	84	+ ``(statistic, p_value)``; negative statistic favours model A.
	85	+ """
	86	+ d = loss_a - loss_b
	87	+ n = len(d)
	88	+ if n < 30:
	89	+ return np.nan, np.nan
	90	+ v = newey_west_var(d, max(h - 1, 0))
	91	+ if v <= 0:
	92	+ return np.nan, np.nan
	93	+ stat = float(d.mean() / np.sqrt(v / n))
	94	+ return stat, float(2 * stats.norm.sf(abs(stat)))
	95	+
	96	+
	97	+# ---------------------------------------------------------------- MCS
	98	+def _stationary_bootstrap_idx(
	99	+ n: int, b: int, avg_block: float, rng: np.random.Generator
	100	+) -> np.ndarray:
	101	+ """Index matrix (b x n) for the Politis-Romano stationary bootstrap."""
	102	+ p = 1.0 / avg_block
	103	+ starts = rng.integers(0, n, size=(b, n))
	104	+ jumps = rng.random((b, n)) < p
	105	+ jumps[:, 0] = True
	106	+ idx = np.zeros((b, n), dtype=np.int64)
	107	+ for j in range(n):
	108	+ if j == 0:
	109	+ idx[:, 0] = starts[:, 0]
	110	+ else:
	111	+ cont = (idx[:, j - 1] + 1) % n
	112	+ idx[:, j] = np.where(jumps[:, j], starts[:, j], cont)
	113	+ return idx
	114	+
	115	+
	116	+def model_confidence_set(
	117	+ losses: pd.DataFrame,
	118	+ alpha: float = 0.10,
	119	+ n_boot: int = 5000,
	120	+ avg_block: float = 10.0,
	121	+ seed: int = 20260810,
	122	+) -> pd.DataFrame:
	123	+ """Model Confidence Set of Hansen, Lunde & Nason (2011), T_max statistic.
	124	+
	125	+ Iteratively eliminates the worst model until the null of equal
	126	+ predictive ability survives at level ``alpha``.
	127	+
	128	+ Parameters
	129	+ ----------
	130	+ losses : pandas.DataFrame
	131	+ Per-observation losses, one column per model (aligned rows).
	132	+ alpha : float
	133	+ MCS level (models kept form the ``1 - alpha`` confidence set).
	134	+ n_boot : int
	135	+ Bootstrap replications.
	136	+ avg_block : float
	137	+ Mean block length of the stationary bootstrap.
	138	+ seed : int
	139	+ RNG seed.
	140	+
	141	+ Returns
	142	+ -------
	143	+ pandas.DataFrame
	144	+ One row per model with its elimination-order ``mcs_pvalue`` and
	145	+ an ``in_mcs`` flag.
	146	+ """
	147	+ losses = losses.dropna()
	148	+ models = list(losses.columns)
	149	+ L = losses.to_numpy()
	150	+ n = len(L)
	151	+ rng = np.random.default_rng(seed)
	152	+ idx = _stationary_bootstrap_idx(n, n_boot, avg_block, rng)
	153	+
	154	+ included = list(range(len(models)))
	155	+ pvals: dict[str, float] = {}
	156	+ running_max_p = 0.0
	157	+ while len(included) > 1:
	158	+ Ls = L[:, included]
	159	+ dbar_i = Ls.mean(axis=0) # mean loss per model
	160	+ d_i_dot = dbar_i - dbar_i.mean() # vs. set average
	161	+ # bootstrap distribution of the max t-stat
	162	+ boot_means = Ls[idx].mean(axis=1) # (b, k)
	163	+ boot_center = boot_means - dbar_i
	164	+ var_i = (boot_center - boot_center.mean(axis=0)).var(axis=0) + EPS
	165	+ t_i = d_i_dot / np.sqrt(var_i)
	166	+ t_max = float(t_i.max())
	167	+ boot_d = boot_center - boot_center.mean(axis=1, keepdims=True)
	168	+ boot_t = boot_d / np.sqrt(var_i)
	169	+ boot_max = boot_t.max(axis=1)
	170	+ p = float(np.mean(boot_max >= t_max))
	171	+ running_max_p = max(running_max_p, p)
	172	+ worst = included[int(np.argmax(t_i))]
	173	+ pvals[models[worst]] = running_max_p
	174	+ if p >= alpha:
	175	+ for i in included:
	176	+ pvals.setdefault(models[i], max(running_max_p, p))
	177	+ break
	178	+ included.remove(worst)
	179	+ if len(included) == 1:
	180	+ pvals[models[included[0]]] = 1.0
	181	+ out = pd.DataFrame(
	182	+ {"model": models, "mcs_pvalue": [pvals[m] for m in models]}
	183	+ )
	184	+ out["in_mcs"] = out["mcs_pvalue"] >= alpha
	185	+ return out.sort_values("mcs_pvalue", ascending=False).reset_index(drop=True)
	186	+
	187	+
	188	+# ---------------------------------------------------------------- MZ / encompassing
	189	+def mincer_zarnowitz(
	190	+ actual: np.ndarray, forecast: np.ndarray, h: int = 1
	191	+) -> dict[str, float]:
	192	+ """Mincer-Zarnowitz levels regression with HAC inference.
	193	+
	194	+ Regresses realized variance on a constant and the forecast and tests
	195	+ the joint null ``(a, b) = (0, 1)``.
	196	+
	197	+ Returns
	198	+ -------
	199	+ dict
	200	+ ``alpha``, ``beta``, ``r2`` and the Wald ``p_joint``.
	201	+ """
	202	+ mask = np.isfinite(actual) & np.isfinite(forecast)
	203	+ y, f = actual[mask], forecast[mask]
	204	+ n = len(y)
	205	+ if n < 30:
	206	+ return {"alpha": np.nan, "beta": np.nan, "r2": np.nan, "p_joint": np.nan}
	207	+ X = np.column_stack([np.ones(n), f])
	208	+ beta, *_ = np.linalg.lstsq(X, y, rcond=None)
	209	+ e = y - X @ beta
	210	+ r2 = 1.0 - e.var() / y.var()
	211	+ # HAC covariance (Newey-West, bandwidth h-1)
	212	+ lags = max(h - 1, 0)
	213	+ Xe = X * e[:, None]
	214	+ S = Xe.T @ Xe / n
	215	+ for k in range(1, min(lags, n - 1) + 1):
	216	+ w = 1.0 - k / (lags + 1.0)
	217	+ G = Xe[k:].T @ Xe[:-k] / n
	218	+ S += w * (G + G.T)
	219	+ XX_inv = np.linalg.inv(X.T @ X / n)
	220	+ V = XX_inv @ S @ XX_inv / n
	221	+ r = beta - np.array([0.0, 1.0])
	222	+ try:
	223	+ wald = float(r @ np.linalg.solve(V, r))
	224	+ p_joint = float(stats.chi2.sf(wald, df=2))
	225	+ except np.linalg.LinAlgError:
	226	+ p_joint = np.nan
	227	+ return {"alpha": float(beta[0]), "beta": float(beta[1]),
	228	+ "r2": float(r2), "p_joint": p_joint}
	229	+
	230	+
	231	+def forecast_encompassing(
	232	+ actual: np.ndarray, f_base: np.ndarray, f_rival: np.ndarray, h: int = 1
	233	+) -> dict[str, float]:
	234	+ """Test whether ``f_base`` encompasses ``f_rival``.
	235	+
	236	+ Runs ``y = a + b1 f_base + b2 f_rival + e`` and reports the HAC
	237	+ t-test of ``b2 = 0``; rejection means the rival adds information.
	238	+
	239	+ Returns
	240	+ -------
	241	+ dict
	242	+ ``b1``, ``b2``, ``t_b2``, ``p_b2``.
	243	+ """
	244	+ mask = np.isfinite(actual) & np.isfinite(f_base) & np.isfinite(f_rival)
	245	+ y, f1, f2 = actual[mask], f_base[mask], f_rival[mask]
	246	+ n = len(y)
	247	+ if n < 30:
	248	+ return {"b1": np.nan, "b2": np.nan, "t_b2": np.nan, "p_b2": np.nan}
	249	+ X = np.column_stack([np.ones(n), f1, f2])
	250	+ beta, *_ = np.linalg.lstsq(X, y, rcond=None)
	251	+ e = y - X @ beta
	252	+ lags = max(h - 1, 0)
	253	+ Xe = X * e[:, None]
	254	+ S = Xe.T @ Xe / n
	255	+ for k in range(1, min(lags, n - 1) + 1):
	256	+ w = 1.0 - k / (lags + 1.0)
	257	+ G = Xe[k:].T @ Xe[:-k] / n
	258	+ S += w * (G + G.T)
	259	+ XX_inv = np.linalg.inv(X.T @ X / n)
	260	+ V = XX_inv @ S @ XX_inv / n
	261	+ t2 = float(beta[2] / np.sqrt(V[2, 2]))
	262	+ return {"b1": float(beta[1]), "b2": float(beta[2]), "t_b2": t2,
	263	+ "p_b2": float(2 * stats.norm.sf(abs(t2)))}
	264	+
	265	+
	266	+# ---------------------------------------------------------------- combinations
	267	+def combine_forecasts(
	268	+ forecasts: pd.DataFrame, actual: pd.Series, window: int = 60
	269	+) -> pd.DataFrame:
	270	+ """Simple and inverse-MSE forecast combinations.
	271	+
	272	+ Parameters
	273	+ ----------
	274	+ forecasts : pandas.DataFrame
	275	+ One column per model, indexed by forecast origin.
	276	+ actual : pandas.Series
	277	+ Realized target aligned on the same index.
	278	+ window : int
	279	+ Trailing window used for the inverse-MSE weights.
	280	+
	281	+ Returns
	282	+ -------
	283	+ pandas.DataFrame
	284	+ Columns ``Comb-Mean`` and ``Comb-InvMSE``.
	285	+ """
	286	+ mean_comb = forecasts.mean(axis=1)
	287	+ sq_err = forecasts.sub(actual, axis=0) ** 2
	288	+ inv_mse = 1.0 / sq_err.rolling(window, min_periods=window // 2).mean().shift(1)
	289	+ weights = inv_mse.div(inv_mse.sum(axis=1), axis=0)
	290	+ inv_comb = (forecasts * weights).sum(axis=1)
	291	+ inv_comb[weights.isna().all(axis=1)] = np.nan
	292	+ return pd.DataFrame({"Comb-Mean": mean_comb, "Comb-InvMSE": inv_comb})

added src/wp12/models_econ.py +379 −0

		`@@ -0,0 +1,379 @@`
	1	+"""
	2	+================================================================
	3	+Auteur : Simon-Pierre Boucher
	4	+Contact : contact@spboucher.ai
	5	+Projet : Prévision de volatilité réalisée multi-actifs
	6	+ (HAR-RV vs GARCH vs Machine Learning)
	7	+Fichier : models_econ.py
	8	+Description : Modèles économétriques — famille HAR (HAR-RV,
	9	+ HAR-J, HAR-CJ, SHAR, HARQ, log-HAR) estimée par
	10	+ projections directes OLS roulantes, et famille
	11	+ GARCH (GARCH, EGARCH, GJR, Realized GARCH) en
	12	+ réestimation roulante.
	13	+================================================================
	14	+"""
	15	+
	16	+from __future__ import annotations
	17	+
	18	+import logging
	19	+import warnings
	20	+
	21	+import numpy as np
	22	+import pandas as pd
	23	+from arch import arch_model
	24	+from scipy import optimize
	25	+
	26	+logger = logging.getLogger(__name__)
	27	+
	28	+EPS = 1e-12
	29	+
	30	+
	31	+# ---------------------------------------------------------------- features
	32	+def har_features(df: pd.DataFrame, measure: str = "rv5ss") -> pd.DataFrame:
	33	+ """Build the HAR information set from a daily realized-measure panel.
	34	+
	35	+ Parameters
	36	+ ----------
	37	+ df : pandas.DataFrame
	38	+ Daily measures of one instrument (output of
	39	+ :func:`wp12.realized_vol.build_daily_rv`).
	40	+ measure : str
	41	+ Column used as the realized-variance measure.
	42	+
	43	+ Returns
	44	+ -------
	45	+ pandas.DataFrame
	46	+ Features at the forecast origin t: daily/weekly/monthly RV
	47	+ averages, continuous and jump components, semivariances,
	48	+ realized quarticity — all lagged consistently (values known
	49	+ at the end of day t).
	50	+ """
	51	+ rv = df[measure].clip(lower=EPS)
	52	+ out = pd.DataFrame(index=df.index)
	53	+ out["rv_d"] = rv
	54	+ out["rv_w"] = rv.rolling(5).mean()
	55	+ out["rv_m"] = rv.rolling(22).mean()
	56	+ out["j_d"] = df["jump"].clip(lower=0.0)
	57	+ out["c_d"] = df["cont"].clip(lower=EPS)
	58	+ out["c_w"] = out["c_d"].rolling(5).mean()
	59	+ out["c_m"] = out["c_d"].rolling(22).mean()
	60	+ out["rsp_d"] = df["rsp"].clip(lower=0.0)
	61	+ out["rsn_d"] = df["rsn"].clip(lower=0.0)
	62	+ out["rvq_d"] = rv * np.sqrt(df["rq"].clip(lower=0.0))
	63	+ return out
	64	+
	65	+
	66	+def horizon_target(rv: pd.Series, h: int) -> pd.Series:
	67	+ """Average realized variance over days t+1 .. t+h, indexed at origin t.
	68	+
	69	+ Parameters
	70	+ ----------
	71	+ rv : pandas.Series
	72	+ Daily realized variance.
	73	+ h : int
	74	+ Forecast horizon in trading days.
	75	+
	76	+ Returns
	77	+ -------
	78	+ pandas.Series
	79	+ ``mean(rv[t+1..t+h])`` aligned on the forecast origin ``t``.
	80	+ """
	81	+ fwd = rv.shift(-1).rolling(h).mean().shift(-(h - 1))
	82	+ fwd.name = f"target_h{h}"
	83	+ return fwd
	84	+
	85	+
	86	+HAR_SPECS: dict[str, list[str]] = {
	87	+ "HAR": ["rv_d", "rv_w", "rv_m"],
	88	+ "HAR-J": ["rv_d", "rv_w", "rv_m", "j_d"],
	89	+ "HAR-CJ": ["c_d", "c_w", "c_m", "j_d"],
	90	+ "SHAR": ["rsp_d", "rsn_d", "rv_w", "rv_m"],
	91	+ "HARQ": ["rv_d", "rvq_d", "rv_w", "rv_m"],
	92	+}
	93	+
	94	+
	95	+def _ols_forecast(
	96	+ x: np.ndarray, y: np.ndarray, x_new: np.ndarray
	97	+) -> tuple[float, np.ndarray]:
	98	+ """Fit OLS with intercept and predict one observation.
	99	+
	100	+ Returns
	101	+ -------
	102	+ tuple
	103	+ ``(prediction, coefficients)``.
	104	+ """
	105	+ xx = np.column_stack([np.ones(len(x)), x])
	106	+ beta, *_ = np.linalg.lstsq(xx, y, rcond=None)
	107	+ return float(np.concatenate([[1.0], x_new]) @ beta), beta
	108	+
	109	+
	110	+def rolling_har_forecasts(
	111	+ df: pd.DataFrame,
	112	+ spec: str,
	113	+ horizons: tuple[int, ...] = (1, 5, 22),
	114	+ measure: str = "rv5ss",
	115	+ window: int = 1000,
	116	+ refit_every: int = 22,
	117	+) -> pd.DataFrame:
	118	+ """Rolling out-of-sample forecasts for one HAR-family specification.
	119	+
	120	+ Direct projections: for each horizon h the average RV over
	121	+ ``t+1..t+h`` is regressed on time-t features; parameters are
	122	+ re-estimated every ``refit_every`` days on the trailing ``window``
	123	+ observations. Forecasts are clipped to the in-sample target range
	124	+ (the "insanity filter" of Bollerslev, Patton & Quaedvlieg 2016).
	125	+
	126	+ Parameters
	127	+ ----------
	128	+ df : pandas.DataFrame
	129	+ Daily realized measures of one instrument.
	130	+ spec : str
	131	+ One of ``HAR``, ``HAR-J``, ``HAR-CJ``, ``SHAR``, ``HARQ``,
	132	+ ``LogHAR``.
	133	+ horizons : tuple of int
	134	+ Forecast horizons (days).
	135	+ measure : str
	136	+ Realized-variance column to forecast.
	137	+ window, refit_every : int
	138	+ Rolling estimation window length and re-estimation frequency.
	139	+
	140	+ Returns
	141	+ -------
	142	+ pandas.DataFrame
	143	+ Long frame with columns ``date`` (forecast origin), ``h``,
	144	+ ``model``, ``forecast`` (average daily variance over the
	145	+ horizon).
	146	+ """
	147	+ is_log = spec == "LogHAR"
	148	+ feats = har_features(df, measure)
	149	+ cols = HAR_SPECS["HAR"] if is_log else HAR_SPECS[spec]
	150	+ x_all = np.log(feats[cols].clip(lower=EPS)) if is_log else feats[cols]
	151	+ rv = df[measure].clip(lower=EPS)
	152	+
	153	+ rows: list[tuple] = []
	154	+ for h in horizons:
	155	+ y_all = horizon_target(rv, h)
	156	+ y_fit = np.log(y_all.clip(lower=EPS)) if is_log else y_all
	157	+ data = pd.concat([x_all, y_fit.rename("y"), y_all.rename("y_raw")], axis=1)
	158	+ valid = data.notna().all(axis=1).to_numpy()
	159	+ beta = None
	160	+ sigma2 = 0.0
	161	+ lo = hi = np.nan
	162	+ for t in range(window + h, len(data)):
	163	+ # fitting sample: origins whose targets are fully realized by t
	164	+ if beta is None or (t - window - h) % refit_every == 0:
	165	+ sl = slice(max(0, t - h - window), t - h)
	166	+ sub = data.iloc[sl]
	167	+ sub = sub[valid[sl]]
	168	+ if len(sub) < 100:
	169	+ continue
	170	+ x = sub[cols].to_numpy()
	171	+ y = sub["y"].to_numpy()
	172	+ xx = np.column_stack([np.ones(len(x)), x])
	173	+ beta, *_ = np.linalg.lstsq(xx, y, rcond=None)
	174	+ resid = y - xx @ beta
	175	+ sigma2 = float(resid.var())
	176	+ lo, hi = float(sub["y_raw"].min()), float(sub["y_raw"].max())
	177	+ if beta is None or not valid[t]:
	178	+ continue
	179	+ x_new = np.concatenate([[1.0], data[cols].iloc[t].to_numpy()])
	180	+ pred = float(x_new @ beta)
	181	+ if is_log:
	182	+ pred = float(np.exp(pred + 0.5 * sigma2)) # lognormal smearing
	183	+ pred = float(np.clip(pred, lo, hi)) # insanity filter
	184	+ rows.append((data.index[t], h, spec, pred))
	185	+ out = pd.DataFrame(rows, columns=["date", "h", "model", "forecast"])
	186	+ logger.info("%s: %d forecasts", spec, len(out))
	187	+ return out
	188	+
	189	+
	190	+# ---------------------------------------------------------------- GARCH family
	191	+def _fit_garch(r_pct: np.ndarray, kind: str):
	192	+ """Fit one GARCH-family model on percentage returns (arch package)."""
	193	+ kw = {"GARCH": dict(vol="GARCH", p=1, q=1),
	194	+ "GJR": dict(vol="GARCH", p=1, o=1, q=1),
	195	+ "EGARCH": dict(vol="EGARCH", p=1, o=1, q=1)}[kind]
	196	+ am = arch_model(r_pct, mean="Constant", dist="normal", rescale=False, **kw)
	197	+ with warnings.catch_warnings():
	198	+ warnings.simplefilter("ignore")
	199	+ return am.fit(disp="off", show_warning=False)
	200	+
	201	+
	202	+def rolling_garch_forecasts(
	203	+ df: pd.DataFrame,
	204	+ kind: str,
	205	+ horizons: tuple[int, ...] = (1, 5, 22),
	206	+ measure: str = "rv5ss",
	207	+ window: int = 1000,
	208	+ refit_every: int = 22,
	209	+ n_sim: int = 500,
	210	+ seed: int = 20260810,
	211	+) -> pd.DataFrame:
	212	+ """Rolling GARCH-family forecasts of average daily realized variance.
	213	+
	214	+ The model is fit on daily close-to-close returns; variance forecasts
	215	+ are mapped into RV units with the in-sample proportionality factor
	216	+ ``c = mean(RV) / mean(eps^2)`` (Hansen & Lunde 2005), which corrects
	217	+ for the overnight-variance wedge in session-limited markets.
	218	+
	219	+ Parameters
	220	+ ----------
	221	+ df : pandas.DataFrame
	222	+ Daily measures with ``ret_cc`` and the RV column.
	223	+ kind : str
	224	+ ``GARCH``, ``GJR`` or ``EGARCH``.
	225	+ horizons, measure, window, refit_every
	226	+ As in :func:`rolling_har_forecasts`.
	227	+ n_sim : int
	228	+ Simulation paths for multi-step EGARCH forecasts.
	229	+ seed : int
	230	+ RNG seed for simulated forecasts.
	231	+
	232	+ Returns
	233	+ -------
	234	+ pandas.DataFrame
	235	+ Long frame ``date, h, model, forecast``.
	236	+ """
	237	+ data = df[["ret_cc", measure]].dropna()
	238	+ r_pct = data["ret_cc"].to_numpy() * 100.0
	239	+ rv = data[measure].to_numpy()
	240	+ idx = data.index
	241	+ hmax = max(horizons)
	242	+
	243	+ rows: list[tuple] = []
	244	+ res = None
	245	+ scale_c = 1.0
	246	+ for t in range(window, len(data)):
	247	+ if res is None or (t - window) % refit_every == 0:
	248	+ r_win = r_pct[t - window:t]
	249	+ try:
	250	+ res = _fit_garch(r_win, kind)
	251	+ except Exception as exc: # noqa: BLE001 — keep last params
	252	+ logger.warning("%s fit failed at %s: %s", kind, idx[t], exc)
	253	+ if res is not None:
	254	+ eps2 = (r_win - r_win.mean()) ** 2 / 1e4
	255	+ scale_c = float(np.mean(rv[t - window:t]) / max(np.mean(eps2), EPS))
	256	+ if res is None:
	257	+ continue
	258	+ try:
	259	+ method = "simulation" if kind == "EGARCH" else "analytic"
	260	+ fc = res.forecast(
	261	+ horizon=hmax, reindex=False, start=None,
	262	+ params=res.params,
	263	+ method=method,
	264	+ simulations=n_sim,
	265	+ x=None,
	266	+ random_state=np.random.RandomState(seed) if method == "simulation" else None,
	267	+ align="origin",
	268	+ )
	269	+ # variance path for the last in-sample date of this window slice
	270	+ var_path = fc.variance.iloc[-1].to_numpy() / 1e4 # back to raw units
	271	+ except Exception as exc: # noqa: BLE001
	272	+ logger.warning("%s forecast failed at %s: %s", kind, idx[t], exc)
	273	+ continue
	274	+ for h in horizons:
	275	+ rows.append((idx[t], h, kind, scale_c * float(np.mean(var_path[:h]))))
	276	+ out = pd.DataFrame(rows, columns=["date", "h", "model", "forecast"])
	277	+ logger.info("%s: %d forecasts", kind, len(out))
	278	+ return out
	279	+
	280	+
	281	+# ---------------------------------------------------------------- Realized GARCH
	282	+def _rg_negll(params: np.ndarray, r: np.ndarray, logx: np.ndarray) -> float:
	283	+ """Negative log-likelihood of the log-linear Realized GARCH(1,1)."""
	284	+ omega, beta, gamma, xi, phi, tau1, tau2, log_su = params
	285	+ su2 = np.exp(2 * log_su)
	286	+ n = len(r)
	287	+ logh = np.empty(n)
	288	+ logh[0] = np.mean(logx)
	289	+ ll = 0.0
	290	+ for t in range(n):
	291	+ if t > 0:
	292	+ logh[t] = omega + beta * logh[t - 1] + gamma * logx[t - 1]
	293	+ h = np.exp(logh[t])
	294	+ z = r[t] / np.sqrt(h)
	295	+ u = logx[t] - xi - phi * logh[t] - tau1 * z - tau2 * (z * z - 1.0)
	296	+ ll += -0.5 * (logh[t] + z * z) - 0.5 * (np.log(su2) + u * u / su2)
	297	+ return -ll
	298	+
	299	+
	300	+def _rg_fit(r: np.ndarray, logx: np.ndarray) -> np.ndarray \| None:
	301	+ """Fit the Realized GARCH by quasi-maximum likelihood (L-BFGS-B)."""
	302	+ x0 = np.array([0.06 * np.mean(logx), 0.55, 0.40, -0.20, 1.00, -0.07, 0.07, np.log(0.4)])
	303	+ bounds = [(-5, 5), (0.0, 0.99), (0.0, 0.99), (-5, 5),
	304	+ (0.5, 1.5), (-1, 1), (-1, 1), (np.log(0.05), np.log(2.0))]
	305	+ with warnings.catch_warnings():
	306	+ warnings.simplefilter("ignore")
	307	+ res = optimize.minimize(
	308	+ _rg_negll, x0, args=(r, logx), method="L-BFGS-B",
	309	+ bounds=bounds, options={"maxiter": 300},
	310	+ )
	311	+ if not np.all(np.isfinite(res.x)) or (res.x[1] + res.x[2] * res.x[4]) >= 0.999:
	312	+ return None
	313	+ return res.x
	314	+
	315	+
	316	+def rolling_realgarch_forecasts(
	317	+ df: pd.DataFrame,
	318	+ horizons: tuple[int, ...] = (1, 5, 22),
	319	+ measure: str = "rv5ss",
	320	+ window: int = 1000,
	321	+ refit_every: int = 22,
	322	+) -> pd.DataFrame:
	323	+ """Rolling forecasts from the log-linear Realized GARCH (Hansen et al. 2012).
	324	+
	325	+ The measurement equation links the latent conditional variance to the
	326	+ realized measure, so the model forecasts RV directly:
	327	+ ``E[x_{t+k}]`` is computed analytically with a lognormal correction.
	328	+
	329	+ Parameters
	330	+ ----------
	331	+ df : pandas.DataFrame
	332	+ Daily measures with ``ret_cc`` and the RV column.
	333	+ horizons, measure, window, refit_every
	334	+ As in :func:`rolling_har_forecasts`.
	335	+
	336	+ Returns
	337	+ -------
	338	+ pandas.DataFrame
	339	+ Long frame ``date, h, model, forecast``.
	340	+ """
	341	+ data = df[["ret_cc", measure]].dropna()
	342	+ r = data["ret_cc"].to_numpy() * 100.0
	343	+ logx = np.log(data[measure].to_numpy() * 1e4).clip(-30, 30) # pct^2 units
	344	+ idx = data.index
	345	+
	346	+ rows: list[tuple] = []
	347	+ params = None
	348	+ for t in range(window, len(data)):
	349	+ if params is None or (t - window) % refit_every == 0:
	350	+ new = _rg_fit(r[t - window:t], logx[t - window:t])
	351	+ if new is not None:
	352	+ params = new
	353	+ elif params is None:
	354	+ continue
	355	+ if params is None:
	356	+ continue
	357	+ omega, beta, gamma, xi, phi, tau1, tau2, log_su = params
	358	+ su2 = np.exp(2 * log_su)
	359	+ vu = tau1*2 + 2 tau2**2 + su2 # variance of the logx innovation
	360	+ rho = beta + gamma * phi # persistence of log h
	361	+ # filter log h_t through the window to get the current state
	362	+ logh = np.mean(logx[t - window:t])
	363	+ for s in range(t - window + 1, t):
	364	+ logh = omega + beta * logh + gamma * logx[s - 1]
	365	+ mu = omega + beta * logh + gamma * logx[t - 1] # E[log h_{t+1}]
	366	+ v = 0.0
	367	+ fpath = []
	368	+ for _k in range(max(horizons)):
	369	+ m_x = xi + phi * mu
	370	+ v_x = phi*2 v + vu
	371	+ fpath.append(np.exp(m_x + 0.5 * v_x) / 1e4) # E[x_{t+k}] raw units
	372	+ mu = omega + gamma * xi + rho * mu
	373	+ v = rho*2 v + gamma*2 vu
	374	+ fpath_arr = np.array(fpath)
	375	+ for h in horizons:
	376	+ rows.append((idx[t], h, "RealGARCH", float(np.mean(fpath_arr[:h]))))
	377	+ out = pd.DataFrame(rows, columns=["date", "h", "model", "forecast"])
	378	+ logger.info("RealGARCH: %d forecasts", len(out))
	379	+ return out

added src/wp12/models_ml.py +619 −0

		`@@ -0,0 +1,619 @@`
	1	+"""
	2	+================================================================
	3	+Auteur : Simon-Pierre Boucher
	4	+Contact : contact@spboucher.ai
	5	+Projet : Prévision de volatilité réalisée multi-actifs
	6	+ (HAR-RV vs GARCH vs Machine Learning)
	7	+Fichier : models_ml.py
	8	+Description : Modèles de machine learning — LASSO/Ridge/Elastic
	9	+ Net, Random Forest, XGBoost, LightGBM (rolling,
	10	+ tuning walk-forward), LSTM et Transformer (multi-
	11	+ horizon), modèle pooled multi-actifs avec
	12	+ identifiants d'actif, SHAP et importance par
	13	+ permutation. Cible : log RV moyen sur l'horizon,
	14	+ re-transformé par lissage de Duan.
	15	+================================================================
	16	+"""
	17	+
	18	+from __future__ import annotations
	19	+
	20	+import logging
	21	+from typing import Any
	22	+
	23	+import numpy as np
	24	+import pandas as pd
	25	+
	26	+from . import config
	27	+from .models_econ import horizon_target
	28	+
	29	+logger = logging.getLogger(__name__)
	30	+
	31	+EPS = 1e-12
	32	+
	33	+HAR_FEATURES = ["lrv_d", "lrv_w", "lrv_m"]
	34	+
	35	+
	36	+# ---------------------------------------------------------------- features
	37	+def ml_features(
	38	+ panel: pd.DataFrame, measure: str = "rv5ss"
	39	+) -> dict[str, pd.DataFrame]:
	40	+ """Build per-ticker ML feature frames from the multi-asset panel.
	41	+
	42	+ Two nested information sets are produced: the HAR set
	43	+ (``lrv_d/w/m``, logs of the daily/weekly/monthly RV averages) and an
	44	+ extended set adding jumps, semivariances, realized quarticity,
	45	+ volume, returns, calendar effects, lagged VIX and cross-asset RV
	46	+ averages. All predictors are known at the end of day t.
	47	+
	48	+ Parameters
	49	+ ----------
	50	+ panel : pandas.DataFrame
	51	+ Long panel of daily realized measures (all tickers).
	52	+ measure : str
	53	+ Realized-variance column.
	54	+
	55	+ Returns
	56	+ -------
	57	+ dict
	58	+ ``ticker -> DataFrame`` of features indexed by date.
	59	+ """
	60	+ # cross-asset predictors on the union calendar
	61	+ rv_wide = panel.pivot_table(index=panel.index, columns="ticker", values=measure)
	62	+ cls_map = panel.groupby("ticker")["cls"].first()
	63	+ cross = {}
	64	+ for cls in ["equity", "fx", "crypto", "futures"]:
	65	+ cols = [c for c in rv_wide.columns if cls_map.get(c) == cls]
	66	+ if cols:
	67	+ cross[f"lrv_cross_{cls}"] = np.log(rv_wide[cols].clip(lower=EPS)).mean(axis=1)
	68	+ cross_df = pd.DataFrame(cross)
	69	+ vix = None
	70	+ if "VIX" in set(panel["ticker"]):
	71	+ vix = np.log(panel.loc[panel["ticker"] == "VIX", "close"]).rename("lvix")
	72	+
	73	+ out: dict[str, pd.DataFrame] = {}
	74	+ for tk, df in panel.groupby("ticker"):
	75	+ if cls_map[tk] == "index":
	76	+ continue
	77	+ df = df.sort_index()
	78	+ rv = df[measure].clip(lower=EPS)
	79	+ f = pd.DataFrame(index=df.index)
	80	+ f["lrv_d"] = np.log(rv)
	81	+ f["lrv_w"] = np.log(rv.rolling(5).mean())
	82	+ f["lrv_m"] = np.log(rv.rolling(22).mean())
	83	+ # extended block
	84	+ f["jump_share"] = (df["jump"] / rv).clip(0, 1)
	85	+ f["lrsp"] = np.log(df["rsp"].clip(lower=EPS))
	86	+ f["lrsn"] = np.log(df["rsn"].clip(lower=EPS))
	87	+ f["sj_norm"] = (df["sj"] / rv).clip(-1, 1)
	88	+ f["lrq"] = np.log(df["rq"].clip(lower=EPS))
	89	+ f["ret_d"] = df["ret_cc"]
	90	+ f["ret_w"] = df["ret_cc"].rolling(5).sum()
	91	+ f["ret_m"] = df["ret_cc"].rolling(22).sum()
	92	+ f["lvol_d"] = np.log(df["volume"].clip(lower=1.0)).diff()
	93	+ f["dow"] = df.index.dayofweek.astype(float)
	94	+ f = f.join(cross_df.reindex(f.index).ffill(limit=5))
	95	+ if vix is not None:
	96	+ f = f.join(vix.reindex(f.index).ffill(limit=5))
	97	+ # target inputs
	98	+ f["_rv"] = rv
	99	+ f["_cls"] = cls_map[tk]
	100	+ out[tk] = f
	101	+ return out
	102	+
	103	+
	104	+def feature_columns(feature_set: str, frame: pd.DataFrame) -> list[str]:
	105	+ """Column list for a named feature set (``har`` or ``extended``)."""
	106	+ if feature_set == "har":
	107	+ return HAR_FEATURES
	108	+ return [c for c in frame.columns if not c.startswith("_")]
	109	+
	110	+
	111	+# ---------------------------------------------------------------- estimators
	112	+def _make_estimator(name: str, params: dict[str, Any], seed: int):
	113	+ """Instantiate a tabular estimator by name with given hyperparameters."""
	114	+ from lightgbm import LGBMRegressor
	115	+ from sklearn.ensemble import RandomForestRegressor
	116	+ from sklearn.linear_model import ElasticNet, Lasso, Ridge
	117	+ from sklearn.pipeline import make_pipeline
	118	+ from sklearn.preprocessing import StandardScaler
	119	+ from xgboost import XGBRegressor
	120	+
	121	+ if name == "Ridge":
	122	+ return make_pipeline(StandardScaler(), Ridge(alpha=params["alpha"]))
	123	+ if name == "LASSO":
	124	+ return make_pipeline(StandardScaler(), Lasso(alpha=params["alpha"], max_iter=5000))
	125	+ if name == "ElasticNet":
	126	+ return make_pipeline(
	127	+ StandardScaler(),
	128	+ ElasticNet(alpha=params["alpha"], l1_ratio=params["l1_ratio"], max_iter=5000),
	129	+ )
	130	+ if name == "RF":
	131	+ return RandomForestRegressor(
	132	+ n_estimators=300, min_samples_leaf=params["min_samples_leaf"],
	133	+ max_features=params["max_features"], random_state=seed, n_jobs=-1,
	134	+ )
	135	+ if name == "XGBoost":
	136	+ return XGBRegressor(
	137	+ n_estimators=400, learning_rate=0.05, max_depth=params["max_depth"],
	138	+ subsample=0.8, colsample_bytree=0.8, random_state=seed,
	139	+ n_jobs=-1, verbosity=0,
	140	+ )
	141	+ if name == "LightGBM":
	142	+ return LGBMRegressor(
	143	+ n_estimators=400, learning_rate=0.05, num_leaves=params["num_leaves"],
	144	+ min_child_samples=20, subsample=0.8, colsample_bytree=0.8,
	145	+ random_state=seed, n_jobs=-1, verbose=-1,
	146	+ )
	147	+ raise ValueError(f"unknown estimator {name}")
	148	+
	149	+
	150	+GRIDS: dict[str, list[dict[str, Any]]] = {
	151	+ "Ridge": [{"alpha": a} for a in (0.1, 1.0, 10.0)],
	152	+ "LASSO": [{"alpha": a} for a in (0.0005, 0.005, 0.05)],
	153	+ "ElasticNet": [{"alpha": a, "l1_ratio": r} for a in (0.001, 0.01) for r in (0.3, 0.7)],
	154	+ "RF": [{"min_samples_leaf": m, "max_features": f}
	155	+ for m in (5, 20) for f in (0.33, "sqrt")],
	156	+ "XGBoost": [{"max_depth": d} for d in (3, 5)],
	157	+ "LightGBM": [{"num_leaves": n} for n in (15, 31)],
	158	+}
	159	+
	160	+
	161	+def _tune(name: str, x: np.ndarray, y: np.ndarray, seed: int) -> dict[str, Any]:
	162	+ """Walk-forward hyperparameter selection inside the training window.
	163	+
	164	+ The last 20% of the window (in time order) serves as validation for
	165	+ a small grid; the winner is refit on the full window by the caller.
	166	+ """
	167	+ split = int(len(x) * 0.8)
	168	+ best, best_loss = GRIDS[name][0], np.inf
	169	+ for params in GRIDS[name]:
	170	+ est = _make_estimator(name, params, seed)
	171	+ est.fit(x[:split], y[:split])
	172	+ pred = est.predict(x[split:])
	173	+ loss = float(np.mean((pred - y[split:]) ** 2))
	174	+ if loss < best_loss:
	175	+ best, best_loss = params, loss
	176	+ return best
	177	+
	178	+
	179	+def rolling_ml_forecasts(
	180	+ frame: pd.DataFrame,
	181	+ model_name: str,
	182	+ feature_set: str = "extended",
	183	+ horizons: tuple[int, ...] = (1, 5, 22),
	184	+ window: int = 1000,
	185	+ refit_every: int = 22,
	186	+ tune_every: int = 12,
	187	+ seed: int \| None = None,
	188	+) -> pd.DataFrame:
	189	+ """Rolling out-of-sample forecasts for one tabular ML model.
	190	+
	191	+ The target is ``log`` average RV over the horizon; predictions are
	192	+ mapped back to variance units with Duan's (1983) smearing factor
	193	+ computed on training residuals, then clipped to the in-sample target
	194	+ range (insanity filter).
	195	+
	196	+ Parameters
	197	+ ----------
	198	+ frame : pandas.DataFrame
	199	+ Feature frame of one ticker (from :func:`ml_features`).
	200	+ model_name : str
	201	+ ``LASSO``, ``Ridge``, ``ElasticNet``, ``RF``, ``XGBoost`` or
	202	+ ``LightGBM``.
	203	+ feature_set : str
	204	+ ``har`` (HAR information set) or ``extended``.
	205	+ horizons, window, refit_every
	206	+ Standard rolling protocol.
	207	+ tune_every : int
	208	+ Hyperparameters re-tuned every ``tune_every`` refits.
	209	+ seed : int, optional
	210	+ RNG seed (defaults to the config seed).
	211	+
	212	+ Returns
	213	+ -------
	214	+ pandas.DataFrame
	215	+ Long frame ``date, h, model, forecast``.
	216	+ """
	217	+ seed = seed if seed is not None else config.load_config()["seed"]
	218	+ cols = feature_columns(feature_set, frame)
	219	+ label = f"{model_name}-{'X' if feature_set == 'extended' else 'H'}"
	220	+ rv = frame["_rv"]
	221	+ rows: list[tuple] = []
	222	+ for h in horizons:
	223	+ y_raw = horizon_target(rv, h)
	224	+ data = pd.concat(
	225	+ [frame[cols], np.log(y_raw.clip(lower=EPS)).rename("y"),
	226	+ y_raw.rename("y_raw")], axis=1,
	227	+ )
	228	+ ok = data.notna().all(axis=1).to_numpy()
	229	+ est = None
	230	+ smear = 1.0
	231	+ lo = hi = np.nan
	232	+ params: dict[str, Any] \| None = None
	233	+ n_refit = 0
	234	+ for t in range(window + h, len(data)):
	235	+ if est is None or (t - window - h) % refit_every == 0:
	236	+ sl = slice(max(0, t - h - window), t - h)
	237	+ sub = data.iloc[sl][ok[sl]]
	238	+ if len(sub) < 200:
	239	+ continue
	240	+ x = sub[cols].to_numpy()
	241	+ y = sub["y"].to_numpy()
	242	+ if params is None or n_refit % tune_every == 0:
	243	+ params = _tune(model_name, x, y, seed)
	244	+ est = _make_estimator(model_name, params, seed)
	245	+ est.fit(x, y)
	246	+ resid = y - est.predict(x)
	247	+ smear = float(np.mean(np.exp(resid)))
	248	+ lo, hi = float(sub["y_raw"].min()), float(sub["y_raw"].max())
	249	+ n_refit += 1
	250	+ if est is None or not ok[t]:
	251	+ continue
	252	+ pred_log = float(est.predict(data[cols].iloc[[t]].to_numpy())[0])
	253	+ pred = float(np.clip(np.exp(pred_log) * smear, lo, hi))
	254	+ rows.append((data.index[t], h, label, pred))
	255	+ out = pd.DataFrame(rows, columns=["date", "h", "model", "forecast"])
	256	+ logger.info("%s: %d forecasts", label, len(out))
	257	+ return out
	258	+
	259	+
	260	+# ---------------------------------------------------------------- pooled model
	261	+def rolling_pooled_lgbm(
	262	+ frames: dict[str, pd.DataFrame],
	263	+ horizons: tuple[int, ...] = (1, 5, 22),
	264	+ window: int = 1000,
	265	+ refit_every: int = 22,
	266	+ seed: int \| None = None,
	267	+ exclude: str \| None = None,
	268	+) -> pd.DataFrame:
	269	+ """Pooled LightGBM trained on all assets jointly (asset identifiers).
	270	+
	271	+ One model per horizon is fit on the stacked window observations of
	272	+ every ticker, with ``ticker`` and asset class as categorical
	273	+ features. Set ``exclude`` to hold one ticker out of training while
	274	+ still predicting it (cross-asset transferability experiment).
	275	+
	276	+ Parameters
	277	+ ----------
	278	+ frames : dict
	279	+ Output of :func:`ml_features`.
	280	+ horizons, window, refit_every
	281	+ Standard rolling protocol.
	282	+ seed : int, optional
	283	+ RNG seed.
	284	+ exclude : str, optional
	285	+ Ticker excluded from the training pool (leave-one-out).
	286	+
	287	+ Returns
	288	+ -------
	289	+ pandas.DataFrame
	290	+ Long frame ``date, ticker, h, model, forecast``.
	291	+ """
	292	+ from lightgbm import LGBMRegressor
	293	+
	294	+ seed = seed if seed is not None else config.load_config()["seed"]
	295	+ label = "Pooled-LGBM" if exclude is None else "Pooled-LGBM-LOO"
	296	+ stacked = []
	297	+ for tk, f in frames.items():
	298	+ cols = feature_columns("extended", f)
	299	+ g = f.copy()
	300	+ g["_ticker"] = tk
	301	+ stacked.append(g)
	302	+ allf = pd.concat(stacked).sort_index()
	303	+ cols = feature_columns("extended", allf.drop(columns=["_ticker"]))
	304	+
	305	+ rows: list[tuple] = []
	306	+ for h in horizons:
	307	+ parts = []
	308	+ for tk, f in frames.items():
	309	+ y_raw = horizon_target(f["_rv"], h)
	310	+ part = f.copy()
	311	+ part["_y"] = np.log(y_raw.clip(lower=EPS))
	312	+ part["_y_raw"] = y_raw
	313	+ part["_ticker"] = tk
	314	+ parts.append(part)
	315	+ data = pd.concat(parts).sort_index()
	316	+ data["_tk_cat"] = data["_ticker"].astype("category")
	317	+ data["_cls_cat"] = data["_cls"].astype("category")
	318	+ xcols = cols + ["_tk_cat", "_cls_cat"]
	319	+ dates = data.index.unique().sort_values()
	320	+ est = None
	321	+ smear = 1.0
	322	+ bounds: dict[str, tuple[float, float]] = {}
	323	+ for i in range(window + h, len(dates)):
	324	+ if est is None or (i - window - h) % refit_every == 0:
	325	+ d_lo, d_hi = dates[max(0, i - h - window)], dates[i - h]
	326	+ sub = data[(data.index >= d_lo) & (data.index < d_hi)].dropna(
	327	+ subset=cols + ["_y"]
	328	+ )
	329	+ if exclude is not None:
	330	+ sub = sub[sub["_ticker"] != exclude]
	331	+ if len(sub) < 1000:
	332	+ continue
	333	+ est = LGBMRegressor(
	334	+ n_estimators=500, learning_rate=0.05, num_leaves=63,
	335	+ min_child_samples=50, subsample=0.8, colsample_bytree=0.8,
	336	+ random_state=seed, n_jobs=-1, verbose=-1,
	337	+ )
	338	+ est.fit(sub[xcols], sub["_y"], categorical_feature=["_tk_cat", "_cls_cat"])
	339	+ resid = sub["_y"] - est.predict(sub[xcols])
	340	+ smear = float(np.mean(np.exp(resid)))
	341	+ bounds = {
	342	+ tk: (float(g["_y_raw"].min()), float(g["_y_raw"].max()))
	343	+ for tk, g in sub.groupby("_ticker")
	344	+ }
	345	+ if est is None:
	346	+ continue
	347	+ today = data[data.index == dates[i]].dropna(subset=cols)
	348	+ if exclude is not None:
	349	+ today = today[today["_ticker"] == exclude]
	350	+ if today.empty:
	351	+ continue
	352	+ preds = np.exp(est.predict(today[xcols])) * smear
	353	+ for (tk, p, yr) in zip(today["_ticker"], preds, today["_y_raw"]):
	354	+ lo, hi = bounds.get(tk, (np.nan, np.nan))
	355	+ if np.isfinite(lo):
	356	+ p = float(np.clip(p, lo, hi))
	357	+ rows.append((dates[i], tk, h, label, float(p)))
	358	+ out = pd.DataFrame(rows, columns=["date", "ticker", "h", "model", "forecast"])
	359	+ logger.info("%s%s: %d forecasts", label, f" (excl {exclude})" if exclude else "", len(out))
	360	+ return out
	361	+
	362	+
	363	+# ---------------------------------------------------------------- deep learning
	364	+def _make_sequences(
	365	+ values: np.ndarray, targets: np.ndarray, seq_len: int
	366	+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
	367	+ """Stack rolling sequences; returns (X, y, end_positions)."""
	368	+ xs, ys, pos = [], [], []
	369	+ for t in range(seq_len - 1, len(values)):
	370	+ if np.isnan(values[t - seq_len + 1: t + 1]).any() or np.isnan(targets[t]).any():
	371	+ continue
	372	+ xs.append(values[t - seq_len + 1: t + 1])
	373	+ ys.append(targets[t])
	374	+ pos.append(t)
	375	+ return np.asarray(xs, dtype=np.float32), np.asarray(ys, dtype=np.float32), np.asarray(pos)
	376	+
	377	+
	378	+def _torch_device():
	379	+ import torch
	380	+
	381	+ if torch.backends.mps.is_available():
	382	+ return torch.device("mps")
	383	+ return torch.device("cpu")
	384	+
	385	+
	386	+def _build_net(arch: str, n_feat: int, n_out: int, cfg: dict[str, Any]):
	387	+ """Construct the LSTM or Transformer forecasting network."""
	388	+ import torch.nn as nn
	389	+
	390	+ if arch == "LSTM":
	391	+ class LSTMNet(nn.Module):
	392	+ def __init__(self) -> None:
	393	+ super().__init__()
	394	+ self.lstm = nn.LSTM(n_feat, cfg["hidden"], cfg["layers"], batch_first=True)
	395	+ self.head = nn.Linear(cfg["hidden"], n_out)
	396	+
	397	+ def forward(self, x):
	398	+ out, _ = self.lstm(x)
	399	+ return self.head(out[:, -1])
	400	+
	401	+ return LSTMNet()
	402	+
	403	+ class TransformerNet(nn.Module):
	404	+ def __init__(self) -> None:
	405	+ super().__init__()
	406	+ d = cfg["d_model"]
	407	+ self.proj = nn.Linear(n_feat, d)
	408	+ self.pos = nn.Parameter(__import__("torch").zeros(1, cfg["seq_len"], d))
	409	+ layer = nn.TransformerEncoderLayer(
	410	+ d_model=d, nhead=cfg["heads"], dim_feedforward=4 * d,
	411	+ batch_first=True, dropout=0.1,
	412	+ )
	413	+ self.enc = nn.TransformerEncoder(layer, cfg["blocks"])
	414	+ self.head = nn.Linear(d, n_out)
	415	+
	416	+ def forward(self, x):
	417	+ z = self.proj(x) + self.pos[:, : x.shape[1]]
	418	+ z = self.enc(z)
	419	+ return self.head(z.mean(dim=1))
	420	+
	421	+ return TransformerNet()
	422	+
	423	+
	424	+def _train_net(net, x_tr, y_tr, x_va, y_va, epochs: int, patience: int, seed: int):
	425	+ """Train with Adam + early stopping on validation MSE; returns the net."""
	426	+ import torch
	427	+
	428	+ torch.manual_seed(seed)
	429	+ dev = _torch_device()
	430	+ net = net.to(dev)
	431	+ opt = torch.optim.Adam(net.parameters(), lr=1e-3)
	432	+ lossf = torch.nn.MSELoss()
	433	+ xt = torch.tensor(x_tr, device=dev)
	434	+ yt = torch.tensor(y_tr, device=dev)
	435	+ xv = torch.tensor(x_va, device=dev)
	436	+ yv = torch.tensor(y_va, device=dev)
	437	+ best_state, best_val, bad = None, np.inf, 0
	438	+ n = len(xt)
	439	+ for _ep in range(epochs):
	440	+ net.train()
	441	+ perm = torch.randperm(n, device=dev)
	442	+ for i in range(0, n, 128):
	443	+ b = perm[i:i + 128]
	444	+ opt.zero_grad()
	445	+ loss = lossf(net(xt[b]), yt[b])
	446	+ loss.backward()
	447	+ opt.step()
	448	+ net.eval()
	449	+ with torch.no_grad():
	450	+ val = float(lossf(net(xv), yv))
	451	+ if val < best_val - 1e-5:
	452	+ best_val, bad = val, 0
	453	+ best_state = {k: v.detach().clone() for k, v in net.state_dict().items()}
	454	+ else:
	455	+ bad += 1
	456	+ if bad >= patience:
	457	+ break
	458	+ if best_state is not None:
	459	+ net.load_state_dict(best_state)
	460	+ net.eval()
	461	+ return net
	462	+
	463	+
	464	+def rolling_deep_forecasts(
	465	+ frame: pd.DataFrame,
	466	+ arch: str,
	467	+ horizons: tuple[int, ...] = (1, 5, 22),
	468	+ window: int = 1000,
	469	+ refit_every: int = 66,
	470	+ seed: int \| None = None,
	471	+) -> pd.DataFrame:
	472	+ """Rolling multi-horizon forecasts from an LSTM or Transformer.
	473	+
	474	+ One network with a 3-dimensional output head (one per horizon) is
	475	+ trained on sequences of standardized features; refits are quarterly
	476	+ (``refit_every=66``) to keep the computational burden manageable.
	477	+ Targets are log average RV; back-transform uses Duan smearing.
	478	+
	479	+ Parameters
	480	+ ----------
	481	+ frame : pandas.DataFrame
	482	+ Feature frame of one ticker (from :func:`ml_features`).
	483	+ arch : str
	484	+ ``LSTM`` or ``Transformer``.
	485	+ horizons, window, refit_every
	486	+ Rolling protocol (quarterly refit by default).
	487	+ seed : int, optional
	488	+ RNG seed.
	489	+
	490	+ Returns
	491	+ -------
	492	+ pandas.DataFrame
	493	+ Long frame ``date, h, model, forecast``.
	494	+ """
	495	+ import torch
	496	+
	497	+ cfg_all = config.load_config()["models"]["ml"]
	498	+ cfg = dict(cfg_all["lstm"] if arch == "LSTM" else cfg_all["transformer"])
	499	+ cfg.setdefault("seq_len", 22)
	500	+ seed = seed if seed is not None else config.load_config()["seed"]
	501	+ seq_len = int(cfg["seq_len"])
	502	+ cols = feature_columns("extended", frame)
	503	+ rv = frame["_rv"]
	504	+ hmax_targets = []
	505	+ for h in horizons:
	506	+ y = horizon_target(rv, h)
	507	+ hmax_targets.append(np.log(y.clip(lower=EPS)))
	508	+ Y = np.column_stack([y.to_numpy() for y in hmax_targets])
	509	+ Yraw = np.column_stack([horizon_target(rv, h).to_numpy() for h in horizons])
	510	+ X = frame[cols].to_numpy(dtype=np.float64)
	511	+ idx = frame.index
	512	+
	513	+ rows: list[tuple] = []
	514	+ net = None
	515	+ mu = sd = None
	516	+ smear = np.ones(len(horizons))
	517	+ lo = np.full(len(horizons), np.nan)
	518	+ hi = np.full(len(horizons), np.nan)
	519	+ hmax = max(horizons)
	520	+ for t in range(window + hmax, len(frame)):
	521	+ if net is None or (t - window - hmax) % refit_every == 0:
	522	+ sl = slice(max(0, t - hmax - window), t - hmax)
	523	+ Xw, Yw = X[sl], Y[sl]
	524	+ mu = np.nanmean(Xw, axis=0)
	525	+ sd = np.nanstd(Xw, axis=0) + 1e-9
	526	+ Xs = (Xw - mu) / sd
	527	+ xs, ys, _pos = _make_sequences(Xs, Yw, seq_len)
	528	+ if len(xs) < 300:
	529	+ continue
	530	+ split = int(len(xs) * 0.85)
	531	+ net = _build_net(arch, xs.shape[2], len(horizons), cfg)
	532	+ net = _train_net(net, xs[:split], ys[:split], xs[split:], ys[split:],
	533	+ int(cfg["epochs"]), int(cfg["patience"]), seed)
	534	+ with torch.no_grad():
	535	+ fit_pred = net(torch.tensor(xs, device=_torch_device())).cpu().numpy()
	536	+ resid = ys - fit_pred
	537	+ smear = np.exp(resid).mean(axis=0)
	538	+ Yraw_w = Yraw[sl]
	539	+ lo = np.nanmin(Yraw_w, axis=0)
	540	+ hi = np.nanmax(Yraw_w, axis=0)
	541	+ if net is None:
	542	+ continue
	543	+ x_seq = (X[t - seq_len + 1: t + 1] - mu) / sd
	544	+ if np.isnan(x_seq).any():
	545	+ continue
	546	+ with torch.no_grad():
	547	+ p_log = net(
	548	+ torch.tensor(x_seq[None].astype(np.float32), device=_torch_device())
	549	+ ).cpu().numpy()[0]
	550	+ for j, h in enumerate(horizons):
	551	+ pred = float(np.clip(np.exp(p_log[j]) * smear[j], lo[j], hi[j]))
	552	+ rows.append((idx[t], h, arch, pred))
	553	+ out = pd.DataFrame(rows, columns=["date", "h", "model", "forecast"])
	554	+ logger.info("%s: %d forecasts", arch, len(out))
	555	+ return out
	556	+
	557	+
	558	+# ---------------------------------------------------------------- interpretation
	559	+def shap_and_permutation(
	560	+ frames: dict[str, pd.DataFrame],
	561	+ h: int = 1,
	562	+ window: int = 1000,
	563	+ seed: int \| None = None,
	564	+) -> tuple[pd.DataFrame, pd.DataFrame]:
	565	+ """Global SHAP values and permutation importance of the pooled LightGBM.
	566	+
	567	+ A pooled LightGBM is fit on the most recent ``window`` days of every
	568	+ ticker; mean absolute SHAP values (TreeExplainer) and permutation
	569	+ importances (on a temporal validation slice) are reported per
	570	+ feature.
	571	+
	572	+ Returns
	573	+ -------
	574	+ tuple of pandas.DataFrame
	575	+ ``(shap_table, permutation_table)`` sorted by importance.
	576	+ """
	577	+ import shap
	578	+ from lightgbm import LGBMRegressor
	579	+ from sklearn.inspection import permutation_importance
	580	+
	581	+ seed = seed if seed is not None else config.load_config()["seed"]
	582	+ parts = []
	583	+ for tk, f in frames.items():
	584	+ y_raw = horizon_target(f["_rv"], h)
	585	+ part = f.copy()
	586	+ part["_y"] = np.log(y_raw.clip(lower=EPS))
	587	+ part["_ticker"] = tk
	588	+ parts.append(part)
	589	+ data = pd.concat(parts).sort_index()
	590	+ cols = feature_columns("extended", data.drop(columns=["_ticker"]))
	591	+ data = data.dropna(subset=cols + ["_y"])
	592	+ dates = data.index.unique().sort_values()
	593	+ data = data[data.index >= dates[-window]]
	594	+ split_date = data.index.unique().sort_values()[int(data.index.nunique() * 0.85)]
	595	+ tr, va = data[data.index < split_date], data[data.index >= split_date]
	596	+
	597	+ est = LGBMRegressor(
	598	+ n_estimators=500, learning_rate=0.05, num_leaves=63, min_child_samples=50,
	599	+ subsample=0.8, colsample_bytree=0.8, random_state=seed, n_jobs=-1, verbose=-1,
	600	+ )
	601	+ est.fit(tr[cols], tr["_y"])
	602	+
	603	+ explainer = shap.TreeExplainer(est)
	604	+ sv = explainer.shap_values(va[cols])
	605	+ shap_tab = (
	606	+ pd.DataFrame({"feature": cols, "mean_abs_shap": np.abs(sv).mean(axis=0)})
	607	+ .sort_values("mean_abs_shap", ascending=False)
	608	+ .reset_index(drop=True)
	609	+ )
	610	+ pi = permutation_importance(
	611	+ est, va[cols], va["_y"], n_repeats=10, random_state=seed, n_jobs=-1
	612	+ )
	613	+ perm_tab = (
	614	+ pd.DataFrame({"feature": cols, "perm_importance": pi.importances_mean,
	615	+ "perm_std": pi.importances_std})
	616	+ .sort_values("perm_importance", ascending=False)
	617	+ .reset_index(drop=True)
	618	+ )
	619	+ return shap_tab, perm_tab

added src/wp12/realized_vol.py +317 −0

		`@@ -0,0 +1,317 @@`
	1	+"""
	2	+================================================================
	3	+Auteur : Simon-Pierre Boucher
	4	+Contact : contact@spboucher.ai
	5	+Projet : Prévision de volatilité réalisée multi-actifs
	6	+ (HAR-RV vs GARCH vs Machine Learning)
	7	+Fichier : realized_vol.py
	8	+Description : Mesures de volatilité réalisée journalières —
	9	+ RV 1-min, RV 5-min sous-échantillonnée, realized
	10	+ kernel (Parzen, BNHLS 2008), bipower variation,
	11	+ test de sauts BNS, semivariances signées et
	12	+ signed jumps. Gestion des sessions par classe.
	13	+================================================================
	14	+"""
	15	+
	16	+from __future__ import annotations
	17	+
	18	+import logging
	19	+from math import gamma as _gamma
	20	+
	21	+import numpy as np
	22	+import pandas as pd
	23	+from scipy import stats
	24	+
	25	+logger = logging.getLogger(__name__)
	26	+
	27	+MU1 = np.sqrt(2.0 / np.pi) # E\|Z\|
	28	+MU43 = 2 ** (2 / 3) * _gamma(7 / 6) / _gamma(1 / 2) # E\|Z\|^{4/3}
	29	+THETA_BNS = (np.pi**2 / 4.0) + np.pi - 5.0 # variance constant of the BNS ratio test
	30	+PARZEN_CSTAR = 3.5134 # optimal bandwidth constant, Parzen kernel
	31	+
	32	+# minimum number of intraday bars for a valid day, per asset class
	33	+MIN_BARS = {"equity": 250, "fx": 700, "crypto": 700, "futures": 700, "index": 250}
	34	+
	35	+# regular trading hours filter (inclusive start, exclusive end), or None for 24h
	36	+SESSION = {
	37	+ "equity": ("09:30", "16:00"),
	38	+ "index": ("09:30", "16:00"),
	39	+ "fx": None,
	40	+ "crypto": None,
	41	+ "futures": None,
	42	+}
	43	+
	44	+
	45	+# ---------------------------------------------------------------- primitives
	46	+def log_returns(prices: np.ndarray) -> np.ndarray:
	47	+ """Log returns of a strictly positive price array.
	48	+
	49	+ Parameters
	50	+ ----------
	51	+ prices : numpy.ndarray
	52	+ Intraday price levels.
	53	+
	54	+ Returns
	55	+ -------
	56	+ numpy.ndarray
	57	+ First differences of log prices (length ``len(prices) - 1``).
	58	+ """
	59	+ return np.diff(np.log(prices))
	60	+
	61	+
	62	+def realized_variance(returns: np.ndarray) -> float:
	63	+ """Plain realized variance: sum of squared intraday returns."""
	64	+ return float(np.sum(returns**2))
	65	+
	66	+
	67	+def rv_subsampled(prices: pd.Series, grid_min: int = 5, offsets: int = 5) -> float:
	68	+ """Subsampled sparse-grid realized variance.
	69	+
	70	+ Averages the realized variance computed on ``offsets`` staggered
	71	+ ``grid_min``-minute grids, following Zhang, Mykland and
	72	+ Aït-Sahalia (2005).
	73	+
	74	+ Parameters
	75	+ ----------
	76	+ prices : pandas.Series
	77	+ 1-minute prices indexed by timestamp (one trading day).
	78	+ grid_min : int
	79	+ Sparse grid spacing in minutes.
	80	+ offsets : int
	81	+ Number of staggered grids averaged.
	82	+
	83	+ Returns
	84	+ -------
	85	+ float
	86	+ Subsampled realized variance.
	87	+ """
	88	+ rvs = []
	89	+ for off in range(offsets):
	90	+ sub = prices.iloc[off::grid_min]
	91	+ if len(sub) > 2:
	92	+ rvs.append(realized_variance(log_returns(sub.to_numpy())))
	93	+ return float(np.mean(rvs)) if rvs else np.nan
	94	+
	95	+
	96	+def bipower_variation(returns: np.ndarray) -> float:
	97	+ """Realized bipower variation (Barndorff-Nielsen & Shephard 2004).
	98	+
	99	+ Robust to jumps; scaled to estimate integrated variance.
	100	+ """
	101	+ n = len(returns)
	102	+ if n < 3:
	103	+ return np.nan
	104	+ absr = np.abs(returns)
	105	+ return float(MU1*-2 (n / (n - 1)) * np.sum(absr[1:] * absr[:-1]))
	106	+
	107	+
	108	+def tripower_quarticity(returns: np.ndarray) -> float:
	109	+ """Tripower quarticity, a jump-robust estimator of integrated quarticity."""
	110	+ n = len(returns)
	111	+ if n < 4:
	112	+ return np.nan
	113	+ a = np.abs(returns) ** (4 / 3)
	114	+ return float(n * MU43*-3 (n / (n - 2)) * np.sum(a[2:] * a[1:-1] * a[:-2]))
	115	+
	116	+
	117	+def bns_jump_test(rv: float, bv: float, tq: float, n: int) -> float:
	118	+ """Ratio-form BNS jump test statistic (Barndorff-Nielsen & Shephard 2006).
	119	+
	120	+ Parameters
	121	+ ----------
	122	+ rv, bv, tq : float
	123	+ Realized variance, bipower variation and tripower quarticity.
	124	+ n : int
	125	+ Number of intraday returns.
	126	+
	127	+ Returns
	128	+ -------
	129	+ float
	130	+ Asymptotically N(0,1) statistic; large positive values indicate a
	131	+ jump day.
	132	+ """
	133	+ if not np.isfinite(rv) or not np.isfinite(bv) or bv <= 0 or n < 4:
	134	+ return np.nan
	135	+ ratio = max(tq / bv**2, 1.0) if np.isfinite(tq) else 1.0
	136	+ denom = np.sqrt(THETA_BNS * ratio / n)
	137	+ return float((1.0 - bv / rv) / denom) if denom > 0 else np.nan
	138	+
	139	+
	140	+def semivariances(returns: np.ndarray) -> tuple[float, float]:
	141	+ """Positive and negative realized semivariance (BNKS 2010).
	142	+
	143	+ Returns
	144	+ -------
	145	+ tuple of float
	146	+ ``(RS+, RS-)`` — sums of squared positive and negative returns.
	147	+ """
	148	+ return (
	149	+ float(np.sum(returns[returns > 0] ** 2)),
	150	+ float(np.sum(returns[returns < 0] ** 2)),
	151	+ )
	152	+
	153	+
	154	+def parzen_kernel(x: np.ndarray) -> np.ndarray:
	155	+ """Parzen kernel weights on [0, 1]."""
	156	+ w = np.zeros_like(x)
	157	+ m1 = x <= 0.5
	158	+ m2 = (x > 0.5) & (x <= 1.0)
	159	+ w[m1] = 1 - 6 * x[m1] ** 2 + 6 * x[m1] ** 3
	160	+ w[m2] = 2 * (1 - x[m2]) ** 3
	161	+ return w
	162	+
	163	+
	164	+def realized_kernel(returns: np.ndarray, iv_proxy: float \| None = None) -> float:
	165	+ """Realized kernel with Parzen weights (Barndorff-Nielsen et al. 2008).
	166	+
	167	+ The bandwidth follows the feasible rule
	168	+ :math:`H^* = c^* \\xi^{4/5} n^{3/5}` with
	169	+ :math:`\\xi^2 = \\hat\\omega^2 / \\sqrt{\\widehat{IQ}}` approximated by
	170	+ :math:`\\hat\\omega^2 / IV` where the noise variance is estimated as
	171	+ ``RV_dense / (2 n)`` and IV by a sparse-grid RV.
	172	+
	173	+ Parameters
	174	+ ----------
	175	+ returns : numpy.ndarray
	176	+ Dense (1-minute) intraday log returns.
	177	+ iv_proxy : float, optional
	178	+ Noise-robust estimate of integrated variance; defaults to the
	179	+ realized variance of the input returns when omitted.
	180	+
	181	+ Returns
	182	+ -------
	183	+ float
	184	+ Realized kernel estimate of integrated variance (non-negative by
	185	+ construction with non-flat-top Parzen weights).
	186	+ """
	187	+ n = len(returns)
	188	+ if n < 10:
	189	+ return np.nan
	190	+ rv_dense = realized_variance(returns)
	191	+ iv = iv_proxy if iv_proxy and np.isfinite(iv_proxy) and iv_proxy > 0 else rv_dense
	192	+ omega2 = rv_dense / (2.0 * n)
	193	+ xi2 = omega2 / iv if iv > 0 else 0.0
	194	+ H = int(np.clip(PARZEN_CSTAR * xi2*0.4 n**0.6, 1, n - 1))
	195	+ gamma0 = float(returns @ returns)
	196	+ k = gamma0
	197	+ weights = parzen_kernel(np.arange(1, H + 1) / (H + 1))
	198	+ for h in range(1, H + 1):
	199	+ gam = float(returns[h:] @ returns[:-h])
	200	+ k += 2.0 * weights[h - 1] * gam
	201	+ return float(max(k, 0.0))
	202	+
	203	+
	204	+# ---------------------------------------------------------------- daily loop
	205	+def _day_measures(day: pd.DataFrame, alpha: float) -> dict[str, float]:
	206	+ """Compute all volatility measures for one trading day.
	207	+
	208	+ Parameters
	209	+ ----------
	210	+ day : pandas.DataFrame
	211	+ 1-minute bars of a single day (columns ``datetime``, ``close``).
	212	+ alpha : float
	213	+ Size of the BNS jump test.
	214	+
	215	+ Returns
	216	+ -------
	217	+ dict
	218	+ All daily realized measures (see module docstring).
	219	+ """
	220	+ prices = day.set_index("datetime")["close"]
	221	+ r1 = log_returns(prices.to_numpy())
	222	+ # sparse 5-min grid (offset 0) for jump-robust quantities
	223	+ p5 = prices.iloc[::5]
	224	+ r5 = log_returns(p5.to_numpy())
	225	+
	226	+ rv1 = realized_variance(r1)
	227	+ rv5 = rv_subsampled(prices, grid_min=5, offsets=5)
	228	+ rk = realized_kernel(r1, iv_proxy=rv5)
	229	+ rv5_plain = realized_variance(r5)
	230	+ bv = bipower_variation(r5)
	231	+ tq = tripower_quarticity(r5)
	232	+ z = bns_jump_test(rv5_plain, bv, tq, len(r5))
	233	+ crit = stats.norm.ppf(1 - alpha)
	234	+ is_jump = bool(np.isfinite(z) and z > crit)
	235	+ jump = max(rv5_plain - bv, 0.0) if is_jump else 0.0
	236	+ cont = rv5_plain - jump
	237	+ rsp, rsn = semivariances(r5)
	238	+ # realized quarticity on the 5-min grid (for HARQ)
	239	+ rq = float(len(r5) / 3.0 * np.sum(r5**4)) if len(r5) > 1 else np.nan
	240	+
	241	+ return {
	242	+ "rv1": rv1,
	243	+ "rv5ss": rv5,
	244	+ "rk": rk,
	245	+ "rv5": rv5_plain,
	246	+ "bv": bv,
	247	+ "tq": tq,
	248	+ "rq": rq,
	249	+ "z_bns": z,
	250	+ "jump": jump,
	251	+ "cont": cont,
	252	+ "rsp": rsp,
	253	+ "rsn": rsn,
	254	+ "sj": rsp - rsn,
	255	+ "n_bars": float(len(prices)),
	256	+ "ret_intraday": float(np.sum(r1)),
	257	+ "close": float(prices.iloc[-1]),
	258	+ "volume": float(day["volume"].sum()),
	259	+ }
	260	+
	261	+
	262	+def build_daily_rv(
	263	+ bars: pd.DataFrame,
	264	+ cls: str,
	265	+ jump_alpha: float = 0.001,
	266	+ min_bars: int \| None = None,
	267	+) -> pd.DataFrame:
	268	+ """Build the daily realized-measure panel for one instrument.
	269	+
	270	+ Parameters
	271	+ ----------
	272	+ bars : pandas.DataFrame
	273	+ 1-minute bars (columns ``datetime``, ``open``, ``high``, ``low``,
	274	+ ``close``, ``volume``) in exchange-local time.
	275	+ cls : str
	276	+ Asset class (``equity``, ``fx``, ``crypto``, ``futures``,
	277	+ ``index``); controls session filtering and the valid-day
	278	+ threshold.
	279	+ jump_alpha : float
	280	+ Size of the BNS jump test used to shrink the jump component.
	281	+ min_bars : int, optional
	282	+ Override the per-class minimum number of intraday bars.
	283	+
	284	+ Returns
	285	+ -------
	286	+ pandas.DataFrame
	287	+ One row per valid trading day, indexed by date, with all
	288	+ realized measures plus the close-to-close daily log return
	289	+ (``ret_cc``, includes overnight for session-limited markets).
	290	+ """
	291	+ if bars.empty:
	292	+ return pd.DataFrame()
	293	+ df = bars.copy()
	294	+ session = SESSION[cls]
	295	+ if session is not None:
	296	+ t = df["datetime"].dt.time
	297	+ lo = pd.Timestamp(f"2000-01-01 {session[0]}").time()
	298	+ hi = pd.Timestamp(f"2000-01-01 {session[1]}").time()
	299	+ df = df[(t >= lo) & (t < hi)]
	300	+ df = df[df["close"] > 0]
	301	+ threshold = min_bars if min_bars is not None else MIN_BARS[cls]
	302	+
	303	+ rows: list[dict[str, float]] = []
	304	+ dates: list[pd.Timestamp] = []
	305	+ for date, day in df.groupby(df["datetime"].dt.normalize()):
	306	+ if len(day) < threshold:
	307	+ continue
	308	+ rows.append(_day_measures(day, jump_alpha))
	309	+ dates.append(date)
	310	+ if not rows:
	311	+ logger.warning("no valid days for class=%s", cls)
	312	+ return pd.DataFrame()
	313	+ out = pd.DataFrame(rows, index=pd.DatetimeIndex(dates, name="date"))
	314	+ out["ret_cc"] = np.log(out["close"]).diff()
	315	+ logger.info("daily RV built: %d valid days (%s -> %s)",
	316	+ len(out), out.index[0].date(), out.index[-1].date())
	317	+ return out

added src/wp12/robustness.py +155 −0

		`@@ -0,0 +1,155 @@`
	1	+"""
	2	+================================================================
	3	+Auteur : Simon-Pierre Boucher
	4	+Contact : contact@spboucher.ai
	5	+Projet : Prévision de volatilité réalisée multi-actifs
	6	+ (HAR-RV vs GARCH vs Machine Learning)
	7	+Fichier : robustness.py
	8	+Description : Analyses de robustesse — sous-périodes, découpage
	9	+ par classe/actif, sensibilité à la fréquence RV et
	10	+ à la taille de fenêtre, test de fluctuation de
	11	+ Giacomini-Rossi (2010).
	12	+================================================================
	13	+"""
	14	+
	15	+from __future__ import annotations
	16	+
	17	+import logging
	18	+
	19	+import numpy as np
	20	+import pandas as pd
	21	+
	22	+from . import evaluation
	23	+
	24	+logger = logging.getLogger(__name__)
	25	+
	26	+# Two-sided 5% critical values of the Giacomini-Rossi (2010) fluctuation
	27	+# test, indexed by the window fraction mu = m / n (their Table 1).
	28	+GR_CRIT_5PCT = {
	29	+ 0.1: 3.393, 0.2: 3.179, 0.3: 3.012, 0.4: 2.890, 0.5: 2.779,
	30	+ 0.6: 2.634, 0.7: 2.560, 0.8: 2.433, 0.9: 2.278, 1.0: 2.062,
	31	+}
	32	+
	33	+
	34	+def loss_by_period(
	35	+ merged: pd.DataFrame,
	36	+ subperiods: dict[str, tuple[str, str]],
	37	+ loss: str = "qlike",
	38	+) -> pd.DataFrame:
	39	+ """Average loss per model within each sub-period.
	40	+
	41	+ Parameters
	42	+ ----------
	43	+ merged : pandas.DataFrame
	44	+ Long forecast table with columns ``date, ticker, h, model,
	45	+ forecast, actual``.
	46	+ subperiods : dict
	47	+ ``name -> (start, end)`` date bounds (inclusive).
	48	+ loss : str
	49	+ ``qlike`` or ``mse``.
	50	+
	51	+ Returns
	52	+ -------
	53	+ pandas.DataFrame
	54	+ Mean loss per (period, h, model).
	55	+ """
	56	+ lf = evaluation.LOSSES[loss]
	57	+ rows = []
	58	+ for name, (lo, hi) in subperiods.items():
	59	+ sub = merged[(merged["date"] >= lo) & (merged["date"] <= hi)]
	60	+ if sub.empty:
	61	+ continue
	62	+ g = sub.groupby(["h", "model"], observed=True).apply(
	63	+ lambda d: float(np.mean(lf(d["actual"].to_numpy(), d["forecast"].to_numpy()))),
	64	+ include_groups=False,
	65	+ )
	66	+ for (h, model), v in g.items():
	67	+ rows.append((name, h, model, v))
	68	+ return pd.DataFrame(rows, columns=["period", "h", "model", loss])
	69	+
	70	+
	71	+def loss_by_group(
	72	+ merged: pd.DataFrame, group_col: str, loss: str = "qlike"
	73	+) -> pd.DataFrame:
	74	+ """Average loss per model within each group (asset class or ticker)."""
	75	+ lf = evaluation.LOSSES[loss]
	76	+ g = merged.groupby([group_col, "h", "model"], observed=True).apply(
	77	+ lambda d: float(np.mean(lf(d["actual"].to_numpy(), d["forecast"].to_numpy()))),
	78	+ include_groups=False,
	79	+ )
	80	+ return g.rename(loss).reset_index()
	81	+
	82	+
	83	+def giacomini_rossi_fluctuation(
	84	+ loss_a: pd.Series,
	85	+ loss_b: pd.Series,
	86	+ h: int = 1,
	87	+ mu: float = 0.3,
	88	+) -> pd.DataFrame:
	89	+ """Giacomini-Rossi (2010) fluctuation test of local relative performance.
	90	+
	91	+ Computes the standardized rolling-mean loss differential of A minus B
	92	+ over a centred window of fraction ``mu`` of the sample, together with
	93	+ the two-sided 5% critical value.
	94	+
	95	+ Parameters
	96	+ ----------
	97	+ loss_a, loss_b : pandas.Series
	98	+ Per-date losses of the two models (aligned on dates).
	99	+ h : int
	100	+ Forecast horizon (HAC bandwidth ``h - 1``).
	101	+ mu : float
	102	+ Rolling window as a fraction of the evaluation sample.
	103	+
	104	+ Returns
	105	+ -------
	106	+ pandas.DataFrame
	107	+ ``date, stat, crit`` — the fluctuation statistic path; values
	108	+ below ``-crit`` mean A significantly better locally.
	109	+ """
	110	+ d = (loss_a - loss_b).dropna()
	111	+ n = len(d)
	112	+ m = max(int(round(mu * n)), 30)
	113	+ if n < m + 10:
	114	+ return pd.DataFrame(columns=["date", "stat", "crit"])
	115	+ # HAC variance estimated on the full sample (GR recommendation)
	116	+ sig2 = evaluation.newey_west_var(d.to_numpy(), max(h - 1, 0))
	117	+ sig = np.sqrt(max(sig2, 1e-30))
	118	+ roll = d.rolling(m).mean()
	119	+ stat = np.sqrt(m) * roll / sig
	120	+ mu_grid = np.array(sorted(GR_CRIT_5PCT))
	121	+ crit = float(np.interp(m / n, mu_grid, [GR_CRIT_5PCT[k] for k in mu_grid]))
	122	+ out = pd.DataFrame({"date": d.index, "stat": stat.to_numpy()})
	123	+ out["crit"] = crit
	124	+ return out.dropna().reset_index(drop=True)
	125	+
	126	+
	127	+def sensitivity_table(
	128	+ results: dict[str, pd.DataFrame], key_name: str, loss: str = "qlike"
	129	+) -> pd.DataFrame:
	130	+ """Stack per-variant loss tables into one comparison frame.
	131	+
	132	+ Parameters
	133	+ ----------
	134	+ results : dict
	135	+ ``variant -> merged forecast table`` (with ``actual`` column).
	136	+ key_name : str
	137	+ Name of the variant dimension (e.g. ``frequency``, ``window``).
	138	+ loss : str
	139	+ Loss function name.
	140	+
	141	+ Returns
	142	+ -------
	143	+ pandas.DataFrame
	144	+ Mean loss per (variant, h, model).
	145	+ """
	146	+ lf = evaluation.LOSSES[loss]
	147	+ rows = []
	148	+ for variant, merged in results.items():
	149	+ g = merged.groupby(["h", "model"], observed=True).apply(
	150	+ lambda d: float(np.mean(lf(d["actual"].to_numpy(), d["forecast"].to_numpy()))),
	151	+ include_groups=False,
	152	+ )
	153	+ for (h, model), v in g.items():
	154	+ rows.append((variant, h, model, v))
	155	+ return pd.DataFrame(rows, columns=[key_name, "h", "model", loss])

added tests/test_evaluation.py +125 −0

		`@@ -0,0 +1,125 @@`
	1	+"""
	2	+================================================================
	3	+Auteur : Simon-Pierre Boucher
	4	+Contact : contact@spboucher.ai
	5	+Projet : Prévision de volatilité réalisée multi-actifs
	6	+ (HAR-RV vs GARCH vs Machine Learning)
	7	+Fichier : test_evaluation.py
	8	+Description : Tests pytest de l'évaluation — propriétés des
	9	+ pertes, taille/puissance du test DM, MCS, Mincer-
	10	+ Zarnowitz, encompassing et combinaisons.
	11	+================================================================
	12	+"""
	13	+
	14	+from __future__ import annotations
	15	+
	16	+import sys
	17	+from pathlib import Path
	18	+
	19	+import numpy as np
	20	+import pandas as pd
	21	+import pytest
	22	+
	23	+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
	24	+
	25	+from wp12 import evaluation as ev # noqa: E402
	26	+
	27	+RNG = np.random.default_rng(20260810)
	28	+
	29	+
	30	+class TestLosses:
	31	+ def test_qlike_zero_at_perfect_forecast(self) -> None:
	32	+ a = RNG.uniform(0.5, 2.0, 100)
	33	+ assert np.allclose(ev.qlike(a, a), 0.0)
	34	+
	35	+ def test_qlike_positive_and_asymmetric(self) -> None:
	36	+ a = np.ones(10)
	37	+ under = ev.qlike(a, 0.5 * a).mean()
	38	+ over = ev.qlike(a, 2.0 * a).mean()
	39	+ assert under > 0 and over > 0
	40	+ assert under > over # QLIKE punishes under-prediction more
	41	+
	42	+ def test_mse_basic(self) -> None:
	43	+ assert ev.mse(np.array([2.0]), np.array([1.0]))[0] == 1.0
	44	+
	45	+
	46	+class TestDieboldMariano:
	47	+ def test_dm_detects_better_model(self) -> None:
	48	+ n = 800
	49	+ actual = np.exp(RNG.normal(0, 0.5, n))
	50	+ good = actual * np.exp(RNG.normal(0, 0.1, n))
	51	+ bad = actual * np.exp(RNG.normal(0.5, 0.5, n))
	52	+ la, lb = ev.qlike(actual, good), ev.qlike(actual, bad)
	53	+ stat, p = ev.diebold_mariano(la, lb, h=1)
	54	+ assert stat < -2 and p < 0.05
	55	+
	56	+ def test_dm_size_under_null(self) -> None:
	57	+ """Two equally noisy forecasts: rejection rate ~ nominal."""
	58	+ rejects = 0
	59	+ for _ in range(200):
	60	+ actual = np.exp(RNG.normal(0, 0.5, 300))
	61	+ f1 = actual * np.exp(RNG.normal(0, 0.3, 300))
	62	+ f2 = actual * np.exp(RNG.normal(0, 0.3, 300))
	63	+ _, p = ev.diebold_mariano(ev.qlike(actual, f1), ev.qlike(actual, f2))
	64	+ rejects += p < 0.05
	65	+ assert rejects / 200 < 0.12
	66	+
	67	+
	68	+class TestMCS:
	69	+ def test_mcs_keeps_good_removes_bad(self) -> None:
	70	+ n = 600
	71	+ actual = np.exp(RNG.normal(0, 0.5, n))
	72	+ losses = pd.DataFrame({
	73	+ "good1": ev.qlike(actual, actual * np.exp(RNG.normal(0, 0.1, n))),
	74	+ "good2": ev.qlike(actual, actual * np.exp(RNG.normal(0, 0.1, n))),
	75	+ "bad": ev.qlike(actual, actual * np.exp(RNG.normal(1.0, 0.5, n))),
	76	+ })
	77	+ out = ev.model_confidence_set(losses, alpha=0.10, n_boot=800)
	78	+ assert not out.loc[out["model"] == "bad", "in_mcs"].iloc[0]
	79	+ assert out["in_mcs"].sum() >= 1
	80	+ good_kept = out.loc[out["model"].str.startswith("good"), "in_mcs"]
	81	+ assert good_kept.any()
	82	+
	83	+
	84	+class TestMZAndEncompassing:
	85	+ def test_mz_unbiased_forecast(self) -> None:
	86	+ n = 1000
	87	+ f = np.exp(RNG.normal(0, 0.5, n))
	88	+ actual = f * np.exp(RNG.normal(-0.005, 0.1, n))
	89	+ res = ev.mincer_zarnowitz(actual, f)
	90	+ assert res["beta"] == pytest.approx(1.0, abs=0.1)
	91	+ assert res["r2"] > 0.8
	92	+
	93	+ def test_mz_biased_forecast_rejected(self) -> None:
	94	+ n = 1000
	95	+ f = np.exp(RNG.normal(0, 0.5, n))
	96	+ actual = 2.0 * f * np.exp(RNG.normal(0, 0.1, n))
	97	+ res = ev.mincer_zarnowitz(actual, f)
	98	+ assert res["p_joint"] < 0.01
	99	+
	100	+ def test_encompassing_detects_added_info(self) -> None:
	101	+ n = 1000
	102	+ true = np.exp(RNG.normal(0, 0.5, n))
	103	+ actual = true * np.exp(RNG.normal(0, 0.2, n))
	104	+ f_base = true * np.exp(RNG.normal(0, 0.4, n)) # noisy version
	105	+ f_rival = 0.5 * f_base + 0.5 * true # contains extra info
	106	+ res = ev.forecast_encompassing(actual, f_base, f_rival)
	107	+ assert res["p_b2"] < 0.05
	108	+ assert res["b2"] > 0
	109	+
	110	+
	111	+class TestCombination:
	112	+ def test_combined_beats_worst(self) -> None:
	113	+ n = 700
	114	+ idx = pd.date_range("2020-01-01", periods=n, freq="B")
	115	+ actual = pd.Series(np.exp(RNG.normal(0, 0.5, n)), index=idx)
	116	+ f = pd.DataFrame({
	117	+ "m1": actual * np.exp(RNG.normal(0, 0.2, n)),
	118	+ "m2": actual * np.exp(RNG.normal(0.3, 0.4, n)),
	119	+ }, index=idx)
	120	+ combs = ev.combine_forecasts(f, actual, window=60)
	121	+ q_mean = ev.qlike(actual.to_numpy(), combs["Comb-Mean"].to_numpy()).mean()
	122	+ q_inv = np.nanmean(ev.qlike(actual.to_numpy(), combs["Comb-InvMSE"].to_numpy()))
	123	+ q_worst = ev.qlike(actual.to_numpy(), f["m2"].to_numpy()).mean()
	124	+ assert q_mean < q_worst
	125	+ assert q_inv < q_worst

added tests/test_realized_vol.py +150 −0

		`@@ -0,0 +1,150 @@`
	1	+"""
	2	+================================================================
	3	+Auteur : Simon-Pierre Boucher
	4	+Contact : contact@spboucher.ai
	5	+Projet : Prévision de volatilité réalisée multi-actifs
	6	+ (HAR-RV vs GARCH vs Machine Learning)
	7	+Fichier : test_realized_vol.py
	8	+Description : Tests pytest des mesures de volatilité réalisée —
	9	+ cohérence RV/BV/kernel sur des diffusions simulées,
	10	+ détection de sauts, semivariances, sessions.
	11	+================================================================
	12	+"""
	13	+
	14	+from __future__ import annotations
	15	+
	16	+import sys
	17	+from pathlib import Path
	18	+
	19	+import numpy as np
	20	+import pandas as pd
	21	+import pytest
	22	+
	23	+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
	24	+
	25	+from wp12 import realized_vol as rvol # noqa: E402
	26	+
	27	+RNG = np.random.default_rng(20260810)
	28	+
	29	+
	30	+def _simulate_day(n: int = 390, sigma_day: float = 0.01, jump: float = 0.0,
	31	+ noise: float = 0.0) -> pd.DataFrame:
	32	+ """Simulate one day of 1-min prices from a Brownian motion (+ optional
	33	+ jump midday and iid microstructure noise)."""
	34	+ r = RNG.normal(0.0, sigma_day / np.sqrt(n), n)
	35	+ if jump:
	36	+ r[n // 2] += jump
	37	+ logp = np.cumsum(r) + np.log(100.0)
	38	+ if noise:
	39	+ logp = logp + RNG.normal(0.0, noise, n)
	40	+ ts = pd.date_range("2024-01-02 09:30", periods=n, freq="1min")
	41	+ return pd.DataFrame({"datetime": ts, "close": np.exp(logp),
	42	+ "open": np.exp(logp), "high": np.exp(logp),
	43	+ "low": np.exp(logp), "volume": 100.0})
	44	+
	45	+
	46	+class TestPrimitives:
	47	+ def test_rv_matches_true_variance(self) -> None:
	48	+ """RV averaged over many simulated days converges to sigma^2."""
	49	+ sigma = 0.012
	50	+ rvs = []
	51	+ for _ in range(300):
	52	+ day = _simulate_day(sigma_day=sigma)
	53	+ rvs.append(rvol.realized_variance(
	54	+ rvol.log_returns(day["close"].to_numpy())))
	55	+ assert np.mean(rvs) == pytest.approx(sigma**2, rel=0.05)
	56	+
	57	+ def test_bipower_robust_to_jump(self) -> None:
	58	+ """BV is nearly unaffected by a jump while RV jumps by ~jump^2."""
	59	+ sigma, jump = 0.01, 0.02
	60	+ rv_j, bv_j = [], []
	61	+ for _ in range(300):
	62	+ day = _simulate_day(sigma_day=sigma, jump=jump)
	63	+ r = rvol.log_returns(day["close"].to_numpy())
	64	+ rv_j.append(rvol.realized_variance(r))
	65	+ bv_j.append(rvol.bipower_variation(r))
	66	+ assert np.mean(rv_j) == pytest.approx(sigma2 + jump2, rel=0.1)
	67	+ # BV has a small finite-sample jump bias (adjacent cross-products,
	68	+ # ~ 2mu1^-2jump*E\|r\|), but must stay FAR below the RV inflation.
	69	+ assert np.mean(bv_j) < sigma*2 + 0.5 jump**2 / 4
	70	+ assert (np.mean(rv_j) - sigma2) / (np.mean(bv_j) - sigma2) > 5
	71	+
	72	+ def test_semivariances_sum_to_rv(self) -> None:
	73	+ day = _simulate_day()
	74	+ r = rvol.log_returns(day["close"].to_numpy())
	75	+ rsp, rsn = rvol.semivariances(r)
	76	+ assert rsp + rsn == pytest.approx(rvol.realized_variance(r), rel=1e-10)
	77	+ assert rsp >= 0 and rsn >= 0
	78	+
	79	+ def test_kernel_removes_noise_bias(self) -> None:
	80	+ """With iid noise, RV explodes with n while the kernel stays close
	81	+ to integrated variance."""
	82	+ sigma, noise = 0.01, 3e-4
	83	+ rv_list, rk_list = [], []
	84	+ for _ in range(200):
	85	+ day = _simulate_day(n=390, sigma_day=sigma, noise=noise)
	86	+ prices = day.set_index("datetime")["close"]
	87	+ r = rvol.log_returns(prices.to_numpy())
	88	+ rv_list.append(rvol.realized_variance(r))
	89	+ rk_list.append(rvol.realized_kernel(
	90	+ r, iv_proxy=rvol.rv_subsampled(prices)))
	91	+ bias_rv = np.mean(rv_list) / sigma**2
	92	+ bias_rk = np.mean(rk_list) / sigma**2
	93	+ assert bias_rv > 1.5 # noise inflates plain RV
	94	+ assert abs(bias_rk - 1) < abs(bias_rv - 1) # kernel shrinks the bias
	95	+ assert bias_rk == pytest.approx(1.0, abs=0.35)
	96	+
	97	+ def test_kernel_non_negative(self) -> None:
	98	+ for _ in range(50):
	99	+ day = _simulate_day(n=60)
	100	+ r = rvol.log_returns(day["close"].to_numpy())
	101	+ assert rvol.realized_kernel(r) >= 0.0
	102	+
	103	+ def test_bns_flags_jump_days(self) -> None:
	104	+ """The BNS statistic is larger on jump days than diffusion days."""
	105	+ z_nojump, z_jump = [], []
	106	+ for _ in range(200):
	107	+ d0 = _simulate_day(sigma_day=0.01)
	108	+ d1 = _simulate_day(sigma_day=0.01, jump=0.02)
	109	+ for day, store in ((d0, z_nojump), (d1, z_jump)):
	110	+ r = rvol.log_returns(day["close"].to_numpy()[::5])
	111	+ z = rvol.bns_jump_test(
	112	+ rvol.realized_variance(r), rvol.bipower_variation(r),
	113	+ rvol.tripower_quarticity(r), len(r))
	114	+ store.append(z)
	115	+ # size: ~1% of diffusion days flagged at the 0.1% level
	116	+ crit = 3.09
	117	+ assert np.mean(np.array(z_nojump) > crit) < 0.05
	118	+ # power: most jump days flagged
	119	+ assert np.mean(np.array(z_jump) > crit) > 0.60
	120	+
	121	+
	122	+class TestDailyPanel:
	123	+ def test_build_daily_rv_sessions(self) -> None:
	124	+ """Equity session filtering keeps only 09:30-16:00 bars."""
	125	+ days = []
	126	+ for d in pd.date_range("2024-01-02", periods=8, freq="B"):
	127	+ ts = pd.date_range(d + pd.Timedelta(hours=4), periods=960, freq="1min")
	128	+ logp = np.cumsum(RNG.normal(0, 3e-4, 960)) + np.log(100)
	129	+ days.append(pd.DataFrame({
	130	+ "datetime": ts, "close": np.exp(logp), "open": np.exp(logp),
	131	+ "high": np.exp(logp), "low": np.exp(logp), "volume": 10.0}))
	132	+ bars = pd.concat(days, ignore_index=True)
	133	+ out = rvol.build_daily_rv(bars, cls="equity")
	134	+ assert len(out) == 8
	135	+ assert (out["n_bars"] == 390).all()
	136	+ # columns present
	137	+ for c in ("rv1", "rv5ss", "rk", "bv", "jump", "cont", "rsp", "rsn", "sj", "rq"):
	138	+ assert c in out.columns
	139	+ # decomposition adds up
	140	+ assert np.allclose(out["cont"] + out["jump"], out["rv5"], rtol=1e-10)
	141	+
	142	+ def test_min_bars_filter(self) -> None:
	143	+ """Days with too few bars are dropped."""
	144	+ ts = pd.date_range("2024-01-02 09:30", periods=100, freq="1min")
	145	+ logp = np.cumsum(RNG.normal(0, 3e-4, 100)) + np.log(100)
	146	+ bars = pd.DataFrame({"datetime": ts, "close": np.exp(logp),
	147	+ "open": np.exp(logp), "high": np.exp(logp),
	148	+ "low": np.exp(logp), "volume": 1.0})
	149	+ out = rvol.build_daily_rv(bars, cls="equity")
	150	+ assert out.empty
	151