650 lines
24 KiB
Python
650 lines
24 KiB
Python
# ============================================================
|
||
# Latent Micro-Regimes in Limit Order Books:
|
||
# Identification and Early Detection
|
||
# ─────────────────────────────────────────
|
||
# Research-grade pipeline — Colab-ready single script
|
||
# Authors: [redacted for blind review]
|
||
# ============================================================
|
||
# Install (uncomment in Colab):
|
||
# !pip install hmmlearn scikit-learn scipy numpy pandas matplotlib
|
||
|
||
import warnings
|
||
warnings.filterwarnings("ignore")
|
||
|
||
import numpy as np
|
||
import pandas as pd
|
||
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 sklearn.metrics import confusion_matrix
|
||
from hmmlearn.hmm import GaussianHMM
|
||
|
||
# ─────────────────────────────────────────
|
||
# 0. Global Configuration
|
||
# ─────────────────────────────────────────
|
||
SEED = 42
|
||
T = 12_000 # total timesteps
|
||
N_REGIMES = 3 # latent states: Normal, Stressed, Crisis
|
||
MAX_LAG = 60 # evaluation window (timesteps)
|
||
FW_WINDOW = 20 # forward window for stress event definition
|
||
N_BOOT = 2_000 # bootstrap replicates
|
||
STRESS_PCT = 95 # percentile threshold for stress definition
|
||
|
||
np.random.seed(SEED)
|
||
|
||
# ─────────────────────────────────────────
|
||
# 1. Structured Latent Regime Data Generator
|
||
# ─────────────────────────────────────────
|
||
|
||
def build_transition_matrix(persist: list[float]) -> np.ndarray:
|
||
"""
|
||
Build a row-stochastic transition matrix from persistence probabilities.
|
||
Off-diagonal mass is split evenly among other states.
|
||
|
||
Parameters
|
||
----------
|
||
persist : list of length N_REGIMES
|
||
Self-transition probability for each state.
|
||
|
||
Returns
|
||
-------
|
||
P : (N_REGIMES, N_REGIMES) ndarray
|
||
"""
|
||
K = len(persist)
|
||
P = np.zeros((K, K))
|
||
for i, p in enumerate(persist):
|
||
P[i, i] = p
|
||
off = (1 - p) / (K - 1)
|
||
for j in range(K):
|
||
if j != i:
|
||
P[i, j] = off
|
||
return P
|
||
|
||
|
||
def simulate_latent_chain(T: int, P: np.ndarray, pi0: np.ndarray,
|
||
rng: np.random.Generator) -> np.ndarray:
|
||
"""
|
||
Simulate a discrete-time Markov chain.
|
||
|
||
Parameters
|
||
----------
|
||
T : int — number of timesteps
|
||
P : (K, K) transition matrix
|
||
pi0 : (K,) initial distribution
|
||
rng : numpy Generator
|
||
|
||
Returns
|
||
-------
|
||
Z : (T,) int array of latent states
|
||
"""
|
||
K = P.shape[0]
|
||
Z = np.empty(T, dtype=int)
|
||
Z[0] = rng.choice(K, p=pi0)
|
||
for t in range(1, T):
|
||
Z[t] = rng.choice(K, p=P[Z[t - 1]])
|
||
return Z
|
||
|
||
|
||
# Regime-specific parameter sets
|
||
# State 0 — Normal: tight spread, deep book, balanced flow
|
||
# State 1 — Stressed: wider spread, shallower book, directional flow
|
||
# State 2 — Crisis: very wide spread, thin book, extreme imbalance
|
||
REGIME_PARAMS = {
|
||
# spread_mu, spread_sig, depth_mu, depth_sig, imb_mu, imb_sig, vol_base
|
||
0: dict(sp_mu=1.5, sp_sig=0.30, dp_mu=120, dp_sig=12, ib_mu=0.00, ib_sig=0.08, vb=0.40),
|
||
1: dict(sp_mu=3.5, sp_sig=0.60, dp_mu= 85, dp_sig=18, ib_mu=0.18, ib_sig=0.14, vb=1.20),
|
||
2: dict(sp_mu=7.0, sp_sig=1.20, dp_mu= 45, dp_sig=22, ib_mu=0.40, ib_sig=0.22, vb=2.80),
|
||
}
|
||
|
||
|
||
def generate_lob_data(T: int, rng: np.random.Generator) -> tuple[np.ndarray, np.ndarray]:
|
||
"""
|
||
Generate synthetic LOB microstructure data driven by a latent Markov chain.
|
||
|
||
Design choices:
|
||
- Regime-conditional mean and volatility for each feature.
|
||
- Spread follows a regime-switching log-normal process (always positive).
|
||
- Depth follows a regime-switching AR(1) process (mean-reverting).
|
||
- Order-flow imbalance is beta-distributed (bounded in [-1, 1]).
|
||
- Cross-sectional noise is injected so no feature alone trivially reveals
|
||
the regime.
|
||
- A Hawkes-inspired self-exciting component adds burst clustering to spread.
|
||
|
||
Returns
|
||
-------
|
||
X : (T, 5) feature matrix [spread, depth, imbalance, roll_vol, ofi]
|
||
Z : (T,) true latent state sequence
|
||
"""
|
||
# 1) Latent chain
|
||
persist = [0.97, 0.93, 0.89] # high persistence → realistic regimes
|
||
P = build_transition_matrix(persist)
|
||
pi0 = np.array([0.80, 0.15, 0.05])
|
||
Z = simulate_latent_chain(T, P, pi0, rng)
|
||
|
||
spread = np.zeros(T)
|
||
depth = np.zeros(T)
|
||
imbalance = np.zeros(T)
|
||
|
||
# ── Spread: log-normal with Hawkes-like self-excitation ──
|
||
hawkes_intensity = np.zeros(T)
|
||
decay = 0.90
|
||
for t in range(T):
|
||
p = REGIME_PARAMS[Z[t]]
|
||
hawkes_intensity[t] = (hawkes_intensity[t - 1] * decay if t > 0 else 0)
|
||
noise = rng.normal(0, p['sp_sig'])
|
||
log_sp = np.log(p['sp_mu']) + 0.15 * hawkes_intensity[t] + noise
|
||
spread[t] = np.exp(log_sp)
|
||
if spread[t] > np.exp(np.log(p['sp_mu']) + p['sp_sig']):
|
||
hawkes_intensity[t] += 0.30 # self-excite on spike
|
||
|
||
# ── Depth: mean-reverting AR(1) per regime ──
|
||
depth[0] = REGIME_PARAMS[Z[0]]['dp_mu']
|
||
phi = 0.92 # AR coefficient
|
||
for t in range(1, T):
|
||
p = REGIME_PARAMS[Z[t]]
|
||
depth[t] = phi * depth[t - 1] + (1 - phi) * p['dp_mu'] + rng.normal(0, p['dp_sig'])
|
||
depth = np.clip(depth, 5, None) # depth always positive
|
||
|
||
# ── Order-flow imbalance: truncated normal ──
|
||
for t in range(T):
|
||
p = REGIME_PARAMS[Z[t]]
|
||
raw = rng.normal(p['ib_mu'], p['ib_sig'])
|
||
imbalance[t] = np.clip(raw, -1, 1)
|
||
|
||
# ── Engineered features ──
|
||
# Rolling 20-step realised volatility of spread
|
||
roll_vol = pd.Series(spread).pct_change().rolling(20).std().fillna(0).values
|
||
|
||
# Order-flow imbalance sign × magnitude (OFI proxy)
|
||
ofi = imbalance * np.abs(np.diff(spread, prepend=spread[0]))
|
||
|
||
X = np.column_stack([spread, depth, imbalance, roll_vol, ofi])
|
||
return X, Z
|
||
|
||
|
||
# ─────────────────────────────────────────
|
||
# 2. Feature Engineering & Normalisation
|
||
# ─────────────────────────────────────────
|
||
|
||
def engineer_features(X_raw: np.ndarray) -> tuple[np.ndarray, StandardScaler]:
|
||
"""
|
||
Augment and normalise the raw feature matrix.
|
||
|
||
Raw columns : [spread, depth, imbalance, roll_vol, ofi]
|
||
Added : [spread/depth ratio, |imbalance|, cumulative ofi (50-step)]
|
||
|
||
Returns
|
||
-------
|
||
X_scaled : (T, n_features) normalised array
|
||
scaler : fitted StandardScaler (for reproducibility)
|
||
"""
|
||
spread = X_raw[:, 0]
|
||
depth = X_raw[:, 1]
|
||
imbalance = X_raw[:, 2]
|
||
roll_vol = X_raw[:, 3]
|
||
ofi = X_raw[:, 4]
|
||
|
||
spread_depth_ratio = spread / (depth + 1e-6)
|
||
abs_imb = np.abs(imbalance)
|
||
cum_ofi = pd.Series(ofi).rolling(50, min_periods=1).mean().values
|
||
|
||
X_full = np.column_stack([
|
||
spread, depth, imbalance, roll_vol, ofi,
|
||
spread_depth_ratio, abs_imb, cum_ofi
|
||
])
|
||
|
||
scaler = StandardScaler()
|
||
X_scaled = scaler.fit_transform(X_full)
|
||
return X_scaled, scaler
|
||
|
||
|
||
# ─────────────────────────────────────────
|
||
# 3. HMM Fitting with Robust Initialisation
|
||
# ─────────────────────────────────────────
|
||
|
||
def fit_hmm(X: np.ndarray, n_components: int = N_REGIMES,
|
||
n_restarts: int = 10, rng_seed: int = SEED) -> GaussianHMM:
|
||
"""
|
||
Fit a Gaussian HMM with multiple random restarts and select
|
||
the run with highest log-likelihood.
|
||
|
||
Parameters
|
||
----------
|
||
X : normalised feature matrix
|
||
n_components: number of latent states
|
||
n_restarts : number of random restarts
|
||
rng_seed : base seed
|
||
|
||
Returns
|
||
-------
|
||
best_model : GaussianHMM with highest converged log-likelihood
|
||
"""
|
||
best_score = -np.inf
|
||
best_model = None
|
||
|
||
for k in range(n_restarts):
|
||
model = GaussianHMM(
|
||
n_components = n_components,
|
||
covariance_type = "full",
|
||
n_iter = 200,
|
||
tol = 1e-5,
|
||
random_state = rng_seed + k,
|
||
init_params = "stmc",
|
||
params = "stmc",
|
||
)
|
||
try:
|
||
model.fit(X)
|
||
score = model.score(X)
|
||
if score > best_score:
|
||
best_score = score
|
||
best_model = model
|
||
except Exception:
|
||
continue
|
||
|
||
if best_model is None:
|
||
raise RuntimeError("HMM fitting failed across all restarts.")
|
||
|
||
return best_model
|
||
|
||
|
||
# ─────────────────────────────────────────
|
||
# 4. Liquidity Stress Event Definition
|
||
# ─────────────────────────────────────────
|
||
|
||
def define_stress_events(X_raw: np.ndarray, fw: int = FW_WINDOW,
|
||
pct: float = STRESS_PCT) -> np.ndarray:
|
||
"""
|
||
Identify liquidity stress events WITHOUT lookahead leakage.
|
||
|
||
A timestep t is a stress event if the mean spread over
|
||
[t+1, t+fw] exceeds the global spread 'pct'-percentile.
|
||
|
||
The threshold is computed on the FULL spread series but
|
||
the forward window ensures the label at t uses only future data.
|
||
|
||
Parameters
|
||
----------
|
||
X_raw : raw feature matrix (column 0 = spread)
|
||
fw : forward window length
|
||
pct : percentile for threshold
|
||
|
||
Returns
|
||
-------
|
||
sigma : 1-D array of stress event timestep indices
|
||
"""
|
||
spread = X_raw[:, 0]
|
||
threshold = np.percentile(spread, pct)
|
||
sigma = []
|
||
for t in range(len(spread) - fw):
|
||
if np.mean(spread[t + 1: t + fw + 1]) > threshold:
|
||
sigma.append(t)
|
||
return np.array(sigma, dtype=int)
|
||
|
||
|
||
# ─────────────────────────────────────────
|
||
# 5. Regime Transition Detection
|
||
# ─────────────────────────────────────────
|
||
|
||
def detect_transitions(Z_hat: np.ndarray) -> np.ndarray:
|
||
"""
|
||
Return indices just *before* each detected regime change.
|
||
|
||
Parameters
|
||
----------
|
||
Z_hat : (T,) array of inferred (or true) state labels
|
||
|
||
Returns
|
||
-------
|
||
tau : 1-D int array of transition indices
|
||
"""
|
||
return np.where(np.diff(Z_hat) != 0)[0]
|
||
|
||
|
||
# ─────────────────────────────────────────
|
||
# 6. Lead-Time Evaluation
|
||
# ─────────────────────────────────────────
|
||
|
||
PENALTY = -MAX_LAG # assigned delta when no stress event found in window
|
||
|
||
|
||
def compute_lead_times(tau: np.ndarray, sigma: np.ndarray,
|
||
max_lag: int = MAX_LAG) -> np.ndarray:
|
||
"""
|
||
For each detected transition τ, find the first stress event σ
|
||
in (τ, τ + max_lag].
|
||
|
||
Returns
|
||
-------
|
||
deltas : (len(tau),) array
|
||
Positive → transition preceded stress (early detection)
|
||
PENALTY → no stress event in window (missed / false alarm)
|
||
"""
|
||
deltas = np.empty(len(tau), dtype=float)
|
||
for i, t in enumerate(tau):
|
||
candidates = sigma[(sigma > t) & (sigma <= t + max_lag)]
|
||
deltas[i] = (candidates[0] - t) if len(candidates) > 0 else PENALTY
|
||
return deltas
|
||
|
||
|
||
def evaluation_metrics(deltas: np.ndarray, max_lag: int = MAX_LAG
|
||
) -> dict:
|
||
"""
|
||
Summarise lead-time performance.
|
||
|
||
Metrics
|
||
-------
|
||
mean_delta : mean lead time (penalised)
|
||
pct_early : fraction of transitions that preceded a stress event
|
||
mean_early : mean lead time conditional on early detection
|
||
std_delta : standard deviation of delta
|
||
n_tau : total number of detected transitions
|
||
n_early : number of early detections
|
||
"""
|
||
valid = deltas > 0
|
||
return dict(
|
||
mean_delta = float(np.mean(deltas)),
|
||
pct_early = float(np.mean(valid)),
|
||
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()),
|
||
)
|
||
|
||
|
||
# ─────────────────────────────────────────
|
||
# 7. Baseline Detectors
|
||
# ─────────────────────────────────────────
|
||
|
||
def imbalance_baseline(X_raw: np.ndarray, pct: float = 90) -> np.ndarray:
|
||
"""
|
||
Trigger when |order-flow imbalance| exceeds a percentile threshold.
|
||
Noise is NOT added here; the baseline uses the same raw features as the
|
||
HMM to ensure a fair comparison.
|
||
"""
|
||
imb = np.abs(X_raw[:, 2])
|
||
threshold = np.percentile(imb, pct)
|
||
return np.where(imb > threshold)[0]
|
||
|
||
|
||
def volatility_baseline(X_raw: np.ndarray, pct: float = 90) -> np.ndarray:
|
||
"""
|
||
Trigger when rolling spread volatility exceeds a percentile threshold.
|
||
"""
|
||
roll_vol = X_raw[:, 3]
|
||
threshold = np.percentile(roll_vol, pct)
|
||
return np.where(roll_vol > threshold)[0]
|
||
|
||
|
||
# ─────────────────────────────────────────
|
||
# 8. Bootstrap Confidence Intervals
|
||
# ─────────────────────────────────────────
|
||
|
||
def bootstrap_ci(deltas: np.ndarray, stat_fn=np.mean,
|
||
n_boot: int = N_BOOT, alpha: float = 0.05,
|
||
seed: int = SEED) -> tuple[float, float]:
|
||
"""
|
||
Percentile bootstrap confidence interval for a scalar statistic.
|
||
|
||
Returns
|
||
-------
|
||
(lower, upper) CI at (1-alpha) level
|
||
"""
|
||
rng = np.random.default_rng(seed)
|
||
boot_stats = np.array([
|
||
stat_fn(rng.choice(deltas, size=len(deltas), replace=True))
|
||
for _ in range(n_boot)
|
||
])
|
||
return (float(np.percentile(boot_stats, 100 * alpha / 2)),
|
||
float(np.percentile(boot_stats, 100 * (1 - alpha / 2))))
|
||
|
||
|
||
def mannwhitney_test(a: np.ndarray, b: np.ndarray) -> tuple[float, float]:
|
||
"""
|
||
Two-sided Mann–Whitney U test (non-parametric, appropriate for
|
||
skewed lead-time distributions).
|
||
|
||
Returns (U-statistic, p-value).
|
||
"""
|
||
return stats.mannwhitneyu(a, b, alternative="two-sided")
|
||
|
||
|
||
# ─────────────────────────────────────────
|
||
# 9. Publication-Quality Visualisation
|
||
# ─────────────────────────────────────────
|
||
|
||
PALETTE = {
|
||
"Model" : "#2C6FAC",
|
||
"Imbalance" : "#D94F3D",
|
||
"Volatility" : "#5AAE61",
|
||
}
|
||
|
||
plt.rcParams.update({
|
||
"font.family" : "serif",
|
||
"font.size" : 11,
|
||
"axes.spines.top" : False,
|
||
"axes.spines.right": False,
|
||
"axes.linewidth" : 0.8,
|
||
"xtick.major.width": 0.8,
|
||
"ytick.major.width": 0.8,
|
||
"figure.dpi" : 150,
|
||
})
|
||
|
||
|
||
def plot_lead_time_densities(delta_dict: dict, max_lag: int = MAX_LAG,
|
||
save_path: str = None):
|
||
"""
|
||
Overlapping KDE plots of lead-time distributions for each detector.
|
||
|
||
Parameters
|
||
----------
|
||
delta_dict : {'Model': deltas, 'Imbalance': deltas, 'Volatility': deltas}
|
||
max_lag : used to shade the penalty region
|
||
save_path : if provided, saves the figure to this path
|
||
"""
|
||
fig, ax = plt.subplots(figsize=(8, 4.5))
|
||
|
||
x_grid = np.linspace(-max_lag - 5, max_lag + 5, 500)
|
||
|
||
for name, deltas in delta_dict.items():
|
||
color = PALETTE[name]
|
||
# KDE on non-penalty values only (for readability)
|
||
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)
|
||
|
||
# Mark penalty mass as a tick at the left edge
|
||
penalty_frac = np.mean(deltas <= PENALTY)
|
||
ax.annotate(
|
||
f" missed={penalty_frac:.1%}",
|
||
xy=(-max_lag, 0),
|
||
xytext=(-max_lag + 2, 0.012 * (list(delta_dict).index(name) + 1)),
|
||
color=color, fontsize=8.5, va="center"
|
||
)
|
||
|
||
ax.axvline(0, color="gray", lw=1.0, ls="--", label="Zero lead-time")
|
||
ax.set_xlabel("Lead time Δ (timesteps)", labelpad=8)
|
||
ax.set_ylabel("Density", labelpad=8)
|
||
ax.set_title("Lead-Time Distribution of Regime Detectors", fontsize=13, pad=10)
|
||
ax.legend(frameon=False, fontsize=10)
|
||
ax.set_xlim(-max_lag - 2, max_lag + 2)
|
||
fig.tight_layout()
|
||
if save_path:
|
||
fig.savefig(save_path, bbox_inches="tight")
|
||
plt.show()
|
||
|
||
|
||
def plot_regime_overlay(X_raw: np.ndarray, Z_true: np.ndarray,
|
||
Z_hat: np.ndarray, sigma: np.ndarray,
|
||
n_show: int = 3000, save_path: str = None):
|
||
"""
|
||
Three-panel figure:
|
||
Top — bid-ask spread with true regime shading
|
||
Middle — inferred HMM state sequence
|
||
Bottom — depth time-series with stress events marked
|
||
"""
|
||
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]
|
||
sigma_visible = sigma[sigma < t_end]
|
||
|
||
regime_colors = {0: "#DDEEFF", 1: "#FFEEDD", 2: "#FFDDDD"}
|
||
|
||
fig, axes = plt.subplots(3, 1, figsize=(11, 7), sharex=True,
|
||
gridspec_kw={"height_ratios": [2, 1, 2]})
|
||
|
||
# ── Panel 1: Spread + true regime shading ──
|
||
ax = axes[0]
|
||
for k, color in regime_colors.items():
|
||
mask = Z_true[:t_end] == k
|
||
ax.fill_between(t_ax, 0, spread.max() * 1.1,
|
||
where=mask, color=color, alpha=0.5,
|
||
label=f"True state {k}")
|
||
ax.plot(t_ax, spread, lw=0.7, color="#1A1A2E")
|
||
ax.set_ylabel("Bid-Ask Spread")
|
||
ax.legend(loc="upper right", fontsize=8, frameon=False, ncol=3)
|
||
|
||
# ── Panel 2: Inferred HMM states ──
|
||
ax = axes[1]
|
||
ax.step(t_ax, Z_hat[:t_end], lw=0.9, color=PALETTE["Model"])
|
||
ax.set_ylabel("HMM State")
|
||
ax.set_yticks([0, 1, 2])
|
||
|
||
# ── Panel 3: Depth + stress events ──
|
||
ax = axes[2]
|
||
ax.plot(t_ax, depth, lw=0.7, color="#2D6A4F")
|
||
ax.vlines(sigma_visible, depth.min(), depth.max(),
|
||
color=PALETTE["Imbalance"], lw=0.6, alpha=0.5, label="Stress event σ")
|
||
ax.set_ylabel("Market Depth")
|
||
ax.set_xlabel("Timestep")
|
||
ax.legend(loc="upper right", fontsize=8, frameon=False)
|
||
|
||
fig.suptitle("Simulated LOB: True Regimes, Inferred States & Stress Events",
|
||
fontsize=13, y=1.01)
|
||
fig.tight_layout()
|
||
if save_path:
|
||
fig.savefig(save_path, bbox_inches="tight")
|
||
plt.show()
|
||
|
||
|
||
def plot_results_table(results_df: pd.DataFrame):
|
||
"""Render results DataFrame as a styled matplotlib table."""
|
||
fig, ax = plt.subplots(figsize=(10, 2.2))
|
||
ax.axis("off")
|
||
col_labels = results_df.columns.tolist()
|
||
rows = results_df.values.tolist()
|
||
tbl = ax.table(cellText=rows, colLabels=col_labels,
|
||
cellLoc="center", loc="center")
|
||
tbl.auto_set_font_size(False)
|
||
tbl.set_fontsize(9.5)
|
||
tbl.scale(1.2, 1.5)
|
||
# Highlight header
|
||
for j in range(len(col_labels)):
|
||
tbl[0, j].set_facecolor("#2C6FAC")
|
||
tbl[0, j].set_text_props(color="white", fontweight="bold")
|
||
fig.tight_layout()
|
||
plt.show()
|
||
|
||
|
||
# ─────────────────────────────────────────
|
||
# 10. Main Experimental Pipeline
|
||
# ─────────────────────────────────────────
|
||
|
||
def run_experiment():
|
||
rng = np.random.default_rng(SEED)
|
||
|
||
# ── 10.1 Data generation ──
|
||
print("=" * 60)
|
||
print(" Step 1 / 6 — Generating structured LOB data …")
|
||
X_raw, Z_true = generate_lob_data(T, rng)
|
||
print(f" Generated {T:,} timesteps | "
|
||
f"Regime distribution: "
|
||
+ " | ".join([f"State {k}: {(Z_true==k).mean():.1%}"
|
||
for k in range(N_REGIMES)]))
|
||
|
||
# ── 10.2 Feature engineering ──
|
||
print("\n Step 2 / 6 — Engineering and normalising features …")
|
||
X_scaled, scaler = engineer_features(X_raw)
|
||
print(f" Feature matrix shape: {X_scaled.shape}")
|
||
|
||
# ── 10.3 HMM fitting ──
|
||
print("\n Step 3 / 6 — Fitting Gaussian HMM (multiple restarts) …")
|
||
model = fit_hmm(X_scaled)
|
||
Z_hat = model.predict(X_scaled)
|
||
ll = model.score(X_scaled)
|
||
print(f" Best log-likelihood: {ll:,.2f} | "
|
||
f"Converged: {model.monitor_.converged}")
|
||
|
||
# ── 10.4 Stress event definition ──
|
||
print("\n Step 4 / 6 — Defining leakage-free stress events …")
|
||
sigma = define_stress_events(X_raw)
|
||
print(f" Stress events identified: {len(sigma):,} "
|
||
f"({len(sigma)/T:.1%} of timesteps)")
|
||
|
||
# ── 10.5 Transition detection for all detectors ──
|
||
print("\n Step 5 / 6 — Evaluating detectors …")
|
||
|
||
tau_model = detect_transitions(Z_hat)
|
||
tau_imb = imbalance_baseline(X_raw)
|
||
tau_vol = volatility_baseline(X_raw)
|
||
|
||
delta_model = compute_lead_times(tau_model, sigma)
|
||
delta_imb = compute_lead_times(tau_imb, sigma)
|
||
delta_vol = compute_lead_times(tau_vol, sigma)
|
||
|
||
delta_dict = {
|
||
"Model" : delta_model,
|
||
"Imbalance" : delta_imb,
|
||
"Volatility": delta_vol,
|
||
}
|
||
|
||
# ── 10.6 Statistical validation ──
|
||
print("\n Step 6 / 6 — Statistical validation (bootstrap + MWU) …")
|
||
|
||
rows = []
|
||
for name, deltas in delta_dict.items():
|
||
m = evaluation_metrics(deltas)
|
||
lo, hi = bootstrap_ci(deltas, stat_fn=np.mean)
|
||
rows.append({
|
||
"Detector" : name,
|
||
"Mean Δ" : f"{m['mean_delta']:+.2f}",
|
||
"95 % CI" : f"[{lo:+.2f}, {hi:+.2f}]",
|
||
"% Early" : f"{m['pct_early']:.1%}",
|
||
"Mean Δ | early" : f"{m['mean_early']:+.2f}",
|
||
"Std Δ" : f"{m['std_delta']:.2f}",
|
||
"N(τ)" : m['n_tau'],
|
||
"N(early)" : m['n_early'],
|
||
})
|
||
|
||
results_df = pd.DataFrame(rows)
|
||
print("\n" + results_df.to_string(index=False))
|
||
|
||
# Pairwise Mann–Whitney tests
|
||
print("\n Pairwise Mann–Whitney U tests (two-sided):")
|
||
pairs = [("Model", "Imbalance"), ("Model", "Volatility"),
|
||
("Imbalance", "Volatility")]
|
||
for a_name, b_name in pairs:
|
||
u, p = mannwhitney_test(delta_dict[a_name], delta_dict[b_name])
|
||
sig = "***" if p < 0.001 else ("**" if p < 0.01 else ("*" if p < 0.05 else "ns"))
|
||
print(f" {a_name:12s} vs {b_name:12s}: U={u:,.0f} p={p:.4f} {sig}")
|
||
|
||
# ── 10.7 Plots ──
|
||
print("\n Rendering publication-quality figures …")
|
||
plot_regime_overlay(X_raw, Z_true, Z_hat, sigma)
|
||
plot_lead_time_densities(delta_dict)
|
||
plot_results_table(results_df)
|
||
|
||
print("\n Experiment complete.")
|
||
return results_df, delta_dict, model, X_raw, Z_true, Z_hat, sigma
|
||
|
||
|
||
# ─────────────────────────────────────────
|
||
# Entry point
|
||
# ─────────────────────────────────────────
|
||
if __name__ == "__main__":
|
||
results_df, delta_dict, model, X_raw, Z_true, Z_hat, sigma = run_experiment()
|