Warrior_EA/research/altdata/screen_ivol.py

"""Instrument-specific implied-volatility screen.

The options market publishes its own forecast of an instrument's future range,
and CBOE's per-instrument vol indices carry 15-25 years of free daily history
on FRED - unlike GEX, which has no free history at all and has to be recorded
forward. If the options market's view of range is worth anything to us, this is
the version we can test TODAY.

Motivation: SP500's `vix`/`vix_chg5` already survived both the family-wise bar
and the incremental (conditional-on-trailing-range) test. Gold, oil, Nasdaq,
Dow and Russell each have their OWN vol index, so the obvious question is
whether the instrument-specific one beats the equity VIX we currently ship.
XAUUSD is the decisive case: it ships vix_chg5 today, and GVZ is gold's own.

Everything is measured against forward RANGE (the proven channel), with the
trailing-range positive control, a noise negative control, a family-wise bar
over the maximum of the null, and the decisive MI|vol column.

Usage:
    python -m altdata.screen_ivol [--perms 499]
"""
import argparse

import numpy as np
import pandas as pd

from .common import DATA_ROOT
from .screen import _tercile, asof_join, cmi_3x3x3, load_bars, mi_3x3, true_range

#--- symbol -> (index series, human label). The cross-pairings are deliberate:
#--- an instrument's own vol index should beat a foreign one, and if it does not,
#--- that is evidence the whole family is just measuring "vol is vol".
PLAN = {
    "XAUUSD": [("gvz", "GVZCLS", "GOLD's own IV"),
               ("vix", "VIXCLS", "equity VIX - what we ship today"),
               ("ovx", "OVXCLS", "OIL IV - cross-check, should be weaker")],
    "SP500":  [("vix", "VIXCLS", "equity VIX - known clean, the benchmark"),
               ("vix3m", "VXVCLS", "3-month VIX"),
               ("vxn", "VXNCLS", "Nasdaq IV - correlated cousin"),
               ("gvz", "GVZCLS", "GOLD IV - cross-check, should be weaker")],
    "EURUSD": [("vix", "VIXCLS", "equity VIX"),
               ("vxd", "VXDCLS", "Dow IV")],
    "USDJPY": [("vix", "VIXCLS", "equity VIX"),
               ("vxd", "VXDCLS", "Dow IV")],
}


def ivol_frame(series_id: str, name: str) -> pd.DataFrame:
    f = pd.read_csv(DATA_ROOT / "fred" / f"{series_id}.csv",
                    parse_dates=["observed", "published"])
    v = f["value"]
    out = pd.DataFrame({"published": f["published"]})
    out[f"{name}"] = v / 100.0                       # level, same a-priori scale as vix
    out[f"{name}_chg5"] = (v - v.shift(5)) / 10.0    # 5-day change
    return out


def main() -> None:
    ap = argparse.ArgumentParser()
    ap.add_argument("--sym", nargs="*", default=list(PLAN))
    ap.add_argument("--perms", type=int, default=499)
    ap.add_argument("--seed", type=int, default=11)
    args = ap.parse_args()
    rng = np.random.default_rng(args.seed)

    for sym in args.sym:
        bars = load_bars(sym)
        tr = true_range(bars)
        atr = tr.ewm(alpha=1 / 14, min_periods=14).mean()
        trailing = tr.rolling(5).sum() / atr
        y_raw = tr.shift(-1).rolling(5).sum().shift(-4) / atr    # forward 5-day range

        feats = pd.DataFrame(index=bars.index)
        for name, sid, _desc in PLAN[sym]:
            j = asof_join(bars, ivol_frame(sid, name))
            for c in j.columns:
                feats[c] = j[c].to_numpy()
        #--- term structure: contango = calm regime, backwardation = stress
        if "vix3m" in feats and "vix" in feats:
            feats["vix_term"] = feats["vix3m"] / feats["vix"].replace(0, np.nan) - 1.0
        feats["CTRL_trailing_range"] = trailing
        feats["CTRL_noise"] = rng.standard_normal(len(bars))

        valid = y_raw.notna() & trailing.notna()
        y_bin = _tercile(y_raw[valid].to_numpy())
        z_bin = _tercile(trailing[valid].to_numpy())
        print(f"\n=== {sym}: {int(valid.sum())} D1 bars, forward-5d range, {args.perms} perms ===")
        for name, sid, desc in PLAN[sym]:
            print(f"    {name:6s} = {sid:8s} {desc}")

        results, null_max = [], np.zeros(args.perms)
        for col in feats.columns:
            x = feats.loc[valid, col].to_numpy(dtype=float)
            ok = ~np.isnan(x)
            if ok.sum() < 500:
                print(f"  {col:22s} SKIPPED - only {ok.sum()} usable rows")
                continue
            xo, yo, zo = x[ok], y_bin[ok], z_bin[ok]
            n = len(xo)
            obs, obs_c = mi_3x3(xo, yo), cmi_3x3x3(xo, yo, zo)
            null = np.empty(args.perms)
            null_c = np.empty(args.perms)
            for k in range(args.perms):
                xs = np.roll(xo, rng.integers(63, n - 63))
                null[k] = mi_3x3(xs, yo)
                null_c[k] = cmi_3x3x3(xs, yo, zo)
            p = (np.sum(null >= obs) + 1) / (args.perms + 1)
            p_c = (np.sum(null_c >= obs_c) + 1) / (args.perms + 1)
            results.append((col, obs, p, obs_c, p_c, n))
            if not col.startswith("CTRL_"):
                null_max = np.maximum(null_max, null)
        fam = np.quantile(null_max, 0.95)
        print(f"  family-wise 5% bar (max of null over the candidates): {fam:.5f}")
        for col, obs, p, obs_c, p_c, n in sorted(results, key=lambda r: -r[3]):
            mark = "FAMILY" if obs > fam else "      "
            inc = "INCREMENTAL" if p_c <= 0.05 else "           "
            print(f"  {col:22s} n={n:5d} MI {obs:.5f} p={p:.4f} {mark} | "
                  f"MI|vol {obs_c:.5f} p={p_c:.4f} {inc}")


if __name__ == "__main__":
    main()