diff --git a/Experiments/v7.py b/Experiments/v7.py new file mode 100644 index 0000000..f25c10f --- /dev/null +++ b/Experiments/v7.py @@ -0,0 +1,1206 @@ +# ============================================================ +# Latent Micro-Regimes in Limit Order Books: +# Identification and Early Detection — v7 +# ───────────────────────────────────────── +# KEY UPGRADES over v6: +# +# 1. THRESHOLD SWEEP — full precision/recall/coverage/Δ +# across signal_pct ∈ [70, 95] +# 2. COVERAGE METRIC — coverage = #early_τ / #σ +# 3. PRECISION–RECALL CURVE — dominance over baselines +# 4. MULTI-REGIME ROBUSTNESS — varied delays, noise, strength +# 5. DETECTION IMPROVEMENT — adaptive + multi-trigger confirmation +# 6. SIGNAL DIAGNOSTICS — which channel triggers earliest +# 7. PUBLICATION-QUALITY FIGS — 8 new figures +# +# Core v6 invariants PRESERVED: +# - causal DGP unchanged +# - evaluation logic unchanged +# - baselines unchanged +# - statistical tests unchanged +# - MAX-trigger + rising-edge architecture unchanged +# ============================================================ + +# !pip install hmmlearn scikit-learn scipy numpy pandas matplotlib +import warnings +warnings.filterwarnings("ignore") + +import numpy as np +import pandas as pd +import matplotlib +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +from scipy import stats +from scipy.stats import gaussian_kde +from sklearn.preprocessing import StandardScaler +from hmmlearn.hmm import GaussianHMM +from itertools import product as iproduct +from collections import defaultdict + +# ───────────────────────────────────────── +# 0. Global Configuration +# ───────────────────────────────────────── +SEED = 42 +T = 14_000 +N_REGIMES = 3 +MAX_LAG = 60 +FW_WINDOW = 20 +STRESS_PCT = 95 +N_BOOT = 2_000 +MIN_GAP = 20 +PENALTY = -MAX_LAG + +# v6 trigger parameters (unchanged) +GAMMA_AMP = 0.35 +DRIFT_RISE_THR = 0.30 +EDGE_DIFF_STEPS = 3 +MIN_LEAD = 5 +SIGNAL_PCT = 85 +SMOOTH_WIN = 5 + +# v7: threshold sweep range +SWEEP_PCTS = np.arange(70, 96, 1) # 70 → 95 inclusive +MULTI_K = 3 # multi-trigger confirmation steps + +np.random.seed(SEED) + +# ───────────────────────────────────────── +# Publication-quality style +# ───────────────────────────────────────── +plt.rcParams.update({ + "font.family" : "serif", + "font.serif" : ["Palatino Linotype", "Palatino", "Georgia", "Times New Roman"], + "font.size" : 11, + "axes.spines.top" : False, + "axes.spines.right": False, + "axes.linewidth" : 0.8, + "axes.titlesize" : 12, + "axes.labelsize" : 11, + "xtick.labelsize" : 9, + "ytick.labelsize" : 9, + "legend.fontsize" : 9, + "figure.dpi" : 150, + "figure.facecolor" : "white", + "axes.facecolor" : "white", + "savefig.facecolor": "white", + "savefig.dpi" : 200, +}) + +PALETTE = { + "Model" : "#1B4F8A", + "Imbalance" : "#B5341B", + "Volatility" : "#2E7D32", + "Adaptive" : "#7B1FA2", + "MultiTrig" : "#E65100", +} +REGIME_FILL = {0: "#D6EAF8", 1: "#FEF9E7", 2: "#FDEDEC"} +CHANNEL_COLORS = { + 'entropy' : "#34495E", + 'prestress' : "#E67E22", + 'drift_spread' : "#27AE60", + 'depth_erosion': "#8E44AD", + 'ofi_momentum' : "#2980B9", +} + +# ───────────────────────────────────────── +# 1. Causal Delayed Stress DGP (UNCHANGED) +# ───────────────────────────────────────── +REGIME_PARAMS = { + 0: dict(sp_mu=1.5, sp_sig=0.20, dp_ar=0.95, dp_mu=120.0, dp_sig=6.0, + ib_mu=0.00, ib_sig=0.06, vol_noise=0.02), + 1: dict(sp_mu=2.4, sp_sig=0.35, dp_ar=0.93, dp_mu=92.0, dp_sig=9.0, + ib_mu=0.12, ib_sig=0.09, vol_noise=0.06), + 2: dict(sp_mu=8.0, sp_sig=1.30, dp_ar=0.88, dp_mu=35.0, dp_sig=18.0, + ib_mu=0.50, ib_sig=0.20, vol_noise=0.40), +} +DELAY_LO, DELAY_HI, BLEND_WIN = 10, 50, 8 + + +def _draw_delay(rng): + return int(rng.integers(DELAY_LO, DELAY_HI + 1)) + +def _draw_crisis_duration(rng): + return int(rng.integers(15, 61)) + +def _draw_stable_duration(rng): + return int(rng.integers(80, 301)) + +def build_regime_sequence(T, rng): + Z = np.zeros(T, dtype=int) + delay_map = {} + t = 0 + while t < T: + dur0 = _draw_stable_duration(rng) + end0 = min(t + dur0, T) + Z[t:end0] = 0 + t = end0 + if t >= T: break + k = _draw_delay(rng) + end1 = min(t + k, T) + Z[t:end1] = 1 + delay_map[t] = k + t = end1 + if t >= T: break + dur2 = _draw_crisis_duration(rng) + end2 = min(t + dur2, T) + Z[t:end2] = 2 + t = end2 + return Z, delay_map + +def _blend(x, Z, win=BLEND_WIN): + out = x.copy() + boundaries = np.where(np.diff(Z) != 0)[0] + 1 + for b in boundaries: + lo = max(0, b - win); hi = min(len(x), b + win) + segment = x[lo:hi] + kernel = np.exp(-0.5 * ((np.arange(len(segment)) - win) / (win / 2))**2) + kernel /= kernel.sum() + out[lo:hi] = np.convolve(segment, kernel, mode='same') + return out + +def generate_lob_data(T, rng, regime_params=None, delay_lo=None, delay_hi=None): + """Parameterised for robustness tests (delay/noise/strength overrides).""" + global DELAY_LO, DELAY_HI + if regime_params is None: + regime_params = REGIME_PARAMS + lo_save, hi_save = DELAY_LO, DELAY_HI + if delay_lo is not None: DELAY_LO = delay_lo + if delay_hi is not None: DELAY_HI = delay_hi + + Z, delay_map = build_regime_sequence(T, rng) + spread = np.zeros(T); depth = np.zeros(T); imbalance = np.zeros(T) + hawkes = 0.0; hawkes_decay = 0.90 + for t in range(T): + p = regime_params[Z[t]] + hawkes *= hawkes_decay + base = np.log(p['sp_mu']) + eps = rng.normal(0, p['sp_sig']) + rng.normal(0, p['vol_noise']) + spread[t] = np.exp(base + 0.10 * hawkes + eps) + if spread[t] > np.exp(base + 0.8 * p['sp_sig']): hawkes += 0.30 + depth[0] = regime_params[Z[0]]['dp_mu'] + for t in range(1, T): + p = regime_params[Z[t]] + depth[t] = (p['dp_ar'] * depth[t-1] + (1 - p['dp_ar']) * p['dp_mu'] + + rng.normal(0, p['dp_sig'])) + depth = np.clip(depth, 5.0, None) + for t in range(T): + p = regime_params[Z[t]] + imbalance[t] = np.clip(rng.normal(p['ib_mu'], p['ib_sig']), -1.0, 1.0) + spread = _blend(spread, Z); depth = _blend(depth, Z); imbalance = _blend(imbalance, Z) + roll_vol = (pd.Series(spread).pct_change().rolling(20, min_periods=1) + .std().fillna(0).values) + ofi = imbalance * np.abs(np.diff(spread, prepend=spread[0])) + X = np.column_stack([spread, depth, imbalance, roll_vol, ofi]) + + DELAY_LO, DELAY_HI = lo_save, hi_save + return X, Z, delay_map + +# ───────────────────────────────────────── +# 2. Feature Engineering (UNCHANGED) +# ───────────────────────────────────────── +def engineer_features(X_raw): + spread = X_raw[:, 0]; depth = X_raw[:, 1] + imbalance = X_raw[:, 2]; roll_vol = X_raw[:, 3]; ofi = X_raw[:, 4] + sd_ratio = spread / (depth + 1e-6) + abs_imb = np.abs(imbalance) + cum_ofi = pd.Series(ofi).rolling(50, min_periods=1).mean().values + roll_depth = pd.Series(depth).rolling(20, min_periods=1).mean().fillna(method='bfill').values + ddepth = -pd.Series(depth).diff(5).fillna(0).values + X_full = np.column_stack([spread, depth, imbalance, roll_vol, ofi, + sd_ratio, abs_imb, cum_ofi, roll_depth, ddepth]) + scaler = StandardScaler() + X_scaled = scaler.fit_transform(X_full) + return X_scaled, scaler + +# ───────────────────────────────────────── +# 3. HMM Fitting (UNCHANGED) +# ───────────────────────────────────────── +def fit_hmm(X, n_components=N_REGIMES, n_restarts=12, rng_seed=SEED): + best_score, best_model = -np.inf, None + for k in range(n_restarts): + model = GaussianHMM(n_components=n_components, covariance_type="full", + n_iter=400, tol=1e-7, random_state=rng_seed + k, + init_params="stmc", params="stmc") + try: + model.fit(X) + sc = model.score(X) + if sc > best_score: + best_score, best_model = sc, model + except Exception: continue + if best_model is None: + raise RuntimeError("HMM fitting failed.") + return best_model + +# ───────────────────────────────────────── +# 4. Stress Events (UNCHANGED) +# ───────────────────────────────────────── +def define_stress_events(X_raw, fw=FW_WINDOW, pct=STRESS_PCT): + spread = X_raw[:, 0] + threshold = np.percentile(spread, pct) + sigma = np.array([t for t in range(len(spread) - fw) + if np.mean(spread[t+1:t+fw+1]) > threshold], dtype=int) + return sigma + +# ───────────────────────────────────────── +# 5. Signal Channels (UNCHANGED) +# ───────────────────────────────────────── +def _norm01(x): + lo, hi = x.min(), x.max() + return (x - lo) / (hi - lo + 1e-12) + +def _causal_rolling(series, window, fn='mean'): + s = pd.Series(series) + if fn == 'mean': return s.rolling(window, min_periods=1).mean().values + elif fn == 'std': return s.rolling(window, min_periods=1).std().fillna(0).values + elif fn == 'sum': return s.rolling(window, min_periods=1).sum().values + +def smooth_posterior(posterior, window=SMOOTH_WIN): + return pd.DataFrame(posterior).rolling(window, min_periods=1).mean().values + +def hmm_entropy_signal(post): + eps = 1e-12 + return -np.sum(post * np.log(post + eps), axis=1) + +def hmm_prestress_signal(post, model): + means_raw = model.means_[:, 0] + state_rank = np.argsort(means_raw) + prestress_id = state_rank[1] + return post[:, prestress_id] + +def spread_drift_signal(spread, short_win=10, long_win=40): + s = pd.Series(spread) + fast_ma = s.rolling(short_win, min_periods=1).mean() + slow_ma = s.rolling(long_win, min_periods=1).mean() + ma_cross = np.clip((fast_ma - slow_ma).values, 0, None) + d_spread = s.diff(1).fillna(0) + spread_mom = np.clip(d_spread.rolling(short_win, min_periods=1).mean().values, 0, None) + cum_drift = d_spread.clip(lower=0).rolling(long_win, min_periods=1).sum().values + return (_norm01(ma_cross) + _norm01(spread_mom) + _norm01(cum_drift)) / 3.0 + +def depth_erosion_signal(depth, win_short=10, win_long=50): + d = pd.Series(depth) + depth_trend = d.rolling(win_short, min_periods=1).mean().values + d_depth = -d.diff(5).fillna(0).values + depth_vel = np.clip(_causal_rolling(d_depth, win_short, fn='mean'), 0, None) + depth_long = d.rolling(win_long, min_periods=1).mean().values + depth_below = np.clip(depth_long - depth_trend, 0, None) + return (_norm01(depth_vel) + _norm01(depth_below)) / 2.0 + +def ofi_momentum_signal(imbalance, ofi, win=30): + abs_imb = np.abs(imbalance) + mom_imb = _causal_rolling(abs_imb, win, fn='mean') + abs_ofi = np.abs(ofi) + mom_ofi = _causal_rolling(abs_ofi, win, fn='mean') + return (_norm01(mom_imb) + _norm01(mom_ofi)) / 2.0 + +# ───────────────────────────────────────── +# 6. Trigger Score (UNCHANGED from v6) +# ───────────────────────────────────────── +def build_trigger_score(post_smooth, model, X_raw): + spread = X_raw[:, 0]; depth = X_raw[:, 1] + imbalance = X_raw[:, 2]; ofi = X_raw[:, 4] + entropy = hmm_entropy_signal(post_smooth) + prestress = hmm_prestress_signal(post_smooth, model) + c_entropy = _norm01(entropy) + c_prestress = _norm01(prestress) + c_drift_sp = _norm01(spread_drift_signal(spread)) + c_depth_det = _norm01(depth_erosion_signal(depth)) + c_ofi_mom = _norm01(ofi_momentum_signal(imbalance, ofi)) + channel_stack = np.column_stack([c_entropy, c_prestress, c_drift_sp, + c_depth_det, c_ofi_mom]) + score_raw = np.max(channel_stack, axis=1) + drift_rising = (c_drift_sp > DRIFT_RISE_THR).astype(float) + score_amp = score_raw * (1.0 + GAMMA_AMP * drift_rising) + d_score = np.zeros_like(score_amp) + k = EDGE_DIFF_STEPS + d_score[k:] = score_amp[k:] - score_amp[:-k] + comps = {'entropy': c_entropy, 'prestress': c_prestress, + 'drift_spread': c_drift_sp, 'depth_erosion': c_depth_det, + 'ofi_momentum': c_ofi_mom} + return score_amp, d_score, comps + +# ───────────────────────────────────────── +# 7. Detection Methods (v7: 4 variants) +# ───────────────────────────────────────── +def deduplicate(indices, min_gap=MIN_GAP): + if len(indices) == 0: + return np.array([], dtype=int) + out = [indices[0]] + for idx in indices[1:]: + if idx - out[-1] >= min_gap: + out.append(idx) + return np.array(out, dtype=int) + +def apply_early_detection_constraint(tau, sigma, min_lead=MIN_LEAD, max_lag=MAX_LAG): + if len(tau) == 0 or len(sigma) == 0: + return tau + kept = [] + for t in tau: + future_sigma = sigma[(sigma > t) & (sigma <= t + max_lag)] + if len(future_sigma) == 0: continue + lead = future_sigma[0] - t + if lead >= min_lead: + kept.append(t) + return np.array(kept, dtype=int) + +def detect_standard(score_amp, d_score, signal_pct, min_gap=MIN_GAP): + """Standard v6: above threshold + rising edge.""" + threshold = np.percentile(score_amp, signal_pct) + above_thr = score_amp > threshold + rising = d_score > 0 + candidates = np.where(above_thr & rising)[0] + return deduplicate(candidates, min_gap=min_gap) + +def detect_adaptive(score_amp, d_score, signal_pct, window=500, min_gap=MIN_GAP): + """ + Adaptive threshold: rolling percentile over trailing window. + More sensitive in quiet periods, stable in volatile ones. + """ + T_ = len(score_amp) + thresh_arr = np.zeros(T_) + for t in range(T_): + lo = max(0, t - window) + thresh_arr[t] = np.percentile(score_amp[lo:t+1], signal_pct) + above_thr = score_amp > thresh_arr + rising = d_score > 0 + candidates = np.where(above_thr & rising)[0] + return deduplicate(candidates, min_gap=min_gap) + +def detect_multitrigger(score_amp, d_score, signal_pct, k_confirm=MULTI_K, min_gap=MIN_GAP): + """ + Multi-trigger confirmation: signal must exceed threshold for k consecutive steps. + Fires at the FIRST step of a confirmed k-step run (reduces false positives). + """ + threshold = np.percentile(score_amp, signal_pct) + above_thr = (score_amp > threshold).astype(int) + # Convolve: position t is True if above_thr[t:t+k] all 1 + confirmed = np.zeros(len(score_amp), dtype=bool) + for t in range(len(score_amp) - k_confirm + 1): + if above_thr[t:t+k_confirm].sum() == k_confirm: + confirmed[t] = True + candidates = np.where(confirmed & (d_score > 0))[0] + return deduplicate(candidates, min_gap=min_gap) + +def full_pipeline(model, X_scaled, X_raw, sigma, + signal_pct=SIGNAL_PCT, smooth_win=SMOOTH_WIN, + min_gap=MIN_GAP, min_lead=MIN_LEAD, + method='standard'): + """Run full detection pipeline for a given method and threshold.""" + posterior = model.predict_proba(X_scaled) + post_smooth = smooth_posterior(posterior, window=smooth_win) + score_amp, d_score, comps = build_trigger_score(post_smooth, model, X_raw) + + if method == 'standard': + tau_raw = detect_standard(score_amp, d_score, signal_pct, min_gap) + elif method == 'adaptive': + tau_raw = detect_adaptive(score_amp, d_score, signal_pct, min_gap=min_gap) + elif method == 'multitrigger': + tau_raw = detect_multitrigger(score_amp, d_score, signal_pct, min_gap=min_gap) + else: + raise ValueError(f"Unknown method: {method}") + + tau = apply_early_detection_constraint(tau_raw, sigma, + min_lead=min_lead, max_lag=MAX_LAG) + return tau, score_amp, d_score, comps, post_smooth, tau_raw + +# ───────────────────────────────────────── +# 8. Baselines (UNCHANGED) +# ───────────────────────────────────────── +def imbalance_baseline(X_raw, pct=90, min_gap=MIN_GAP): + imb = np.abs(X_raw[:, 2]) + raw = np.where(imb > np.percentile(imb, pct))[0] + return deduplicate(raw, min_gap=min_gap) + +def volatility_baseline(X_raw, pct=90, min_gap=MIN_GAP): + rv = X_raw[:, 3] + raw = np.where(rv > np.percentile(rv, pct))[0] + return deduplicate(raw, min_gap=min_gap) + +# ───────────────────────────────────────── +# 9. Evaluation (UNCHANGED + coverage) +# ───────────────────────────────────────── +def compute_lead_times(tau, sigma, max_lag=MAX_LAG): + deltas = np.empty(len(tau), dtype=float) + for i, t in enumerate(tau): + cands = sigma[(sigma > t) & (sigma <= t + max_lag)] + deltas[i] = (cands[0] - t) if len(cands) > 0 else PENALTY + return deltas + +def compute_coverage(tau, sigma, max_lag=MAX_LAG): + """Fraction of stress events covered by at least one early signal.""" + if len(tau) == 0 or len(sigma) == 0: + return 0.0 + covered = 0 + for s in sigma: + # Is there any τ in (s - max_lag, s)? + preceding = tau[(tau < s) & (tau >= s - max_lag)] + if len(preceding) > 0: + covered += 1 + return covered / len(sigma) + +def evaluation_metrics(deltas, sigma, tau, max_lag=MAX_LAG): + valid = deltas > 0 + cov = compute_coverage(tau, sigma, max_lag) + return dict( + mean_delta = float(np.mean(deltas)), + pct_early = float(np.mean(valid)), # Precision + recall = float(cov), # Coverage = Recall + mean_early = float(np.mean(deltas[valid])) if valid.any() else 0.0, + std_delta = float(np.std(deltas)), + n_tau = int(len(deltas)), + n_early = int(valid.sum()), + coverage = float(cov), + ) + +def bootstrap_ci(deltas, stat_fn=np.mean, n_boot=N_BOOT, alpha=0.05, seed=SEED): + rng = np.random.default_rng(seed) + boot = np.array([ + stat_fn(rng.choice(deltas, size=len(deltas), replace=True)) + for _ in range(n_boot) + ]) + return (float(np.percentile(boot, 100*alpha/2)), + float(np.percentile(boot, 100*(1-alpha/2)))) + +def mannwhitney_test(a, b): + return stats.mannwhitneyu(a, b, alternative="two-sided") + +# ───────────────────────────────────────── +# 10. Threshold Sweep (NEW v7) +# ───────────────────────────────────────── +def threshold_sweep(model, X_scaled, X_raw, sigma, + pct_range=SWEEP_PCTS, + smooth_win=SMOOTH_WIN, + min_gap=MIN_GAP, + min_lead=MIN_LEAD, + method='standard'): + """ + Sweep signal_pct and record Precision, Recall (Coverage), Mean Δ, N(τ). + Returns a DataFrame indexed by pct. + """ + rows = [] + posterior = model.predict_proba(X_scaled) + post_smooth = smooth_posterior(posterior, window=smooth_win) + score_amp, d_score, comps = build_trigger_score(post_smooth, model, X_raw) + + for pct in pct_range: + if method == 'standard': + tau_raw = detect_standard(score_amp, d_score, pct, min_gap) + elif method == 'adaptive': + tau_raw = detect_adaptive(score_amp, d_score, pct, min_gap=min_gap) + elif method == 'multitrigger': + tau_raw = detect_multitrigger(score_amp, d_score, pct, min_gap=min_gap) + + tau = apply_early_detection_constraint(tau_raw, sigma, + min_lead=min_lead, max_lag=MAX_LAG) + if len(tau) == 0: + rows.append({'pct': pct, 'mean_delta': np.nan, + 'precision': 0.0, 'recall': 0.0, + 'n_tau': 0, 'n_early': 0}) + continue + + deltas = compute_lead_times(tau, sigma) + m = evaluation_metrics(deltas, sigma, tau) + rows.append({ + 'pct' : pct, + 'mean_delta': m['mean_delta'], + 'precision' : m['pct_early'], + 'recall' : m['coverage'], + 'n_tau' : m['n_tau'], + 'n_early' : m['n_early'], + }) + return pd.DataFrame(rows) + +def baseline_pr_point(tau_bl, sigma): + """Single precision/recall point for a baseline.""" + if len(tau_bl) == 0: + return 0.0, 0.0 + deltas = compute_lead_times(tau_bl, sigma) + prec = float(np.mean(deltas > 0)) + rec = compute_coverage(tau_bl, sigma) + return prec, rec + +# ───────────────────────────────────────── +# 11. Multi-Regime Robustness Tests (NEW) +# ───────────────────────────────────────── +def build_robust_regime_params(noise_scale=1.0, strength_scale=1.0): + """ + Perturb base regime parameters for robustness testing. + noise_scale: multiplies vol_noise and signal sigmas + strength_scale: multiplies Regime-2 spread mean offset above Regime-0 + """ + base_sp = REGIME_PARAMS[2]['sp_mu'] + r2_sp = REGIME_PARAMS[0]['sp_mu'] + (base_sp - REGIME_PARAMS[0]['sp_mu']) * strength_scale + params = { + 0: dict(**REGIME_PARAMS[0]), + 1: dict(**REGIME_PARAMS[1]), + 2: dict(**REGIME_PARAMS[2]), + } + for k in [0, 1, 2]: + params[k]['sp_sig'] *= noise_scale + params[k]['vol_noise'] *= noise_scale + params[2]['sp_mu'] = max(r2_sp, REGIME_PARAMS[1]['sp_mu'] + 0.5) + return params + +def run_robustness_grid(model_base, scaler_base, + sigma_base, X_scaled_base, X_raw_base, + n_reps=3): + """ + Tests across: + - delay regimes: short (10-20), default (10-50), long (30-60) + - noise levels: low (0.6×), medium (1.0×), high (1.5×) + - strength levels: weak (0.7×), medium (1.0×), strong (1.3×) + For each, fits a new HMM and runs detection. + Returns a tidy DataFrame of results. + """ + delay_configs = [('Short', 10, 20), ('Default', 10, 50), ('Long', 30, 60)] + noise_configs = [('Low', 0.6), ('Medium', 1.0), ('High', 1.5)] + strength_configs = [('Weak', 0.7), ('Medium', 1.0), ('Strong', 1.3)] + + rows = [] + rng = np.random.default_rng(SEED + 99) + + for (dname, dlo, dhi), (nname, nsc), (sname, ssc) in iproduct( + delay_configs, noise_configs, strength_configs): + regime_params = build_robust_regime_params(noise_scale=nsc, strength_scale=ssc) + for rep in range(n_reps): + try: + X_r, Z_r, dm_r = generate_lob_data( + T, rng, regime_params=regime_params, + delay_lo=dlo, delay_hi=dhi) + X_sc_r, _ = engineer_features(X_r) + mdl_r = fit_hmm(X_sc_r, n_restarts=6) + sig_r = define_stress_events(X_r) + tau_r, sc_r, ds_r, _, _, _ = full_pipeline( + mdl_r, X_sc_r, X_r, sig_r, + signal_pct=SIGNAL_PCT, method='standard') + if len(tau_r) == 0: + rows.append({'delay': dname, 'noise': nname, 'strength': sname, + 'rep': rep, 'mean_delta': np.nan, + 'precision': 0.0, 'recall': 0.0, + 'n_tau': 0, 'n_sigma': len(sig_r)}) + continue + deltas_r = compute_lead_times(tau_r, sig_r) + m = evaluation_metrics(deltas_r, sig_r, tau_r) + rows.append({'delay': dname, 'noise': nname, 'strength': sname, + 'rep': rep, + 'mean_delta': m['mean_delta'], + 'precision' : m['pct_early'], + 'recall' : m['coverage'], + 'n_tau' : m['n_tau'], + 'n_sigma' : len(sig_r)}) + except Exception as e: + print(f" [robustness] skipped ({dname},{nname},{sname},rep{rep}): {e}") + return pd.DataFrame(rows) + +# ───────────────────────────────────────── +# 12. Signal Diagnostics (NEW) +# ───────────────────────────────────────── +def channel_lead_time_analysis(comps, sigma, Z_true, max_lag=MAX_LAG, + pct=85, min_gap=MIN_GAP, min_lead=MIN_LEAD): + """ + For each channel, treat its values as a standalone detector: + threshold at pct, find candidate rises, compute lead times. + Returns dict: channel → mean lead time at stress events. + """ + results = {} + for key, ch in comps.items(): + thr = np.percentile(ch, pct) + above = ch > thr + # rising edge + d_ch = np.zeros_like(ch) + d_ch[EDGE_DIFF_STEPS:] = ch[EDGE_DIFF_STEPS:] - ch[:-EDGE_DIFF_STEPS] + cands = np.where(above & (d_ch > 0))[0] + tau = deduplicate(cands, min_gap) + tau = apply_early_detection_constraint(tau, sigma, + min_lead=min_lead, max_lag=max_lag) + if len(tau) == 0: + results[key] = dict(mean_delta=np.nan, pct_early=0.0, + coverage=0.0, n_tau=0) + continue + deltas = compute_lead_times(tau, sigma, max_lag) + m = evaluation_metrics(deltas, sigma, tau, max_lag) + results[key] = dict( + mean_delta = m['mean_delta'], + pct_early = m['pct_early'], + coverage = m['coverage'], + n_tau = m['n_tau'], + ) + return results + +def channel_earliest_trigger_analysis(comps, sigma, max_lag=MAX_LAG): + """ + For each stress event σ, find which channel provides the earliest + (longest lead time) signal. Returns frequency counts per channel. + """ + channel_keys = list(comps.keys()) + channel_pcts = {k: np.percentile(comps[k], SIGNAL_PCT) for k in channel_keys} + winner_counts = defaultdict(int) + lead_by_channel = defaultdict(list) + + for s in sigma: + lo = max(0, s - max_lag) + best_lead = -1; best_ch = None + for key in channel_keys: + ch = comps[key] + thr = channel_pcts[key] + # Find earliest crossing in (s-max_lag, s) + hits = np.where((ch[lo:s] > thr))[0] + lo + if len(hits) > 0: + earliest = hits[0] + lead = s - earliest + if lead > best_lead: + best_lead = lead + best_ch = key + if best_ch is not None: + winner_counts[best_ch] += 1 + lead_by_channel[best_ch].append(best_lead) + + return winner_counts, lead_by_channel + +# ───────────────────────────────────────── +# 13. Sanity Check (UNCHANGED) +# ───────────────────────────────────────── +def check_baseline_blindness(X_raw, Z_true): + imb = np.abs(X_raw[:, 2]) + roll_vol = X_raw[:, 3] + imb_thr = np.percentile(imb, 90) + vol_thr = np.percentile(roll_vol, 90) + print(" ── Baseline Blindness Sanity Check ──────────────────────") + for k in range(3): + mask = Z_true == k + print(f" Regime {k} | mean |imb| = {imb[mask].mean():.4f} " + f"(frac > thr: {(imb[mask] > imb_thr).mean():.2%}) | " + f"mean rv = {roll_vol[mask].mean():.4f} " + f"(frac > thr: {(roll_vol[mask] > vol_thr).mean():.2%})") + print() + +# ───────────────────────────────────────── +# 14. Publication-Quality Visualizations +# ───────────────────────────────────────── + +def _shade_regimes(ax, Z, t_end, y_lo, y_hi): + for k, c in REGIME_FILL.items(): + ax.fill_between(np.arange(t_end), y_lo, y_hi, + where=Z[:t_end] == k, color=c, alpha=0.55) + +def plot_dgp_causal_structure(X_raw, Z_true, delay_map, n_show=2500): + t_end = min(n_show, len(Z_true)); t_ax = np.arange(t_end) + spread = X_raw[:t_end, 0]; depth = X_raw[:t_end, 1]; imb = X_raw[:t_end, 2] + fig, axes = plt.subplots(3, 1, figsize=(13, 7), sharex=True) + fig.suptitle("Causal DGP: Hidden Build-Up (Regime 1) → Delayed Stress (Regime 2)", + fontsize=13, fontweight="bold", y=1.01) + for ax, y, ylabel in zip(axes, [spread, depth, imb], + ["Bid–Ask Spread", "Market Depth", "Order Imbalance"]): + _shade_regimes(ax, Z_true, t_end, y.min(), y.max()) + ax.plot(t_ax, y, lw=0.65, color="#1A1A2E"); ax.set_ylabel(ylabel) + shown = 0 + for t_entry, k in sorted(delay_map.items()): + if t_entry >= t_end: break + t_stress = min(t_entry + k, t_end - 1) + y_ann = spread[t_entry] * 1.08 + axes[0].annotate("", xy=(t_stress, y_ann * 1.06), xytext=(t_entry, y_ann), + arrowprops=dict(arrowstyle="->", color="#CC6600", lw=1.3)) + axes[0].text(t_entry, y_ann * 1.02, f"k={k}", fontsize=7, + color="#CC6600", ha="left") + shown += 1 + if shown >= 5: break + from matplotlib.patches import Patch + legend_elems = [Patch(fc=REGIME_FILL[k], label=f"Regime {k}") for k in range(3)] + axes[0].legend(handles=legend_elems, loc="upper right", + fontsize=8, frameon=False, ncol=3) + axes[2].set_xlabel("Timestep") + fig.tight_layout(); plt.savefig("fig1_dgp_structure.pdf", bbox_inches="tight"); plt.show() + +def plot_threshold_sweep(sweep_std, sweep_adp, sweep_mtr): + """ + Three panels: (a) Mean Δ vs threshold, (b) Coverage vs threshold, + (c) Precision vs threshold. Three method lines per panel. + """ + fig, axes = plt.subplots(1, 3, figsize=(14, 4.5)) + fig.suptitle("Threshold Sweep — Mean Lead Time, Coverage, and Precision\n" + "across Detection Methods [v7]", + fontsize=12, fontweight="bold") + methods = [ + (sweep_std, 'Standard', PALETTE['Model'], '-'), + (sweep_adp, 'Adaptive', PALETTE['Adaptive'], '--'), + (sweep_mtr, 'Multi-Trig.', PALETTE['MultiTrig'], ':'), + ] + ylabels = ["Mean Lead Time Δ (steps)", "Coverage (Recall)", "Precision (% Early)"] + cols = ['mean_delta', 'recall', 'precision'] + for ax, col, ylabel in zip(axes, cols, ylabels): + for df, label, color, ls in methods: + valid = df.dropna(subset=[col]) + ax.plot(valid['pct'], valid[col], lw=2.0, color=color, + ls=ls, label=label, marker='o', markersize=3.5) + ax.axhline(0, color='gray', lw=0.8, ls='--') + ax.set_xlabel("Signal Threshold Percentile") + ax.set_ylabel(ylabel) + ax.legend(frameon=False) + ax.set_xlim(SWEEP_PCTS[0] - 0.5, SWEEP_PCTS[-1] + 0.5) + fig.tight_layout() + plt.savefig("fig2_threshold_sweep.pdf", bbox_inches="tight") + plt.show() + +def plot_precision_recall(sweep_std, sweep_adp, sweep_mtr, tau_imb, tau_vol, sigma): + """ + Precision–Recall curve. Model methods form a curve; baselines are single points. + """ + fig, ax = plt.subplots(figsize=(7.5, 5.5)) + methods = [ + (sweep_std, 'Standard', PALETTE['Model'], '-', 'o'), + (sweep_adp, 'Adaptive', PALETTE['Adaptive'], '--', 's'), + (sweep_mtr, 'Multi-Trig.', PALETTE['MultiTrig'], ':', '^'), + ] + for df, label, color, ls, mk in methods: + valid = df.dropna(subset=['precision', 'recall']) + ax.plot(valid['recall'], valid['precision'], lw=2.2, color=color, + ls=ls, label=label, marker=mk, markersize=4.5, zorder=3) + # annotate a few percentile points + for _, row in valid.iloc[::6].iterrows(): + ax.annotate(f"{int(row['pct'])}%", + (row['recall'], row['precision']), + fontsize=6.5, color=color, + xytext=(4, 2), textcoords='offset points') + + # Baselines as single points + for tau_bl, name, color in [(tau_imb, 'Imbalance', PALETTE['Imbalance']), + (tau_vol, 'Volatility', PALETTE['Volatility'])]: + prec, rec = baseline_pr_point(tau_bl, sigma) + ax.scatter(rec, prec, s=90, color=color, zorder=5, + label=f"{name} baseline", marker='D', edgecolors='white', linewidth=0.8) + ax.annotate(name, (rec, prec), fontsize=9, color=color, + xytext=(6, 3), textcoords='offset points', fontweight='bold') + + ax.set_xlabel("Recall (Coverage = # stress events covered / # total)") + ax.set_ylabel("Precision (% Early detections)") + ax.set_title("Precision–Recall Trade-Off:\nModel Methods vs Baselines [v7]", + fontsize=12, fontweight='bold') + ax.set_xlim(-0.02, 1.05); ax.set_ylim(-0.02, 1.05) + ax.axline((0, 0), slope=1, color='lightgray', lw=0.8, ls='--', zorder=0) + ax.legend(frameon=False, fontsize=9, loc='lower left') + fig.tight_layout() + plt.savefig("fig3_precision_recall.pdf", bbox_inches="tight") + plt.show() + +def plot_robustness_heatmap(rob_df): + """ + Heatmap: mean Δ and Recall across delay × noise configurations. + Averages over strength and rep. + """ + if rob_df.empty: + print(" [Warning] Robustness grid is empty — skipping heatmap.") + return + + fig, axes = plt.subplots(1, 2, figsize=(13, 4)) + fig.suptitle("Robustness: Mean Lead Time Δ and Coverage across\n" + "Delay Regime × Noise Level [v7]", + fontsize=12, fontweight="bold") + + for ax, metric, title in zip(axes, + ['mean_delta', 'recall'], + ['Mean Lead Time Δ', 'Coverage (Recall)']): + pivot = (rob_df.groupby(['delay', 'noise'])[metric] + .mean() + .unstack('noise')) + # reorder + order_delay = ['Short', 'Default', 'Long'] + order_noise = ['Low', 'Medium', 'High'] + pivot = pivot.reindex(index=[d for d in order_delay if d in pivot.index], + columns=[n for n in order_noise if n in pivot.columns]) + im = ax.imshow(pivot.values, cmap='RdYlGn', aspect='auto', + vmin=pivot.values[~np.isnan(pivot.values)].min() if not np.all(np.isnan(pivot.values)) else 0, + vmax=pivot.values[~np.isnan(pivot.values)].max() if not np.all(np.isnan(pivot.values)) else 1) + ax.set_xticks(range(len(pivot.columns))); ax.set_xticklabels(pivot.columns) + ax.set_yticks(range(len(pivot.index))); ax.set_yticklabels(pivot.index) + ax.set_xlabel("Noise Level"); ax.set_ylabel("Delay Regime") + ax.set_title(title) + plt.colorbar(im, ax=ax, shrink=0.85) + for i in range(len(pivot.index)): + for j in range(len(pivot.columns)): + val = pivot.values[i, j] + if not np.isnan(val): + ax.text(j, i, f"{val:.2f}", ha='center', va='center', + fontsize=9, fontweight='bold', + color='white' if abs(val) > 0.5 * pivot.values.max() else 'black') + fig.tight_layout() + plt.savefig("fig4_robustness_heatmap.pdf", bbox_inches="tight") + plt.show() + +def plot_channel_diagnostics(comps, diag_results, winner_counts, lead_by_channel): + """ + Two panels: + (a) Channel standalone performance (mean Δ, precision, recall bars) + (b) Frequency of "earliest trigger" per channel + """ + ch_labels = { + 'entropy' : "HMM\nEntropy", + 'prestress' : "HMM\nPre-Stress", + 'drift_spread' : "Spread\nDrift", + 'depth_erosion': "Depth\nErosion", + 'ofi_momentum' : "OFI\nMomentum", + } + keys = list(ch_labels.keys()) + colors = [CHANNEL_COLORS[k] for k in keys] + + fig, axes = plt.subplots(1, 2, figsize=(13, 5)) + fig.suptitle("Signal Channel Diagnostics [v7]\n" + "Per-channel detection performance and earliest-trigger frequency", + fontsize=12, fontweight="bold") + + # Panel A: grouped bar + metrics = ['mean_delta', 'pct_early', 'coverage'] + m_labels = ['Mean Δ (scaled)', 'Precision', 'Coverage'] + x = np.arange(len(keys)) + width = 0.22 + ax = axes[0] + for mi, (met, mlab) in enumerate(zip(metrics, m_labels)): + vals = [] + for k in keys: + v = diag_results.get(k, {}).get(met, np.nan) + if met == 'mean_delta' and not np.isnan(v): + v = v / MAX_LAG # normalize to [-1,1] + vals.append(v if not np.isnan(v) else 0.0) + ax.bar(x + mi * width, vals, width, label=mlab, + alpha=0.82, edgecolor='white') + ax.axhline(0, color='gray', lw=0.7, ls='--') + ax.set_xticks(x + width); ax.set_xticklabels([ch_labels[k] for k in keys], fontsize=9) + ax.set_ylabel("Metric Value (Mean Δ normalized to [−1, 1])") + ax.set_title("Standalone Channel Performance") + ax.legend(frameon=False, fontsize=9) + + # Panel B: earliest-trigger pie + ax = axes[1] + total = sum(winner_counts.values()) + if total > 0: + sizes = [winner_counts.get(k, 0) for k in keys] + labels = [f"{ch_labels[k]}\n({winner_counts.get(k,0)})" for k in keys] + wedges, texts, autotexts = ax.pie( + sizes, labels=labels, colors=colors, + autopct=lambda p: f'{p:.1f}%' if p > 4 else '', + startangle=140, pctdistance=0.78, + wedgeprops=dict(edgecolor='white', linewidth=1.5)) + for at in autotexts: + at.set_fontsize(8) + ax.set_title("Earliest-Trigger Frequency per Channel\n" + "(Which channel fires first before each stress event?)") + else: + ax.text(0.5, 0.5, "No winners found", transform=ax.transAxes, + ha='center', va='center') + + fig.tight_layout() + plt.savefig("fig5_channel_diagnostics.pdf", bbox_inches="tight") + plt.show() + +def plot_composite_signal(X_raw, Z_true, score_amp, tau_model, sigma, + method_label="Standard", n_show=3000): + t_end = min(n_show, len(Z_true)); t_ax = np.arange(t_end) + tau_vis = tau_model[tau_model < t_end]; sigma_vis = sigma[sigma < t_end] + sc = score_amp[:t_end] + thresh = np.percentile(score_amp, SIGNAL_PCT) + spread = X_raw[:t_end, 0]; depth = X_raw[:t_end, 1]; imb = X_raw[:t_end, 2] + + fig, axes = plt.subplots(4, 1, figsize=(13, 10), sharex=True, + gridspec_kw={"height_ratios": [2, 2.5, 1.5, 1.5]}) + fig.suptitle(f"Trigger-Based Instability Detector — v7 [{method_label}]\n" + "MAX-trigger + Rising-Edge + Early-Detection Constraint", + fontsize=13, fontweight="bold") + + ax = axes[0] + _shade_regimes(ax, Z_true, t_end, 0, spread.max()*1.1) + ax.plot(t_ax, spread, lw=0.65, color="#1A1A2E") + ax.set_ylabel("Spread") + from matplotlib.patches import Patch + ax.legend(handles=[Patch(fc=REGIME_FILL[k], label=f"Regime {k}") for k in range(3)], + loc="upper right", fontsize=8, frameon=False, ncol=3) + + ax = axes[1] + ax.plot(t_ax, sc, lw=0.8, color="#333333", alpha=0.85, label="MAX score") + ax.axhline(thresh, color="#E67E22", lw=1.2, ls="--", + label=f"{SIGNAL_PCT}th pct threshold") + ax.fill_between(t_ax, thresh, sc, where=sc > thresh, + color=PALETTE["Model"], alpha=0.18) + ax.vlines(tau_vis, sc.min(), sc.max(), + color=PALETTE["Model"], lw=1.0, alpha=0.75, label="Signal τ") + ax.vlines(sigma_vis, sc.min(), sc.max(), + color=PALETTE["Imbalance"], lw=0.6, ls=":", alpha=0.45, + label="Stress event σ") + ax.set_ylabel("Instability Score") + ax.legend(loc="upper right", fontsize=8, frameon=False, ncol=2) + + axes[2].plot(t_ax, depth, lw=0.7, color="#2D6A4F"); axes[2].set_ylabel("Depth") + axes[3].plot(t_ax, imb, lw=0.7, color="#6A3D9A", alpha=0.85) + axes[3].set_ylabel("Imbalance"); axes[3].set_xlabel("Timestep") + fig.tight_layout() + plt.savefig(f"fig6_composite_{method_label.lower()}.pdf", bbox_inches="tight") + plt.show() + +def plot_lead_time_densities(delta_dict): + fig, ax = plt.subplots(figsize=(9, 5)) + x_grid = np.linspace(-MAX_LAG - 5, MAX_LAG + 5, 800) + ordered = [('Model', PALETTE['Model']), + ('Adaptive', PALETTE['Adaptive']), + ('Multi-Trig.',PALETTE['MultiTrig']), + ('Imbalance', PALETTE['Imbalance']), + ('Volatility', PALETTE['Volatility'])] + for i, (name, color) in enumerate(ordered): + if name not in delta_dict: continue + deltas = delta_dict[name] + valid = deltas[deltas > PENALTY] + if len(valid) > 5: + kde = gaussian_kde(valid, bw_method='scott') + ax.plot(x_grid, kde(x_grid), lw=2.2, color=color, label=name) + ax.fill_between(x_grid, kde(x_grid), alpha=0.12, color=color) + missed = np.mean(deltas <= PENALTY) + mean_v = np.mean(deltas[deltas > 0]) if (deltas > 0).any() else 0 + ax.annotate(f"{name} missed={missed:.1%} E[Δ|early]={mean_v:+.1f}", + xy=(-MAX_LAG + 1, 0.006 * (i + 1)), + color=color, fontsize=8.5, fontweight="bold") + ax.axvline(0, color='gray', lw=1.2, ls='--', label='Zero lead-time') + ax.set_xlabel("Lead Time Δ (timesteps before stress event)") + ax.set_ylabel("Density") + ax.set_title("Lead-Time Distributions: All Detectors [v7]", + fontsize=13, pad=10) + ax.legend(frameon=False, fontsize=9); ax.set_xlim(-MAX_LAG - 2, MAX_LAG + 2) + fig.tight_layout() + plt.savefig("fig7_lead_time_densities.pdf", bbox_inches="tight") + plt.show() + +def plot_results_table(results_df): + fig, ax = plt.subplots(figsize=(16, 3.0)) + ax.axis("off") + tbl = ax.table(cellText=results_df.values, colLabels=results_df.columns, + cellLoc="center", loc="center") + tbl.auto_set_font_size(False); tbl.set_fontsize(9.2); tbl.scale(1.15, 1.8) + header_color = "#1B4F8A" + row_colors = ["#EDF4FF", "#FFFFFF"] + for j in range(len(results_df.columns)): + tbl[0, j].set_facecolor(header_color) + tbl[0, j].set_text_props(color="white", fontweight="bold") + for i in range(1, len(results_df) + 1): + for j in range(len(results_df.columns)): + tbl[i, j].set_facecolor(row_colors[(i - 1) % 2]) + fig.suptitle("Detection Performance Summary — v7\n" + "(Threshold-Based Early Detection with Coverage and Precision–Recall)", + fontsize=10, y=1.04, fontweight="bold") + fig.tight_layout() + plt.savefig("fig8_results_table.pdf", bbox_inches="tight") + plt.show() + +def plot_sweep_coverage_delta(sweep_std, sweep_adp, sweep_mtr): + """ + Dual-axis plot: Mean Δ (left) and Coverage (right) vs threshold pct. + Shows trade-off visually in one panel. + """ + fig, ax1 = plt.subplots(figsize=(9, 5)) + ax2 = ax1.twinx() + methods = [ + (sweep_std, 'Standard', PALETTE['Model'], '-'), + (sweep_adp, 'Adaptive', PALETTE['Adaptive'], '--'), + (sweep_mtr, 'Multi-Trig.', PALETTE['MultiTrig'], ':'), + ] + for df, label, color, ls in methods: + valid = df.dropna(subset=['mean_delta', 'recall']) + ax1.plot(valid['pct'], valid['mean_delta'], lw=2.0, color=color, + ls=ls, label=f"{label} – Mean Δ", marker='o', markersize=3.5) + ax2.plot(valid['pct'], valid['recall'], lw=1.5, color=color, + ls='--', alpha=0.55, label=f"{label} – Coverage") + ax1.axhline(0, color='gray', lw=0.8, ls='--') + ax1.set_xlabel("Signal Threshold Percentile") + ax1.set_ylabel("Mean Lead Time Δ (steps)") + ax2.set_ylabel("Coverage (Recall)", color='#444444') + ax2.tick_params(axis='y', colors='#444444') + ax1.set_title("Early Detection vs Coverage Trade-Off [v7]", + fontsize=12, fontweight='bold') + lines1, labels1 = ax1.get_legend_handles_labels() + lines2, labels2 = ax2.get_legend_handles_labels() + ax1.legend(lines1 + lines2, labels1 + labels2, frameon=False, fontsize=8.5, + loc='upper left', ncol=2) + fig.tight_layout() + plt.savefig("fig9_sweep_coverage_delta.pdf", bbox_inches="tight") + plt.show() + +# ───────────────────────────────────────── +# 15. Main Pipeline +# ───────────────────────────────────────── +def run_experiment(): + rng = np.random.default_rng(SEED) + + print("=" * 70) + print(" LOB Micro-Regime Detection — v7") + print(" Threshold Sweep · Coverage · Precision–Recall · Robustness") + print("=" * 70) + + # ── 1. Data ───────────────────────────────────────────────── + print("\n Step 1 / 7 — Generating causal LOB data ...") + X_raw, Z_true, delay_map = generate_lob_data(T, rng) + rd = " | ".join([f"Regime {k}: {(Z_true==k).mean():.1%}" for k in range(N_REGIMES)]) + print(f" {T:,} timesteps | {rd}") + print(f" Regime-1 episodes: {len(delay_map)} | " + f"Mean delay: {np.mean(list(delay_map.values())):.1f} steps") + + # ── 2. Features ───────────────────────────────────────────── + print("\n Step 2 / 7 — Feature engineering ...") + X_scaled, scaler = engineer_features(X_raw) + print(f" Feature matrix: {X_scaled.shape}") + + # ── 3. HMM ────────────────────────────────────────────────── + print("\n Step 3 / 7 — Fitting HMM (12 restarts) ...") + model = fit_hmm(X_scaled) + ll = model.score(X_scaled); conv = model.monitor_.converged + print(f" Best log-likelihood: {ll:,.2f} | Converged: {conv}") + + # ── 4. Stress events ───────────────────────────────────────── + print("\n Step 4 / 7 — Stress event definition ...") + sigma = define_stress_events(X_raw) + print(f" Stress events: {len(sigma):,} ({len(sigma)/T:.1%})") + check_baseline_blindness(X_raw, Z_true) + + # ── 5. Detection (three methods) ───────────────────────────── + print(" Step 5 / 7 — Running detection methods ...") + tau_std, score_amp, d_score, comps, post_smooth, tau_raw_std = full_pipeline( + model, X_scaled, X_raw, sigma, signal_pct=SIGNAL_PCT, method='standard') + tau_adp, _, _, _, _, _ = full_pipeline( + model, X_scaled, X_raw, sigma, signal_pct=SIGNAL_PCT, method='adaptive') + tau_mtr, _, _, _, _, _ = full_pipeline( + model, X_scaled, X_raw, sigma, signal_pct=SIGNAL_PCT, method='multitrigger') + tau_imb = imbalance_baseline(X_raw) + tau_vol = volatility_baseline(X_raw) + + print(f" Detections — Standard: {len(tau_std)} | Adaptive: {len(tau_adp)} | " + f"Multi-Trig: {len(tau_mtr)} | Imbalance: {len(tau_imb)} | " + f"Volatility: {len(tau_vol)}") + + # ── 5b. Evaluation ───────────────────────────────────────── + print("\n Step 5b — Evaluation ...") + det_methods = { + 'Model' : tau_std, + 'Adaptive' : tau_adp, + 'Multi-Trig.': tau_mtr, + 'Imbalance' : tau_imb, + 'Volatility' : tau_vol, + } + delta_dict = {nm: compute_lead_times(tau, sigma) + for nm, tau in det_methods.items()} + + rows = [] + for name, (tau, deltas) in zip(det_methods.keys(), + zip(det_methods.values(), delta_dict.values())): + m = evaluation_metrics(deltas, sigma, tau) + lo, hi = bootstrap_ci(deltas) + rows.append({ + "Detector" : name, + "Mean Δ" : f"{m['mean_delta']:+.2f}", + "95% CI" : f"[{lo:+.2f}, {hi:+.2f}]", + "Precision" : f"{m['pct_early']:.1%}", + "Coverage" : f"{m['coverage']:.1%}", + "Mean Δ|early" : f"{m['mean_early']:+.2f}", + "N(τ)" : m['n_tau'], + "N(early)" : m['n_early'], + }) + results_df = pd.DataFrame(rows) + print("\n" + results_df.to_string(index=False)) + + print("\n Pairwise Mann–Whitney U tests (two-sided):") + pairs = [("Model","Imbalance"), ("Model","Volatility"), + ("Adaptive","Imbalance"), ("Multi-Trig.","Imbalance"), + ("Imbalance","Volatility")] + for a, b in pairs: + u, p = mannwhitney_test(delta_dict[a], delta_dict[b]) + sig = ("***" if p < 0.001 else "**" if p < 0.01 else + "*" if p < 0.05 else "ns") + print(f" {a:14s} vs {b:14s}: U={u:,.0f} p={p:.4f} {sig}") + + # ── 6. Threshold sweep ─────────────────────────────────────── + print("\n Step 6 / 7 — Threshold sweep across methods ...") + sweep_std = threshold_sweep(model, X_scaled, X_raw, sigma, + pct_range=SWEEP_PCTS, method='standard') + sweep_adp = threshold_sweep(model, X_scaled, X_raw, sigma, + pct_range=SWEEP_PCTS, method='adaptive') + sweep_mtr = threshold_sweep(model, X_scaled, X_raw, sigma, + pct_range=SWEEP_PCTS, method='multitrigger') + print(f" Sweep complete. Standard: {sweep_std['pct_early'].notna().sum()} " + f"valid thresholds / {len(SWEEP_PCTS)}") + + # ── 7. Signal diagnostics ──────────────────────────────────── + print("\n Step 6b — Signal diagnostics ...") + diag_results = channel_lead_time_analysis(comps, sigma, Z_true) + winner_counts, lead_by_ch = channel_earliest_trigger_analysis(comps, sigma) + print(" Channel standalone performance:") + for key, res in diag_results.items(): + print(f" {key:15s}: mean Δ={res['mean_delta']:+.2f} " + f"precision={res['pct_early']:.1%} coverage={res['coverage']:.1%} " + f"n_tau={res['n_tau']}") + print(" Earliest-trigger counts per channel:") + for key, cnt in sorted(winner_counts.items(), key=lambda x: -x[1]): + ml = np.mean(lead_by_ch[key]) if lead_by_ch[key] else 0 + print(f" {key:15s}: {cnt} events mean earliest lead = {ml:.1f} steps") + + # ── 7. Robustness ──────────────────────────────────────────── + print("\n Step 7 / 7 — Robustness grid (3 delay × 3 noise × 3 strength × 3 reps) ...") + print(" [This may take several minutes]") + rob_df = run_robustness_grid(model, scaler, sigma, X_scaled, X_raw, n_reps=3) + if not rob_df.empty: + rob_summary = (rob_df.groupby(['delay','noise']) + [['mean_delta','precision','recall']] + .mean().round(3)) + print("\n Robustness summary (mean over strength and reps):") + print(rob_summary.to_string()) + + # ── Summary ────────────────────────────────────────────────── + m_std = evaluation_metrics(delta_dict['Model'], sigma, tau_std) + print(f"\n ── FINAL SUMMARY ─────────────────────────────────────────") + print(f" Model Mean Δ : {m_std['mean_delta']:+.2f} steps") + print(f" Model Precision : {m_std['pct_early']:.1%}") + print(f" Model Coverage : {m_std['coverage']:.1%}") + print(f" Model Mean Δ|early: {m_std['mean_early']:+.2f} steps") + checks = [ + (m_std['mean_delta'] > 0, "Positive mean lead-time"), + (m_std['pct_early'] > 0.60, "> 60% precision"), + (m_std['coverage'] > 0.30, "> 30% coverage"), + ] + for ok, msg in checks: + print(f" {'✓' if ok else '✗'} {msg}") + print() + + # ── Figures ────────────────────────────────────────────────── + print(" Rendering figures (9 publication-quality plots) ...") + plot_dgp_causal_structure(X_raw, Z_true, delay_map) + plot_threshold_sweep(sweep_std, sweep_adp, sweep_mtr) + plot_precision_recall(sweep_std, sweep_adp, sweep_mtr, tau_imb, tau_vol, sigma) + plot_robustness_heatmap(rob_df) + plot_channel_diagnostics(comps, diag_results, winner_counts, lead_by_ch) + plot_composite_signal(X_raw, Z_true, score_amp, tau_std, sigma, + method_label="Standard") + plot_lead_time_densities(delta_dict) + plot_results_table(results_df) + plot_sweep_coverage_delta(sweep_std, sweep_adp, sweep_mtr) + + print("\n Experiment complete. Outputs saved as fig1_*.pdf … fig9_*.pdf") + + return dict( + results_df = results_df, + delta_dict = delta_dict, + model = model, + X_raw = X_raw, + Z_true = Z_true, + sigma = sigma, + delay_map = delay_map, + score_amp = score_amp, + d_score = d_score, + comps = comps, + post_smooth = post_smooth, + tau_std = tau_std, + tau_adp = tau_adp, + tau_mtr = tau_mtr, + sweep_std = sweep_std, + sweep_adp = sweep_adp, + sweep_mtr = sweep_mtr, + rob_df = rob_df, + diag_results = diag_results, + winner_counts = winner_counts, + lead_by_ch = lead_by_ch, + ) + + +if __name__ == "__main__": + results = run_experiment()