#!/usr/bin/env python3 # MIT License # # Copyright (c) 2026 Options Lab (ashares options learning site) # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE # SOFTWARE. """Is the day IV hits an extreme the best day to sell options? Reproduce the headline numbers. For learning only. Past results do not guarantee future ones; not investment advice. What it does ------------ 1. Downloads public data with AKShare: daily closes of six A-share ETFs and their tracking indices (Sina), and each ETF option's QVIX, a VIX-style 30-day implied-vol index. 2. Builds a daily panel: every QVIX print date of every ETF whose next 30 calendar days have closed. Per row: the ATM-IV estimate (QVIX plus a fixed per-ETF offset, see ATM_OFFSET below), the realised vol of the next 30 days, and point-in-time features of the IV series: its percentile over the last 365 days, the trading days since its 60-day peak, and how fast it rose. 3. The premium an option seller collects = IV - the realised vol that followed, in vol points. Reported: * by group (IV high and fresh / high and settled / middle / low), with 95 % confidence intervals from a 3-month block bootstrap; * by trading days since the 60-day IV peak, on days with IV high; * on "extreme" episodes (QVIX in its own top 2 %): selling on the peak, the first day above the line, or the first day above a line known on that day, against waiting 5 / 10 / 20 trading days; * the forecast form: next-month realised vol = IV - c against k x IV, out of sample, with Diebold-Mariano t; * the "just spiked" penalty: fitted on 2015-2020, checked on 2021 on. Run --- pip install -r requirements.txt python iv_timing.py # download, then print the numbers (~1 minute) python iv_timing.py --cache-dir data # keep the download as CSV and reuse it next time python iv_timing.py --end none # use every month available, not just the published run The CSVs are the same as vol_forecast.py's, so one --cache-dir serves both. Python 3.10+. Vol is in decimals inside the code (0.20 = 20 %); the printout is in vol points (1 pt = 0.01). """ from __future__ import annotations import argparse import bisect import math import statistics from datetime import date, timedelta from pathlib import Path import numpy as np import pandas as pd # --------------------------------------------------------------------------- # Data # --------------------------------------------------------------------------- # code: (name, QVIX series, AKShare QVIX function, Sina ETF symbol, tracking index, Sina index symbol) UNDERLYINGS = { "510050": ("50ETF", "50etf", "index_option_50etf_qvix", "sh510050", "000016", "sh000016"), "510300": ("300ETF", "300etf", "index_option_300etf_qvix", "sh510300", "000300", "sh000300"), "510500": ("500ETF", "500etf", "index_option_500etf_qvix", "sh510500", "000905", "sh000905"), "159915": ("ChiNext ETF", "cyb", "index_option_cyb_qvix", "sz159915", "399006", "sz399006"), "588000": ("STAR 50 ETF", "kcb", "index_option_kcb_qvix", "sh588000", "000688", "sh000688"), "159901": ("SZSE 100 ETF", "100etf", "index_option_100etf_qvix", "sz159901", "399330", "sz399330"), } CODES = tuple(UNDERLYINGS) # QVIX includes the skew and sits above at-the-money (ATM) IV. The study turns # QVIX into an ATM-IV estimate with one constant per ETF: the median gap # (ATM IV - QVIX) between the site's own ATM IV, read from the option chains, # and QVIX, 2025-10-01 to 2026-09-29. That one input is not public, so it is # fixed here (decimal; -0.0179 = -1.79 vol points). ATM_OFFSET = { "510050": -0.0178583, "510300": -0.0218805, "510500": -0.0052150, "159915": -0.0114000, "588000": -0.0202000, "159901": -0.0143403, } PUBLISHED_END = "2026-09-30" # the data the published numbers were computed on def fetch(cache_dir: str | Path | None = None) -> dict[str, pd.DataFrame]: """Download (or reload) the three tables. Returns long frames: ``qvix`` (series, date, close), ``etf`` and ``idx`` (code, date, close). ``close`` is in index points for QVIX (18.5 = 18.5 %). """ names = ("qvix", "etf", "idx") if cache_dir is not None: paths = {n: Path(cache_dir) / f"{n}.csv" for n in names} if all(p.exists() for p in paths.values()): return {n: pd.read_csv(p, parse_dates=["date"], dtype={"code": str, "series": str}) for n, p in paths.items()} import akshare as ak # imported here so compute() works without it qv, etf, idx = [], [], [] for code, (_, series, qfunc, etf_sym, ix_code, ix_sym) in UNDERLYINGS.items(): q = getattr(ak, qfunc)().rename(columns=str.lower)[["date", "close"]] qv.append(q.assign(series=series)) e = ak.fund_etf_hist_sina(symbol=etf_sym)[["date", "close"]] etf.append(e.assign(code=code)) i = ak.stock_zh_index_daily(symbol=ix_sym)[["date", "close"]] idx.append(i.assign(code=ix_code)) out = { "qvix": pd.concat(qv, ignore_index=True), "etf": pd.concat(etf, ignore_index=True), "idx": pd.concat(idx, ignore_index=True), } for df in out.values(): df["date"] = pd.to_datetime(df["date"]) if cache_dir is not None: Path(cache_dir).mkdir(parents=True, exist_ok=True) for n, df in out.items(): df.to_csv(Path(cache_dir) / f"{n}.csv", index=False) return out def cut(frames: dict[str, pd.DataFrame], end: str | None) -> dict[str, pd.DataFrame]: if end is None: return frames return {k: v[v["date"] <= pd.Timestamp(end)] for k, v in frames.items()} def _series(df: pd.DataFrame, col: str, key: str) -> pd.Series: s = df[df[col] == key].set_index("date")["close"].astype(float).sort_index() return s[~s.index.duplicated(keep="last")] def qvix_series(frames, code: str) -> pd.Series: s = frames["qvix"][frames["qvix"]["series"] == UNDERLYINGS[code][1]].set_index("date")["close"].astype(float) return (s[s > 0] / 100.0).sort_index() # QVIX is quoted in %; zero rows are gaps def price_series(frames, code: str) -> pd.Series: return _series(frames["etf"], "code", code) def index_series(frames, code: str) -> pd.Series: return _series(frames["idx"], "code", UNDERLYINGS[code][4]) # A return is a split / consolidation, not a price move, when the ETF moved # > 10 % while its index stayed calm (< 3 %), the gap is > 5 %, and the move # is not a +-10 % daily limit move. CA_GAP, CA_MIN_MOVE, CA_INDEX_CALM, LIMIT_TOL = 0.05, 0.10, 0.03, 0.001 def corporate_action_dates(frames, code: str) -> list[date]: df = pd.DataFrame({"px": price_series(frames, code), "ix": index_series(frames, code)}).dropna() r_px, r_ix = np.log(df["px"]).diff(), np.log(df["ix"]).diff() ok = r_px.notna() & r_ix.notna() r_px, r_ix = r_px[ok], r_ix[ok] at_limit = np.zeros(len(r_px), dtype=bool) for lim in (np.log(0.9), np.log(1.1)): at_limit |= ((r_px - lim).abs() < LIMIT_TOL).to_numpy() ca = ((r_px - r_ix).abs() > CA_GAP) & (r_ix.abs() < CA_INDEX_CALM) & (r_px.abs() > CA_MIN_MOVE) & ~at_limit return [ts.date() for ts in r_px.index[ca.to_numpy()]] # --------------------------------------------------------------------------- # Realised vol: the next 30 calendar days, every trading-day return # --------------------------------------------------------------------------- FORWARD_DAYS = 30 MIN_RETURNS = 10 def forward_rv(points, origin: date, drops=()) -> float | None: """Annualised close-to-close vol of the returns in (origin, origin + 30 days]. Log returns, population std (divide by n), x sqrt(252); a return that lands on a corporate-action date is dropped. """ as_of = origin + timedelta(days=FORWARD_DAYS) window = [(d, p) for d, p in points if origin <= d <= as_of and p > 0] if len(window) < MIN_RETURNS + 1: return None drop = set(drops) rets = [math.log(p1 / p0) for (_, p0), (d1, p1) in zip(window, window[1:]) if d1 not in drop] if len(rets) < MIN_RETURNS: return None m = sum(rets) / len(rets) return math.sqrt(sum((r - m) ** 2 for r in rets) / len(rets) * 252) # --------------------------------------------------------------------------- # Point-in-time features of the IV series (only prints dated <= t) # --------------------------------------------------------------------------- PCT_DAYS = 365 # IV percentile window (calendar days), as on the site PEAK_PRINTS = 60 # trailing window for the IV peak D5_TOP_Q, D5_TOP_MIN = 0.90, 120 LM_LAG, LM_LAGS = 21, 6 HIGH, LOW = 0.70, 0.20 def iv_percentile(hist: list[float], current: float) -> float | None: """The site's IV percentile: share of the last year's prints <= today, after a 2-sd filter.""" if len(hist) < 20: return None mean, sd = statistics.mean(hist), statistics.stdev(hist) filt = [v for v in hist if mean - 2 * sd <= v <= mean + 2 * sd] if len(filt) < 5: filt = hist return sum(1 for v in filt if v <= current) / len(filt) def iv_features(iv: pd.Series) -> pd.DataFrame: """Per print date: IV percentile, prints since the 60-print peak, fast-rise flag, slow-belief gap.""" v, dates = iv.to_numpy(), iv.index n = len(v) lo365 = dates.searchsorted(dates - pd.Timedelta(days=PCT_DAYS), side="left") pct, dsp, lm = (np.full(n, np.nan) for _ in range(3)) w = 0.5 ** np.arange(LM_LAGS) w = w / w.sum() for i in range(n): p = iv_percentile(v[lo365[i]: i + 1].tolist(), v[i]) pct[i] = np.nan if p is None else p if i + 1 >= PEAK_PRINTS: dsp[i] = int(np.argmax(v[i + 1 - PEAK_PRINTS: i + 1][::-1])) # 0 = today is the peak (latest on ties) if i >= LM_LAG * LM_LAGS: # IV minus a slowly updated belief (Lochstoer-Muir) lm[i] = v[i] - float(w @ v[[i - LM_LAG * k for k in range(1, LM_LAGS + 1)]]) d5 = pd.Series(v, index=dates).diff(5).to_numpy() d5_top = np.zeros(n, dtype=bool) # 5-print rise in the top 10 % of the last year's rises for i in range(n): h = d5[lo365[i]: i + 1] h = h[~np.isnan(h)] if len(h) >= D5_TOP_MIN and not np.isnan(d5[i]): d5_top[i] = d5[i] >= np.quantile(h, D5_TOP_Q) return pd.DataFrame({"iv": v, "iv_pct": pct, "days_since_peak": dsp, "d5_top": d5_top, "lm_gap": 100 * lm}, index=dates) def group(iv_pct: float, dsp: float, d5_top: bool) -> str | None: if np.isnan(iv_pct) or np.isnan(dsp): return None if iv_pct >= HIGH: if dsp <= 5 or d5_top: return "high & fresh" return "high & settled" if dsp >= 10 else "high, 6-9 d" return "low" if iv_pct <= LOW else "middle" # --------------------------------------------------------------------------- # The daily panel # --------------------------------------------------------------------------- GRID_START = date(2015, 5, 11) N_GRIDS = 21 TRAIN_END, TEST_START = pd.Timestamp("2020-12-31"), pd.Timestamp("2021-01-01") def origin_grid(calendar: list[date]) -> list[date]: """The vol premium study's grid: from 2015-05-11, each date the first trading day >= previous + 30 days.""" i = bisect.bisect_left(calendar, GRID_START) out = [] while i < len(calendar): out.append(calendar[i]) i = bisect.bisect_left(calendar, calendar[i] + timedelta(days=FORWARD_DAYS)) return out def build_panel(frames) -> pd.DataFrame: calendar = [ts.date() for ts in price_series(frames, "510050").index] phase0 = set(origin_grid(calendar)) grids = {d: i % N_GRIDS for i, d in enumerate(calendar)} # 21 grids, every 21st trading day parts = [] for code in CODES: feats = iv_features(qvix_series(frames, code)) px = price_series(frames, code) px_dates = [ts.date() for ts in px.index] pts = list(zip(px_dates, px.to_numpy(dtype=float))) drops = corporate_action_dates(frames, code) rows = [] for ts in feats.index: t = ts.date() if t + timedelta(days=FORWARD_DAYS) > px_dates[-1] or ts not in px.index: continue fwd = pts[bisect.bisect_left(px_dates, t): bisect.bisect_right(px_dates, t + timedelta(days=FORWARD_DAYS))] rv = forward_rv(fwd, t, drops) if rv is not None: rows.append((ts, rv, t in phase0, grids.get(t, -1))) out = feats.loc[[r[0] for r in rows]].copy() out["rv_fwd"] = [r[1] for r in rows] out["phase0"] = [r[2] for r in rows] out["tgrid"] = [r[3] for r in rows] out["iv_atm"] = out["iv"] + ATM_OFFSET[code] out = out.rename_axis("origin").reset_index() out.insert(0, "underlying", code) parts.append(out) df = pd.concat(parts, ignore_index=True) df["group"] = [group(p, d, t) for p, d, t in zip(df["iv_pct"], df["days_since_peak"], df["d5_top"])] df["prem"] = 100 * (df["iv_atm"] - df["rv_fwd"]) # what the seller collects, vol pts df["month"] = df["origin"].dt.to_period("M") end = df["origin"] + pd.Timedelta(days=FORWARD_DAYS) df["split"] = np.where(end <= TRAIN_END, "train", np.where(df["origin"] >= TEST_START, "test", None)) return df # --------------------------------------------------------------------------- # Statistics: describe, and a moving-block bootstrap over calendar months # --------------------------------------------------------------------------- # Windows overlap (a new sample every day, each looking 30 days ahead) and the # six ETFs move together, so a plain sd / sqrt(n) would be far too narrow. # Resample blocks of 3 consecutive origin months, all ETFs together. N_BOOT, BLOCK, SEED = 2000, 3, 11 def describe(x) -> dict: x = pd.Series(x, dtype=float).dropna() p5 = x.quantile(0.05) sd = x.std(ddof=1) return dict(n=len(x), mean=x.mean(), sd=sd, mean_sd=x.mean() / sd, share_pos=(x > 0).mean(), worst5=x[x <= p5].mean(), worst=x.min()) def block_weights(n_months: int) -> np.ndarray: """(draws, months) multiplicities of a moving-block bootstrap; each row sums to n_months.""" rng = np.random.default_rng(SEED) b = min(BLOCK, n_months) n_blocks = math.ceil(n_months / b) starts = rng.integers(0, n_months - b + 1, size=(N_BOOT, n_blocks)) picks = (starts[:, :, None] + np.arange(b)).reshape(N_BOOT, -1)[:, :n_months] w = np.zeros((N_BOOT, n_months)) np.add.at(w, (np.repeat(np.arange(N_BOOT), n_months), picks.ravel()), 1.0) return w class MonthBootstrap: """Statistics rebuilt from per-month sums, so every draw keeps whole months together.""" def __init__(self, months: pd.Series): cats = pd.PeriodIndex(sorted(months.unique())) self.month_of_row = cats.get_indexer(months) self.n_months = len(cats) self.w = block_weights(self.n_months) def mean_and_ratio(self, x, mask): x, mask = np.asarray(x, dtype=float), np.asarray(mask, dtype=bool) m = mask & ~np.isnan(x) k = self.month_of_row[m] n = self.w @ np.bincount(k, minlength=self.n_months) s1 = self.w @ np.bincount(k, weights=x[m], minlength=self.n_months) s2 = self.w @ np.bincount(k, weights=x[m] ** 2, minlength=self.n_months) with np.errstate(invalid="ignore", divide="ignore"): mean = s1 / n var = (s2 / n - mean ** 2) * n / (n - 1) return mean, mean / np.sqrt(var) def ols(self, X: np.ndarray, y: np.ndarray) -> np.ndarray: p = X.shape[1] xtx, xty = np.zeros((self.n_months, p * p)), np.zeros((self.n_months, p)) np.add.at(xtx, self.month_of_row, (X[:, :, None] * X[:, None, :]).reshape(len(X), -1)) np.add.at(xty, self.month_of_row, X * y[:, None]) A, B = (self.w @ xtx).reshape(-1, p, p), self.w @ xty return np.linalg.solve(A, B[:, :, None])[:, :, 0] def ci(draws: np.ndarray) -> tuple[float, float]: d = draws[np.isfinite(draws)] return float(np.quantile(d, 0.025)), float(np.quantile(d, 0.975)) def newey_west_t(x: np.ndarray, lags: int = 2) -> float: x = np.asarray(x, dtype=float) n = len(x) u = x - x.mean() s = u @ u / n for lag in range(1, lags + 1): s += 2 * (1 - lag / (lags + 1)) * (u[lag:] @ u[:-lag]) / n return float(x.mean() / math.sqrt(s / n)) # --------------------------------------------------------------------------- # Groups and time since the peak # --------------------------------------------------------------------------- HIGH_GROUPS = ("high & fresh", "high & settled", "high, 6-9 d") GROUPS = ("high & fresh", "high & settled", "all high", "middle", "low") DSP_BUCKETS = [(0, 0), (1, 2), (3, 5), (6, 9), (10, 19), (20, 39), (40, 59)] def group_masks(df: pd.DataFrame) -> dict[str, np.ndarray]: m = {g: (df["group"] == g).to_numpy() for g in ("high & fresh", "high & settled", "middle", "low")} m["all high"] = df["group"].isin(HIGH_GROUPS).to_numpy() return m def group_stats(sample: pd.DataFrame) -> dict: mb, masks, x = MonthBootstrap(sample["month"]), group_masks(sample), sample["prem"].to_numpy() out, draws = {}, {} for g in GROUPS: d = describe(x[masks[g]]) mean_d, ratio_d = mb.mean_and_ratio(x, masks[g]) d["dates"] = int(sample.loc[masks[g], "origin"].nunique()) d["mean_ci"], d["ratio_ci"] = ci(mean_d), ci(ratio_d) out[g], draws[g] = d, (mean_d, ratio_d) for a, b in (("all high", "middle"), ("high & fresh", "high & settled")): out[f"{a} - {b}"] = dict( mean=out[a]["mean"] - out[b]["mean"], mean_ci=ci(draws[a][0] - draws[b][0]), ratio=out[a]["mean_sd"] - out[b]["mean_sd"], ratio_ci=ci(draws[a][1] - draws[b][1]), share_le_0=float(np.mean((draws[a][0] - draws[b][0]) <= 0))) return out def since_peak(sample: pd.DataFrame) -> dict: mb, x = MonthBootstrap(sample["month"]), sample["prem"].to_numpy() high = (sample["iv_pct"] >= HIGH).to_numpy() dsp = sample["days_since_peak"].to_numpy() out = {} for lo, hi in DSP_BUCKETS: m = high & (dsp >= lo) & (dsp <= hi) d = describe(x[m]) mean_d, _ = mb.mean_and_ratio(x, m) d["dates"], d["mean_ci"] = int(sample.loc[m, "origin"].nunique()), ci(mean_d) out[f"{lo}" if lo == hi else f"{lo}-{hi}"] = d return out # --------------------------------------------------------------------------- # Event study: IV in its own top 2 % # --------------------------------------------------------------------------- TOP_Q, MAX_GAP, MERGE_GAP, EXPANDING_MIN = 0.98, 20, 20, 250 LAGS = (0, 5, 10, 20) def find_episodes(iv: pd.Series, thr: pd.Series) -> list[dict]: """Runs of prints with IV >= the line; a gap of more than 20 prints starts a new episode.""" v, t = iv.to_numpy(), thr.reindex(iv.index).to_numpy() hits = np.flatnonzero(np.nan_to_num(v >= t, nan=False).astype(bool) & ~np.isnan(t)) if hits.size == 0: return [] out = [] for chunk in np.split(hits, np.flatnonzero(np.diff(hits) > MAX_GAP + 1) + 1): lo, hi = int(chunk[0]), int(chunk[-1]) peak = lo + int(np.argmax(v[lo: hi + 1])) out.append(dict(first_pos=lo, peak_pos=peak, peak=iv.index[peak], peak_iv=float(v[peak]))) return out def episodes(frames, expanding: bool) -> pd.DataFrame: rows = [] for code in CODES: iv = qvix_series(frames, code) thr = (iv.expanding(min_periods=EXPANDING_MIN).quantile(TOP_Q) if expanding else pd.Series(float(iv.quantile(TOP_Q)), index=iv.index)) rows += [dict(underlying=code, **e) for e in find_episodes(iv, thr)] ep = pd.DataFrame(rows) # ETFs whose peaks fall within 20 trading days of a cluster's first peak are one independent episode. calendar = price_series(frames, "510050").index pos = pd.Series(calendar.searchsorted(ep["peak"]), index=ep.index).sort_values(kind="stable") ep["event"], current, anchor = 0, -1, None for i, p in pos.items(): if anchor is None or p - anchor > MERGE_GAP: current, anchor = current + 1, p ep.loc[i, "event"] = current return ep def event_ladder(frames, panel: pd.DataFrame, expanding: bool, anchor: str) -> dict: """Sell on day 0 (the peak or the first day above the line) or 5 / 10 / 20 prints later.""" ep = episodes(frames, expanding) prem = panel.set_index(["underlying", "origin"])["prem"] rows = [] for _, e in ep.iterrows(): iv = qvix_series(frames, e["underlying"]) for lag in LAGS: pos = e[anchor] + lag day = iv.index[pos] if pos < len(iv) else None rows.append(dict(event=e["event"], lag=lag, prem=prem.get((e["underlying"], day), np.nan))) per_event = pd.DataFrame(rows).groupby(["event", "lag"])["prem"].mean().unstack("lag") # ETFs averaged first peaks = ep.groupby("event")["peak"].min() out = dict(episodes=len(per_event), peaks=[str(p.date()) for p in peaks], per_event=per_event) for lag in LAGS: x = per_event[lag].dropna() d = (per_event[lag] - per_event[0]).dropna() out[lag] = dict(mean=x.mean(), worst=x.min(), waited_better=int((d > 0).sum()) if lag else None) return out # --------------------------------------------------------------------------- # Forecast form: IV - c against k x IV, out of sample # --------------------------------------------------------------------------- FLOOR, OOS_MIN = 0.01, 24 LOSSES = { "squared": lambda t, f: 1e4 * (t - f) ** 2, "log": lambda t, f: np.log(t / f) ** 2, "qlike": lambda t, f: (t / f) ** 2 - np.log((t / f) ** 2) - 1, } def oos_forms(df: pd.DataFrame) -> pd.DataFrame: """At each date, fit on samples whose 30-day window had closed (pooled); floor at 1 %. F1: IV - c, c = mean(IV - RV). F2: k x IV, least squares. F2': k x IV, k = mean RV / mean IV. """ df = df.dropna(subset=["iv_atm", "rv_fwd"]).sort_values(["origin", "underlying"]).reset_index(drop=True) iv, rv = df["iv_atm"].to_numpy(), df["rv_fwd"].to_numpy() out = {k: np.full(len(df), np.nan) for k in ("f1", "f2", "f2r")} for t, idx in df.groupby("origin").groups.items(): past = (df["origin"] + pd.Timedelta(days=FORWARD_DAYS) <= t).to_numpy() if past.sum() < OOS_MIN: continue pos = df.index.get_indexer(idx) out["f1"][pos] = iv[pos] - np.mean(iv[past] - rv[past]) out["f2"][pos] = (iv[past] @ rv[past] / (iv[past] @ iv[past])) * iv[pos] out["f2r"][pos] = (rv[past].sum() / iv[past].sum()) * iv[pos] for k, v in out.items(): df[k] = np.maximum(v, FLOOR) return df def dm_t(df: pd.DataFrame, a: str, b: str, loss: str) -> float: """Diebold-Mariano: loss(a) - loss(b) averaged per date, Newey-West t; t > 0 = b is better.""" m = df[[a, b, "rv_fwd"]].notna().all(axis=1) sub = df[m] d = LOSSES[loss](sub["rv_fwd"], sub[a]) - LOSSES[loss](sub["rv_fwd"], sub[b]) return newey_west_t(d.groupby(sub["origin"]).mean().sort_index().to_numpy()) def forecast_form(panel: pd.DataFrame) -> dict: vrp0 = panel[panel["phase0"]] f = oos_forms(vrp0).dropna(subset=["f1"]) t = f["rv_fwd"] table = {} for key in ("f1", "f2", "f2r"): table[key] = dict(n=len(f), bias=100 * (t - f[key]).mean(), rmse=math.sqrt(LOSSES["squared"](t, f[key]).mean()), log_rmse=math.sqrt(LOSSES["log"](t, f[key]).mean()), qlike=LOSSES["qlike"](t, f[key]).mean(), dm=None if key == "f1" else {loss: dm_t(f, "f1", key, loss) for loss in LOSSES}) k = lambda d: float(d["rv_fwd"].sum() / d["iv_atm"].sum()) # noqa: E731 better = dict(rmse=0, log_rmse=0) for g in range(N_GRIDS): # the same out-of-sample race on 21 date grids that never overlap fg = oos_forms(panel[panel["tgrid"] == g]) e1, e2 = fg["rv_fwd"] - fg["f1"], fg["rv_fwd"] - fg["f2r"] better["rmse"] += int((e2 ** 2).mean() < (e1 ** 2).mean()) better["log_rmse"] += int((np.log(fg["rv_fwd"] / fg["f2r"]) ** 2).mean() < (np.log(fg["rv_fwd"] / fg["f1"]) ** 2).mean()) return dict(samples=len(vrp0), dates=int(vrp0["origin"].nunique()), table=table, c_full_pts=100 * float((vrp0["iv_atm"] - vrp0["rv_fwd"]).mean()), k_full=k(vrp0), k_train=k(vrp0[vrp0["split"] == "train"]), k_test=k(vrp0[vrp0["split"] == "test"]), grids_better=better) # --------------------------------------------------------------------------- # The "just spiked" penalty: RV = a + b x IV + c x trigger, train -> test # --------------------------------------------------------------------------- def triggers(df: pd.DataFrame) -> dict[str, np.ndarray]: """The six candidates fixed before the test, all on days with IV percentile >= 0.70.""" high, dsp = (df["iv_pct"] >= HIGH).to_numpy(), df["days_since_peak"].to_numpy() out = {f"peak <= {n}": high & (dsp <= n) for n in (5, 10, 15, 20)} out["fast rise"] = high & df["d5_top"].to_numpy(dtype=bool) out["above slow belief"] = high & (df["lm_gap"].to_numpy() > 0) return out def trigger_fit(df: pd.DataFrame, trig: np.ndarray) -> dict: X = np.column_stack([np.ones(len(df)), df["iv_atm"].to_numpy(), trig.astype(float)]) y = df["rv_fwd"].to_numpy() beta = np.linalg.lstsq(X, y, rcond=None)[0] draws = MonthBootstrap(df["month"]).ols(X, y)[:, 2] lo, hi = ci(draws) return dict(a=beta[0], b=beta[1], c=100 * beta[2], c_ci=(100 * lo, 100 * hi), t=float(beta[2] / np.nanstd(draws, ddof=1))) def spike_penalty(panel: pd.DataFrame) -> dict: df = panel.dropna(subset=["iv_atm", "iv_pct", "days_since_peak"]).reset_index(drop=True) trig = triggers(df) train, test = (df["split"] == "train").to_numpy(), (df["split"] == "test").to_numpy() fits = {name: trigger_fit(df[train].reset_index(drop=True), t[train]) for name, t in trig.items()} chosen = max((n for n in fits if fits[n]["c"] > 0), key=lambda n: fits[n]["t"]) fit, check = fits[chosen], trigger_fit(df[test].reset_index(drop=True), trig[chosen][test]) # Test on the out-of-sample years, one date every ~30 days: train-fitted a + b x IV with and without c. b, a = np.polyfit(df.loc[train, "iv_atm"], df.loc[train, "rv_fwd"], 1) m = test & df["phase0"].to_numpy() chk = df[m].copy() chk["base"] = a + b * chk["iv_atm"] chk["penalised"] = fit["a"] + fit["b"] * chk["iv_atm"] + fit["c"] / 100 * trig[chosen][m] dm = dm_t(chk, "base", "penalised", "squared") kept = fit["c_ci"][0] > 0 and check["c"] > 0 and check["c"] >= 0.5 * fit["c"] and dm >= 0 return dict(chosen=chosen, train=fit, test=check, dm_t=dm, kept=kept, full=trigger_fit(df, trig[chosen])) # --------------------------------------------------------------------------- def compute(frames: dict[str, pd.DataFrame]) -> dict: panel = build_panel(frames) sample = panel[panel["group"].notna()].reset_index(drop=True) return dict( panel=panel, rows=len(panel), dates=int(panel["origin"].nunique()), sample_rows=len(sample), groups=group_stats(sample), since_peak=since_peak(sample), ladder={ "peak": event_ladder(frames, panel, expanding=False, anchor="peak_pos"), "first": event_ladder(frames, panel, expanding=False, anchor="first_pos"), "first, line known that day": event_ladder(frames, panel, expanding=True, anchor="first_pos"), }, forms=forecast_form(panel), spike=spike_penalty(panel), ) def _ci(c) -> str: return f"[{c[0]:+.1f}, {c[1]:+.1f}]" def main(argv=None) -> int: ap = argparse.ArgumentParser(description="IV timing study (public data)") ap.add_argument("--cache-dir", type=Path, help="save / reuse the downloaded CSVs here") ap.add_argument("--end", default=PUBLISHED_END, help="ignore data after this date ('none' = use all)") args = ap.parse_args(argv) r = compute(cut(fetch(args.cache_dir), None if args.end == "none" else args.end)) print(f"daily panel: {r['rows']} rows on {r['dates']} dates; {r['sample_rows']} with a group\n") print("1. The look-ahead ladder: mean IV - RV over the next 30 days, vol pts (ATM IV estimate)") print(f" {'sell on':<30}{'episodes':>9}{'day 0':>8}{'+5':>7}{'+10':>7}{'+20':>7} waiting won (+5/+10/+20) day-0 worst") for name, s in r["ladder"].items(): won = " / ".join(str(s[lag]["waited_better"]) for lag in LAGS[1:]) print(f" {name:<30}{s['episodes']:>9}" + "".join(f"{s[lag]['mean']:>+7.1f}" for lag in LAGS) + f" {won:>10} of {s['episodes']:<12}{s[0]['worst']:+.1f}") pe = r["ladder"]["peak"]["per_event"] print(f" 2015: peak 2015-06-29 day 0 {pe.loc[0, 0]:+.1f}, +20 {pe.loc[0, 20]:+.1f}; " f"peak 2015-08-26 day 0 {pe.loc[1, 0]:+.1f}") print(f" independent episodes (peaks): {', '.join(r['ladder']['peak']['peaks'])}") print("\n2. Days with IV high (percentile >= 0.70), by trading days since the 60-day IV peak") for k, d in r["since_peak"].items(): print(f" {k:>6}: n {d['n']:>4} on {d['dates']:>3} dates, mean {d['mean']:+.1f} {_ci(d['mean_ci'])}, " f"mean/sd {d['mean_sd']:+.2f}, worst 5 % mean {d['worst5']:+.1f}") print("\n3. Groups (every trading day; 95 % CIs, 3-month block bootstrap)") for g in GROUPS: d = r["groups"][g] print(f" {g:<16} {d['n']:>5} / {d['dates']:>4} dates mean {d['mean']:+.1f} {_ci(d['mean_ci'])} " f"mean/sd {d['mean_sd']:+.2f} [{d['ratio_ci'][0]:+.2f}, {d['ratio_ci'][1]:+.2f}] sd {d['sd']:.1f}") for key in ("all high - middle", "high & fresh - high & settled"): d = r["groups"][key] print(f" {key}: mean {d['mean']:+.1f} {_ci(d['mean_ci'])}, mean/sd {d['ratio']:+.2f} " f"[{d['ratio_ci'][0]:+.2f}, {d['ratio_ci'][1]:+.2f}], {d['share_le_0']:.0%} of draws <= 0") fm = r["forms"] print(f"\n4. Forecast form, out of sample ({fm['table']['f1']['n']} forecasts; " f"{fm['samples']} samples on {fm['dates']} dates)") names = {"f1": "IV - c", "f2": "k x IV, least squares", "f2r": "k x IV, k = mean RV / mean IV"} for key, row in fm["table"].items(): dm = "" if row["dm"] is None else " DM t " + " / ".join(f"{v:+.2f}" for v in row["dm"].values()) print(f" {names[key]:<31} bias {row['bias']:+.2f} RMSE {row['rmse']:.2f} log-RMSE {row['log_rmse']:.3f}" f" QLIKE {row['qlike']:.3f}{dm}") print(f" full sample: c = {fm['c_full_pts']:.2f} pts, k = {fm['k_full']:.3f}; " f"k to 2020 {fm['k_train']:.3f}, from 2021 {fm['k_test']:.3f}") print(f" k x IV better than IV - c on {fm['grids_better']['rmse']} of {N_GRIDS} grids (RMSE), " f"{fm['grids_better']['log_rmse']} of {N_GRIDS} (log-RMSE)") sp = r["spike"] print(f"\n5. 'Just spiked' penalty: chosen on 2015-2020: '{sp['chosen']}' (IV percentile >= 0.70)") print(f" c (pts of extra realised vol) train {sp['train']['c']:+.1f} {_ci(sp['train']['c_ci'])}, " f"test {sp['test']['c']:+.1f} {_ci(sp['test']['c_ci'])}, full {sp['full']['c']:+.1f} {_ci(sp['full']['c_ci'])}") print(f" DM t on the test years {sp['dm_t']:+.2f} (> 0 = the penalty helps) -> {'kept' if sp['kept'] else 'not kept'}") return 0 if __name__ == "__main__": raise SystemExit(main())