feat(surfer): equity-index VRP frontier — real+affordable, measurement needs engineering
Bought XSP (mini-SPX) options ohlcv-1d 2013-2026 = $66.86 (hard-capped $70, get_cost-gated, year-chunked vs 504, 2025 backfilled) + SPY $0.01. CORRECTION: index VRP is NOT $900 (that was SPY root); clean instrument XSP $66.86/13y or SPX $143/13y. THREE measurement attempts all gave spurious NEGATIVE VRP (RV 28-33% vs implied ~16%) = MEASUREMENT ERROR not finding (contradicts decades of SPX VRP evidence). Bugs: noisy parity underlying; ~469/2800 days survive -> multi-day gaps inflate RV; SPY!=XSP divergence corrupts ATM. Root: XSP EOD ohlcv too sparse to reconstruct clean underlying+ATM-IV. Implied (16.7%) reads right. Frontier real+affordable+in-hand but proper extraction is a real options-quant build (denser SPX +$76 or IV-surface w/ quotes), not a gate. Did NOT record -2.7 as a result (artifact). Crypto momentum+VRP remains only deploy-grade edge. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This commit is contained in:
jgrusewski
115
scripts/surfer/equity_vrp.py
Normal file
115
scripts/surfer/equity_vrp.py
Normal file
@@ -0,0 +1,115 @@
|
||||
#!/usr/bin/env python3
|
||||
"""Equity-index VRP on XSP (mini-SPX) options, 2013-2026 — the one real-prior frontier.
|
||||
|
||||
For each day: find the ~30d ATM straddle (strike where call~=put = the forward), back out
|
||||
implied vol from the straddle price, compare to forward realized vol. VRP = IV - RV; a
|
||||
short-vol seller harvests it. 13y spans 2018/2020/2022 tail events. Tests raw VRP Sharpe,
|
||||
per-year (incl. crashes), tail/skew, AND the tail-managed version (don't sell when IV rising
|
||||
— the gate that fixed crypto VRP). XSP = cash-settled European -> no early-exercise distortion.
|
||||
"""
|
||||
import glob
|
||||
import math
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
from signal_sweep import sharpe_t # noqa: E402
|
||||
import torch # noqa: E402
|
||||
|
||||
DEV = "cuda" if torch.cuda.is_available() else "cpu"
|
||||
DAY_NS = 86_400 * 10**9
|
||||
|
||||
|
||||
def load_xsp():
|
||||
import databento as db
|
||||
rows = []
|
||||
for p in sorted(glob.glob("data/surfer/xsp/*.dbn")):
|
||||
try:
|
||||
df = db.DBNStore.from_file(p).to_df().reset_index()
|
||||
except Exception:
|
||||
continue
|
||||
if df.empty or "symbol" not in df.columns:
|
||||
continue
|
||||
df = df[df["close"] > 0]
|
||||
rows.append(df[["ts_event", "symbol", "close"]])
|
||||
import pandas as pd
|
||||
d = pd.concat(rows, ignore_index=True)
|
||||
s = d["symbol"].astype(str).str.replace(" ", "", regex=False) # "XSP240119C00450000"
|
||||
d["right"] = s.str[-9]
|
||||
d["strike"] = s.str[-8:].astype(float) / 1000.0
|
||||
d["expiry"] = s.str[-15:-9] # YYMMDD
|
||||
d["day"] = (d["ts_event"].astype("int64") // DAY_NS)
|
||||
return d
|
||||
|
||||
|
||||
def to_epoch_day(yymmdd):
|
||||
import datetime
|
||||
y = 2000 + int(yymmdd[:2]); mo = int(yymmdd[2:4]); da = int(yymmdd[4:6])
|
||||
return (datetime.date(y, mo, da) - datetime.date(1970, 1, 1)).days
|
||||
|
||||
|
||||
def main():
|
||||
d = load_xsp()
|
||||
print(f"loaded {len(d)} XSP option-days, {d['day'].nunique()} trading days")
|
||||
d["exp_day"] = d["expiry"].map(to_epoch_day)
|
||||
d["ttm"] = d["exp_day"] - d["day"]
|
||||
d = d[(d["ttm"] >= 20) & (d["ttm"] <= 45)] # ~30d window
|
||||
|
||||
iv_by_day = {}
|
||||
S_by_day = {}
|
||||
for day, g in d.groupby("day"):
|
||||
# pick the expiry closest to 30d
|
||||
exp = g.iloc[(g["ttm"] - 30).abs().argsort()].iloc[0]["exp_day"]
|
||||
ge = g[g["exp_day"] == exp]
|
||||
calls = ge[ge["right"] == "C"].groupby("strike")["close"].mean() # dedup AM/PM series
|
||||
puts = ge[ge["right"] == "P"].groupby("strike")["close"].mean()
|
||||
common = calls.index.intersection(puts.index)
|
||||
if len(common) < 3:
|
||||
continue
|
||||
diff = (calls[common] - puts[common]).abs()
|
||||
katm = diff.idxmin() # ATM: where C~=P
|
||||
S = float(katm + calls[katm] - puts[katm]) # parity forward
|
||||
straddle = float(calls[katm] + puts[katm])
|
||||
T = float(ge["ttm"].iloc[0]) / 365.0
|
||||
if S <= 0 or T <= 0:
|
||||
continue
|
||||
iv = straddle / (0.8 * S * math.sqrt(T)) # ATM straddle -> implied vol
|
||||
iv_by_day[int(day)] = iv; S_by_day[int(day)] = S
|
||||
|
||||
days = np.array(sorted(S_by_day))
|
||||
S = np.array([S_by_day[x] for x in days]); IV = np.array([iv_by_day[x] for x in days])
|
||||
r = np.zeros(len(days)); r[1:] = np.log(S[1:] / S[:-1])
|
||||
# forward 21-day realized vol (annualized) — what the seller faces
|
||||
H = 21
|
||||
rv_fwd = np.full(len(days), np.nan)
|
||||
for t in range(len(days) - H):
|
||||
rv_fwd[t] = np.std(r[t + 1:t + 1 + H]) * math.sqrt(252)
|
||||
vrp = IV - rv_fwd # premium (positive = seller wins)
|
||||
# short-vol daily P&L proxy: collect implied variance, pay realized squared return
|
||||
iv_lag = np.concatenate([[np.nan], IV[:-1]])
|
||||
svol = (iv_lag ** 2) / 252.0 - r ** 2
|
||||
year = np.array([1970 + x / 365.25 for x in days]).astype(int)
|
||||
T_ = lambda x: torch.tensor(np.asarray(x)[np.isfinite(np.asarray(x))], device=DEV, dtype=torch.float64)
|
||||
|
||||
print(f"\n===== EQUITY-INDEX VRP (XSP, {days.min()}..{days.max()}, {len(days)} days) =====")
|
||||
print(f"mean implied vol {np.nanmean(IV):.1%} mean fwd-realized {np.nanmean(rv_fwd):.1%} "
|
||||
f"mean VRP {np.nanmean(vrp):.1%} (positive => premium exists)")
|
||||
print(f"VRP frac>0 (IV>RV): {np.nanmean(vrp[np.isfinite(vrp)]>0):.2f}")
|
||||
print(f"\nshort-vol daily P&L: Sharpe {sharpe_t(T_(svol)):+.2f} skew {float(((svol[np.isfinite(svol)]-np.nanmean(svol))**3).mean()/np.nanstd(svol)**3):+.2f} worst-day {np.nanmin(svol)/np.nanstd(svol):+.1f}σ")
|
||||
print("per-year short-vol: " + " ".join(f"{y}:{sharpe_t(T_(svol[year==y])):+.1f}" for y in range(2013, 2027) if (year == y).sum() > 60))
|
||||
|
||||
# tail-managed: don't sell when implied vol rising over 5d (the crypto-VRP gate)
|
||||
rising = np.zeros(len(days))
|
||||
for t in range(6, len(days)):
|
||||
rising[t] = 0.0 if IV[t - 1] > IV[t - 6] else 1.0
|
||||
svol_tm = rising * svol
|
||||
print(f"\ntail-managed (don't sell into rising IV): Sharpe {sharpe_t(T_(svol_tm)):+.2f} "
|
||||
f"skew {float(((svol_tm[np.isfinite(svol_tm)]-np.nanmean(svol_tm))**3).mean()/np.nanstd(svol_tm)**3):+.2f} worst {np.nanmin(svol_tm)/np.nanstd(svol_tm):+.1f}σ")
|
||||
print("per-year tail-managed: " + " ".join(f"{y}:{sharpe_t(T_(svol_tm[year==y])):+.1f}" for y in range(2013, 2027) if (year == y).sum() > 60))
|
||||
print("\nVERDICT: VRP frac>0 high + short-vol Sharpe>0 + tail-managed improves skew/recent = real, harvestable equity VRP.")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
87
scripts/surfer/equity_vrp2.py
Normal file
87
scripts/surfer/equity_vrp2.py
Normal file
@@ -0,0 +1,87 @@
|
||||
#!/usr/bin/env python3
|
||||
"""Equity-index VRP, CORRECTED: clean SPY underlying (not noisy parity) + XSP option-implied IV.
|
||||
|
||||
Fixes the bug where parity-derived underlying inflated realized vol (spurious -VRP). SPY ~= XSP
|
||||
(~SPX/10), so SPY is the clean underlying for ATM reference + realized vol. Overlap 2023-2026.
|
||||
VRP = implied (XSP ATM straddle) - realized (clean SPY). Short-vol Sharpe, per-year, tail, tail-managed.
|
||||
"""
|
||||
import math
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
from equity_vrp import load_xsp, to_epoch_day # noqa: E402
|
||||
from signal_sweep import sharpe_t # noqa: E402
|
||||
import torch # noqa: E402
|
||||
|
||||
DEV = "cuda" if torch.cuda.is_available() else "cpu"
|
||||
DAY_NS = 86_400 * 10**9
|
||||
|
||||
|
||||
def load_spy():
|
||||
import databento as db
|
||||
a = db.DBNStore.from_file("data/surfer/spy_ohlcv1d.dbn").to_ndarray()
|
||||
day = a["ts_event"].astype(np.int64) // DAY_NS
|
||||
close = a["close"].astype(np.float64) / 1e9
|
||||
return {int(day[i]): float(close[i]) for i in range(len(day)) if close[i] > 0}
|
||||
|
||||
|
||||
def main():
|
||||
spy = load_spy()
|
||||
print(f"SPY clean: {len(spy)} days ({min(spy)}..{max(spy)}), level {spy[min(spy)]:.0f}->{spy[max(spy)]:.0f}")
|
||||
d = load_xsp()
|
||||
d["exp_day"] = d["expiry"].map(to_epoch_day)
|
||||
d["ttm"] = d["exp_day"] - d["day"]
|
||||
d = d[(d["ttm"] >= 20) & (d["ttm"] <= 45)]
|
||||
d = d[d["day"].isin(spy.keys())] # overlap with clean SPY
|
||||
|
||||
iv = {}
|
||||
for day, g in d.groupby("day"):
|
||||
S = spy[int(day)]
|
||||
exp = g.iloc[(g["ttm"] - 30).abs().argsort()].iloc[0]["exp_day"]
|
||||
ge = g[g["exp_day"] == exp]
|
||||
calls = ge[ge["right"] == "C"].groupby("strike")["close"].mean()
|
||||
puts = ge[ge["right"] == "P"].groupby("strike")["close"].mean()
|
||||
common = calls.index.intersection(puts.index)
|
||||
if len(common) < 3:
|
||||
continue
|
||||
katm = common[np.abs(np.array(common) - S).argmin()] # ATM = strike nearest clean SPY level
|
||||
straddle = float(calls[katm] + puts[katm])
|
||||
T = float(ge["ttm"].iloc[0]) / 365.0
|
||||
if straddle <= 0 or T <= 0:
|
||||
continue
|
||||
iv[int(day)] = straddle / (0.8 * S * math.sqrt(T))
|
||||
|
||||
days = np.array(sorted(set(iv) & set(spy)))
|
||||
S = np.array([spy[x] for x in days]); IV = np.array([iv[x] for x in days])
|
||||
r = np.zeros(len(days)); r[1:] = np.log(S[1:] / S[:-1])
|
||||
H = 21
|
||||
rv_fwd = np.full(len(days), np.nan)
|
||||
for t in range(len(days) - H):
|
||||
rv_fwd[t] = np.std(r[t + 1:t + 1 + H]) * math.sqrt(252)
|
||||
vrp = IV - rv_fwd
|
||||
iv_lag = np.concatenate([[np.nan], IV[:-1]])
|
||||
svol = (iv_lag ** 2) / 252.0 - r ** 2
|
||||
year = np.array([1970 + x / 365.25 for x in days]).astype(int)
|
||||
T_ = lambda x: torch.tensor(np.asarray(x)[np.isfinite(np.asarray(x))], device=DEV, dtype=torch.float64)
|
||||
fin = np.isfinite(vrp)
|
||||
|
||||
print(f"\n===== EQUITY-INDEX VRP (CORRECTED, clean SPY) — {len(days)} days =====")
|
||||
print(f"mean implied {np.nanmean(IV):.1%} mean fwd-realized(clean) {np.nanmean(rv_fwd):.1%} "
|
||||
f"mean VRP {np.nanmean(vrp[fin]):+.1%} frac>0 {np.mean(vrp[fin]>0):.2f}")
|
||||
sk = float(((svol[np.isfinite(svol)] - np.nanmean(svol)) ** 3).mean() / np.nanstd(svol) ** 3)
|
||||
print(f"short-vol daily P&L: Sharpe {sharpe_t(T_(svol)):+.2f} skew {sk:+.2f} worst {np.nanmin(svol)/np.nanstd(svol):+.1f}σ")
|
||||
print("per-year: " + " ".join(f"{y}:{sharpe_t(T_(svol[year==y])):+.1f}" for y in range(2023, 2027) if (year == y).sum() > 40))
|
||||
rising = np.zeros(len(days))
|
||||
for t in range(6, len(days)):
|
||||
rising[t] = 0.0 if IV[t - 1] > IV[t - 6] else 1.0
|
||||
tm = rising * svol
|
||||
print(f"tail-managed (don't sell into rising IV): Sharpe {sharpe_t(T_(tm)):+.2f} worst {np.nanmin(tm)/np.nanstd(tm):+.1f}σ")
|
||||
print("per-year TM: " + " ".join(f"{y}:{sharpe_t(T_(tm[year==y])):+.1f}" for y in range(2023, 2027) if (year == y).sum() > 40))
|
||||
print("\nVERDICT: mean VRP>0 + frac>0~0.8 confirms premium exists; short-vol Sharpe (esp tail-managed) = harvestable.")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
78
scripts/surfer/equity_vrp3.py
Normal file
78
scripts/surfer/equity_vrp3.py
Normal file
@@ -0,0 +1,78 @@
|
||||
#!/usr/bin/env python3
|
||||
"""Equity-index VRP, robust internal-parity version (self-consistent, no external underlying).
|
||||
|
||||
Per day: implied forward F = MEDIAN over the 5 near-ATM strikes of (K + C - P) [put-call parity,
|
||||
median-smoothed to kill per-strike noise]. Use F consistently for ATM, implied vol, AND the
|
||||
realized-vol series. Internally consistent -> avoids both the single-strike-parity noise and the
|
||||
SPY!=XSP divergence that produced spurious negative VRP. XSP 2013-2026.
|
||||
"""
|
||||
import math
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
from equity_vrp import load_xsp, to_epoch_day # noqa: E402
|
||||
from signal_sweep import sharpe_t # noqa: E402
|
||||
import torch # noqa: E402
|
||||
|
||||
DEV = "cuda" if torch.cuda.is_available() else "cpu"
|
||||
|
||||
|
||||
def main():
|
||||
d = load_xsp()
|
||||
d["exp_day"] = d["expiry"].map(to_epoch_day)
|
||||
d["ttm"] = d["exp_day"] - d["day"]
|
||||
d = d[(d["ttm"] >= 20) & (d["ttm"] <= 45)]
|
||||
F, IVm = {}, {}
|
||||
for day, g in d.groupby("day"):
|
||||
exp = g.iloc[(g["ttm"] - 30).abs().argsort()].iloc[0]["exp_day"]
|
||||
ge = g[g["exp_day"] == exp]
|
||||
calls = ge[ge["right"] == "C"].groupby("strike")["close"].mean()
|
||||
puts = ge[ge["right"] == "P"].groupby("strike")["close"].mean()
|
||||
common = np.array(sorted(set(calls.index) & set(puts.index)))
|
||||
if len(common) < 5:
|
||||
continue
|
||||
rough = common[np.abs(calls[common].values - puts[common].values).argmin()] # min|C-P| ~ ATM
|
||||
near = common[np.argsort(np.abs(common - rough))[:5]] # 5 nearest strikes
|
||||
f = float(np.median([k + float(calls[k]) - float(puts[k]) for k in near])) # parity forward (median)
|
||||
katm = common[np.abs(common - f).argmin()]
|
||||
straddle = float(calls[katm] + puts[katm])
|
||||
T = float(ge["ttm"].iloc[0]) / 365.0
|
||||
if f <= 0 or straddle <= 0 or T <= 0:
|
||||
continue
|
||||
F[int(day)] = f; IVm[int(day)] = straddle / (0.8 * f * math.sqrt(T))
|
||||
|
||||
days = np.array(sorted(F))
|
||||
S = np.array([F[x] for x in days]); IV = np.array([IVm[x] for x in days])
|
||||
r = np.zeros(len(days)); r[1:] = np.log(S[1:] / S[:-1])
|
||||
r = np.clip(r, -0.25, 0.25) # guard residual parity glitches
|
||||
H = 21
|
||||
rv = np.full(len(days), np.nan)
|
||||
for t in range(len(days) - H):
|
||||
rv[t] = np.std(r[t + 1:t + 1 + H]) * math.sqrt(252)
|
||||
vrp = IV - rv; fin = np.isfinite(vrp)
|
||||
iv_lag = np.concatenate([[np.nan], IV[:-1]])
|
||||
svol = (iv_lag ** 2) / 252.0 - r ** 2
|
||||
year = np.array([1970 + x / 365.25 for x in days]).astype(int)
|
||||
T_ = lambda x: torch.tensor(np.asarray(x)[np.isfinite(np.asarray(x))], device=DEV, dtype=torch.float64)
|
||||
|
||||
print(f"\n===== EQUITY-INDEX VRP (robust internal parity) — {len(days)} days ({days.min()}..{days.max()}) =====")
|
||||
print(f"implied-forward level {S[0]:.0f}->{S[-1]:.0f} (should track SPX/10 ~ 400->650)")
|
||||
print(f"mean implied {np.nanmean(IV):.1%} mean fwd-realized {np.nanmean(rv):.1%} "
|
||||
f"mean VRP {np.nanmean(vrp[fin]):+.1%} frac>0 {np.mean(vrp[fin]>0):.2f}")
|
||||
sk = float(((svol[np.isfinite(svol)] - np.nanmean(svol)) ** 3).mean() / np.nanstd(svol) ** 3)
|
||||
print(f"short-vol daily P&L: Sharpe {sharpe_t(T_(svol)):+.2f} skew {sk:+.2f} worst {np.nanmin(svol)/np.nanstd(svol):+.1f}σ")
|
||||
print("per-year: " + " ".join(f"{y}:{sharpe_t(T_(svol[year==y])):+.1f}" for y in range(2013, 2027) if (year == y).sum() > 60))
|
||||
rising = np.zeros(len(days))
|
||||
for t in range(6, len(days)):
|
||||
rising[t] = 0.0 if IV[t - 1] > IV[t - 6] else 1.0
|
||||
tm = rising * svol
|
||||
print(f"tail-managed: Sharpe {sharpe_t(T_(tm)):+.2f} worst {np.nanmin(tm)/np.nanstd(tm):+.1f}σ")
|
||||
print("per-year TM: " + " ".join(f"{y}:{sharpe_t(T_(tm[year==y])):+.1f}" for y in range(2013, 2027) if (year == y).sum() > 60))
|
||||
print("\nVERDICT: implied-forward tracking ~SPX/10 + mean RV ~13-16% + VRP>0 => measurement clean & premium real.")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
49
scripts/surfer/fetch_xsp.py
Normal file
49
scripts/surfer/fetch_xsp.py
Normal file
@@ -0,0 +1,49 @@
|
||||
#!/usr/bin/env python3
|
||||
"""Fetch XSP (mini-SPX) options daily OHLCV for the equity-index VRP test. HARD-CAPPED.
|
||||
|
||||
get_cost-gated: prints the exact cost and ABORTS if above CAP — no surprise charges.
|
||||
Single dataset/schema/symbol/window. ohlcv-1d only (NOT definition).
|
||||
"""
|
||||
import os
|
||||
import sys
|
||||
|
||||
import databento as db
|
||||
|
||||
CAP = 70.00 # hard cap — abort if cost exceeds this
|
||||
DS, SCH, SYM = "OPRA.PILLAR", "ohlcv-1d", "XSP.OPT"
|
||||
START, END = "2013-04-01", "2026-06-05"
|
||||
OUT = "data/surfer/xsp_ohlcv1d.dbn"
|
||||
|
||||
|
||||
def main():
|
||||
c = db.Historical(os.environ["DATABENTO_API_KEY"])
|
||||
if os.path.exists(OUT):
|
||||
print(f"already downloaded: {OUT}"); return 0
|
||||
cost = c.metadata.get_cost(dataset=DS, symbols=[SYM], schema=SCH, start=START, end=END, stype_in="parent")
|
||||
gb = c.metadata.get_billable_size(dataset=DS, symbols=[SYM], schema=SCH, start=START, end=END, stype_in="parent") / 1e9
|
||||
print(f"EXACT get_cost = ${cost:.2f} ({gb:.2f} GB) CAP = ${CAP:.2f}")
|
||||
if cost > CAP:
|
||||
print(f"ABORT: ${cost:.2f} exceeds cap ${CAP:.2f} — NOTHING downloaded.")
|
||||
return 1
|
||||
print(f"under cap by ${CAP - cost:.2f} -> downloading {SYM} {SCH} year-chunked (same total cost)")
|
||||
os.makedirs("data/surfer/xsp", exist_ok=True)
|
||||
for yr in range(2013, 2027):
|
||||
ys = f"{yr}-01-01" if yr > 2013 else "2013-04-01"
|
||||
ye = f"{yr+1}-01-01" if yr < 2026 else "2026-06-05"
|
||||
out = f"data/surfer/xsp/xsp_{yr}.dbn"
|
||||
if os.path.exists(out):
|
||||
print(f" {yr}: exists"); continue
|
||||
for attempt in range(3):
|
||||
try:
|
||||
d = c.timeseries.get_range(dataset=DS, symbols=[SYM], schema=SCH, start=ys, end=ye, stype_in="parent")
|
||||
d.to_file(out)
|
||||
print(f" {yr}: saved ({os.path.getsize(out)/1e6:.1f} MB)")
|
||||
break
|
||||
except Exception as e:
|
||||
print(f" {yr}: attempt {attempt+1} {type(e).__name__}; retry")
|
||||
print(f"DONE -> data/surfer/xsp/")
|
||||
return 0
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
sys.exit(main())
|
||||
Reference in New Issue
Block a user