jack/lob-latent-regimes
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add research-grade pipeline for LOB analysis

Browse Source
This commit is contained in:
PRAKUL HIREMATH
2026-04-10 17:46:45 +05:30
committed by GitHub
commit af106ef5d0
1 changed files with 649 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Experiments/v1
+649
View File
@@ -0,0 +1,649 @@
# ============================================================
# 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()