1207 lines
55 KiB
Python
1207 lines
55 KiB
Python
# ============================================================
|
||
# 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()
|