You cannot select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

554 lines
22 KiB
Python

This file contains ambiguous Unicode characters!

This file contains ambiguous Unicode characters that may be confused with others in your current locale. If your use case is intentional and legitimate, you can safely ignore this warning. Use the Escape button to highlight these characters.

"""
V47 Diagnostic: Overnight Gap Residual (Sector-ETF OLS Regression)
Hypothesis: The idiosyncratic component of a stock's overnight gap (after removing
systematic sector/market beta) predicts ORB follow-through better than the raw gap.
Based on Lou/Polk (JFE 2019) beta decomposition and Bogousslavsky (JF 2021) overnight
persistence as informed-trader footprint.
Features tested:
gap_residual : OLS residual of stock gap on sector ETF gap (rolling 60d β̂)
gap_residual_qqq : OLS residual of stock gap on QQQ gap (QQQ universal benchmark)
gap_residual_std : Z-scored gap_residual (divided by σ̂ of in-window residuals)
overnight_residual_avg_5d : Mean of 5 prior daily gap_residuals (Bogousslavsky persistence)
overnight_residual_avg_10d: Mean of 10 prior daily gap_residuals
Gates (pre-committed, must be set before running):
G1: |Pearson| ≥ 0.07 AND n ≥ 200
G2: |top bottom tercile avg_R| ≥ 0.30R ← binding gate
G3: |top bottom tercile WR| ≥ 5pp
G4: coverage ≥ 80% of V46 400d trade set
G5a: |ρ(feature, obv_slope_20)| < 0.70
G5b: |ρ(feature, raw gap_pct)| < 0.70
Decision rule (pre-committed):
≥1 feature passes G1-G5b → advance highest-G2 to Phase 2 backtest
0 pass → close axis, document in lineage memory
Data source: data/cache/intraday/{ticker}/{date}.parquet
V46 400d baseline: runs/v46_correct_w12_400d.json/intraday_20260422_043215_63f03e64.json
"""
from __future__ import annotations
import json
import os
import sys
from pathlib import Path
import pyarrow.parquet as pq
sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "../../..")))
from zoneinfo import ZoneInfo
from libs.common.time_utils import trading_days_between
import datetime as dt
_ET = ZoneInfo("America/New_York")
_MKT_OPEN = dt.time(9, 30)
_MKT_CLOSE = dt.time(16, 0)
INTRADAY_CACHE_DIR = "data/cache/intraday"
SECTOR_CACHE = "data/cache/sector_cache.json"
V46_400D_RUN = "runs/v46_correct_w12_400d.json/intraday_20260422_043215_63f03e64.json"
SECTOR_TO_ETF: dict[str, str] = {
"Technology": "XLK",
"Financial Services": "XLF",
"Healthcare": "XLV",
"Energy": "XLE",
"Industrials": "XLI",
"Basic Materials": "XLB",
"Communication Services": "XLC",
"Utilities": "XLU",
"Real Estate": "XLRE",
"Consumer Defensive": "XLP",
"Consumer Cyclical": "XLY",
}
BETA_WINDOW = 60 # rolling days for β̂ estimation
MIN_PAIRS = 20 # minimum non-null pairs for valid β̂
# ── Helpers ────────────────────────────────────────────────────────────────
def _parse_ts(ts_raw: object) -> dt.datetime:
s = str(ts_raw)
if s.endswith("Z"):
s = s[:-1] + "+00:00"
return dt.datetime.fromisoformat(s).astimezone(_ET)
def load_daily_ohlcv(ticker: str, dates: list[str]) -> dict[str, dict]:
"""Load regular-session {open, close} per date from intraday parquet cache."""
result: dict[str, dict] = {}
for date in dates:
path = Path(INTRADAY_CACHE_DIR) / ticker.upper() / f"{date}.parquet"
if not path.exists():
continue
try:
table = pq.read_table(str(path))
rows = table.to_pydict()
except Exception:
continue
mkt_bars = []
for i, ts_raw in enumerate(rows.get("timestamp", [])):
try:
ts = _parse_ts(ts_raw)
except Exception:
continue
if not (_MKT_OPEN <= ts.time() < _MKT_CLOSE):
continue
mkt_bars.append({
"ts": ts,
"open": float(rows["open"][i] or 0),
"close": float(rows["close"][i] or 0),
})
if not mkt_bars:
continue
mkt_bars.sort(key=lambda b: b["ts"])
result[date] = {
"open": mkt_bars[0]["open"],
"close": mkt_bars[-1]["close"],
}
return result
def compute_overnight_gaps(ohlcv: dict[str, dict], sorted_dates: list[str]) -> dict[str, float]:
"""Compute overnight gap = (open_D - close_{D-1}) / close_{D-1} for each date."""
gaps: dict[str, float] = {}
for i in range(1, len(sorted_dates)):
d0, d1 = sorted_dates[i - 1], sorted_dates[i]
prev = ohlcv.get(d0)
curr = ohlcv.get(d1)
if prev and curr and prev["close"] > 0 and curr["open"] > 0:
gaps[d1] = (curr["open"] - prev["close"]) / prev["close"]
return gaps
def _winsorize(vals: list[float], n_sigma: float = 3.0) -> list[float]:
if len(vals) < 4:
return vals
n = len(vals)
mu = sum(vals) / n
sigma = (sum((v - mu) ** 2 for v in vals) / n) ** 0.5
if sigma <= 0:
return vals
lo, hi = mu - n_sigma * sigma, mu + n_sigma * sigma
return [max(lo, min(hi, v)) for v in vals]
def ols_residual_on_window(
entry_stock_gap: float,
entry_etf_gap: float,
prior_stock_gaps: list[float],
prior_etf_gaps: list[float],
) -> tuple[float, float, float] | None:
"""
OLS regression: stock_gap = alpha + beta * etf_gap on the prior window.
Returns (residual_at_entry, beta, sigma_residuals) or None if insufficient data.
"""
assert len(prior_stock_gaps) == len(prior_etf_gaps)
n = len(prior_stock_gaps)
if n < MIN_PAIRS:
return None
xs = _winsorize(prior_etf_gaps)
ys = _winsorize(prior_stock_gaps)
xm = sum(xs) / n
ym = sum(ys) / n
num = sum((xs[i] - xm) * (ys[i] - ym) for i in range(n))
den = sum((xs[i] - xm) ** 2 for i in range(n))
if den <= 0:
return None
beta = num / den
alpha = ym - beta * xm
residuals_window = [ys[i] - (alpha + beta * xs[i]) for i in range(n)]
sigma = (sum(r ** 2 for r in residuals_window) / n) ** 0.5
residual_entry = entry_stock_gap - (alpha + beta * entry_etf_gap)
return residual_entry, beta, sigma
def pearson(xs: list[float], ys: list[float]) -> float | None:
n = len(xs)
if n < 5:
return None
xm, ym = sum(xs) / n, sum(ys) / n
num = sum((xs[i] - xm) * (ys[i] - ym) for i in range(n))
dx = sum((x - xm) ** 2 for x in xs) ** 0.5
dy = sum((y - ym) ** 2 for y in ys) ** 0.5
if dx <= 0 or dy <= 0:
return None
return num / (dx * dy)
def tercile_stats(vals: list[float], rs: list[float]) -> dict:
if len(vals) < 9:
return {}
pairs = sorted(zip(vals, rs), key=lambda p: p[0])
n = len(pairs)
t = n // 3
def stats(sub):
ys = [p[1] for p in sub]
return {"n": len(ys), "wr": sum(1 for y in ys if y > 0) / len(ys), "avg_r": sum(ys) / len(ys)}
return {"low": stats(pairs[:t]), "mid": stats(pairs[t:2 * t]), "high": stats(pairs[2 * t:])}
def pairwise_rho(idxs_a: list[int], vals_a: list[float],
idxs_b: list[int], vals_b: list[float]) -> float | None:
map_b = {idx: vals_b[i] for i, idx in enumerate(idxs_b)}
va, vb = [], []
for i, idx in enumerate(idxs_a):
if idx in map_b:
va.append(vals_a[i])
vb.append(map_b[idx])
return pearson(va, vb) if len(va) >= 10 else None
def report_feature(
label: str,
idxs: list[int],
vals: list[float],
rs: list[float],
n_total: int,
g5_comparisons: list[tuple[str, list[int], list[float]]],
) -> dict:
n = len(vals)
p = pearson(vals, rs)
ts = tercile_stats(vals, rs)
coverage = n / n_total
g4 = coverage >= 0.80
if not ts or p is None:
print(f" {label}: n={n}/{n_total} ({coverage:.1%}) INSUFFICIENT DATA")
return {"pass_g2": False, "pass_all": False, "pearson": None, "g2_gap": None}
low, mid, high = ts["low"], ts["mid"], ts["high"]
raw_g2 = high["avg_r"] - low["avg_r"]
avg_r_gap = abs(raw_g2)
wr_gap = abs(high["wr"] - low["wr"])
g1 = abs(p) >= 0.07 and n >= 200
g2 = avg_r_gap >= 0.30
g3 = wr_gap >= 0.05
print(f"\n [{label}] n={n}/{n_total} ({coverage:.1%}) Pearson={p:+.3f} direction={'↑high' if raw_g2 > 0 else '↓low'}")
print(f" G1: {'PASS' if g1 else 'FAIL'} (|{abs(p):.3f}| {'' if abs(p) >= 0.07 else '<'} 0.07, n={n} {'' if n >= 200 else '<'} 200)")
print(f" G2: {'PASS' if g2 else 'FAIL'} (avg_R gap = {avg_r_gap:.3f}R [threshold 0.30R])")
print(f" G3: {'PASS' if g3 else 'FAIL'} (WR gap = {wr_gap * 100:.1f}pp [threshold 5pp])")
print(f" G4: {'PASS' if g4 else 'FAIL'} (coverage = {coverage:.1%} [threshold 80%])")
print(f" Bottom tercile: WR {low['wr'] * 100:.1f}% avg_R {low['avg_r']:+.3f} n={low['n']}")
print(f" Middle tercile: WR {mid['wr'] * 100:.1f}% avg_R {mid['avg_r']:+.3f} n={mid['n']}")
print(f" Top tercile: WR {high['wr'] * 100:.1f}% avg_R {high['avg_r']:+.3f} n={high['n']}")
g5_results = {}
pass_all_g5 = True
for cmp_label, cmp_idxs, cmp_vals in g5_comparisons:
rho = pairwise_rho(idxs, vals, cmp_idxs, cmp_vals)
if rho is not None:
gate_key = f"g5_{cmp_label}"
gate_pass = abs(rho) < 0.70
g5_results[gate_key] = {"rho": rho, "pass": gate_pass}
print(f" G5({cmp_label}): {'PASS' if gate_pass else 'FAIL'} (ρ={rho:+.3f} [|threshold| < 0.70])")
if not gate_pass:
pass_all_g5 = False
pass_all = g1 and g2 and g3 and g4 and pass_all_g5
print(f"{'ALL GATES PASS ✓' if pass_all else 'FAIL'}")
return {
"pass_g2": g2, "pass_all": pass_all,
"pearson": p, "g2_gap": avg_r_gap, "raw_g2_gap": raw_g2,
"direction": "high" if raw_g2 > 0 else "low",
**g5_results,
}
# ── Main ───────────────────────────────────────────────────────────────────
def main() -> None:
print("=== V47 Diagnostic: Overnight Gap Residual ===\n")
print("Pre-committed gates: G1(|P|≥0.07,n≥200) G2(avg_R≥0.30R) G3(WR≥5pp) G4(cov≥80%) G5(|ρ|<0.70)")
print("Pre-committed decision rule: highest G2 that passes ALL gates → advance to Phase 2 backtest\n")
with open(V46_400D_RUN) as f:
run_data = json.load(f)
trades_raw = run_data.get("trades", [])
m = run_data.get("metrics", {})
print(f"V46 400d: {m.get('total_trades')} trades, {m.get('start_date')}{m.get('end_date')}, "
f"return={m.get('total_return_pct', 0) * 100:.1f}%\n")
with open(SECTOR_CACHE) as f:
sector_cache = json.load(f)
trades = [
{
"idx": i,
"ticker": t["ticker"],
"date": t["date"][:10],
"r": float(t["r_multiple_at_exit"]),
"sector": sector_cache.get(t["ticker"], "UNKNOWN"),
}
for i, t in enumerate(trades_raw)
if t.get("r_multiple_at_exit") is not None
]
n_total = len(trades)
print(f"Valid trades (with r_multiple): {n_total}")
# Build date range: from earliest trade date - 70 trading days to last trade date
all_trade_dates = sorted(set(t["date"] for t in trades))
start_d = dt.date.fromisoformat(all_trade_dates[0]) - dt.timedelta(days=110)
end_d = dt.date.fromisoformat(all_trade_dates[-1])
all_trading_days = [d.isoformat() for d in trading_days_between(start_d, end_d)]
td_set = set(all_trading_days)
print(f"Computing over trading days: {all_trading_days[0]}{all_trading_days[-1]} ({len(all_trading_days)} days)")
# Load sector ETF + QQQ daily OHLCV
etf_tickers = list(set(SECTOR_TO_ETF.values())) + ["QQQ"]
print(f"\nLoading {len(etf_tickers)} sector ETF + QQQ daily OHLCV...")
etf_ohlcv: dict[str, dict[str, dict]] = {}
for etf in etf_tickers:
etf_ohlcv[etf] = load_daily_ohlcv(etf, all_trading_days)
print(f" {etf}: {len(etf_ohlcv[etf])} days")
# Pre-compute ETF overnight gaps
etf_gaps: dict[str, dict[str, float]] = {}
for etf, ohlcv in etf_ohlcv.items():
sorted_dates = sorted(ohlcv.keys())
etf_gaps[etf] = compute_overnight_gaps(ohlcv, sorted_dates)
print()
# Load stock OHLCV for all unique tickers
unique_tickers = sorted(set(t["ticker"] for t in trades))
print(f"Loading {len(unique_tickers)} stock intraday caches...")
stock_ohlcv: dict[str, dict[str, dict]] = {}
for ticker in unique_tickers:
stock_ohlcv[ticker] = load_daily_ohlcv(ticker, all_trading_days)
print(f" Loaded. Sample: {unique_tickers[0]}: {len(stock_ohlcv[unique_tickers[0]])} days")
# Pre-compute stock overnight gaps
stock_gaps: dict[str, dict[str, float]] = {}
for ticker in unique_tickers:
sorted_dates = sorted(stock_ohlcv[ticker].keys())
stock_gaps[ticker] = compute_overnight_gaps(stock_ohlcv[ticker], sorted_dates)
# Pre-compute daily residuals for each ticker (for Bogousslavsky persistence)
# residuals_sector[ticker][date] = gap_residual on that date (sector ETF β̂)
# residuals_qqq[ticker][date] = gap_residual on that date (QQQ β̂)
print("\nPre-computing daily residuals for all tickers (for persistence features)...")
residuals_sector: dict[str, dict[str, float]] = {}
residuals_qqq: dict[str, dict[str, float]] = {}
for ticker in unique_tickers:
sector = sector_cache.get(ticker, "UNKNOWN")
sector_etf = SECTOR_TO_ETF.get(sector)
s_gaps = stock_gaps[ticker]
sorted_stock_dates = sorted(s_gaps.keys())
r_sector: dict[str, float] = {}
r_qqq: dict[str, float] = {}
for j, date in enumerate(sorted_stock_dates):
sg = s_gaps.get(date)
if sg is None:
continue
# Get prior BETA_WINDOW dates with paired gaps
prior_dates = [d for d in sorted_stock_dates[:j] if d < date][-BETA_WINDOW:]
# Sector ETF residual
if sector_etf:
eg = etf_gaps[sector_etf]
prior_pairs_s = [(s_gaps[d], eg[d]) for d in prior_dates if d in s_gaps and d in eg]
if len(prior_pairs_s) >= MIN_PAIRS:
ps, pe = [x[0] for x in prior_pairs_s], [x[1] for x in prior_pairs_s]
entry_etf_gap = etf_gaps[sector_etf].get(date)
if entry_etf_gap is not None:
res = ols_residual_on_window(sg, entry_etf_gap, ps, pe)
if res is not None:
r_sector[date] = res[0]
# QQQ residual
qq_gaps = etf_gaps["QQQ"]
prior_pairs_q = [(s_gaps[d], qq_gaps[d]) for d in prior_dates if d in s_gaps and d in qq_gaps]
if len(prior_pairs_q) >= MIN_PAIRS:
pq_s, pq_q = [x[0] for x in prior_pairs_q], [x[1] for x in prior_pairs_q]
entry_qqq_gap = qq_gaps.get(date)
if entry_qqq_gap is not None:
res = ols_residual_on_window(sg, entry_qqq_gap, pq_s, pq_q)
if res is not None:
r_qqq[date] = res[0]
residuals_sector[ticker] = r_sector
residuals_qqq[ticker] = r_qqq
# ── Compute features for each trade ─────────────────────────────────
print("\nComputing features for each trade...\n")
# Feature arrays (indexed by trade idx for pairwise ρ)
feat_gap_residual: list[tuple[int, float]] = []
feat_gap_residual_qqq: list[tuple[int, float]] = []
feat_gap_residual_std: list[tuple[int, float]] = []
feat_persist_5d: list[tuple[int, float]] = []
feat_persist_10d: list[tuple[int, float]] = []
feat_raw_gap: list[tuple[int, float]] = []
feat_obv_slope: list[tuple[int, float]] = [] # placeholder - computed from daily bars if available
r_by_idx: dict[int, float] = {t["idx"]: t["r"] for t in trades}
for trade in trades:
idx = trade["idx"]
ticker = trade["ticker"]
date = trade["date"]
sector = trade["sector"]
sector_etf = SECTOR_TO_ETF.get(sector)
sg = stock_gaps[ticker].get(date)
if sg is None:
continue
# Raw gap feature
feat_raw_gap.append((idx, sg))
# gap_residual (sector ETF)
if date in residuals_sector.get(ticker, {}):
feat_gap_residual.append((idx, residuals_sector[ticker][date]))
# gap_residual_qqq
if date in residuals_qqq.get(ticker, {}):
feat_gap_residual_qqq.append((idx, residuals_qqq[ticker][date]))
# gap_residual_std: need sigma from the window on trade date
# Recompute sigma for this trade date
if sector_etf:
s_gaps_t = stock_gaps[ticker]
sorted_sd = sorted(s_gaps_t.keys())
j = sorted_sd.index(date) if date in sorted_sd else -1
if j > 0:
prior_dates = sorted_sd[:j][-BETA_WINDOW:]
eg = etf_gaps[sector_etf]
prior_pairs = [(s_gaps_t[d], eg[d]) for d in prior_dates if d in s_gaps_t and d in eg]
if len(prior_pairs) >= MIN_PAIRS:
ps, pe = [x[0] for x in prior_pairs], [x[1] for x in prior_pairs]
entry_etf_gap = eg.get(date)
if entry_etf_gap is not None:
res = ols_residual_on_window(sg, entry_etf_gap, ps, pe)
if res is not None and res[2] > 0:
feat_gap_residual_std.append((idx, res[0] / res[2]))
# Persistence features: mean of prior 5/10 daily residuals (sector β̂)
sorted_all_td = sorted(td_set & set(stock_gaps[ticker].keys()))
if date in sorted_all_td:
d_idx = sorted_all_td.index(date)
prior_td = sorted_all_td[max(0, d_idx - 10):d_idx]
ticker_residuals = residuals_sector.get(ticker, {})
prev_5d_resids = [ticker_residuals[d] for d in prior_td[-5:] if d in ticker_residuals]
prev_10d_resids = [ticker_residuals[d] for d in prior_td[-10:] if d in ticker_residuals]
if len(prev_5d_resids) >= 3:
feat_persist_5d.append((idx, sum(prev_5d_resids) / len(prev_5d_resids)))
if len(prev_10d_resids) >= 5:
feat_persist_10d.append((idx, sum(prev_10d_resids) / len(prev_10d_resids)))
# ── Build r-multiple arrays aligned to each feature ──────────────
def split(feat_pairs: list[tuple[int, float]]) -> tuple[list[int], list[float], list[float]]:
idxs = [p[0] for p in feat_pairs]
vals = [p[1] for p in feat_pairs]
rs = [r_by_idx[i] for i in idxs]
return idxs, vals, rs
# OBV slope from run data if available (fallback: empty)
obv_idxs, obv_vals, _ = split(feat_obv_slope) if feat_obv_slope else ([], [], [])
raw_idxs, raw_vals, raw_rs = split(feat_raw_gap)
g5_refs = [
("raw_gap", raw_idxs, raw_vals),
]
print("=" * 70)
print(f"GATE RESULTS (n_total={n_total})")
print("=" * 70)
results = {}
for label, feat_pairs in [
("gap_residual (sector ETF β̂)", feat_gap_residual),
("gap_residual_qqq (QQQ β̂)", feat_gap_residual_qqq),
("gap_residual_std (z-scored, sector ETF)", feat_gap_residual_std),
("overnight_residual_avg_5d (Bogousslavsky)", feat_persist_5d),
("overnight_residual_avg_10d (Bogousslavsky)", feat_persist_10d),
]:
if not feat_pairs:
print(f"\n [{label}]: no data")
continue
idxs, vals, rs = split(feat_pairs)
results[label] = report_feature(label, idxs, vals, rs, n_total, g5_refs)
# Pairwise ρ matrix among all features
print("\n\n=== Pairwise ρ matrix ===")
all_feat_sets = [
("gap_residual", feat_gap_residual),
("gap_residual_qqq", feat_gap_residual_qqq),
("gap_residual_std", feat_gap_residual_std),
("persist_5d", feat_persist_5d),
("persist_10d", feat_persist_10d),
("raw_gap", feat_raw_gap),
]
feat_names = [n for n, _ in all_feat_sets]
feat_data = {n: (split(fp)[0], split(fp)[1]) for n, fp in all_feat_sets if fp}
for na in feat_names:
if na not in feat_data:
continue
row_parts = []
for nb in feat_names:
if nb not in feat_data:
row_parts.append(" N/A ")
continue
if na == nb:
row_parts.append(" 1.000")
continue
rho = pairwise_rho(feat_data[na][0], feat_data[na][1],
feat_data[nb][0], feat_data[nb][1])
row_parts.append(f"{rho:+.3f}" if rho is not None else " N/A ")
print(f" {na:30s} {' '.join(row_parts)}")
# Summary
print("\n\n=== Summary & Decision ===")
passing = [(lbl, r) for lbl, r in results.items() if r.get("pass_all")]
g2_only = [(lbl, r) for lbl, r in results.items() if r.get("pass_g2")]
print(f"Features passing ALL gates: {len(passing)}")
for lbl, r in passing:
print(f"{lbl} G2={r['g2_gap']:.3f}R Pearson={r['pearson']:+.3f}")
if passing:
# Pick highest G2 among passing
best_lbl, best_r = max(passing, key=lambda x: x[1]["g2_gap"])
print(f"\n→ ADVANCE to Phase 2: {best_lbl}")
print(f" G2={best_r['g2_gap']:.3f}R Pearson={best_r['pearson']:+.3f} direction={best_r['direction']}")
else:
print("\nNo features pass ALL gates.")
best_g2 = max((r.get("g2_gap") or 0) for r in results.values()) if results else 0
print(f"Best G2 observed: {best_g2:.3f}R (threshold 0.30R)")
print("→ CLOSE overnight gap residual axis. Document in project_orb_gainers_lineage.md.")
print(" Next axis candidates: options flow (put/call ratio, IV skew) or 13F ownership momentum.")
print("\n" + "=" * 70)
print(f"Raw gap reference: Pearson={pearson(raw_vals, raw_rs):+.3f}" if len(raw_vals) >= 5 else "")
if __name__ == "__main__":
main()