# ============================================================ # 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)), precision = 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['precision'], '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['precision'], '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, precision=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'], precision = m['precision'], 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', 'precision', '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['precision']:.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['precision'].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['precision']:.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['precision']:.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['precision'] > 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()