Repository Public Research preview
1"""Rolling-origin forecasting driver (one job = one instrument x one model x one innovation distribution).23Protocol (configs/forecast_protocol.yaml, frozen by docs/preregistered_protocol.md):4 * estimation window: rolling, W observations ending at the origin t (expanding = robustness variant);5 * re-estimation every `refit_every` trading days; between refits the conditional variance is FILTERED daily with the6 last estimated parameters over the current window (no parameter look-ahead, no future data);7 * at every origin t in the out-of-sample range, per-step variance forecasts for k = 1..22 are produced from8 information up to and including day t (returns r_1..r_t and, for Tier B, realized measures x_1..x_t observed at the9 end of day t); cumulative forecasts for h = 1, 5, 22 are sums of per-step forecasts;10 * failure policy: a refit that raises or does not converge after the bounded retry sequence is logged as a failure;11 the job then keeps the LAST successful parameter vector (flag stale_fit = True) if one exists, otherwise the forecast12 is NaN. Inadmissible optima are kept and flagged (admissible = False). Nothing is ever replaced by another model.13Outputs: forecasts table (one row per origin) and fits table (one row per refit attempt), both as pandas DataFrames.14"""15from __future__ import annotations1617import time18import traceback19from dataclasses import dataclass2021import numpy as np22import pandas as pd2324from wp46.models.arch_models import SPECS as ARCH_SPECS, ArchModel25from wp46.models.baselines import HARRV, HistVar, RMRandomWalk26from wp46.models.custom import MODELS as CUSTOM_MODELS27from wp46.models.extensions import EXTENSIONS28from wp46.scoring.proxies import proxy_series2930HMAX = 22313233@dataclass34class JobSpec:35 symbol: str36 model_id: str37 dist: str38 tier: str # 'A' daily-only, 'B' realized-measure39 window: int = 100040 refit_every: int = 541 expanding: bool = False42 rm_proxy: str = "rv5ss_on" # realized measure fed to Tier B models43 oos_start: str = "2012-01-01"44 oos_end: str = "2026-10-30"45 seed: int = 4646 max_restarts: int = 34748 @property49 def spec_id(self) -> str:50 w = "exp" if self.expanding else f"w{self.window}"51 return f"{self.model_id}__{self.dist}__{w}_r{self.refit_every}"525354def make_model(model_id: str):55 if model_id in ARCH_SPECS:56 return ArchModel(model_id)57 if model_id in CUSTOM_MODELS:58 return CUSTOM_MODELS[model_id]()59 if model_id in EXTENSIONS:60 return EXTENSIONS[model_id]()61 if model_id.startswith("histvar"):62 return HistVar(int(model_id[7:]))63 if model_id == "har_rv":64 return HARRV()65 if model_id == "loghar_rv":66 return HARRV(log=True)67 if model_id == "rm_rw":68 return RMRandomWalk()69 raise KeyError(model_id)707172def _empirical_tail(model, fit, r: np.ndarray, x) -> dict:73 """Empirical 1 % / 5 % quantiles and tail means of the standardized in-window residuals z = (r - mu)/sqrt(sigma2)74 (used for the VaR/ES of distribution-free models and as a diagnostic for parametric ones)."""75 try:76 s2 = np.asarray(model.filter(fit, r, x), float)77 mu = float(fit.params[0]) if len(fit.params) else float(np.mean(r))78 z = (np.asarray(r, float) - mu) / np.sqrt(np.maximum(s2, 1e-12))79 z = z[np.isfinite(z)]80 out = {}81 for p in (0.01, 0.05):82 q = float(np.quantile(z, p))83 out[f"emp_q{int(p*100):02d}"] = q84 out[f"emp_es{int(p*100):02d}"] = float(np.mean(z[z <= q])) if np.any(z <= q) else q85 out["z_kurtosis"] = float(np.mean(z ** 4) / np.mean(z ** 2) ** 2)86 return out87 except Exception: # noqa: BLE00188 return {}899091def run_job(panel: pd.DataFrame, job: JobSpec, progress: bool = False) -> tuple[pd.DataFrame, pd.DataFrame]:92 model = make_model(job.model_id)93 needs_rm = bool(getattr(model, "needs_rm", False))94 r_all = panel["r"].values.astype(float)95 dates = pd.to_datetime(panel["session"]).values96 x_all = proxy_series(panel, job.rm_proxy) if needs_rm else None97 if needs_rm:98 # Tier B input: missing realized measures (invalid sessions) are replaced by the squared return of that day99 # (documented imputation, flagged in the panel) so that the recursion stays defined100 x_all = np.where(np.isfinite(x_all), x_all, r_all ** 2)101 x_all = np.maximum(x_all, 1e-8)102 T = len(r_all)103 ok = np.isfinite(r_all)104 oos_mask = (dates >= np.datetime64(job.oos_start)) & (dates <= np.datetime64(job.oos_end))105 origins = np.where(oos_mask)[0]106 origins = origins[origins >= job.window]107 fits_rows, fc_rows = [], []108 fit = None109 fit_id = -1110 n_since = job.refit_every111 stale = False112 rng_seed = job.seed113 for i, t in enumerate(origins):114 lo = 0 if job.expanding else t + 1 - job.window115 r = r_all[lo:t + 1]116 x = x_all[lo:t + 1] if needs_rm else None117 if not np.all(ok[lo:t + 1]): # missing returns inside the window: skip origin (recorded)118 fc_rows.append(dict(origin=dates[t], fit_id=fit_id, status="missing_input"))119 continue120 if n_since >= job.refit_every or fit is None:121 fit_id += 1122 t0 = time.perf_counter()123 try:124 if job.model_id == "garch_midas":125 months = (dates[lo:t + 1].astype("datetime64[M]")).astype(int)126 new_fit = model.fit(r, x, dist=job.dist, max_restarts=job.max_restarts, months=months)127 else:128 new_fit = model.fit(r, x, dist=job.dist, max_restarts=job.max_restarts)129 err = ""130 except Exception as e: # noqa: BLE001131 new_fit, err = None, f"{type(e).__name__}: {str(e)[:120]}"132 dt = time.perf_counter() - t0133 emp = {}134 if new_fit is not None and new_fit.converged:135 fit, stale = new_fit, False136 status = "ok" if new_fit.admissible else "inadmissible"137 emp = _empirical_tail(model, new_fit, r, x)138 else:139 status = "failed"140 stale = fit is not None141 fits_rows.append(dict(fit_id=fit_id, origin=dates[t], status=status, error=err or (new_fit.notes if new_fit else ""), **emp,142 loglik=new_fit.loglik if new_fit else np.nan, n_iter=new_fit.n_iter if new_fit else 0,143 n_restarts=new_fit.n_restarts if new_fit else 0, time_s=dt, n_obs=len(r),144 admissible=bool(new_fit.admissible) if new_fit else False,145 params=list(map(float, new_fit.params)) if new_fit else None,146 param_names=list(new_fit.param_names) if new_fit else None))147 n_since = 0148 n_since += 1149 if fit is None:150 fc_rows.append(dict(origin=dates[t], fit_id=fit_id, status="no_fit"))151 continue152 t0 = time.perf_counter()153 try:154 fc = model.forecast(fit, r, x, HMAX, seed=rng_seed + t)155 sv = np.asarray(fc.step_var, float)156 if not np.all(np.isfinite(sv)) or np.any(sv <= 0):157 raise FloatingPointError("non-positive or non-finite variance forecast")158 row = dict(origin=dates[t], fit_id=fit_id, status="stale_fit" if stale else "ok", method=fc.method,159 h1=float(sv[0]), h5=float(np.sum(sv[:5])), h22=float(np.sum(sv[:22])), sigma2_t=fc.sigma2_t,160 mu=fc.mean, time_fc_s=time.perf_counter() - t0, admissible=bool(fit.admissible))161 dp = np.asarray(fc.dist_params, float)162 row["dp1"] = float(dp[0]) if len(dp) > 0 else np.nan163 row["dp2"] = float(dp[1]) if len(dp) > 1 else np.nan164 except Exception as e: # noqa: BLE001165 row = dict(origin=dates[t], fit_id=fit_id, status="forecast_failed", error=f"{type(e).__name__}: {str(e)[:120]}",166 time_fc_s=time.perf_counter() - t0)167 fc_rows.append(row)168 if progress and i % 250 == 0:169 print(f" {job.symbol} {job.spec_id} {i}/{len(origins)} {str(dates[t])[:10]}", flush=True)170 F = pd.DataFrame(fc_rows)171 G = pd.DataFrame(fits_rows)172 for df in (F, G):173 df["symbol"] = job.symbol174 df["model_id"] = job.model_id175 df["dist"] = job.dist176 df["spec_id"] = job.spec_id177 df["tier"] = job.tier178 df["window"] = job.window179 df["refit_every"] = job.refit_every180 df["expanding"] = job.expanding181 return F, G182