Repository Public Research preview
1"""Parameter-free / least-squares comparators: rolling historical variance, HAR-RV (Corsi 2009) direct projections,2log-HAR with smearing correction, and the realized-measure random walk. These have no innovation distribution of their3own; for VaR they are paired with the empirical distribution of standardized returns in the estimation window (documented)."""4from __future__ import annotations56import time78import numpy as np910from .base import FitResult, Forecast111213class HistVar:14 """Rolling sample variance of the last `window` daily returns (percent^2), flat multi-step forecast."""15 family = "baseline"1617 def __init__(self, window: int):18 self.window = window19 self.model_id = f"histvar{window}"2021 def fit(self, r, x=None, dist="empirical", **kw):22 t0 = time.perf_counter()23 return FitResult(self.model_id, "empirical", ["mu"], np.array([float(np.mean(r))]), float("nan"), True, len(r),24 time_s=time.perf_counter() - t0, aux={})2526 def filter(self, fit, r, x=None):27 r = np.asarray(r, float)28 out = np.full(len(r), np.nan)29 for t in range(1, len(r)):30 w = r[max(0, t - self.window):t]31 out[t] = float(np.var(w, ddof=0)) if len(w) >= 2 else np.nan32 out[0] = out[1]33 return out3435 def forecast(self, fit, r, x, h, **kw):36 r = np.asarray(r, float)37 v = float(np.var(r[-self.window:], ddof=0))38 return Forecast(self.model_id, "empirical", np.full(h, v), float(np.mean(r[-self.window:])), np.empty(0), "direct", v)394041class RMRandomWalk:42 """Realized-measure random walk: forecast of each future day's variance = last observed full-day proxy."""43 family = "baseline"44 model_id = "rm_rw"45 needs_rm = True4647 def fit(self, r, x=None, dist="empirical", **kw):48 return FitResult(self.model_id, "empirical", ["mu"], np.array([float(np.mean(r))]), float("nan"), True, len(r), aux={})4950 def filter(self, fit, r, x=None):51 x = np.asarray(x, float)52 out = np.empty(len(x))53 out[0] = x[0]54 out[1:] = x[:-1]55 return out5657 def forecast(self, fit, r, x, h, **kw):58 v = float(x[-1])59 return Forecast(self.model_id, "empirical", np.full(h, v), float(np.mean(r[-250:])), np.empty(0), "direct", v)606162def _har_X(y: np.ndarray) -> np.ndarray:63 """Regressors at origin t (row t): [1, y_t, mean(y_{t-4..t}), mean(y_{t-21..t})]; NaN for t < 21."""64 T = len(y)65 X = np.full((T, 4), np.nan)66 c = np.cumsum(np.insert(y, 0, 0.0))67 for t in range(21, T):68 X[t] = [1.0, y[t], (c[t + 1] - c[t - 4]) / 5.0, (c[t + 1] - c[t - 21]) / 22.0]69 return X707172class HARRV:73 """HAR-RV (Corsi 2009) on the full-day variance proxy y_t (percent^2): for horizon h the h-day cumulative target74 sum_{k=1..h} y_{t+k} is projected directly on [1, y_t, y_t^(w), y_t^(m)] by OLS (one regression per horizon, 'direct'75 multi-step). log=True estimates the regression in logs of the per-day average and back-transforms with Duan's76 smearing factor (mean of exp(residuals)). Forecasts are floored at 1e-8 percent^2."""77 family = "realized"78 needs_rm = True7980 def __init__(self, horizons=(1, 5, 22), log: bool = False):81 self.horizons = tuple(horizons)82 self.log = log83 self.model_id = "loghar_rv" if log else "har_rv"8485 def fit(self, r, x=None, dist="empirical", **kw):86 t0 = time.perf_counter()87 y = np.asarray(x, float)88 X = _har_X(np.log(np.maximum(y, 1e-8)) if self.log else y)89 T = len(y)90 c = np.cumsum(np.insert(y, 0, 0.0))91 coefs, smear = {}, {}92 ok = True93 for h in self.horizons:94 tgt = np.full(T, np.nan)95 for t in range(T - h):96 tgt[t] = (c[t + 1 + h] - c[t + 1]) / h # average daily variance over t+1..t+h97 if self.log:98 tgt = np.log(np.maximum(tgt, 1e-8))99 m = np.isfinite(tgt) & np.all(np.isfinite(X), axis=1)100 if m.sum() < 50:101 ok = False102 coefs[h] = np.array([np.nan] * 4)103 continue104 b, *_ = np.linalg.lstsq(X[m], tgt[m], rcond=None)105 coefs[h] = b106 if self.log:107 res = tgt[m] - X[m] @ b108 smear[h] = float(np.mean(np.exp(res)))109 params = np.concatenate([coefs[h] for h in self.horizons])110 names = [f"b{h}_{k}" for h in self.horizons for k in ("c", "d", "w", "m")]111 return FitResult(self.model_id, "empirical", names, params, float("nan"), ok, int(T), time_s=time.perf_counter() - t0,112 admissible=ok, aux={"coefs": coefs, "smear": smear})113114 def filter(self, fit, r, x=None):115 y = np.asarray(x, float)116 X = _har_X(np.log(np.maximum(y, 1e-8)) if self.log else y)117 b = fit.aux["coefs"][1]118 f = X @ b119 if self.log:120 f = np.exp(f) * fit.aux["smear"].get(1, 1.0)121 out = np.full(len(y), np.nan)122 out[1:] = f[:-1] # forecast for day t made at t-1123 return np.maximum(out, 1e-8)124125 def forecast(self, fit, r, x, h, **kw):126 y = np.asarray(x, float)127 X = _har_X(np.log(np.maximum(y, 1e-8)) if self.log else y)[-1]128 step = np.full(h, np.nan)129 # fill per-step variances consistent with the direct forecasts of the cumulative targets130 prev_h, prev_cum = 0, 0.0131 avail = sorted(hh for hh in self.horizons if hh <= h)132 for hh in avail:133 f = float(X @ fit.aux["coefs"][hh])134 if self.log:135 f = float(np.exp(f) * fit.aux["smear"].get(hh, 1.0))136 cum = max(f, 1e-8) * hh137 inc = max(cum - prev_cum, 1e-8 * (hh - prev_h))138 step[prev_h:hh] = inc / (hh - prev_h)139 prev_h, prev_cum = hh, prev_cum + inc140 if prev_h < h:141 step[prev_h:] = step[prev_h - 1]142 return Forecast(self.model_id, "empirical", step, float(np.mean(r[-250:])), np.empty(0), "direct", float(y[-1]))143