"""GEX/DEX aggregation + key levels.""" import logging from datetime import date, datetime, time, timedelta from typing import Dict, List, Optional, Tuple from zoneinfo import ZoneInfo import numpy as np import pandas as pd try: import pandas_market_calendars as mcal _NYSE = mcal.get_calendar("NYSE") except Exception: # pragma: no cover _NYSE = None from .config import GexConfig from .greeks import bs_gamma_vec logger = logging.getLogger(__name__) _ET = ZoneInfo("America/New_York") # trading session length (09:30-16:00 ET) _SESSION = timedelta(hours=6.5) def business_days_to_expiry(expiry: date, today: date) -> int: """Count NYSE business days from today (inclusive) to expiry. FIX 48: a PAST expiry returns -1 (not 0). Returning 0 let an already-expired contract look like a 0DTE contract and survive the DTE<0 filter, after which its floored T made Black-Scholes gamma explode as spot approached its strike. """ if expiry < today: return -1 if _NYSE is not None: try: start = pd.Timestamp(today) end = pd.Timestamp(expiry) sched = _NYSE.schedule(start_date=start, end_date=end) return max(len(sched) - 1, 0) # exclude today itself except Exception: pass # fallback: naive weekday count days = np.busday_count(today, expiry) return max(int(days), 0) def _settlement_time(am_settled: bool, instrument_class: str = "equity_etf") -> time: """FIX 84: the settlement clock depends on instrument class. - AM-settled index expiries (third-Friday): 09:30 ET (one day earlier). - PM-settled INDEX options (instrument_class == "index", not AM-settled): 16:15 ET. Index options settle at the close + 15 min, NOT 16:00. - equity/ETF options: 16:00 ET. Using 16:00 for a PM-settled index understates time-to-settlement by 15 min — for a 0DTE NDX captured at 15:54 that is 6 min vs the true 21 min, which pushes gamma ∝ 1/sqrt(T) ~1.9x too high and destabilises reconciliation. """ if am_settled: return time(9, 30) if instrument_class == "index": return time(16, 15) return time(16, 0) def time_to_expiry_years(expiry: date, now_et: datetime, cfg: GexConfig, am_settled: bool = False, instrument_class: str = "equity_etf") -> float: """FIX 65: calendar-time T for ALL expiries. T = (minutes from snapshot to settlement) / (365 * 24 * 60) Settlement clock is instrument-class-aware (FIX 84): 16:15 ET for PM-settled index options, 09:30 ET (one day earlier) for AM-settled index expiries, 16:00 ET for equity/ETF. Replaces the business-day path entirely. Cboe's reported gamma reflects actual calendar time remaining, not trading sessions. The old business-day path (full_days/252) diverged from calendar time by sqrt(252/365 * 7/5) ≈ 1.28 in gamma for multi-day expiries — the dominant residual error after FIX 62 fixed the 0DTE branch. SMH 2026-07-31: T_business=0.0190 vs T_calendar=0.0115, ratio 1.64, sqrt=1.28 — matching the observed reported/recomputed=1.235. `now_et` must be the SNAPSHOT timestamp converted to ET (not wall-clock now) so that --from-cache reproduces exactly. FIX 45: for AM-settled index options (third-Friday expiries), settlement is at the Thursday CLOSE, not Friday 16:00 — so the effective expiry is one day earlier. Pass am_settled=True to apply that adjustment. """ eff_expiry = expiry - timedelta(days=1) if am_settled else expiry settle = _settlement_time(am_settled, instrument_class) settle_dt = datetime.combine(eff_expiry, settle, tzinfo=now_et.tzinfo) minutes = max((settle_dt - now_et).total_seconds() / 60.0, 0.0) T = minutes / (365.0 * 24.0 * 60.0) return max(T, 1.0 / (252.0 * 13.0)) # floor ≈ 2.67 calendar hours (gamma cap) def minutes_to_settlement(expiry: date, now_et: datetime, am_settled: bool = False, instrument_class: str = "equity_etf") -> float: """FIX 80: RAW (unfloored) minutes from the snapshot to settlement for one expiry. Mirrors time_to_expiry_years exactly (same instrument-class-aware settlement clock, FIX 84) but returns the raw minute count instead of the floored year fraction. T is floored at ~160 min so it cannot be reversed to recover the true minutes; reconciliation's near-settlement guard needs the real value to tell a 5-min expiry from a 160-min one. """ eff_expiry = expiry - timedelta(days=1) if am_settled else expiry settle = _settlement_time(am_settled, instrument_class) settle_dt = datetime.combine(eff_expiry, settle, tzinfo=now_et.tzinfo) return max((settle_dt - now_et).total_seconds() / 60.0, 0.0) def gamma_precision(df: pd.DataFrame) -> Tuple[int, bool, str]: """FIX 66 / FIX 70: detect Cboe's gamma publication precision for this chain. Cboe publishes gamma to 4 decimal places. For underlyings priced in the tens of thousands (NDX ~28 000), the true ATM gamma is ~0.00014 — one significant figure at 4dp. Consecutive strikes all report exactly 0.0001, making the bars coarse and the reconciliation gap a DATA limit, not a code bug. Returns (sig_figs, low_precision, label): sig_figs: significant figures in the median non-zero reported gamma. low_precision: True when sig_figs == 1 (coarse bars; flip bars_reliable). label: FIX 70 three-state — "high" (3+ digits) | "adequate" (2) | "coarse" (1). The coarse-bars warning shows only for "coarse"; "adequate" (2 sig figs, ~1.3% granularity) is corroborated by SMH's 0.99% reconciliation and does not warrant flipping bars_reliable. """ gammas = df["gamma"].dropna() gammas = gammas[gammas > 0] if gammas.empty: return 4, False, "high" med = float(gammas.median()) if med <= 0: return 4, False, "high" import math # significant figures: count digits from first non-zero digit # e.g. 0.0001 -> 1 sig fig; 0.00014 -> 2; 0.00123 -> 3 exponent = math.floor(math.log10(med)) # round to 4dp (Cboe's publication precision) and count non-zero digits rounded = round(med, 4) if rounded == 0: return 1, True, "coarse" s = f"{rounded:.10f}".rstrip("0") # strip leading "0." and count remaining digits digits = s.split(".")[-1].lstrip("0") if "." in s else s.lstrip("0") sig_figs = max(len(digits), 1) if sig_figs >= 3: label = "high" elif sig_figs == 2: label = "adequate" else: label = "coarse" return sig_figs, sig_figs <= 1, label def filter_contracts_full(contracts: List[Dict], spot: float, cfg: GexConfig, now_et, instrument_class: str = "equity_etf") -> pd.DataFrame: """Apply §2 rules 1-4 ONLY (no strike band). Returns the FULL chain frame. This frame feeds BOTH profile curves — the far OTM open interest shapes the profile wings and creates the zero-crossing, so it must not be truncated. `now_et` may be a datetime (preferred — the snapshot timestamp in ET, used for the continuous time-to-expiry) or a bare date (back-compat; treated as 15:35 ET). FIX 45: for instrument_class "index", third-Friday expiries are AM-settled (settlement at the Thursday close), so their T ends one day earlier. FIX 48: compare RAW CALENDAR DATES before any business-day maths — a contract whose expiry_date < snapshot_date is already expired and must be dropped. The count is stored in df.attrs["expired_contracts_dropped"]. """ if isinstance(now_et, datetime): today = now_et.date() else: # bare date — default to a mid-afternoon snapshot time today = now_et now_et = datetime.combine(today, time(15, 35), tzinfo=_ET) rows = [] expired_dropped = 0 for c in contracts: # FIX 48: raw calendar date guard BEFORE any business-day maths if c["expiry"] < today: expired_dropped += 1 continue dte = (c["expiry"] - today).days if dte < 0: # rule 1 (redundant with FIX 48 guard, kept for safety) continue if dte > cfg.dte_max: # rule 2 continue if c["oi"] <= 0: # rule 3 continue if c["iv"] <= 0 or c["gamma"] == 0: # rule 4 continue am_settled = (instrument_class == "index" and is_third_friday(c["expiry"])) T = time_to_expiry_years(c["expiry"], now_et, cfg, am_settled=am_settled, instrument_class=instrument_class) rows.append({**c, "dte": dte, "T": T}) df = pd.DataFrame(rows) df.attrs["expired_contracts_dropped"] = expired_dropped if expired_dropped > 0: logger.warning("FIX 48: dropped %d expired contracts (expiry < %s) for chain.", expired_dropped, today) return df def filter_contracts_band(df_full: pd.DataFrame, spot: float, cfg: GexConfig) -> pd.DataFrame: """Apply §2 rule 5 (±12% strike band) to the full frame — for the BARS only.""" if df_full.empty: return df_full lo = spot * (1 - cfg.strike_band) hi = spot * (1 + cfg.strike_band) return df_full[(df_full["strike"] >= lo) & (df_full["strike"] <= hi)].copy() def filter_contracts(contracts: List[Dict], spot: float, cfg: GexConfig, today: date) -> pd.DataFrame: """Back-compat convenience: full filter then band (rules 1-5). Bars use this.""" return filter_contracts_band(filter_contracts_full(contracts, spot, cfg, today), spot, cfg) def aggregate(df: pd.DataFrame, spot: float, cfg: GexConfig) -> pd.DataFrame: """ Per-strike GEX/DEX aggregation across all expirations. Returns DataFrame indexed by strike with columns: gex_call, gex_put, net_gex, dex, oi_call, oi_put """ M = cfg.contract_multiplier S2 = spot * spot df = df.copy() sign = np.where(df["cp"] == "C", 1.0, -1.0) df["gex"] = sign * df["gamma"].abs() * df["oi"] * M * S2 * 0.01 df["dex_c"] = df["delta"] * df["oi"] * M * spot g = df.groupby("strike") call_mask = df["cp"] == "C" put_mask = df["cp"] == "P" out = pd.DataFrame(index=g.groups.keys()) out["gex_call"] = df.loc[call_mask].groupby("strike")["gex"].sum() out["gex_put"] = df.loc[put_mask].groupby("strike")["gex"].sum() out["gex_call"] = out["gex_call"].fillna(0.0) out["gex_put"] = out["gex_put"].fillna(0.0) out["net_gex"] = out["gex_call"] + out["gex_put"] out["dex"] = df.groupby("strike")["dex_c"].sum() out["oi_call"] = df.loc[call_mask].groupby("strike")["oi"].sum().reindex(out.index).fillna(0) out["oi_put"] = df.loc[put_mask].groupby("strike")["oi"].sum().reindex(out.index).fillna(0) return out.sort_index() def total_net_gex_from_contracts(df: pd.DataFrame, spot: float, cfg: GexConfig) -> float: """Total net GEX summed over a contract frame (FIX 15). Uses the SAME formula as the profile (sign*|gamma|*OI*M*S²*0.01) so that the bars and the profile are directly comparable. Pass the FULL chain frame to get total_net_gex_full (the quantity the profile reproduces at s=spot). """ if df.empty: return 0.0 M = cfg.contract_multiplier S2 = spot * spot sign = np.where(df["cp"].to_numpy() == "C", 1.0, -1.0) gex = sign * df["gamma"].abs().to_numpy() * df["oi"].to_numpy() * M * S2 * 0.01 return float(gex.sum()) def reconcile(total_net_gex: float, grid: np.ndarray, gex_prof: np.ndarray, spot: float, threshold: float = 0.10) -> Dict: """FIX 15: bars and profiles must agree at spot. The simulated GEX profile evaluated at s = spot should approximately equal the total net GEX summed over the SAME contract set the profile uses (the full chain). A large relative error indicates the two paths (reported gamma vs recomputed Black-Scholes gamma, or a time-to-expiry mismatch) have diverged. """ prof_at_spot = float(np.interp(spot, grid, gex_prof)) denom = max(abs(total_net_gex), 1.0) rel_err = abs(prof_at_spot - total_net_gex) / denom return { "total_net_gex": float(total_net_gex), "profile_at_spot": prof_at_spot, "rel_err": float(rel_err), "rel_err_denominator": float(denom), # FIX 73: what the signed error is scored against "pass": bool(rel_err < threshold), } # FIX 60: bar rendering targets a bar COUNT, not a fixed spacing. NDX lists # strikes every 10 points while OI clusters at round numbers, so increment-sized # bars render as thin spikes with gaps. Bucketing for RENDERING ONLY aggregates # strikes into ~target_bars bins; key levels stay at true strike resolution. _BUCKET_LADDER = [0.5, 1, 2.5, 5, 10, 25, 50, 100, 250, 500] def render_bucket(window_span: float, strike_increment: float, target_bars: int = 40) -> float: """FIX 69: bucket = ladder value CLOSEST to raw = window_span / target_bars, then clamped to at least strike_increment. raw = visible_window_span / 40 bucket = closest ladder value [0.5,1,2.5,5,10,25,50,100,250,500] to raw bucket = max(bucket, strike_increment) SMH: span ~98 -> raw 2.45 -> closest 2.5 -> max(2.5, 2.5) = 2.5 (39 bars). NDX: span ~2310 -> raw 57.75 -> closest 50 -> max(50, 10) = 50 (46 bars). """ if window_span <= 0: return float(max(_BUCKET_LADDER[0], strike_increment)) raw = window_span / target_bars bucket = min(_BUCKET_LADDER, key=lambda v: abs(v - raw)) return float(max(bucket, strike_increment)) def reconcile_by_expiry(df_full: pd.DataFrame, spot: float, cfg: GexConfig, worst_n: int = 5) -> Dict: """FIX 45 / FIX 63 / FIX 64: per-expiry reconciliation breakdown. For EVERY expiry, compare Σ(reported-gamma GEX) with Σ(recomputed-BS-gamma GEX) at spot. FIX 63: rank by DOLLAR gap (abs(reported - recomputed)) descending, not rel_err — rel_err surfaces tiny expiries (a -2.8M expiry at 36% rel_err) while missing the dominant one. Publish the top `worst_n`, each with its share of the total dollar gap so it is obvious what to chase. rel_err is kept per entry. FIX 64: also compute an UNSIGNED reconciliation — Σ|per-expiry gap| / Σ|reported GEX| — published as rel_err_unsigned. The signed headline (profile@spot vs total) can pass on cancelling errors; the unsigned figure cannot, so it exposes a historical "pass" that was really offsetting mistakes. """ from scipy.stats import norm if df_full.empty: return {"reconciliation_worst_expiries": [], "zero_greek_contracts_dropped": 0, "rel_err_unsigned": None, "rel_err_unsigned_denominator": 0.0} M = cfg.contract_multiplier S2 = spot * spot r, q = cfg.risk_free_rate, cfg.dividend_yield zero_greek = 0 # per-contract reported-gamma GEX sign = np.where(df_full["cp"].to_numpy() == "C", 1.0, -1.0) reported_gex = sign * df_full["gamma"].abs().to_numpy() * df_full["oi"].to_numpy() * M * S2 * 0.01 df_tmp = df_full.copy() df_tmp["_reported_gex"] = reported_gex rows = [] total_abs_gap = 0.0 total_abs_reported = 0.0 for exp, sub in df_tmp.groupby("expiry"): rep = float(sub["_reported_gex"].sum()) # recomputed BS gamma at spot K = sub["strike"].to_numpy(dtype=float) T = sub["T"].to_numpy(dtype=float) iv = sub["iv"].to_numpy(dtype=float) oi = sub["oi"].to_numpy(dtype=float) sgn = np.where(sub["cp"].to_numpy() == "C", 1.0, -1.0) valid = (T > 0) & (iv > 0) if not valid.any(): rec_gex = 0.0 else: Kv, Tv, ivv, oiv, sgnv = K[valid], T[valid], iv[valid], oi[valid], sgn[valid] sqrt_T = np.sqrt(Tv) d1 = (np.log(spot / Kv) + (r - q + 0.5 * ivv ** 2) * Tv) / (ivv * sqrt_T) gamma_bs = np.exp(-q * Tv) * norm.pdf(d1) / (spot * ivv * sqrt_T) rec_gex = float((sgnv * gamma_bs * oiv * M * S2 * 0.01).sum()) dollar_gap = abs(rec_gex - rep) denom = max(abs(rep), 1.0) rel = dollar_gap / denom total_abs_gap += dollar_gap total_abs_reported += abs(rep) rows.append({ "expiry": str(exp), "reported_gex": rep, "recomputed_gex": rec_gex, "dollar_gap": dollar_gap, "rel_err": float(rel), }) # FIX 63: rank by dollar gap; share of the total dollar gap per entry. rows.sort(key=lambda r: r["dollar_gap"], reverse=True) for r in rows: r["gap_share"] = float(r["dollar_gap"] / total_abs_gap) if total_abs_gap > 0 else 0.0 top = rows[:worst_n] # FIX 64: unsigned reconciliation (immune to signed cancellation). rel_err_unsigned = float(total_abs_gap / total_abs_reported) if total_abs_reported > 0 else None return { "reconciliation_worst_expiries": top, "zero_greek_contracts_dropped": zero_greek, "rel_err_unsigned": rel_err_unsigned, # FIX 73: Σ|reported GEX| — what the unsigned error is scored against. A # variant that drops the largest expiry is scored against a smaller book. "rel_err_unsigned_denominator": float(total_abs_reported), } def reconciliation_floor_unsigned(df_full: pd.DataFrame, spot: float, cfg: GexConfig) -> Dict: """FIX 87: the irreducible unsigned reconciliation error from PUBLICATION ROUNDING, measured by the grid_rounding method. Method (grid_rounding): take the recomputed Black-Scholes gamma as the full-precision TRUTH, round it to the observed publication grid (Cboe publishes gamma to 4 decimal places, so the grid interval is ±0.5 × 10^-4), and measure the aggregate unsigned error of the ROUNDED-vs-UNROUNDED difference. FIX 93: this is DETERMINISTIC. Rounding to a fixed grid is not a random process — each gamma rounds to exactly one published value — so the floor is a single scalar, not a Monte Carlo distribution. The previous version simulated 2000 uniform draws within the grid interval and reported a median/p95 spread, which implied a variability that does not exist. The deterministic floor is computed once: round every recomputed gamma to 4dp, aggregate by expiry, and measure the unsigned error of rounded-vs-unrounded. Why the actual unsigned error can sit BELOW this floor: the floor measures the quantisation noise from rounding the RECOMPUTED gamma to the publication grid. The actual unsigned error measures the net difference between the REPORTED (Cboe-rounded) gamma and the recomputed gamma. The reported gamma is itself rounded to the same 4dp grid, so the two rounding errors (reported vs recomputed) can partially cancel. The floor is an upper bound on the rounding contribution to the unsigned error, not a prediction of the actual error. When the actual error (NDX: 0.10345) sits below the floor (0.13606), it means the model error and the rounding error are partially offsetting — the reconciliation is better than the rounding noise alone would suggest. Why this replaces the FIX 85 estimator: the old method perturbed reported_gamma (ALREADY rounded) by ±0.5 × 10^-4 and compared it against the recomputed gamma. That double-counts rounding (the reported value is rounded, then perturbed again) AND folds in the genuine model error, so it overstates the floor. The grid_rounding method isolates pure rounding: the truth and the rounded value are BOTH derived from the recomputed gamma, so the model error cancels and only the quantisation remains. If the floor brackets the unsigned gate, the gate is dominated by rounding noise (see FIX 88 for the capped-gate / indeterminate handling). """ from scipy.stats import norm if df_full.empty: return {"method": "grid_rounding", "rounding_interval": None, "floor": None, "brackets_gate": False} M = cfg.contract_multiplier S2 = spot * spot r, q = cfg.risk_free_rate, cfg.dividend_yield # Publication grid: Cboe publishes gamma to 4dp -> half a unit in the last place. rounding_interval = 0.5e-4 # Recompute the full-precision BS gamma per contract — this is the TRUTH that # rounding is applied to. n = len(df_full) gamma_bs = np.zeros(n) K = df_full["strike"].to_numpy(dtype=float) T = df_full["T"].to_numpy(dtype=float) iv = df_full["iv"].to_numpy(dtype=float) valid = (T > 0) & (iv > 0) if valid.any(): Kv, Tv, ivv = K[valid], T[valid], iv[valid] sqrt_T = np.sqrt(Tv) d1 = (np.log(spot / Kv) + (r - q + 0.5 * ivv ** 2) * Tv) / (ivv * sqrt_T) gamma_bs[valid] = np.exp(-q * Tv) * norm.pdf(d1) / (spot * ivv * sqrt_T) sign = np.where(df_full["cp"].to_numpy() == "C", 1.0, -1.0) oi = df_full["oi"].to_numpy(dtype=float) # Map each contract to its expiry ordinal so we can aggregate per expiry. exp_list = list(df_full["expiry"].unique()) exp_ordinal = {e: i for i, e in enumerate(exp_list)} contract_exp = np.array([exp_ordinal[e] for e in df_full["expiry"]]) n_exp = len(exp_list) # Per-expiry TRUTH GEX (unrounded recomputed gamma) and ROUNDED GEX (gamma # rounded to the 4dp publication grid). The difference is the pure quantisation. truth_gex_contract = sign * gamma_bs * oi * M * S2 * 0.01 gamma_rounded = np.round(gamma_bs, 4) rounded_gex_contract = sign * gamma_rounded * oi * M * S2 * 0.01 truth_by_exp = np.zeros(n_exp) rounded_by_exp = np.zeros(n_exp) for i in range(n_exp): mask = contract_exp == i truth_by_exp[i] = truth_gex_contract[mask].sum() rounded_by_exp[i] = rounded_gex_contract[mask].sum() total_abs_truth = float(np.abs(truth_by_exp).sum()) total_abs_gap = float(np.abs(rounded_by_exp - truth_by_exp).sum()) floor = (total_abs_gap / total_abs_truth) if total_abs_truth > 0 else 0.0 return { "method": "grid_rounding", "rounding_interval": rounding_interval, "floor": floor, "brackets_gate": bool(floor >= cfg.unsigned_gate), } def _profile_outlier_guard(contrib: np.ndarray, meta: List[Dict], grid: np.ndarray, spot: float, ratio_threshold: float = 5.0, width_threshold_pct: float = 5.0) -> Tuple[np.ndarray, List[Dict]]: """FIX 48 / v1.6.2: drop a contract that is BOTH dominant AND a narrow spike. The expired-contract explosion this targets has a distinctive signature: a floored T (1/(252*13)) makes Black-Scholes gamma collapse to a delta function, so the contract's contribution is concentrated in a narrow price band around its strike while dwarfing everything else there. A legitimate dominant contract — the ATM 0DTE (~50% of others, ratio ~0.5) or a far-OTM wing-shaper (locally huge but spread across its whole wing) — is never both dominant AND narrow. v1.6.2: the width test is now in PRICE terms (% of spot) rather than a fraction of grid points, so one threshold means the same thing on every grid (GEX, the wider HVL grid, the DEX grid) and every chain regardless of profile_band or grid resolution. Thresholds (calibrated on three known cases): - ratio_threshold = 5.0: peak |contribution| > 5x the sum of all others. Expired explosion: 406x. Far-OTM wing-shaper: 568x. ATM 0DTE (realistic): ~0.5x (kept by this test alone). - width_threshold_pct = 5.0: half-peak width < 5% of spot. Expired explosion: 1.1% (narrow → dropped). Far-OTM wing-shaper: 12.0% (wide → kept). ATM 0DTE: 0.96% (narrow, but ratio 0.5 < 5.0 → kept). Margins: expired 1.1% vs 5% threshold (3.9pp headroom); far-OTM 12.0% vs 5% (7.0pp headroom). Since FIX 48 drops expired contracts by raw calendar date before they reach the profile, this guard should be unreachable in normal operation. Any firing is a fault, logged at ERROR. `contrib` is (n_contracts, n_grid); `meta[i]` carries strike/expiry for contract i. `grid` is the price grid; `spot` is the current spot price. Returns (profile, outliers) where profile excludes the dropped contracts and outliers is a list of {strike, expiry, value, price_width_pct}. """ n, G = contrib.shape if n < 2: return contrib.sum(axis=0), [] abs_c = np.abs(contrib) total_abs = abs_c.sum(axis=0) # (G,) peak_idx = abs_c.argmax(axis=1) # (n,) grid idx of each contract's peak peak_val = abs_c[np.arange(n), peak_idx] # (n,) others_at_peak = total_abs[peak_idx] - peak_val # (n,) peak_ratio = peak_val / (others_at_peak + 1e-30) # (n,) # spike width in PRICE terms: span of grid points where |contrib| >= half peak, # expressed as % of spot. Grid-independent: same contract measures the same on # any grid resolution or band width. half = 0.5 * peak_val[:, None] above = abs_c >= half # (n, G) grid_arr = np.asarray(grid, dtype=float) price_width_pct = np.zeros(n) for i in range(n): idx = np.where(above[i])[0] if len(idx) > 0: price_width_pct[i] = (grid_arr[idx[-1]] - grid_arr[idx[0]]) / spot * 100.0 # an expired/floored-T contract is dominant AND narrow; legitimate structure is # either not dominant (ATM 0DTE, ratio ~0.5) or not narrow (far-OTM, width 12%) is_outlier = (peak_ratio > ratio_threshold) & (price_width_pct < width_threshold_pct) outliers = [] if is_outlier.any(): for i in np.where(is_outlier)[0]: outliers.append({ "strike": float(meta[i]["strike"]), "expiry": str(meta[i]["expiry"]), "value": float(peak_val[i]), "price_width_pct": round(float(price_width_pct[i]), 3), }) logger.error("FIX 48 FAULT: profile outlier guard fired — %d contract(s) " "dropped (dominant narrow spike; should be unreachable after " "FIX 48 calendar-date filter): %s", len(outliers), outliers[:5]) profile = contrib[~is_outlier].sum(axis=0) return profile, outliers def gex_profile(contracts_df: pd.DataFrame, spot: float, cfg: GexConfig, grid: np.ndarray) -> Tuple[np.ndarray, List[Dict]]: """§4b GEX Profile: re-evaluate total net GEX as if spot were at each grid level. Uses reported iv held constant, recomputed BS gamma. FIX 48 / v1.6.2: returns (profile, outliers). A per-contract outlier guard drops any contract that is BOTH dominant (peak |contribution| > 5x the sum of all others at that grid point) AND a narrow spike (half-peak width < 5% of spot in price terms). This is the expired-contract gamma-explosion signature: a floored T collapses BS gamma to a delta function. Legitimate dominant structure (the ATM 0DTE, ratio ~0.5; far-OTM wing-shapers, width ~12%) is kept. Since FIX 48 drops expired contracts by raw calendar date before they reach the profile, this guard should be unreachable in normal operation; any firing is a fault (logged at ERROR, surfaced on the chart). Outliers are reported as a list of {strike, expiry, value, price_width_pct}. (FIX 14: the unused `agg` argument was removed from the signature.) """ from scipy.stats import norm M = cfg.contract_multiplier r, q = cfg.risk_free_rate, cfg.dividend_yield s = np.asarray(grid, dtype=float) if contracts_df.empty or len(s) == 0: return np.zeros_like(s, dtype=float), [] K = contracts_df["strike"].to_numpy(dtype=float) T = contracts_df["T"].to_numpy(dtype=float) iv = contracts_df["iv"].to_numpy(dtype=float) oi = contracts_df["oi"].to_numpy(dtype=float) sign = np.where(contracts_df["cp"].to_numpy() == "C", 1.0, -1.0) expiry = contracts_df["expiry"].to_numpy() # valid mask per contract (T>0, iv>0); s>0 guaranteed for price grid valid = (T > 0) & (iv > 0) K, T, iv, oi, sign, expiry = K[valid], T[valid], iv[valid], oi[valid], sign[valid], expiry[valid] if len(K) == 0: return np.zeros_like(s, dtype=float), [] sqrt_T = np.sqrt(T)[:, None] # (n,1) s_b = s[None, :] # (1,G) d1 = (np.log(s_b / K[:, None]) + (r - q + 0.5 * iv[:, None] ** 2) * T[:, None]) \ / (iv[:, None] * sqrt_T) gamma = np.exp(-q * T[:, None]) * norm.pdf(d1) / (s_b * iv[:, None] * sqrt_T) contrib = sign[:, None] * gamma * oi[:, None] * M * (s_b ** 2) * 0.01 # (n, G) meta = [{"strike": K[i], "expiry": expiry[i]} for i in range(len(K))] return _profile_outlier_guard(contrib, meta, s, spot) def dex_profile(contracts_full: pd.DataFrame, cfg: GexConfig, grid: np.ndarray, spot: float) -> Tuple[np.ndarray, List[Dict]]: """Total dealer delta exposure re-evaluated as if spot were at each grid level. Per contract: bs_delta(s, K, T, iv, cp, r, q) * OI * M * s. Call delta positive, put delta already negative — NO extra sign applied. Returns (profile, outliers) — same length as grid, in dollars. Shape (FIX 12): generally rising in spot, but with a V-shaped minimum where deep-ITM put delta dominates. NOT monotonic — has an interior minimum (see dex_min_price). Window-independent. FIX 48: applies the same per-contract dominance guard as gex_profile. """ from .greeks import bs_delta_2d M = cfg.contract_multiplier r, q = cfg.risk_free_rate, cfg.dividend_yield s = np.asarray(grid, dtype=float) if contracts_full.empty or len(s) == 0: return np.zeros_like(s, dtype=float), [] K = contracts_full["strike"].to_numpy(dtype=float) T = contracts_full["T"].to_numpy(dtype=float) iv = contracts_full["iv"].to_numpy(dtype=float) oi = contracts_full["oi"].to_numpy(dtype=float) cp_sign = np.where(contracts_full["cp"].to_numpy() == "C", 1.0, -1.0) expiry = contracts_full["expiry"].to_numpy() valid = (T > 0) & (iv > 0) K, T, iv, oi, cp_sign, expiry = K[valid], T[valid], iv[valid], oi[valid], cp_sign[valid], expiry[valid] if len(K) == 0: return np.zeros_like(s, dtype=float), [] delta = bs_delta_2d(s[None, :], K[:, None], T[:, None], iv[:, None], cp_sign[:, None], r, q) # (n, G) contrib = delta * oi[:, None] * M * s[None, :] # (n, G) meta = [{"strike": K[i], "expiry": expiry[i]} for i in range(len(K))] return _profile_outlier_guard(contrib, meta, s, spot) def find_delta_neutral(grid: np.ndarray, dex_prof: np.ndarray, spot: float) -> Tuple[Optional[float], list]: """DEX-profile zero-crossings (FIX 12). A V-shaped DEX profile can cross zero twice. Returns (nearest, crossings): - crossings: ALL grid prices where the profile crosses zero (linear interp). - nearest: the crossing nearest spot (same rule as HVL), or None if no crossing. """ crossings = [] for i in range(1, len(dex_prof)): y0, y1 = dex_prof[i - 1], dex_prof[i] if (y0 < 0 <= y1) or (y0 > 0 >= y1): x0, x1 = grid[i - 1], grid[i] cross = x0 if y1 == y0 else x0 - y0 * (x1 - x0) / (y1 - y0) crossings.append(float(cross)) if not crossings: return None, [] nearest = min(crossings, key=lambda c: abs(c - spot)) return nearest, crossings def dex_min_price(grid: np.ndarray, dex_prof: np.ndarray) -> Tuple[Optional[float], str]: """Price at the V-shaped minimum of the DEX profile (FIX 12). FIX 26: routed through interior_extremum — if the minimum sits within 2 grid steps of either edge it is a boundary artifact (the true V-minimum lies beyond the grid), so return (None, "at_grid_boundary") rather than a misleading number. The DEX grid is widened to ±40% (profile_band_dex) so the true minimum can be located; if it is still at the edge there, it is genuinely off-grid. """ return interior_extremum(dex_prof, grid, which="min", edge_tol=2) def detect_increment(strikes: np.ndarray) -> float: """Detect dominant strike increment as mode of diffs.""" s = np.sort(np.unique(strikes)) if len(s) < 2: return 1.0 diffs = np.round(np.diff(s), 6) diffs = diffs[diffs > 0] if len(diffs) == 0: return 1.0 vals, counts = np.unique(np.round(diffs, 4), return_counts=True) return float(vals[np.argmax(counts)]) def detect_render_spacing(strikes: np.ndarray, net_gex: np.ndarray, fallback_increment: float) -> float: """FIX 44: bar height should follow POPULATED strike spacing, not listed spacing. detect_increment returns the mode of ALL listed diffs (e.g. 25 for NDX), but strikes carrying non-trivial OI may sit ~100 apart, giving tiny bars with gaps on a wide axis. This returns the median spacing among strikes that actually carry a bar (|net_gex| >= 1% of the max), so bars render contiguous. Falls back to fallback_increment if fewer than 3 populated strikes. """ if len(strikes) == 0: return fallback_increment max_abs = float(np.max(np.abs(net_gex))) if len(net_gex) else 0.0 if max_abs <= 0: return fallback_increment sig_mask = np.abs(net_gex) >= 0.01 * max_abs sig = np.sort(strikes[sig_mask]) if len(sig) < 3: return fallback_increment diffs = np.diff(sig) diffs = diffs[diffs > 0] if len(diffs) == 0: return fallback_increment return float(np.median(diffs)) def _hvl_zero_cross(grid: np.ndarray, profile: np.ndarray, spot: float) -> Tuple[Optional[float], list]: """Sign flip (neg->pos) of the simulated GEX profile. Returns (nearest_to_spot, all_crossings). This is the ONLY HVL definition as of v1.5.0 (FIX 29): a level that separates a positive-gamma regime above from a negative-gamma regime below IS by definition the sign change of the gamma profile. Any other construction cannot partition the price axis that way. """ crossings = [] for i in range(1, len(profile)): if profile[i - 1] < 0 and profile[i] >= 0: x0, x1 = grid[i - 1], grid[i] y0, y1 = profile[i - 1], profile[i] cross = x0 if y1 == y0 else x0 - y0 * (x1 - x0) / (y1 - y0) crossings.append(float(cross)) if not crossings: return None, [] nearest = min(crossings, key=lambda c: abs(c - spot)) return nearest, crossings def interior_extremum(arr: np.ndarray, grid: np.ndarray, which: str = "min", edge_tol: int = 2) -> Tuple[Optional[float], str]: """FIX 26: return (price, status) for an extremum, guarding grid-boundary artifacts. If the argmin/argmax falls within `edge_tol` grid steps of either edge, the extremum is a boundary artifact (the true turning point lies beyond the grid), so return (None, "at_grid_boundary"). Otherwise return (price, "interior"). """ if len(arr) == 0: return None, "empty" idx = int(np.argmin(arr) if which == "min" else np.argmax(arr)) if idx < edge_tol or idx >= len(arr) - edge_tol: return None, "at_grid_boundary" return float(grid[idx]), "interior" def _hvl_bar_flip(agg: pd.DataFrame, spot: float) -> Optional[float]: """Linear interpolation of the strike where the PER-STRIKE net_gex bar series changes sign neg->pos, crossing nearest spot. NOTE (FIX 30): this is NOT an HVL candidate. It measures where the per-strike call/put OI composition flips locally; the HVL zero crossing measures where the whole portfolio's gamma flips sign. Different quantities. Promoted to its own named level `gex_transition` (see compute_gex_transition). """ if agg.empty: return None strikes = agg.index.to_numpy(dtype=float) net = agg["net_gex"].to_numpy(dtype=float) order = np.argsort(strikes) strikes, net = strikes[order], net[order] crossings = [] for i in range(1, len(net)): if net[i - 1] < 0 and net[i] >= 0: x0, x1 = strikes[i - 1], strikes[i] y0, y1 = net[i - 1], net[i] cross = x0 if y1 == y0 else x0 - y0 * (x1 - x0) / (y1 - y0) crossings.append(float(cross)) if not crossings: return None return min(crossings, key=lambda c: abs(c - spot)) def compute_gex_transition(agg: pd.DataFrame, spot: float, persistence: int = 3) -> Dict: """FIX 30: GEX Transition — where strike-level net gamma changes sign locally. This is the per-strike bar sign flip (formerly "bar_flip"), promoted to its own named level. It is NOT an HVL candidate: HVL marks where TOTAL portfolio gamma changes sign; GEX Transition marks where strike-level net gamma changes sign locally. They coincide only in balanced chains; a large concentrated wall separates them. Hardened against noise: the sign must hold for at least `persistence` consecutive populated strikes on EACH side of the candidate. If no candidate satisfies that, publish null with gex_transition_status = "no_persistent_flip". """ result = { "gex_transition": None, "gex_transition_status": "no_persistent_flip", "gex_transition_distance_pct": None, } if agg.empty: return result strikes = agg.index.to_numpy(dtype=float) net = agg["net_gex"].to_numpy(dtype=float) order = np.argsort(strikes) strikes, net = strikes[order], net[order] # candidate crossings neg->pos crossings = [] for i in range(1, len(net)): if net[i - 1] < 0 and net[i] >= 0: x0, x1 = strikes[i - 1], strikes[i] y0, y1 = net[i - 1], net[i] cross = x0 if y1 == y0 else x0 - y0 * (x1 - x0) / (y1 - y0) crossings.append((float(cross), i)) if not crossings: return result # persistence check: `persistence` consecutive strikes negative BEFORE the flip # and `persistence` consecutive strikes positive AFTER it. persistent = [] for cross, i in crossings: before = net[max(0, i - persistence):i] after = net[i:i + persistence] if len(before) >= persistence and len(after) >= persistence \ and (before < 0).all() and (after > 0).all(): persistent.append(cross) if not persistent: result["gex_transition_status"] = "no_persistent_flip" return result nearest = min(persistent, key=lambda c: abs(c - spot)) result.update({ "gex_transition": float(nearest), "gex_transition_status": "ok", "gex_transition_distance_pct": float((nearest - spot) / spot), }) return result def compute_hvl(grid: np.ndarray, profile: np.ndarray, spot: float, cfg: GexConfig) -> Dict: """FIX 29: HVL = the gamma-profile zero crossing nearest spot. Always defined. No "indeterminate". Distance is information, not a defect — publish hvl_distance_pct and a hvl_regime_note. If no crossing exists in the grid, return hvl=None with hvl_status="no_flip_in_range" (the caller widens the grid to ±40% and retries; if still none, that is a real market state, not an error). """ nearest, crossings = _hvl_zero_cross(grid, profile, spot) if nearest is None: return { "hvl": None, "hvl_status": "no_flip_in_range", "hvl_distance_pct": None, "hvl_regime_note": None, "hvl_crossings": [], } distance_pct = (nearest - spot) / spot d = abs(distance_pct) sign = "positive" if profile_at_spot_sign(profile, grid, spot) > 0 else "negative" if d <= 0.03: note = "near spot — regime flip in play" elif d <= 0.08: note = "moderately distant" else: note = f"far from spot — no nearby regime flip; sustained {sign} gamma" return { "hvl": float(nearest), "hvl_status": "ok", "hvl_distance_pct": float(distance_pct), "hvl_regime_note": note, "hvl_crossings": [float(c) for c in crossings], } def profile_at_spot_sign(profile: np.ndarray, grid: np.ndarray, spot: float) -> float: """Sign of the interpolated GEX profile at spot (helper for regime note).""" return float(np.interp(spot, grid, profile)) def compute_hvl_candidates(grid: np.ndarray, profile: np.ndarray, agg: pd.DataFrame, spot: float, cfg: GexConfig, increment: float = 0.0) -> Dict: """DEPRECATED (v1.5.0): kept for back-compat only. Use compute_hvl + compute_gex_transition. The old three-candidate system (inflection / zero_cross / bar_flip) with indeterminacy gating is retired. HVL is now always the zero crossing (FIX 29); bar_flip is promoted to gex_transition (FIX 30); inflection is deleted (FIX 28). """ hvl_info = compute_hvl(grid, profile, spot, cfg) gt_info = compute_gex_transition(agg, spot, cfg.gex_transition_persistence) return { **hvl_info, **gt_info, "hvl_rule_used": "zero_cross", "hvl_confidence": "high" if hvl_info["hvl"] is not None else "n/a", "hvl_candidates": {"zero_cross": hvl_info["hvl"]}, "hvl_spread_pct": 0.0, "hvl_range": None, "hvl_inflection_masked_strikes": [], "hvl_inflection_status": "retired_v1.5.0", } def find_hvl(grid: np.ndarray, profile: np.ndarray, spot: float) -> Tuple[float, str, str, list]: """Back-compat wrapper (pre-FIX-16 signature). Returns (hvl, rule, confidence, crossings) using the zero-crossing rule. New code should use compute_hvl. FIX 28: inflection fallback deleted — it structurally re-finds the dominant put wall (curvature peaks where gamma concentrates). """ crossings = [] for i in range(1, len(profile)): if profile[i - 1] < 0 and profile[i] >= 0: x0, x1 = grid[i - 1], grid[i] y0, y1 = profile[i - 1], profile[i] cross = x0 if y1 == y0 else x0 - y0 * (x1 - x0) / (y1 - y0) crossings.append(float(cross)) if crossings: best = min(crossings, key=lambda c: abs(c - spot)) return best, "zero_crossing", "high", crossings # No crossing: return midpoint with low confidence (no inflection fallback) mid = float(grid[len(grid) // 2]) return mid, "no_crossing_midpoint", "low", crossings def compute_levels(agg: pd.DataFrame, grid: np.ndarray, profile: np.ndarray, spot: float, increment: float, cfg: GexConfig = None, profile_at_spot: float = None) -> Dict: """Compute all §5/§6 key levels + diagnostics. FIX 29: HVL = gamma-profile zero crossing nearest spot. Always a single number. Distance is information, not a defect — publish hvl_distance_pct + hvl_regime_note. FIX 28: inflection rule retired (it structurally re-finds the dominant put wall). FIX 30: bar_flip promoted to gex_transition (separate named level, persistence-hardened). FIX 22: gamma_condition = sign of the simulated GEX profile at spot (unchanged). `cfg` is optional for back-compat; when omitted a default config is used. """ if cfg is None: cfg = GexConfig() call_res = float(agg["gex_call"].idxmax()) put_sup = float(agg["gex_put"].idxmin()) # FIX 29: HVL = zero crossing nearest spot hvl_info = compute_hvl(grid, profile, spot, cfg) hvl_raw = hvl_info["hvl"] # snap HVL to nearest strike increment if hvl_raw is not None: hvl = round(hvl_raw / increment) * increment if increment > 0 else hvl_raw else: hvl = None # FIX 30: GEX Transition (formerly bar_flip, now its own named level) gt_info = compute_gex_transition(agg, spot, cfg.gex_transition_persistence) total_net = float(agg["net_gex"].sum()) # FIX 22: gamma_condition = sign of the simulated profile at spot (direct measure) if profile_at_spot is None: profile_at_spot = float(np.interp(spot, grid, profile)) gamma_condition = "POSITIVE" if profile_at_spot > 0 else "NEGATIVE" # distance to the regime flip distance_to_flip_pct = hvl_info["hvl_distance_pct"] sum_gex_put = float(agg["gex_put"].sum()) sum_gex_call = float(agg["gex_call"].sum()) gex_pc_ratio = abs(sum_gex_put) / sum_gex_call if sum_gex_call != 0 else float("inf") sum_oi_put = float(agg["oi_put"].sum()) sum_oi_call = float(agg["oi_call"].sum()) oi_pc_ratio = sum_oi_put / sum_oi_call if sum_oi_call != 0 else float("inf") return { "call_resistance": call_res, "put_support": put_sup, "hvl": hvl, "hvl_raw": hvl_raw, "hvl_status": hvl_info["hvl_status"], "hvl_distance_pct": hvl_info["hvl_distance_pct"], "hvl_regime_note": hvl_info["hvl_regime_note"], "hvl_crossings": hvl_info["hvl_crossings"], "hvl_rule": "zero_cross", "hvl_rule_used": "zero_cross", "hvl_confidence": "high" if hvl is not None else "n/a", "hvl_candidates": {"zero_cross": hvl_raw}, "hvl_spread_pct": 0.0, "hvl_range": None, "hvl_inflection_masked_strikes": [], # FIX 30: GEX Transition "gex_transition": gt_info["gex_transition"], "gex_transition_status": gt_info["gex_transition_status"], "gex_transition_distance_pct": gt_info["gex_transition_distance_pct"], # FIX 22: gamma condition "spot": float(spot), "total_net_gex": total_net, "gamma_condition": gamma_condition, "gamma_condition_basis": "sign of simulated GEX profile at spot", "net_gex_at_spot": float(profile_at_spot), "distance_to_flip_pct": distance_to_flip_pct, "gex_put_call_ratio": gex_pc_ratio, "oi_put_call_ratio": oi_pc_ratio, } def _atm_iv_for_expiry(ref: pd.DataFrame, spot: float, contract_oi_floor: float = 250.0 ) -> Tuple[Optional[float], list, str]: """ATM IV from a single expiry frame (FIX 24; gates made relative in FIX 43). Requires BOTH bracketing strikes to pass (OI >= contract_oi_floor, two-sided market bid > 0, 0.05 < iv < 2.0); averages call & put IV at each strike, interpolates in strike to spot. Returns (atm_iv, source_strikes, status). If only one strike passes, still computes but status = "single_strike_no_interpolation" — never "ok". """ below = ref[ref["strike"] <= spot] above = ref[ref["strike"] >= spot] bracket_strikes = [] if not below.empty: bracket_strikes.append(float(below["strike"].max())) if not above.empty: k = float(above["strike"].min()) if k not in bracket_strikes: bracket_strikes.append(k) if not bracket_strikes: return None, [], "no_bracket" pts = [] # (strike, iv) for k in bracket_strikes: legs = ref[ref["strike"] == k] good = legs[(legs["oi"] >= contract_oi_floor) & (legs["bid"] > 0) & (legs["iv"] > 0.05) & (legs["iv"] < 2.0)] if good.empty: continue pts.append((k, float(good["iv"].mean()))) if not pts: return None, [], "no_strike_passes" if len(pts) == 1: return pts[0][1], [p[0] for p in pts], "single_strike_no_interpolation" pts.sort() (k0, iv0), (k1, iv1) = pts[0], pts[1] if k1 == k0: atm_iv = 0.5 * (iv0 + iv1) else: w = (spot - k0) / (k1 - k0) atm_iv = iv0 + w * (iv1 - iv0) return atm_iv, [p[0] for p in pts], "ok" def atm_expected_move(contracts_df: pd.DataFrame, spot: float) -> Dict: """1d expected move from ATM IV (FIX 17 picker, hardened in FIX 24). Reference expiry = nearest expiry with DTE >= 5 AND total expiry OI >= max(5000, 0.02 * chain_total_oi) — a liquidity gate so we never source IV from a near-dead expiry (FIX 24a). On the two strikes bracketing spot, require BOTH to pass (OI >= 250, two-sided market bid > 0, 0.05 < iv < 2.0); average call & put IV at each strike, interpolate in strike to spot. Cross-check (FIX 24d): compute atm_iv from the two nearest qualifying expiries; if they differ by more than 0.15 absolute, the term structure is unstable and the status says so. SANITY GATE: if atm_iv is outside [0.05, 1.50] or no strike passes, publish nulls and atm_iv_status = "rejected: ". Never publish a number we cannot defend, and never report "ok" unless two strikes on a liquid expiry were interpolated. """ result = { "atm_iv": None, "exp_move_pct": None, "min_price": None, "max_price": None, "atm_iv_expiry": None, "atm_iv_dte": None, "atm_iv_source_strikes": None, "atm_iv_expiry_oi": None, "atm_iv_expiry_rank": None, "atm_iv_alt": None, "atm_iv_term_spread": None, "atm_iv_status": None, } if contracts_df.empty or "dte" not in contracts_df.columns: result["atm_iv_status"] = "rejected: empty or missing DTE" return result chain_total_oi = float(contracts_df["oi"].sum()) # FIX 43: chain-relative gates (were absolute 5000 / 250, SMH-scaled). oi_floor = max(0.02 * chain_total_oi, 200.0) contract_oi_floor = max(0.0002 * chain_total_oi, 10.0) # FIX 24a: expiry-level liquidity gate (DTE >= 5 AND expiry OI >= floor) elig = contracts_df[contracts_df["dte"] >= 5] if elig.empty: result["atm_iv_status"] = "rejected: no expiry with DTE >= 5" return result exp_oi = elig.groupby("expiry")["oi"].sum() liquid = exp_oi[exp_oi >= oi_floor] if liquid.empty: result["atm_iv_status"] = ( f"rejected: no DTE>=5 expiry holds OI >= {oi_floor:.0f} " f"(max(200, 2% of chain {chain_total_oi:.0f}))") return result # rank qualifying expiries by DTE (nearest first) elig_dte = elig.groupby("expiry")["dte"].min().sort_values() qual_exps = [e for e in elig_dte.index if e in liquid.index] ref_exp = qual_exps[0] ref = contracts_df[contracts_df["expiry"] == ref_exp].copy() ref_dte = int(ref["dte"].iloc[0]) ref_oi = int(exp_oi[ref_exp]) ref_rank = int(list(elig_dte.index).index(ref_exp)) + 1 # FIX 43: try the nearest qualifying expiry, then widen to the three nearest # before falling back — a low-OI chain may have a thin front expiry but a # perfectly good second or third. atm_iv = None src_strikes = [] status = "no_strike_passes" for candidate_exp in qual_exps[:3]: cand_ref = contracts_df[contracts_df["expiry"] == candidate_exp].copy() cand_iv, cand_strikes, cand_status = _atm_iv_for_expiry( cand_ref, spot, contract_oi_floor) if cand_iv is not None: atm_iv = cand_iv src_strikes = cand_strikes status = cand_status # re-point the published metadata at the expiry that actually resolved ref_exp = candidate_exp ref_dte = int(cand_ref["dte"].iloc[0]) ref_oi = int(exp_oi[candidate_exp]) ref_rank = int(list(elig_dte.index).index(candidate_exp)) + 1 break if atm_iv is None: result["atm_iv_status"] = f"rejected: {status} on expiry {ref_exp}" result["atm_iv_expiry"] = str(ref_exp) result["atm_iv_dte"] = ref_dte result["atm_iv_expiry_oi"] = ref_oi result["atm_iv_expiry_rank"] = ref_rank return result # sanity gate on the primary pick if not (0.05 <= atm_iv <= 1.50): result["atm_iv_status"] = f"rejected: atm_iv {atm_iv:.3f} outside [0.05, 1.50]" result["atm_iv_source_strikes"] = src_strikes result["atm_iv_expiry"] = str(ref_exp) result["atm_iv_dte"] = ref_dte result["atm_iv_expiry_oi"] = ref_oi result["atm_iv_expiry_rank"] = ref_rank return result # FIX 24d: term-structure cross-check against the nearest OTHER qualifying # expiry (FIX 43: the resolved expiry may not be qual_exps[0] any more). atm_iv_alt = None term_spread = None alt_candidates = [e for e in qual_exps if e != ref_exp] if alt_candidates: alt_exp = alt_candidates[0] alt_ref = contracts_df[contracts_df["expiry"] == alt_exp].copy() alt_iv, _, alt_status = _atm_iv_for_expiry(alt_ref, spot, contract_oi_floor) if alt_iv is not None and 0.05 <= alt_iv <= 1.50: atm_iv_alt = float(alt_iv) term_spread = float(abs(atm_iv - alt_iv)) if term_spread > 0.15 and status == "ok": status = "term_structure_unstable" exp_move_pct = atm_iv * np.sqrt(1 / 252.0) result.update({ "atm_iv": float(atm_iv), "exp_move_pct": float(exp_move_pct), "min_price": float(spot * (1 - exp_move_pct)), "max_price": float(spot * (1 + exp_move_pct)), "atm_iv_expiry": str(ref_exp), "atm_iv_dte": ref_dte, "atm_iv_source_strikes": src_strikes, "atm_iv_expiry_oi": ref_oi, "atm_iv_expiry_rank": ref_rank, "atm_iv_alt": atm_iv_alt, "atm_iv_term_spread": term_spread, "atm_iv_status": status, }) return result def build_outlier_report(agg: pd.DataFrame, df_band: pd.DataFrame, cfg: GexConfig, spot: float, top_n: int = 5) -> Dict: """FIX 13: top strikes by |net_gex| with OI + per-expiry breakdown. Lets a reader judge whether a dominant bar is real (e.g. concentrated LEAPS put OI) or a parsing artifact, instead of just clipping it off the chart. FIX 91: the per-expiry breakdown uses spot^2, the SAME price term as the bars it explains (aggregate()). It previously used each strike's own K^2, which made the breakdown disagree with its parent bar by (K/S)^2 and broke the invariant Sigma by_expiry[].gex == net_gex. spot is passed in explicitly so the decomposition cannot drift from the aggregate it documents. """ M = cfg.contract_multiplier S2 = spot * spot if agg.empty: return {"top_strikes": []} order = agg["net_gex"].abs().sort_values(ascending=False).index top = [] for k in order[:top_n]: row = agg.loc[k] sub = df_band[df_band["strike"] == k] by_exp = [] for expiry, g in sub.groupby("expiry"): gc = g[g["cp"] == "C"] gp = g[g["cp"] == "P"] # per-expiry gex contribution (signed): calls +, puts - gex_exp = (gc["gamma"].abs() * gc["oi"]).sum() - (gp["gamma"].abs() * gp["oi"]).sum() by_exp.append({ "expiry": str(expiry), "oi_call": int(gc["oi"].sum()), "oi_put": int(gp["oi"].sum()), "gex": float(gex_exp * M * S2 * 0.01), }) top.append({ "strike": float(k), "net_gex": float(row["net_gex"]), "oi_call": int(row["oi_call"]), "oi_put": int(row["oi_put"]), "by_expiry": by_exp, }) return {"top_strikes": top} def level_front_expiry_pct(df_full: pd.DataFrame, strike: float, spot: float, cfg: GexConfig) -> Dict: """FIX 27 / FIX 46: fraction of a key level's GEX attributable to the FRONT expiry. Put Support / Call Resistance can be almost entirely 0DTE open interest that ceases to exist at today's close — an ephemeral level. FIX 46: the old denominator |net gex at strike| could be smaller than the numerator when later expiries offset with opposite sign, yielding a "127%" footnote that reads as a bug. Now publish TWO figures: - front_expiry_abs_share = |front gex| / Σ|gex across all expiries| (<= 100%) - front_expiry_net_ratio = |front gex| / |net gex at strike| (diagnostic, may exceed 1.0) Returns a dict; empty dict if the strike has no GEX. """ M = cfg.contract_multiplier S2 = spot * spot if df_full.empty: return {} sub = df_full[df_full["strike"] == strike] if sub.empty: return {} sign = np.where(sub["cp"].to_numpy() == "C", 1.0, -1.0) gex = sign * sub["gamma"].abs().to_numpy() * sub["oi"].to_numpy() * M * S2 * 0.01 abs_total = float(np.abs(gex).sum()) net_total = abs(float(gex.sum())) if abs_total == 0: return {} front_exp = df_full["expiry"].min() front_mask = (sub["expiry"] == front_exp).to_numpy() front_gex = abs(float(gex[front_mask].sum())) result = { "front_expiry": str(front_exp), "front_expiry_abs_share": front_gex / abs_total, } if net_total > 0: result["front_expiry_net_ratio"] = front_gex / net_total return result def build_oi_totals(df_full: pd.DataFrame) -> Dict: """FIX 14: OI sanity-check totals + DTE-bucket breakdown (full chain).""" if df_full.empty: return {"call_oi": 0, "put_oi": 0, "n_contracts": 0, "n_expiries": 0, "oi_by_dte_bucket": {}} calls = df_full[df_full["cp"] == "C"] puts = df_full[df_full["cp"] == "P"] buckets = {"0-7": 0, "8-30": 0, "31-90": 0, "91-365": 0} for _, r in df_full.iterrows(): oi = int(r["oi"]) dte = int(r["dte"]) if dte <= 7: buckets["0-7"] += oi elif dte <= 30: buckets["8-30"] += oi elif dte <= 90: buckets["31-90"] += oi else: buckets["91-365"] += oi return { "call_oi": int(calls["oi"].sum()), "put_oi": int(puts["oi"].sum()), "n_contracts": int(len(df_full)), "n_expiries": int(df_full["expiry"].nunique()), "oi_by_dte_bucket": buckets, } def build_gex_by_expiry(df_full: pd.DataFrame, spot: float, cfg: GexConfig) -> Dict: """FIX 18: per-expiry GEX concentration. Returns front_expiry_share = |net_gex from the nearest expiry| / Σ|net_gex| over all strikes, plus gex_by_expiry = [{expiry, dte, net_gex, share}] for every expiry. FIX 94: the field was previously named sum_abs_gex but it is |Σ GEX| (the absolute value of the net per-expiry GEX), NOT Σ|GEX| (the sum of absolute per-contract GEX). Renamed to net_gex to match what it actually measures. This exposes 0DTE domination: when one expiry drives most of the total |GEX|, the chart is really a single-expiry picture and should be read as such. FIX 25: the concentration warning must key on the DOMINANT expiry (argmax over gex_by_expiry), not the front one — in the exfront variant the front expiry is a tiny remnant while a later expiry holds the concentration. Also publishes max_expiry_share / max_expiry / max_expiry_dte and top3_expiry_share (sum of the three largest), so the warning fires on whichever expiry actually dominates. """ M = cfg.contract_multiplier S2 = spot * spot if df_full.empty: return {"front_expiry_share": 0.0, "front_expiry": None, "front_expiry_dte": None, "gex_by_expiry": [], "max_expiry_share": 0.0, "max_expiry": None, "max_expiry_dte": None, "top3_expiry_share": 0.0} df = df_full.copy() sign = np.where(df["cp"].to_numpy() == "C", 1.0, -1.0) df["gex"] = sign * df["gamma"].abs() * df["oi"] * M * S2 * 0.01 by_exp = df.groupby("expiry")["gex"].sum() per_strike_net = df.groupby("strike")["gex"].sum() denom = float(per_strike_net.abs().sum()) front_exp = df["expiry"].min() front_dte = int(df[df["expiry"] == front_exp]["dte"].min()) rows = [] for exp, g in by_exp.items(): dte = int(df[df["expiry"] == exp]["dte"].min()) abs_g = abs(float(g)) rows.append({"expiry": str(exp), "dte": dte, "net_gex": abs_g, "share": abs_g / denom if denom > 0 else 0.0}) rows.sort(key=lambda r: r["dte"]) front_share = abs(float(by_exp.get(front_exp, 0.0))) / denom if denom > 0 else 0.0 # FIX 25: dominant expiry (argmax by share) and top-3 concentration by_share = sorted(rows, key=lambda r: r["share"], reverse=True) max_row = by_share[0] if by_share else None top3_share = float(sum(r["share"] for r in by_share[:3])) return { "front_expiry_share": float(front_share), "front_expiry": str(front_exp), "front_expiry_dte": front_dte, "gex_by_expiry": rows, "max_expiry_share": float(max_row["share"]) if max_row else 0.0, "max_expiry": max_row["expiry"] if max_row else None, "max_expiry_dte": max_row["dte"] if max_row else None, "top3_expiry_share": top3_share, } def detect_spread_candidates(df_full: pd.DataFrame, spot: float, cfg: GexConfig, increment: float) -> Dict: """FIX 31 (revised FIX 33b, FIX 37): detect probable vertical-spread structures. Within each expiry and right (P/C), find strike pairs where: min(oi_a, oi_b) / max(oi_a, oi_b) >= cfg.spread_ratio_min, both OI >= max(cfg.spread_oi_min_abs, cfg.spread_oi_min_frac * chain_oi), |strike_a - strike_b| <= cfg.spread_max_width * increment. FIX 33b: greedy-dedupe — sort pairs by combined_abs_gex descending, accept only if NEITHER strike is already used (within that expiry+right). This prevents 527.5/530, 527.5/532.5, 530/532.5 all appearing and inflating the share. FIX 33a: sensitivity uses a SIGN FLIP of the smaller leg (gross scores it -g*OI; netted is +g*OI, so the adjustment is 2x), not removal. Returns spread_candidates, spread_flagged_share, and (if flagged) a sensitivity figure recomputing total_net_gex and HVL with the smaller leg sign-flipped. """ M = cfg.contract_multiplier S2 = spot * spot result = { "spread_candidates": [], "spread_flagged_share": 0.0, "sensitivity_smaller_leg_sign_flipped": None, } if df_full.empty or increment <= 0: return result # per-expiry, per-right OI by strike df = df_full.copy() sign = np.where(df["cp"].to_numpy() == "C", 1.0, -1.0) df["gex"] = sign * df["gamma"].abs() * df["oi"] * M * S2 * 0.01 total_abs_gex = float(df["gex"].abs().sum()) if total_abs_gex == 0: return result # FIX 37: relative OI floor chain_oi = float(df["oi"].sum()) oi_floor = max(cfg.spread_oi_min_abs, cfg.spread_oi_min_frac * chain_oi) candidates = [] max_width = cfg.spread_max_width * increment for (expiry, right), grp in df.groupby(["expiry", "cp"]): by_strike = grp.groupby("strike").agg( oi=("oi", "sum"), gex=("gex", "sum")).reset_index() by_strike = by_strike.sort_values("strike") strikes = by_strike["strike"].to_numpy() ois = by_strike["oi"].to_numpy() gexs = by_strike["gex"].to_numpy() for i in range(len(strikes)): for j in range(i + 1, len(strikes)): if strikes[j] - strikes[i] > max_width: break oi_a, oi_b = ois[i], ois[j] if oi_a < oi_floor or oi_b < oi_floor: continue ratio = min(oi_a, oi_b) / max(oi_a, oi_b) if ratio < cfg.spread_ratio_min: continue combined_abs = abs(gexs[i]) + abs(gexs[j]) candidates.append({ "expiry": str(expiry), "right": right, "strike_low": float(strikes[i]), "strike_high": float(strikes[j]), "oi_low": int(oi_a), "oi_high": int(oi_b), "combined_abs_gex": float(combined_abs), "share_of_total_abs_gex": float(combined_abs / total_abs_gex), }) # FIX 33b: greedy dedupe — sort by combined_abs_gex desc, accept only if # neither strike is already used within that expiry+right. candidates.sort(key=lambda c: c["combined_abs_gex"], reverse=True) deduped = [] used = {} # (expiry, right) -> set of used strikes for c in candidates: key = (c["expiry"], c["right"]) used_set = used.setdefault(key, set()) if c["strike_low"] in used_set or c["strike_high"] in used_set: continue used_set.add(c["strike_low"]) used_set.add(c["strike_high"]) deduped.append(c) result["spread_candidates"] = deduped flagged_share = sum(c["share_of_total_abs_gex"] for c in deduped) result["spread_flagged_share"] = float(flagged_share) if flagged_share > cfg.spread_flag_threshold and deduped: logger.warning( "SPREAD DETECTION: %d probable vertical-spread structures = %.0f%% of |GEX|; " "gross-OI proxy overstates net dealer gamma here.", len(deduped), flagged_share * 100) # FIX 33a: SENSITIVITY — sign-flip the smaller leg (adjustment is 2x removal). # Gross scores a put leg as -g*OI; netted (dealer-long) is +g*OI. # So the adjustment = +2 * |gex_smaller_leg| for puts, -2 * |gex_smaller_leg| for calls. df_netted = df.copy() for c in deduped: if c["oi_low"] <= c["oi_high"]: flip_strike, flip_oi = c["strike_low"], c["oi_low"] else: flip_strike, flip_oi = c["strike_high"], c["oi_high"] mask = ((df_netted["expiry"].astype(str) == c["expiry"]) & (df_netted["cp"] == c["right"]) & (df_netted["strike"] == flip_strike)) idx = df_netted.index[mask] if len(idx) > 0: total_oi_at_strike = df_netted.loc[idx, "oi"].sum() if total_oi_at_strike > 0: frac = min(flip_oi / total_oi_at_strike, 1.0) # Sign flip: negate the gex contribution of the flipped fraction. # Original gex = sign * |gamma| * oi * M * S2 * 0.01 # Flipped gex = -sign * |gamma| * oi * M * S2 * 0.01 # Adjustment = -2 * original gex for the flipped fraction orig_gex = df_netted.loc[idx, "gex"].sum() df_netted.loc[idx, "gex"] -= 2.0 * orig_gex * frac netted_total = float(df_netted["gex"].sum()) # recompute HVL on a quick grid from .config import GexConfig as _Cfg _cfg = _Cfg() n_pts = cfg.profile_grid_points grid_n = np.linspace(spot * (1 - _cfg.profile_band_hvl), spot * (1 + _cfg.profile_band_hvl), n_pts) gp_n, _ = gex_profile(df_netted, spot, cfg, grid_n) hvl_n_info = compute_hvl(grid_n, gp_n, spot, cfg) result["sensitivity_smaller_leg_sign_flipped"] = { "total_net_gex": netted_total, "hvl": hvl_n_info["hvl"], "hvl_distance_pct": hvl_n_info["hvl_distance_pct"], "note": "illustrative bound — smaller leg of each flagged spread sign-flipped; " "NOT the headline number", } return result def compute_realised_vol(closes: list, days: int = 20) -> Optional[float]: """FIX 32: annualised realised volatility from a list of daily closes. Uses log returns over the last `days` observations. Returns None if fewer than 5 data points. """ if len(closes) < 5: return None arr = np.array(closes[-days - 1:], dtype=float) # need days+1 for days returns if len(arr) < 2: return None log_ret = np.diff(np.log(arr)) if len(log_ret) == 0: return None return float(np.std(log_ret, ddof=1) * np.sqrt(252)) def atm_iv_cross_check(atm_iv: Optional[float], realised_vol: Optional[float], cfg: GexConfig) -> Dict: """FIX 32: cross-check ATM IV against realised volatility. Publishes iv_hv_ratio and vol_regime. If the ratio is an outlier (> 2.5 or < 0.4), flags atm_iv_status accordingly. """ result = { "realised_vol_20d": realised_vol, "iv_hv_ratio": None, "vol_regime": None, } if atm_iv is None or realised_vol is None or realised_vol <= 0: return result ratio = atm_iv / realised_vol result["iv_hv_ratio"] = float(ratio) result["vol_regime"] = "IV > HV" if ratio > 1.0 else "IV < HV" if ratio > cfg.iv_hv_outlier_hi or ratio < cfg.iv_hv_outlier_lo: result["iv_hv_outlier"] = True else: result["iv_hv_outlier"] = False return result def compute_bands(atm_iv: Optional[float], cfg: GexConfig, realised_vol: Optional[float] = None, atm_iv_status: Optional[str] = None) -> Dict: """FIX 35: derive volatility-scaled bands from the chain's own ATM IV. Two-pass: PASS 1 gets ATM IV (already computed by the caller). DERIVE bands. PASS 2 runs the pipeline with those bands. sigma_30d = atm_iv * sqrt(30/365) strike_band = clip(band_strike_mult * sigma_30d, limits) plot_band = clip(band_plot_mult * sigma_30d, limits) profile_band = clip(band_profile_mult * sigma_30d, limits) dex_band = clip(band_dex_mult * sigma_30d, limits) FIX 43: if atm_iv is None, fall back to realised vol (clip(rv*1.1, 0.10, 1.00)) instead of a hardcoded 0.30, and set band_basis="fallback_from_rv". Only if RV is also unavailable do we use cfg.atm_iv_fallback with band_basis="fallback". FIX 67: when atm_iv resolved from a single strike (no interpolation), set band_basis="atm_iv_single_strike" so the audit trail is honest about the input. """ if atm_iv is not None and 0.01 < atm_iv < 3.0: iv_used = atm_iv basis = "atm_iv_single_strike" if atm_iv_status == "single_strike_no_interpolation" else "atm_iv" elif realised_vol is not None and realised_vol > 0: iv_used = float(np.clip(realised_vol * 1.1, 0.10, 1.00)) basis = "fallback_from_rv" else: iv_used = cfg.atm_iv_fallback basis = "fallback" sigma_30d = iv_used * np.sqrt(30.0 / 365.0) lim = cfg.band_limits def _clip(val, key): lo, hi = lim[key] return float(np.clip(val, lo, hi)) strike_band = _clip(cfg.band_strike_mult * sigma_30d, "strike_band") plot_band = _clip(cfg.band_plot_mult * sigma_30d, "plot_band") profile_band = _clip(cfg.band_profile_mult * sigma_30d, "profile_band") dex_band = _clip(cfg.band_dex_mult * sigma_30d, "dex_band") return { "atm_iv_used": float(iv_used), "sigma_30d": float(sigma_30d), "strike_band": strike_band, "plot_band": plot_band, "profile_band": profile_band, "dex_band": dex_band, "band_basis": basis, } def snap_to_increment(value: Optional[float], increment: float) -> Optional[float]: """FIX 36: round a published level to the detected strike increment.""" if value is None or increment <= 0: return value return round(value / increment) * increment def is_third_friday(d: date) -> bool: """FIX 37: check if a date is the third Friday of its month (monthly index expiry).""" if d.weekday() != 4: # not a Friday return False # Third Friday: day-of-month is 15-21 return 15 <= d.day <= 21