Python 75.9%
TeX 24%
1"""2================================================================3Auteur : Simon-Pierre Boucher4Contact : contact@spboucher.ai5Projet : Prévision de volatilité réalisée multi-actifs6 (HAR-RV vs GARCH vs Machine Learning)7Fichier : realized_vol.py8Description : Mesures de volatilité réalisée journalières —9 RV 1-min, RV 5-min sous-échantillonnée, realized10 kernel (Parzen, BNHLS 2008), bipower variation,11 test de sauts BNS, semivariances signées et12 signed jumps. Gestion des sessions par classe.13================================================================14"""1516from __future__ import annotations1718import logging19from math import gamma as _gamma2021import numpy as np22import pandas as pd23from scipy import stats2425logger = logging.getLogger(__name__)2627MU1 = np.sqrt(2.0 / np.pi) # E|Z|28MU43 = 2 ** (2 / 3) * _gamma(7 / 6) / _gamma(1 / 2) # E|Z|^{4/3}29THETA_BNS = (np.pi**2 / 4.0) + np.pi - 5.0 # variance constant of the BNS ratio test30PARZEN_CSTAR = 3.5134 # optimal bandwidth constant, Parzen kernel3132# minimum number of intraday bars for a valid day, per asset class33MIN_BARS = {"equity": 250, "fx": 700, "crypto": 700, "futures": 700, "index": 250}3435# regular trading hours filter (inclusive start, exclusive end), or None for 24h36SESSION = {37 "equity": ("09:30", "16:00"),38 "index": ("09:30", "16:00"),39 "fx": None,40 "crypto": None,41 "futures": None,42}434445# ---------------------------------------------------------------- primitives46def log_returns(prices: np.ndarray) -> np.ndarray:47 """Log returns of a strictly positive price array.4849 Parameters50 ----------51 prices : numpy.ndarray52 Intraday price levels.5354 Returns55 -------56 numpy.ndarray57 First differences of log prices (length ``len(prices) - 1``).58 """59 return np.diff(np.log(prices))606162def realized_variance(returns: np.ndarray) -> float:63 """Plain realized variance: sum of squared intraday returns."""64 return float(np.sum(returns**2))656667def rv_subsampled(prices: pd.Series, grid_min: int = 5, offsets: int = 5) -> float:68 """Subsampled sparse-grid realized variance.6970 Averages the realized variance computed on ``offsets`` staggered71 ``grid_min``-minute grids, following Zhang, Mykland and72 Aït-Sahalia (2005).7374 Parameters75 ----------76 prices : pandas.Series77 1-minute prices indexed by timestamp (one trading day).78 grid_min : int79 Sparse grid spacing in minutes.80 offsets : int81 Number of staggered grids averaged.8283 Returns84 -------85 float86 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.nan949596def bipower_variation(returns: np.ndarray) -> float:97 """Realized bipower variation (Barndorff-Nielsen & Shephard 2004).9899 Robust to jumps; scaled to estimate integrated variance.100 """101 n = len(returns)102 if n < 3:103 return np.nan104 absr = np.abs(returns)105 return float(MU1**-2 * (n / (n - 1)) * np.sum(absr[1:] * absr[:-1]))106107108def 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.nan113 a = np.abs(returns) ** (4 / 3)114 return float(n * MU43**-3 * (n / (n - 2)) * np.sum(a[2:] * a[1:-1] * a[:-2]))115116117def bns_jump_test(rv: float, bv: float, tq: float, n: int) -> float:118 """Ratio-form BNS jump test statistic (Barndorff-Nielsen & Shephard 2006).119120 Parameters121 ----------122 rv, bv, tq : float123 Realized variance, bipower variation and tripower quarticity.124 n : int125 Number of intraday returns.126127 Returns128 -------129 float130 Asymptotically N(0,1) statistic; large positive values indicate a131 jump day.132 """133 if not np.isfinite(rv) or not np.isfinite(bv) or bv <= 0 or n < 4:134 return np.nan135 ratio = max(tq / bv**2, 1.0) if np.isfinite(tq) else 1.0136 denom = np.sqrt(THETA_BNS * ratio / n)137 return float((1.0 - bv / rv) / denom) if denom > 0 else np.nan138139140def semivariances(returns: np.ndarray) -> tuple[float, float]:141 """Positive and negative realized semivariance (BNKS 2010).142143 Returns144 -------145 tuple of float146 ``(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 )152153154def parzen_kernel(x: np.ndarray) -> np.ndarray:155 """Parzen kernel weights on [0, 1]."""156 w = np.zeros_like(x)157 m1 = x <= 0.5158 m2 = (x > 0.5) & (x <= 1.0)159 w[m1] = 1 - 6 * x[m1] ** 2 + 6 * x[m1] ** 3160 w[m2] = 2 * (1 - x[m2]) ** 3161 return w162163164def realized_kernel(returns: np.ndarray, iv_proxy: float | None = None) -> float:165 """Realized kernel with Parzen weights (Barndorff-Nielsen et al. 2008).166167 The bandwidth follows the feasible rule168 :math:`H^* = c^* \\xi^{4/5} n^{3/5}` with169 :math:`\\xi^2 = \\hat\\omega^2 / \\sqrt{\\widehat{IQ}}` approximated by170 :math:`\\hat\\omega^2 / IV` where the noise variance is estimated as171 ``RV_dense / (2 n)`` and *IV* by a sparse-grid RV.172173 Parameters174 ----------175 returns : numpy.ndarray176 Dense (1-minute) intraday log returns.177 iv_proxy : float, optional178 Noise-robust estimate of integrated variance; defaults to the179 realized variance of the input returns when omitted.180181 Returns182 -------183 float184 Realized kernel estimate of integrated variance (non-negative by185 construction with non-flat-top Parzen weights).186 """187 n = len(returns)188 if n < 10:189 return np.nan190 rv_dense = realized_variance(returns)191 iv = iv_proxy if iv_proxy and np.isfinite(iv_proxy) and iv_proxy > 0 else rv_dense192 omega2 = rv_dense / (2.0 * n)193 xi2 = omega2 / iv if iv > 0 else 0.0194 H = int(np.clip(PARZEN_CSTAR * xi2**0.4 * n**0.6, 1, n - 1))195 gamma0 = float(returns @ returns)196 k = gamma0197 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] * gam201 return float(max(k, 0.0))202203204# ---------------------------------------------------------------- daily loop205def _day_measures(day: pd.DataFrame, alpha: float) -> dict[str, float]:206 """Compute all volatility measures for one trading day.207208 Parameters209 ----------210 day : pandas.DataFrame211 1-minute bars of a single day (columns ``datetime``, ``close``).212 alpha : float213 Size of the BNS jump test.214215 Returns216 -------217 dict218 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 quantities223 p5 = prices.iloc[::5]224 r5 = log_returns(p5.to_numpy())225226 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.0236 cont = rv5_plain - jump237 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.nan240241 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 }260261262def 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.269270 Parameters271 ----------272 bars : pandas.DataFrame273 1-minute bars (columns ``datetime``, ``open``, ``high``, ``low``,274 ``close``, ``volume``) in exchange-local time.275 cls : str276 Asset class (``equity``, ``fx``, ``crypto``, ``futures``,277 ``index``); controls session filtering and the valid-day278 threshold.279 jump_alpha : float280 Size of the BNS jump test used to shrink the jump component.281 min_bars : int, optional282 Override the per-class minimum number of intraday bars.283284 Returns285 -------286 pandas.DataFrame287 One row per valid trading day, indexed by date, with all288 realized measures plus the close-to-close daily log return289 (``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.time297 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]302303 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 continue308 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 out318