#!/usr/bin/env python3
"""
MM Negative Gamma Magnet Effect Research
Tests whether ES price is drawn toward or repelled from strikes where MMs are most short gamma.
"""

import pandas as pd
import numpy as np
import json
import glob
import os
import warnings
from datetime import timedelta
from collections import defaultdict

warnings.filterwarnings('ignore')

DATA_DIR = '/Users/daniel/.openclaw/workspace/data'
TRACE_DIR = os.path.join(DATA_DIR, 'trace_api')
OUTPUT_JSON = os.path.join(DATA_DIR, 'mm_gamma_magnet_research.json')
OUTPUT_MD = os.path.join(DATA_DIR, 'mm_gamma_magnet_research.md')

# RTH hours in ET
RTH_START_HOUR = 9
RTH_START_MIN = 30
RTH_END_HOUR = 16
RTH_END_MIN = 0

###############################################################################
# 1. Load ES price data
###############################################################################
print("Loading ES 1-min bars...")
es = pd.read_csv(os.path.join(DATA_DIR, 'es_1min_bars.csv'), parse_dates=['ts_event'])
es = es.rename(columns={'ts_event': 'timestamp'})
# Convert to ET
es['timestamp'] = es['timestamp'].dt.tz_convert('US/Eastern') if es['timestamp'].dt.tz is not None else es['timestamp'].dt.tz_localize('UTC').dt.tz_convert('US/Eastern')
es = es.set_index('timestamp').sort_index()
# Keep only RTH
es_rth = es.between_time('09:30', '16:00')
print(f"  ES bars: {len(es):,} total, {len(es_rth):,} RTH, date range: {es.index[0].date()} to {es.index[-1].date()}")

###############################################################################
# 2. Load all TRACE intradayStrikeGEX files
###############################################################################
print("Loading TRACE intradayStrikeGEX files...")
trace_files = sorted(glob.glob(os.path.join(TRACE_DIR, 'intradayStrikeGEX_*.csv')))
print(f"  Found {len(trace_files)} files")

all_trace = []
for f in trace_files:
    df = pd.read_csv(f, parse_dates=['timestamp'])
    all_trace.append(df)
    
trace = pd.concat(all_trace, ignore_index=True)
trace['timestamp'] = pd.to_datetime(trace['timestamp'], utc=True).dt.tz_convert('US/Eastern')
trace['date'] = trace['timestamp'].dt.date
print(f"  TRACE rows: {len(trace):,}, dates: {trace['date'].min()} to {trace['date'].max()}")
print(f"  Unique dates: {trace['date'].nunique()}")

###############################################################################
# 3. Build ES price lookup (round to nearest minute)
###############################################################################
print("Building ES price lookup...")
# Create a minute-level close price series
es_close = es['close'].copy()
es_close.index = es_close.index.floor('min')
es_close = es_close[~es_close.index.duplicated(keep='last')]

# Also load intraday_gex_5min for spot reference
gex5 = pd.read_csv(os.path.join(DATA_DIR, 'intraday_gex_5min.csv'), parse_dates=['timestamp'])
gex5['timestamp'] = pd.to_datetime(gex5['timestamp'], utc=True).dt.tz_convert('US/Eastern')
gex5 = gex5.set_index('timestamp').sort_index()
gex5_spot = gex5['spot'].copy()
gex5_spot.index = gex5_spot.index.floor('min')
gex5_spot = gex5_spot[~gex5_spot.index.duplicated(keep='last')]

def get_es_price(ts):
    """Get ES price at timestamp, trying es_1min first, then gex5min spot."""
    ts_floor = ts.floor('min')
    if ts_floor in es_close.index:
        return es_close[ts_floor]
    # Try within 1 minute
    mask = es_close.index[(es_close.index >= ts_floor - pd.Timedelta(minutes=1)) & 
                           (es_close.index <= ts_floor + pd.Timedelta(minutes=1))]
    if len(mask) > 0:
        return es_close[mask[-1]]
    # Try gex5
    if ts_floor in gex5_spot.index:
        return gex5_spot[ts_floor]
    mask2 = gex5_spot.index[(gex5_spot.index >= ts_floor - pd.Timedelta(minutes=5)) & 
                             (gex5_spot.index <= ts_floor + pd.Timedelta(minutes=5))]
    if len(mask2) > 0:
        return gex5_spot[mask2[-1]]
    return np.nan

###############################################################################
# 4. For each TRACE timestamp, compute signals and forward returns
###############################################################################
print("Computing signals for each TRACE snapshot...")

# Get unique timestamps during RTH
trace_rth = trace[(trace['timestamp'].dt.hour > RTH_START_HOUR) | 
                  ((trace['timestamp'].dt.hour == RTH_START_HOUR) & (trace['timestamp'].dt.minute >= RTH_START_MIN))]
trace_rth = trace_rth[(trace_rth['timestamp'].dt.hour < RTH_END_HOUR) | 
                       ((trace_rth['timestamp'].dt.hour == RTH_END_HOUR) & (trace_rth['timestamp'].dt.minute == 0))]

unique_ts = sorted(trace_rth['timestamp'].unique())
print(f"  RTH TRACE timestamps: {len(unique_ts):,}")

# Build EOD price lookup per date
eod_prices = {}
for d in trace_rth['date'].unique():
    d_pd = pd.Timestamp(d)
    eod_ts = pd.Timestamp(f"{d} 16:00:00", tz='US/Eastern')
    # Find nearest ES price to 16:00
    mask = es_close.index[(es_close.index >= eod_ts - pd.Timedelta(minutes=5)) & 
                           (es_close.index <= eod_ts + pd.Timedelta(minutes=5))]
    if len(mask) > 0:
        eod_prices[d] = es_close[mask[-1]]
    else:
        # Try 15:59
        mask2 = es_close.index[(es_close.index >= eod_ts - pd.Timedelta(minutes=10)) & 
                                (es_close.index <= eod_ts)]
        if len(mask2) > 0:
            eod_prices[d] = es_close[mask2[-1]]

# Pre-build a dict of TRACE data by timestamp for fast lookup
trace_by_ts = {}
for ts, grp in trace_rth.groupby('timestamp'):
    trace_by_ts[ts] = grp

# Process each timestamp
results = []
processed = 0
skipped = 0

for ts in unique_ts:
    grp = trace_by_ts[ts]
    d = ts.date() if hasattr(ts, 'date') else pd.Timestamp(ts).date()
    
    # Get current ES price
    current_price = get_es_price(pd.Timestamp(ts))
    if np.isnan(current_price) if isinstance(current_price, float) else pd.isna(current_price):
        skipped += 1
        continue
    
    # Get forward prices
    ts_pd = pd.Timestamp(ts)
    
    # 1H forward
    ts_1h = ts_pd + pd.Timedelta(hours=1)
    price_1h = get_es_price(ts_1h)
    
    # 3H forward
    ts_3h = ts_pd + pd.Timedelta(hours=3)
    price_3h = get_es_price(ts_3h)
    
    # EOD
    price_eod = eod_prices.get(d, np.nan)
    
    # Also get high/low between now and EOD for touch detection
    eod_ts = pd.Timestamp(f"{d} 16:00:00", tz='US/Eastern')
    future_bars = es_rth.loc[ts_pd:eod_ts]
    future_high = future_bars['high'].max() if len(future_bars) > 0 else np.nan
    future_low = future_bars['low'].min() if len(future_bars) > 0 else np.nan
    
    # === TEST 1: Biggest negative MM gamma strike ===
    neg_gamma = grp[grp['mm_gamma'] < 0].copy()
    if len(neg_gamma) == 0:
        skipped += 1
        continue
    
    most_neg_idx = neg_gamma['mm_gamma'].idxmin()
    most_neg_strike = neg_gamma.loc[most_neg_idx, 'strike_price']
    most_neg_gamma = neg_gamma.loc[most_neg_idx, 'mm_gamma']
    
    dist_to_most_neg = most_neg_strike - current_price  # positive = above
    abs_dist = abs(dist_to_most_neg)
    
    # Direction toward magnet
    magnet_direction = 1 if most_neg_strike > current_price else -1  # 1=need to go up, -1=need to go down
    
    # Forward moves
    move_1h = (price_1h - current_price) if not pd.isna(price_1h) else np.nan
    move_3h = (price_3h - current_price) if not pd.isna(price_3h) else np.nan
    move_eod = (price_eod - current_price) if not pd.isna(price_eod) else np.nan
    
    # Did price move toward the magnet?
    toward_1h = (move_1h * magnet_direction > 0) if not pd.isna(move_1h) else np.nan
    toward_3h = (move_3h * magnet_direction > 0) if not pd.isna(move_3h) else np.nan
    toward_eod = (move_eod * magnet_direction > 0) if not pd.isna(move_eod) else np.nan
    
    # Did price touch the magnet strike (within 5 pts)?
    touched = False
    if not pd.isna(future_high) and not pd.isna(future_low):
        touched = (future_low <= most_neg_strike + 5) and (future_high >= most_neg_strike - 5)
    
    # Distance reduction
    dist_1h = abs(price_1h - most_neg_strike) if not pd.isna(price_1h) else np.nan
    dist_3h = abs(price_3h - most_neg_strike) if not pd.isna(price_3h) else np.nan
    dist_eod = abs(price_eod - most_neg_strike) if not pd.isna(price_eod) else np.nan
    
    closer_1h = (dist_1h < abs_dist) if not pd.isna(dist_1h) else np.nan
    closer_3h = (dist_3h < abs_dist) if not pd.isna(dist_3h) else np.nan
    closer_eod = (dist_eod < abs_dist) if not pd.isna(dist_eod) else np.nan
    
    # === TEST 2: Top 5 negative gamma cluster ===
    top5_neg = neg_gamma.nsmallest(5, 'mm_gamma')
    top5_strikes = top5_neg['strike_price'].values
    top5_gammas = top5_neg['mm_gamma'].values
    
    # Gravity center (weighted by magnitude of negative gamma)
    weights = np.abs(top5_gammas)
    gravity_center = np.average(top5_strikes, weights=weights)
    
    # Cluster spread
    cluster_spread = top5_strikes.max() - top5_strikes.min()
    is_clustered = cluster_spread <= 25
    
    # Direction toward gravity center
    gravity_dist = gravity_center - current_price
    gravity_direction = 1 if gravity_center > current_price else -1
    
    toward_gravity_1h = (move_1h * gravity_direction > 0) if not pd.isna(move_1h) else np.nan
    toward_gravity_3h = (move_3h * gravity_direction > 0) if not pd.isna(move_3h) else np.nan
    toward_gravity_eod = (move_eod * gravity_direction > 0) if not pd.isna(move_eod) else np.nan
    
    closer_gravity_1h = (abs(price_1h - gravity_center) < abs(gravity_dist)) if not pd.isna(price_1h) else np.nan
    closer_gravity_3h = (abs(price_3h - gravity_center) < abs(gravity_dist)) if not pd.isna(price_3h) else np.nan
    closer_gravity_eod = (abs(price_eod - gravity_center) < abs(gravity_dist)) if not pd.isna(price_eod) else np.nan
    
    # Touch gravity center
    touched_gravity = False
    if not pd.isna(future_high) and not pd.isna(future_low):
        touched_gravity = (future_low <= gravity_center + 5) and (future_high >= gravity_center - 5)
    
    # === TEST 3: 0DTE overlay ===
    neg_gamma_0dte = grp[grp['mm_gamma_0'] < 0].copy()
    has_0dte_overlap = False
    most_neg_0dte_strike = np.nan
    if len(neg_gamma_0dte) > 0:
        most_neg_0dte_idx = neg_gamma_0dte['mm_gamma_0'].idxmin()
        most_neg_0dte_strike = neg_gamma_0dte.loc[most_neg_0dte_idx, 'strike_price']
        # Check if total mm_gamma is also negative at this strike
        row_at_strike = grp[grp['strike_price'] == most_neg_0dte_strike]
        if len(row_at_strike) > 0 and row_at_strike['mm_gamma'].values[0] < 0:
            has_0dte_overlap = True
    
    # 0DTE magnet
    if not pd.isna(most_neg_0dte_strike):
        dte0_direction = 1 if most_neg_0dte_strike > current_price else -1
        toward_0dte_1h = (move_1h * dte0_direction > 0) if not pd.isna(move_1h) else np.nan
        toward_0dte_eod = (move_eod * dte0_direction > 0) if not pd.isna(move_eod) else np.nan
        dte0_dist = abs(most_neg_0dte_strike - current_price)
        closer_0dte_eod = (abs(price_eod - most_neg_0dte_strike) < dte0_dist) if not pd.isna(price_eod) else np.nan
    else:
        toward_0dte_1h = np.nan
        toward_0dte_eod = np.nan
        closer_0dte_eod = np.nan
    
    # === TEST 4: Nearest negative gamma wall ===
    nearby_neg = neg_gamma[(neg_gamma['strike_price'] >= current_price - 50) & 
                            (neg_gamma['strike_price'] <= current_price + 50)].copy()
    nearest_above = np.nan
    nearest_below = np.nan
    nearest_above_gamma = 0
    nearest_below_gamma = 0
    
    if len(nearby_neg) > 0:
        above = nearby_neg[nearby_neg['strike_price'] > current_price]
        below = nearby_neg[nearby_neg['strike_price'] < current_price]
        
        if len(above) > 0:
            # Nearest above with significant negative gamma
            above_sorted = above.sort_values('strike_price')
            nearest_above = above_sorted.iloc[0]['strike_price']
            nearest_above_gamma = above_sorted.iloc[0]['mm_gamma']
        
        if len(below) > 0:
            below_sorted = below.sort_values('strike_price', ascending=False)
            nearest_below = below_sorted.iloc[0]['strike_price']
            nearest_below_gamma = below_sorted.iloc[0]['mm_gamma']
    
    # Which nearest wall is stronger?
    nearest_wall_direction = np.nan
    if not pd.isna(nearest_above) and not pd.isna(nearest_below):
        # More negative = stronger magnet
        if nearest_above_gamma < nearest_below_gamma:
            nearest_wall_direction = 1  # above is stronger
        else:
            nearest_wall_direction = -1  # below is stronger
    elif not pd.isna(nearest_above):
        nearest_wall_direction = 1
    elif not pd.isna(nearest_below):
        nearest_wall_direction = -1
    
    toward_nearest_wall_1h = np.nan
    toward_nearest_wall_eod = np.nan
    if not pd.isna(nearest_wall_direction):
        toward_nearest_wall_1h = (move_1h * nearest_wall_direction > 0) if not pd.isna(move_1h) else np.nan
        toward_nearest_wall_eod = (move_eod * nearest_wall_direction > 0) if not pd.isna(move_eod) else np.nan
    
    # === TEST 5: Gamma gravity score ===
    # Sum of mm_gamma / distance^2 for all strikes with negative gamma
    gravity_score = 0
    for _, row in neg_gamma.iterrows():
        d_strike = row['strike_price'] - current_price
        if abs(d_strike) > 1:  # avoid division by zero
            # Negative gamma * sign(direction) / distance^2
            # If strike is above price and gamma is negative: pulls UP (positive contribution)
            gravity_score += row['mm_gamma'] / (d_strike * abs(d_strike))  # preserves sign of direction, scales by dist^2
    
    # Reinterpret: gravity_score > 0 means net pull DOWN (negative gamma above repels? no...)
    # Actually: mm_gamma is negative. If strike is above (d_strike > 0): 
    #   contribution = negative / positive = negative → suggests pull down? 
    # Hmm, let me think about this differently.
    # The magnet theory says negative gamma ATTRACTS. So strikes with neg gamma above should pull UP.
    # Let's define: magnet_pull = sum of |mm_gamma| * sign(strike - price) / distance^2
    # Positive magnet_pull = net attraction upward
    
    magnet_pull = 0
    for _, row in neg_gamma.iterrows():
        d_strike = row['strike_price'] - current_price
        if abs(d_strike) > 1:
            magnet_pull += abs(row['mm_gamma']) * np.sign(d_strike) / (d_strike ** 2)
    
    # magnet_pull > 0 → net pull UP, < 0 → net pull DOWN
    gravity_predicts_up = magnet_pull > 0
    gravity_correct_1h = ((move_1h > 0) == gravity_predicts_up) if not pd.isna(move_1h) else np.nan
    gravity_correct_3h = ((move_3h > 0) == gravity_predicts_up) if not pd.isna(move_3h) else np.nan
    gravity_correct_eod = ((move_eod > 0) == gravity_predicts_up) if not pd.isna(move_eod) else np.nan
    
    # === TEST 7: Participant overlay ===
    # Firm AND MM both negative at same strike
    firm_mm_neg = grp[(grp['mm_gamma'] < 0) & (grp['firm_gamma'] < 0)]
    has_firm_mm_overlap = len(firm_mm_neg) > 0
    
    # Customer positive where MM negative (more OI = stronger magnet)
    cust_pos_mm_neg = grp[(grp['mm_gamma'] < 0) & (grp['cust_gamma'] > 0)]
    has_cust_overlay = len(cust_pos_mm_neg) > 0
    
    # Time of day bucket
    hour = pd.Timestamp(ts).hour
    minute = pd.Timestamp(ts).minute
    time_minutes = hour * 60 + minute
    if time_minutes < 11 * 60:  # before 11am ET
        tod = 'morning'
    elif time_minutes < 14 * 60:  # before 2pm ET
        tod = 'midday'
    else:
        tod = 'afternoon'
    
    results.append({
        'timestamp': str(ts),
        'date': str(d),
        'current_price': float(current_price),
        'time_of_day': tod,
        'hour': hour,
        
        # Test 1
        'most_neg_strike': float(most_neg_strike),
        'most_neg_gamma': float(most_neg_gamma),
        'dist_to_most_neg': float(dist_to_most_neg),
        'abs_dist': float(abs_dist),
        'magnet_above': most_neg_strike > current_price,
        'toward_1h': toward_1h,
        'toward_3h': toward_3h,
        'toward_eod': toward_eod,
        'closer_1h': closer_1h,
        'closer_3h': closer_3h,
        'closer_eod': closer_eod,
        'touched_eod': touched,
        'move_1h': float(move_1h) if not pd.isna(move_1h) else None,
        'move_3h': float(move_3h) if not pd.isna(move_3h) else None,
        'move_eod': float(move_eod) if not pd.isna(move_eod) else None,
        
        # Test 2
        'gravity_center': float(gravity_center),
        'cluster_spread': float(cluster_spread),
        'is_clustered': bool(is_clustered),
        'toward_gravity_1h': toward_gravity_1h,
        'toward_gravity_3h': toward_gravity_3h,
        'toward_gravity_eod': toward_gravity_eod,
        'closer_gravity_1h': closer_gravity_1h,
        'closer_gravity_3h': closer_gravity_3h,
        'closer_gravity_eod': closer_gravity_eod,
        'touched_gravity_eod': touched_gravity,
        
        # Test 3
        'has_0dte_overlap': has_0dte_overlap,
        'most_neg_0dte_strike': float(most_neg_0dte_strike) if not pd.isna(most_neg_0dte_strike) else None,
        'toward_0dte_1h': toward_0dte_1h,
        'toward_0dte_eod': toward_0dte_eod,
        'closer_0dte_eod': closer_0dte_eod,
        
        # Test 4
        'nearest_above': float(nearest_above) if not pd.isna(nearest_above) else None,
        'nearest_below': float(nearest_below) if not pd.isna(nearest_below) else None,
        'nearest_wall_direction': float(nearest_wall_direction) if not pd.isna(nearest_wall_direction) else None,
        'toward_nearest_wall_1h': toward_nearest_wall_1h,
        'toward_nearest_wall_eod': toward_nearest_wall_eod,
        
        # Test 5
        'magnet_pull': float(magnet_pull),
        'gravity_predicts_up': bool(gravity_predicts_up),
        'gravity_correct_1h': gravity_correct_1h,
        'gravity_correct_3h': gravity_correct_3h,
        'gravity_correct_eod': gravity_correct_eod,
        
        # Test 7
        'has_firm_mm_overlap': has_firm_mm_overlap,
        'has_cust_overlay': has_cust_overlay,
    })
    
    processed += 1
    if processed % 500 == 0:
        print(f"  Processed {processed}/{len(unique_ts)} timestamps...")

print(f"  Done. Processed: {processed}, Skipped: {skipped}")

###############################################################################
# 5. Analyze results
###############################################################################
print("\n" + "="*80)
print("ANALYZING RESULTS")
print("="*80)

df = pd.DataFrame(results)
df['timestamp'] = pd.to_datetime(df['timestamp'])

# Convert boolean columns properly (they may have np.nan mixed in)
bool_cols = ['toward_1h', 'toward_3h', 'toward_eod', 'closer_1h', 'closer_3h', 'closer_eod',
             'toward_gravity_1h', 'toward_gravity_3h', 'toward_gravity_eod',
             'closer_gravity_1h', 'closer_gravity_3h', 'closer_gravity_eod',
             'toward_0dte_1h', 'toward_0dte_eod', 'closer_0dte_eod',
             'toward_nearest_wall_1h', 'toward_nearest_wall_eod',
             'gravity_correct_1h', 'gravity_correct_3h', 'gravity_correct_eod']

# IS/OOS split by date
dates = sorted(df['date'].unique())
n_dates = len(dates)
is_cutoff = dates[int(n_dates * 2/3)]  # ~120 IS, ~60 OOS
print(f"\nTotal dates: {n_dates}, IS cutoff: {is_cutoff}")
print(f"IS: {dates[0]} to {is_cutoff}, OOS: next to {dates[-1]}")

df['is_oos'] = df['date'].apply(lambda x: 'IS' if x <= is_cutoff else 'OOS')

analysis = {}

def safe_mean(series):
    """Compute mean of boolean series ignoring NaN."""
    s = series.dropna()
    if len(s) == 0:
        return np.nan
    return float(s.astype(float).mean())

def safe_count(series):
    return int(series.dropna().sum()) if len(series.dropna()) > 0 else 0

def safe_n(series):
    return int(series.dropna().count())

def compute_tstat(series, null=0.5):
    """t-stat vs null hypothesis of 50% hit rate."""
    s = series.dropna().astype(float)
    n = len(s)
    if n < 2:
        return 0
    p = s.mean()
    se = np.sqrt(p * (1-p) / n) if p > 0 and p < 1 else 1e-6
    return (p - null) / se

###############################################################################
# TEST 1: MM Negative Gamma Magnet
###############################################################################
print("\n--- TEST 1: MM Negative Gamma Magnet ---")

for split in ['IS', 'OOS', 'ALL']:
    sub = df if split == 'ALL' else df[df['is_oos'] == split]
    n = len(sub)
    
    t1 = {
        'n': n,
        'toward_1h_rate': safe_mean(sub['toward_1h']),
        'toward_3h_rate': safe_mean(sub['toward_3h']),
        'toward_eod_rate': safe_mean(sub['toward_eod']),
        'closer_1h_rate': safe_mean(sub['closer_1h']),
        'closer_3h_rate': safe_mean(sub['closer_3h']),
        'closer_eod_rate': safe_mean(sub['closer_eod']),
        'touch_rate_eod': safe_mean(sub['touched_eod']),
        'toward_1h_tstat': compute_tstat(sub['toward_1h']),
        'toward_eod_tstat': compute_tstat(sub['toward_eod']),
        'closer_eod_tstat': compute_tstat(sub['closer_eod']),
        'avg_dist_pts': float(sub['abs_dist'].mean()),
        'median_dist_pts': float(sub['abs_dist'].median()),
        'pct_magnet_above': float(sub['magnet_above'].mean()),
    }
    
    print(f"\n  [{split}] N={n}")
    print(f"    Toward magnet: 1H={t1['toward_1h_rate']:.1%}, 3H={t1['toward_3h_rate']:.1%}, EOD={t1['toward_eod_rate']:.1%}")
    print(f"    Closer to magnet: 1H={t1['closer_1h_rate']:.1%}, 3H={t1['closer_3h_rate']:.1%}, EOD={t1['closer_eod_rate']:.1%}")
    print(f"    Touch rate (EOD): {t1['touch_rate_eod']:.1%}")
    print(f"    t-stat toward EOD: {t1['toward_eod_tstat']:.2f}, closer EOD: {t1['closer_eod_tstat']:.2f}")
    print(f"    Avg distance: {t1['avg_dist_pts']:.1f} pts, Magnet above: {t1['pct_magnet_above']:.1%}")
    
    analysis[f'test1_{split}'] = t1

###############################################################################
# TEST 2: Gravity Center (Cluster)
###############################################################################
print("\n--- TEST 2: Gravity Center (Cluster Analysis) ---")

for split in ['IS', 'OOS', 'ALL']:
    sub = df if split == 'ALL' else df[df['is_oos'] == split]
    
    t2 = {
        'toward_gravity_1h': safe_mean(sub['toward_gravity_1h']),
        'toward_gravity_3h': safe_mean(sub['toward_gravity_3h']),
        'toward_gravity_eod': safe_mean(sub['toward_gravity_eod']),
        'closer_gravity_eod': safe_mean(sub['closer_gravity_eod']),
        'touch_gravity_eod': safe_mean(sub['touched_gravity_eod']),
        'toward_gravity_eod_tstat': compute_tstat(sub['toward_gravity_eod']),
    }
    
    # Clustered vs dispersed
    clustered = sub[sub['is_clustered'] == True]
    dispersed = sub[sub['is_clustered'] == False]
    t2['clustered_toward_eod'] = safe_mean(clustered['toward_gravity_eod'])
    t2['dispersed_toward_eod'] = safe_mean(dispersed['toward_gravity_eod'])
    t2['clustered_closer_eod'] = safe_mean(clustered['closer_gravity_eod'])
    t2['dispersed_closer_eod'] = safe_mean(dispersed['closer_gravity_eod'])
    t2['pct_clustered'] = float(sub['is_clustered'].mean())
    t2['n_clustered'] = int(clustered.shape[0])
    t2['n_dispersed'] = int(dispersed.shape[0])
    
    print(f"\n  [{split}] N={len(sub)}")
    print(f"    Toward gravity: 1H={t2['toward_gravity_1h']:.1%}, EOD={t2['toward_gravity_eod']:.1%} (t={t2['toward_gravity_eod_tstat']:.2f})")
    print(f"    Closer to gravity EOD: {t2['closer_gravity_eod']:.1%}")
    print(f"    Touch rate: {t2['touch_gravity_eod']:.1%}")
    print(f"    Clustered ({t2['n_clustered']}): toward={t2['clustered_toward_eod']:.1%}, closer={t2['clustered_closer_eod']:.1%}")
    print(f"    Dispersed ({t2['n_dispersed']}): toward={t2['dispersed_toward_eod']:.1%}, closer={t2['dispersed_closer_eod']:.1%}")
    
    analysis[f'test2_{split}'] = t2

###############################################################################
# TEST 3: 0DTE Overlay
###############################################################################
print("\n--- TEST 3: 0DTE MM Gamma Overlay ---")

for split in ['IS', 'OOS', 'ALL']:
    sub = df if split == 'ALL' else df[df['is_oos'] == split]
    
    with_0dte = sub[sub['has_0dte_overlap'] == True]
    without_0dte = sub[sub['has_0dte_overlap'] == False]
    
    t3 = {
        'pct_with_0dte_overlap': float(sub['has_0dte_overlap'].mean()),
        'with_0dte_toward_1h': safe_mean(with_0dte['toward_0dte_1h']),
        'with_0dte_toward_eod': safe_mean(with_0dte['toward_0dte_eod']),
        'with_0dte_closer_eod': safe_mean(with_0dte['closer_0dte_eod']),
        'without_0dte_toward_eod': safe_mean(without_0dte['toward_eod']),
        'n_with': int(len(with_0dte)),
        'n_without': int(len(without_0dte)),
        # Compare: when 0DTE overlaps, does the total mm_gamma magnet work better?
        'overlap_total_toward_eod': safe_mean(with_0dte['toward_eod']),
        'no_overlap_total_toward_eod': safe_mean(without_0dte['toward_eod']),
    }
    
    print(f"\n  [{split}] N={len(sub)}, {t3['pct_with_0dte_overlap']:.1%} with 0DTE overlap")
    print(f"    With 0DTE overlap ({t3['n_with']}): toward_0dte={t3['with_0dte_toward_eod']:.1%}, closer={t3['with_0dte_closer_eod']:.1%}")
    print(f"    With overlap → total magnet toward EOD: {t3['overlap_total_toward_eod']:.1%}")
    print(f"    Without overlap → total magnet toward EOD: {t3['no_overlap_total_toward_eod']:.1%}")
    
    analysis[f'test3_{split}'] = t3

###############################################################################
# TEST 4: Nearest Negative Gamma Wall
###############################################################################
print("\n--- TEST 4: Nearest Negative Gamma Wall ---")

for split in ['IS', 'OOS', 'ALL']:
    sub = df if split == 'ALL' else df[df['is_oos'] == split]
    has_wall = sub[sub['nearest_wall_direction'].notna()]
    
    t4 = {
        'n_with_nearby_wall': int(len(has_wall)),
        'toward_nearest_1h': safe_mean(has_wall['toward_nearest_wall_1h']),
        'toward_nearest_eod': safe_mean(has_wall['toward_nearest_wall_eod']),
        'toward_nearest_eod_tstat': compute_tstat(has_wall['toward_nearest_wall_eod']),
    }
    
    # Above vs below
    wall_up = has_wall[has_wall['nearest_wall_direction'] == 1]
    wall_down = has_wall[has_wall['nearest_wall_direction'] == -1]
    t4['wall_above_toward_eod'] = safe_mean(wall_up['toward_nearest_wall_eod'])
    t4['wall_below_toward_eod'] = safe_mean(wall_down['toward_nearest_wall_eod'])
    t4['n_wall_above'] = int(len(wall_up))
    t4['n_wall_below'] = int(len(wall_down))
    
    print(f"\n  [{split}] N with nearby wall={t4['n_with_nearby_wall']}")
    print(f"    Toward nearest wall: 1H={t4['toward_nearest_1h']:.1%}, EOD={t4['toward_nearest_eod']:.1%} (t={t4['toward_nearest_eod_tstat']:.2f})")
    print(f"    Wall above ({t4['n_wall_above']}): toward EOD={t4['wall_above_toward_eod']:.1%}")
    print(f"    Wall below ({t4['n_wall_below']}): toward EOD={t4['wall_below_toward_eod']:.1%}")
    
    analysis[f'test4_{split}'] = t4

###############################################################################
# TEST 5: Gamma Gravity Score
###############################################################################
print("\n--- TEST 5: Gamma Gravity Score ---")

for split in ['IS', 'OOS', 'ALL']:
    sub = df if split == 'ALL' else df[df['is_oos'] == split]
    
    t5 = {
        'gravity_correct_1h': safe_mean(sub['gravity_correct_1h']),
        'gravity_correct_3h': safe_mean(sub['gravity_correct_3h']),
        'gravity_correct_eod': safe_mean(sub['gravity_correct_eod']),
        'gravity_correct_eod_tstat': compute_tstat(sub['gravity_correct_eod']),
    }
    
    # Quintile analysis
    sub_valid = sub[sub['magnet_pull'].notna() & sub['move_eod'].notna()].copy()
    if len(sub_valid) > 50:
        sub_valid['pull_quintile'] = pd.qcut(sub_valid['magnet_pull'], 5, labels=[1,2,3,4,5], duplicates='drop')
        q_means = sub_valid.groupby('pull_quintile')['move_eod'].mean()
        t5['quintile_returns_eod'] = {str(k): float(v) for k, v in q_means.items()}
        # Monotonicity: Q5 (strong up pull) should have highest return
        q_vals = [q_means.get(i, 0) for i in range(1, 6)]
        if len(q_vals) == 5:
            mono = sum(1 for i in range(4) if q_vals[i+1] > q_vals[i])
            t5['monotonicity'] = mono / 4  # 1.0 = perfect
        else:
            t5['monotonicity'] = None
    else:
        t5['quintile_returns_eod'] = {}
        t5['monotonicity'] = None
    
    print(f"\n  [{split}] N={len(sub)}")
    print(f"    Gravity correct: 1H={t5['gravity_correct_1h']:.1%}, 3H={t5['gravity_correct_3h']:.1%}, EOD={t5['gravity_correct_eod']:.1%} (t={t5['gravity_correct_eod_tstat']:.2f})")
    if t5['quintile_returns_eod']:
        print(f"    Quintile EOD returns: {t5['quintile_returns_eod']}")
        print(f"    Monotonicity: {t5['monotonicity']}")
    
    analysis[f'test5_{split}'] = t5

###############################################################################
# TEST 6: Time of Day
###############################################################################
print("\n--- TEST 6: Time of Day Effect ---")

for split in ['IS', 'OOS', 'ALL']:
    sub = df if split == 'ALL' else df[df['is_oos'] == split]
    
    t6 = {}
    for tod in ['morning', 'midday', 'afternoon']:
        tod_sub = sub[sub['time_of_day'] == tod]
        t6[tod] = {
            'n': int(len(tod_sub)),
            'toward_1h': safe_mean(tod_sub['toward_1h']),
            'toward_eod': safe_mean(tod_sub['toward_eod']),
            'closer_eod': safe_mean(tod_sub['closer_eod']),
            'touch_eod': safe_mean(tod_sub['touched_eod']),
            'gravity_correct_1h': safe_mean(tod_sub['gravity_correct_1h']),
            'gravity_correct_eod': safe_mean(tod_sub['gravity_correct_eod']),
            'toward_eod_tstat': compute_tstat(tod_sub['toward_eod']),
        }
        print(f"\n  [{split}] {tod} (N={t6[tod]['n']})")
        print(f"    Toward magnet: 1H={t6[tod]['toward_1h']:.1%}, EOD={t6[tod]['toward_eod']:.1%} (t={t6[tod]['toward_eod_tstat']:.2f})")
        print(f"    Closer EOD: {t6[tod]['closer_eod']:.1%}, Touch: {t6[tod]['touch_eod']:.1%}")
        print(f"    Gravity correct: 1H={t6[tod]['gravity_correct_1h']:.1%}, EOD={t6[tod]['gravity_correct_eod']:.1%}")
    
    analysis[f'test6_{split}'] = t6

###############################################################################
# TEST 7: Participant Overlay
###############################################################################
print("\n--- TEST 7: Participant Overlay ---")

for split in ['IS', 'OOS', 'ALL']:
    sub = df if split == 'ALL' else df[df['is_oos'] == split]
    
    firm_mm = sub[sub['has_firm_mm_overlap'] == True]
    no_firm = sub[sub['has_firm_mm_overlap'] == False]
    cust_overlay = sub[sub['has_cust_overlay'] == True]
    no_cust = sub[sub['has_cust_overlay'] == False]
    
    t7 = {
        'pct_firm_mm_overlap': float(sub['has_firm_mm_overlap'].mean()),
        'firm_mm_toward_eod': safe_mean(firm_mm['toward_eod']),
        'no_firm_toward_eod': safe_mean(no_firm['toward_eod']),
        'firm_mm_closer_eod': safe_mean(firm_mm['closer_eod']),
        'no_firm_closer_eod': safe_mean(no_firm['closer_eod']),
        'pct_cust_overlay': float(sub['has_cust_overlay'].mean()),
        'cust_overlay_toward_eod': safe_mean(cust_overlay['toward_eod']),
        'no_cust_toward_eod': safe_mean(no_cust['toward_eod']),
        'n_firm_mm': int(len(firm_mm)),
        'n_cust_overlay': int(len(cust_overlay)),
    }
    
    print(f"\n  [{split}]")
    print(f"    Firm+MM overlap ({t7['n_firm_mm']}, {t7['pct_firm_mm_overlap']:.1%}): toward EOD={t7['firm_mm_toward_eod']:.1%}, closer={t7['firm_mm_closer_eod']:.1%}")
    print(f"    No firm overlap: toward EOD={t7['no_firm_toward_eod']:.1%}")
    print(f"    Cust overlay ({t7['n_cust_overlay']}, {t7['pct_cust_overlay']:.1%}): toward EOD={t7['cust_overlay_toward_eod']:.1%}")
    print(f"    No cust overlay: toward EOD={t7['no_cust_toward_eod']:.1%}")
    
    analysis[f'test7_{split}'] = t7

###############################################################################
# Additional: Distance-based analysis
###############################################################################
print("\n--- BONUS: Distance-Based Analysis ---")

for split in ['IS', 'OOS', 'ALL']:
    sub = df if split == 'ALL' else df[df['is_oos'] == split]
    
    # Near (< 25 pts) vs Far (> 50 pts)
    near = sub[sub['abs_dist'] < 25]
    mid = sub[(sub['abs_dist'] >= 25) & (sub['abs_dist'] <= 50)]
    far = sub[sub['abs_dist'] > 50]
    
    print(f"\n  [{split}]")
    print(f"    Near (<25 pts, N={len(near)}): toward EOD={safe_mean(near['toward_eod']):.1%}, touch={safe_mean(near['touched_eod']):.1%}")
    print(f"    Mid (25-50, N={len(mid)}): toward EOD={safe_mean(mid['toward_eod']):.1%}, touch={safe_mean(mid['touched_eod']):.1%}")
    print(f"    Far (>50, N={len(far)}): toward EOD={safe_mean(far['toward_eod']):.1%}, touch={safe_mean(far['touched_eod']):.1%}")
    
    analysis[f'distance_{split}'] = {
        'near_toward_eod': safe_mean(near['toward_eod']),
        'near_touch_eod': safe_mean(near['touched_eod']),
        'near_n': int(len(near)),
        'mid_toward_eod': safe_mean(mid['toward_eod']),
        'mid_touch_eod': safe_mean(mid['touched_eod']),
        'mid_n': int(len(mid)),
        'far_toward_eod': safe_mean(far['toward_eod']),
        'far_touch_eod': safe_mean(far['touched_eod']),
        'far_n': int(len(far)),
    }

###############################################################################
# 3-day manual validation
###############################################################################
print("\n--- 3-DAY MANUAL VALIDATION ---")
sample_dates = sorted(df['date'].unique())[:3]
for d in sample_dates:
    day_df = df[df['date'] == d]
    morning = day_df[day_df['time_of_day'] == 'morning'].iloc[0] if len(day_df[day_df['time_of_day'] == 'morning']) > 0 else None
    if morning is not None:
        print(f"\n  {d}: ES={morning['current_price']:.0f}, Magnet={morning['most_neg_strike']:.0f} ({morning['most_neg_gamma']:.0f}), "
              f"Dist={morning['abs_dist']:.0f}pts, Above={morning['magnet_above']}")
        print(f"    Move EOD={morning['move_eod']}, Toward={morning['toward_eod']}, Touched={morning['touched_eod']}")
        print(f"    Gravity={morning['gravity_center']:.0f}, Gravity correct EOD={morning['gravity_correct_eod']}")

###############################################################################
# Save JSON
###############################################################################
print("\nSaving results...")

# Convert any remaining numpy types
def convert(obj):
    if isinstance(obj, (np.bool_, bool)):
        return bool(obj)
    if isinstance(obj, (np.integer,)):
        return int(obj)
    if isinstance(obj, (np.floating,)):
        return float(obj) if not np.isnan(obj) else None
    if isinstance(obj, dict):
        return {k: convert(v) for k, v in obj.items()}
    if isinstance(obj, list):
        return [convert(v) for v in obj]
    return obj

analysis_clean = convert(analysis)

with open(OUTPUT_JSON, 'w') as f:
    json.dump(analysis_clean, f, indent=2)

###############################################################################
# Generate markdown report
###############################################################################
a = analysis_clean

md = """# MM Negative Gamma Magnet Effect Research

## Executive Summary

"""

# Determine overall verdict
t1_all = a['test1_ALL']
t1_oos = a['test1_OOS']

is_toward = t1_all['toward_eod_rate']
oos_toward = t1_oos['toward_eod_rate']
is_closer = t1_all['closer_eod_rate']
oos_closer = t1_oos['closer_eod_rate']

if oos_toward and oos_toward > 0.53 and t1_oos['toward_eod_tstat'] > 2.0:
    verdict = "YES — Price shows statistically significant tendency to move TOWARD MM short-gamma strikes."
elif oos_toward and oos_toward < 0.47 and t1_oos['toward_eod_tstat'] < -2.0:
    verdict = "REVERSE — Price is actually REPELLED from MM short-gamma strikes (wall effect, not magnet)."
elif oos_toward and abs(oos_toward - 0.5) < 0.03:
    verdict = "NO — No meaningful magnet or wall effect detected. Essentially random (coin flip)."
else:
    if oos_toward and oos_toward > 0.5:
        verdict = f"WEAK YES — Slight tendency toward magnet ({oos_toward:.1%}) but not statistically significant (t={t1_oos['toward_eod_tstat']:.2f})."
    else:
        verdict = f"WEAK NO — Slight tendency away from magnet ({oos_toward:.1%}) but not statistically significant (t={t1_oos['toward_eod_tstat']:.2f})."

md += f"**Does price get pulled toward MM short-gamma strikes?**\n\n> {verdict}\n\n"

md += f"""## Data Coverage
- **TRACE intradayStrikeGEX**: 180 days ({dates[0]} to {dates[-1]})
- **ES 1-min bars**: {es.index[0].date()} to {es.index[-1].date()}
- **Total snapshots analyzed**: {len(df):,}
- **IS/OOS split**: {int(n_dates*2/3)} days IS / {n_dates - int(n_dates*2/3)} days OOS (cutoff: {is_cutoff})

## Test 1: Single Biggest MM Negative Gamma Strike (Magnet)

| Metric | IS | OOS |
|--------|-----|------|
| N snapshots | {a['test1_IS']['n']:,} | {a['test1_OOS']['n']:,} |
| Move toward 1H | {a['test1_IS']['toward_1h_rate']:.1%} | {a['test1_OOS']['toward_1h_rate']:.1%} |
| Move toward 3H | {a['test1_IS']['toward_3h_rate']:.1%} | {a['test1_OOS']['toward_3h_rate']:.1%} |
| Move toward EOD | {a['test1_IS']['toward_eod_rate']:.1%} | {a['test1_OOS']['toward_eod_rate']:.1%} |
| **Closer to magnet EOD** | **{a['test1_IS']['closer_eod_rate']:.1%}** | **{a['test1_OOS']['closer_eod_rate']:.1%}** |
| Touch rate EOD (±5pts) | {a['test1_IS']['touch_rate_eod']:.1%} | {a['test1_OOS']['touch_rate_eod']:.1%} |
| t-stat (toward EOD) | {a['test1_IS']['toward_eod_tstat']:.2f} | {a['test1_OOS']['toward_eod_tstat']:.2f} |
| t-stat (closer EOD) | {a['test1_IS']['closer_eod_tstat']:.2f} | {a['test1_OOS']['closer_eod_tstat']:.2f} |
| Avg distance (pts) | {a['test1_IS']['avg_dist_pts']:.1f} | {a['test1_OOS']['avg_dist_pts']:.1f} |
| % magnet above price | {a['test1_IS']['pct_magnet_above']:.1%} | {a['test1_OOS']['pct_magnet_above']:.1%} |

**Interpretation**: "Move toward" = price moves in the direction of the magnet. "Closer" = absolute distance to magnet decreases. 50% = random.

## Test 2: Gravity Center (Top 5 Negative Gamma Cluster)

| Metric | IS | OOS |
|--------|-----|------|
| Toward gravity EOD | {a['test2_IS']['toward_gravity_eod']:.1%} | {a['test2_OOS']['toward_gravity_eod']:.1%} |
| Closer to gravity EOD | {a['test2_IS']['closer_gravity_eod']:.1%} | {a['test2_OOS']['closer_gravity_eod']:.1%} |
| Touch gravity (±5pts) | {a['test2_IS']['touch_gravity_eod']:.1%} | {a['test2_OOS']['touch_gravity_eod']:.1%} |
| t-stat | {a['test2_IS']['toward_gravity_eod_tstat']:.2f} | {a['test2_OOS']['toward_gravity_eod_tstat']:.2f} |
| **Clustered (≤25pts) toward EOD** | **{a['test2_IS']['clustered_toward_eod']:.1%}** (N={a['test2_IS']['n_clustered']}) | **{a['test2_OOS']['clustered_toward_eod']:.1%}** (N={a['test2_OOS']['n_clustered']}) |
| Dispersed (>25pts) toward EOD | {a['test2_IS']['dispersed_toward_eod']:.1%} (N={a['test2_IS']['n_dispersed']}) | {a['test2_OOS']['dispersed_toward_eod']:.1%} (N={a['test2_OOS']['n_dispersed']}) |

## Test 3: 0DTE MM Gamma Overlay

| Metric | IS | OOS |
|--------|-----|------|
| % with 0DTE overlap | {a['test3_IS']['pct_with_0dte_overlap']:.1%} | {a['test3_OOS']['pct_with_0dte_overlap']:.1%} |
| With 0DTE → total magnet toward EOD | {a['test3_IS']['overlap_total_toward_eod']:.1%} | {a['test3_OOS']['overlap_total_toward_eod']:.1%} |
| Without 0DTE → total magnet toward EOD | {a['test3_IS']['no_overlap_total_toward_eod']:.1%} | {a['test3_OOS']['no_overlap_total_toward_eod']:.1%} |
| With 0DTE → closer EOD | {a['test3_IS']['with_0dte_closer_eod']:.1%} | {a['test3_OOS']['with_0dte_closer_eod']:.1%} |

## Test 4: Nearest Negative Gamma Wall

| Metric | IS | OOS |
|--------|-----|------|
| Toward nearest wall 1H | {a['test4_IS']['toward_nearest_1h']:.1%} | {a['test4_OOS']['toward_nearest_1h']:.1%} |
| **Toward nearest wall EOD** | **{a['test4_IS']['toward_nearest_eod']:.1%}** | **{a['test4_OOS']['toward_nearest_eod']:.1%}** |
| t-stat | {a['test4_IS']['toward_nearest_eod_tstat']:.2f} | {a['test4_OOS']['toward_nearest_eod_tstat']:.2f} |
| Wall above → toward EOD | {a['test4_IS']['wall_above_toward_eod']:.1%} (N={a['test4_IS']['n_wall_above']}) | {a['test4_OOS']['wall_above_toward_eod']:.1%} (N={a['test4_OOS']['n_wall_above']}) |
| Wall below → toward EOD | {a['test4_IS']['wall_below_toward_eod']:.1%} (N={a['test4_IS']['n_wall_below']}) | {a['test4_OOS']['wall_below_toward_eod']:.1%} (N={a['test4_OOS']['n_wall_below']}) |

## Test 5: Gamma Gravity Score (Composite Directional Predictor)

| Metric | IS | OOS |
|--------|-----|------|
| Correct direction 1H | {a['test5_IS']['gravity_correct_1h']:.1%} | {a['test5_OOS']['gravity_correct_1h']:.1%} |
| Correct direction 3H | {a['test5_IS']['gravity_correct_3h']:.1%} | {a['test5_OOS']['gravity_correct_3h']:.1%} |
| **Correct direction EOD** | **{a['test5_IS']['gravity_correct_eod']:.1%}** | **{a['test5_OOS']['gravity_correct_eod']:.1%}** |
| t-stat | {a['test5_IS']['gravity_correct_eod_tstat']:.2f} | {a['test5_OOS']['gravity_correct_eod_tstat']:.2f} |
"""

if a['test5_ALL'].get('quintile_returns_eod'):
    md += f"| Quintile returns (1=down pull, 5=up pull) | {a['test5_IS'].get('quintile_returns_eod', {})} | {a['test5_OOS'].get('quintile_returns_eod', {})} |\n"
    md += f"| Monotonicity | {a['test5_IS'].get('monotonicity', 'N/A')} | {a['test5_OOS'].get('monotonicity', 'N/A')} |\n"

md += f"""
## Test 6: Time of Day Effect

| Period | Toward EOD (IS) | Toward EOD (OOS) | Gravity Correct EOD (IS) | Gravity Correct EOD (OOS) |
|--------|----------------|------------------|--------------------------|---------------------------|
| Morning (9:30-11:00) | {a['test6_IS']['morning']['toward_eod']:.1%} (t={a['test6_IS']['morning']['toward_eod_tstat']:.2f}) | {a['test6_OOS']['morning']['toward_eod']:.1%} (t={a['test6_OOS']['morning']['toward_eod_tstat']:.2f}) | {a['test6_IS']['morning']['gravity_correct_eod']:.1%} | {a['test6_OOS']['morning']['gravity_correct_eod']:.1%} |
| Midday (11:00-14:00) | {a['test6_IS']['midday']['toward_eod']:.1%} (t={a['test6_IS']['midday']['toward_eod_tstat']:.2f}) | {a['test6_OOS']['midday']['toward_eod']:.1%} (t={a['test6_OOS']['midday']['toward_eod_tstat']:.2f}) | {a['test6_IS']['midday']['gravity_correct_eod']:.1%} | {a['test6_OOS']['midday']['gravity_correct_eod']:.1%} |
| Afternoon (14:00-16:00) | {a['test6_IS']['afternoon']['toward_eod']:.1%} (t={a['test6_IS']['afternoon']['toward_eod_tstat']:.2f}) | {a['test6_OOS']['afternoon']['toward_eod']:.1%} (t={a['test6_OOS']['afternoon']['toward_eod_tstat']:.2f}) | {a['test6_IS']['afternoon']['gravity_correct_eod']:.1%} | {a['test6_OOS']['afternoon']['gravity_correct_eod']:.1%} |

## Test 7: Participant Overlay

| Metric | IS | OOS |
|--------|-----|------|
| Firm+MM both negative → toward EOD | {a['test7_IS']['firm_mm_toward_eod']:.1%} (N={a['test7_IS']['n_firm_mm']}) | {a['test7_OOS']['firm_mm_toward_eod']:.1%} (N={a['test7_OOS']['n_firm_mm']}) |
| No firm overlap → toward EOD | {a['test7_IS']['no_firm_toward_eod']:.1%} | {a['test7_OOS']['no_firm_toward_eod']:.1%} |
| Cust positive + MM negative → toward EOD | {a['test7_IS']['cust_overlay_toward_eod']:.1%} (N={a['test7_IS']['n_cust_overlay']}) | {a['test7_OOS']['cust_overlay_toward_eod']:.1%} (N={a['test7_OOS']['n_cust_overlay']}) |
| No cust overlay → toward EOD | {a['test7_IS']['no_cust_toward_eod']:.1%} | {a['test7_OOS']['no_cust_toward_eod']:.1%} |

## Distance-Based Analysis

| Distance to Magnet | Toward EOD (IS) | Touch Rate (IS) | Toward EOD (OOS) | Touch Rate (OOS) |
|---------------------|-----------------|------------------|-------------------|-------------------|
| Near (<25 pts) | {a['distance_IS']['near_toward_eod']:.1%} (N={a['distance_IS']['near_n']}) | {a['distance_IS']['near_touch_eod']:.1%} | {a['distance_OOS']['near_toward_eod']:.1%} (N={a['distance_OOS']['near_n']}) | {a['distance_OOS']['near_touch_eod']:.1%} |
| Mid (25-50 pts) | {a['distance_IS']['mid_toward_eod']:.1%} (N={a['distance_IS']['mid_n']}) | {a['distance_IS']['mid_touch_eod']:.1%} | {a['distance_OOS']['mid_toward_eod']:.1%} (N={a['distance_OOS']['mid_n']}) | {a['distance_OOS']['mid_touch_eod']:.1%} |
| Far (>50 pts) | {a['distance_IS']['far_toward_eod']:.1%} (N={a['distance_IS']['far_n']}) | {a['distance_IS']['far_touch_eod']:.1%} | {a['distance_OOS']['far_toward_eod']:.1%} (N={a['distance_OOS']['far_n']}) | {a['distance_OOS']['far_touch_eod']:.1%} |

## Magnet vs Wall: Which Is It?

"""

# Final assessment
if oos_toward and oos_toward > 0.53:
    md += "**MAGNET effect confirmed** — price is drawn toward MM short-gamma strikes. The dealer hedging amplification theory appears valid.\n\n"
elif oos_toward and oos_toward < 0.47:
    md += "**WALL effect observed** — price is repelled from MM short-gamma strikes. This contradicts the SpotGamma magnet theory; these strikes act as barriers, not attractors.\n\n"
else:
    md += "**Neither clear magnet nor wall** — the effect is too weak to be actionable. MM negative gamma positioning alone does not reliably predict ES price direction.\n\n"

md += """## Actionable Trading Rules

"""

# Generate rules based on what works
best_signals = []

# Check if any test shows OOS > 55% with t > 2
for test_name, metric_key, desc in [
    ('test1_OOS', 'toward_eod_rate', 'Trade toward biggest MM neg-gamma strike'),
    ('test1_OOS', 'closer_eod_rate', 'Trade expecting price to get closer to biggest MM neg-gamma'),
    ('test2_OOS', 'toward_gravity_eod', 'Trade toward gravity center of top-5 neg-gamma cluster'),
    ('test4_OOS', 'toward_nearest_eod', 'Trade toward nearest MM neg-gamma wall'),
    ('test5_OOS', 'gravity_correct_eod', 'Trade in direction of gamma gravity score'),
]:
    val = a.get(test_name, {}).get(metric_key)
    if val and val > 0.55:
        best_signals.append((desc, val))

if best_signals:
    md += "Based on OOS results (>55% hit rate):\n\n"
    for desc, rate in sorted(best_signals, key=lambda x: -x[1]):
        md += f"1. **{desc}** — {rate:.1%} OOS hit rate\n"
else:
    md += "**No signal exceeded 55% OOS hit rate.** The MM negative gamma magnet effect, as tested here, does not produce actionable standalone trading signals.\n\n"
    md += "Consider:\n"
    md += "- Combining with other signals (GEX regime, volatility regime, trend)\n"
    md += "- Using as a filter rather than a primary signal\n"
    md += "- The effect may be real but too small to trade profitably after costs\n"

md += """
## Methodology Notes

- **No lookahead**: all signals computed from TRACE data at time T, forward returns from ES bars after T
- **Timezone handling**: TRACE timestamps are ET-aware; ES bars converted from UTC to US/Eastern  
- **Forward horizons**: 1H, 3H, and EOD (16:00 ET) measured from signal time
- **"Toward"**: price moves in the direction of the magnet (up if magnet above, down if below)
- **"Closer"**: absolute distance to magnet decreases
- **Touch**: ES high/low between signal and EOD comes within ±5 pts of magnet strike
- **IS/OOS**: First 2/3 of dates = IS, last 1/3 = OOS (time-ordered)
- **t-stat**: vs null hypothesis of 50% (coin flip)
"""

with open(OUTPUT_MD, 'w') as f:
    f.write(md)

print(f"\nSaved: {OUTPUT_JSON}")
print(f"Saved: {OUTPUT_MD}")
print("DONE.")