%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
class LimitOrderBook:
"""Price-time priority LOB matching engine."""
def __init__(self, tick_size=0.01):
self.tick_size = tick_size
self.bids = {} # price -> volume
self.asks = {} # price -> volume
self.best_bid = None
self.best_ask = None
self.mid_price = None
def _snap(self, p):
return round(round(p / self.tick_size) * self.tick_size, 8)
def _refresh(self):
self.bids = {p: v for p, v in self.bids.items() if v > 1e-9}
self.asks = {p: v for p, v in self.asks.items() if v > 1e-9}
self.best_bid = max(self.bids) if self.bids else None
self.best_ask = min(self.asks) if self.asks else None
if self.best_bid is not None and self.best_ask is not None:
self.mid_price = (self.best_bid + self.best_ask) / 2.0
def add_limit(self, side, price, volume):
price = self._snap(price)
if side == 'buy':
while volume > 1e-9 and self.asks and min(self.asks) <= price:
best = min(self.asks)
trade = min(self.asks[best], volume)
self.asks[best] -= trade
volume -= trade
self._refresh()
if volume > 1e-9:
self.bids[price] = self.bids.get(price, 0) + volume
else:
while volume > 1e-9 and self.bids and max(self.bids) >= price:
best = max(self.bids)
trade = min(self.bids[best], volume)
self.bids[best] -= trade
volume -= trade
self._refresh()
if volume > 1e-9:
self.asks[price] = self.asks.get(price, 0) + volume
self._refresh()
def add_market(self, side, volume):
filled = 0
if side == 'buy':
while volume > 1e-9 and self.asks:
best = min(self.asks)
trade = min(self.asks[best], volume)
self.asks[best] -= trade
volume -= trade
filled += trade
else:
while volume > 1e-9 and self.bids:
best = max(self.bids)
trade = min(self.bids[best], volume)
self.bids[best] -= trade
volume -= trade
filled += trade
self._refresh()
return filled
def cancel(self, side, price, vol=None):
price = self._snap(price)
book = self.bids if side == 'buy' else self.asks
if price in book:
if vol is None or vol >= book[price]:
del book[price]
else:
book[price] = max(0.0, book[price] - vol)
self._refresh()
def snapshot(self, n=12):
bids = sorted([(p, v) for p, v in self.bids.items()], reverse=True)[:n]
asks = sorted([(p, v) for p, v in self.asks.items()])[:n]
return bids, asks
def plot_lob_dynamics(snapshots, title='LOB Dynamics'):
times = [s['t'] for s in snapshots]
mids = [s['mid'] for s in snapshots]
t_b, p_b, v_b = [], [], []
t_a, p_a, v_a = [], [], []
for s in snapshots:
for price, vol in s['bids']:
t_b.append(s['t']); p_b.append(price); v_b.append(vol)
for price, vol in s['asks']:
t_a.append(s['t']); p_a.append(price); v_a.append(vol)
t_b = np.array(t_b); p_b = np.array(p_b); v_b = np.array(v_b)
t_a = np.array(t_a); p_a = np.array(p_a); v_a = np.array(v_a)
all_v = np.concatenate([v_b, v_a]) if (len(v_b) and len(v_a)) else np.array([1.0])
vref = max(np.percentile(all_v, 90), 1e-9)
def sz(v): return np.clip(v / vref, 0.05, 1.5) * 50 + 5
fig, (ax1, ax2) = plt.subplots(
2, 1, figsize=(12, 8),
gridspec_kw={'height_ratios': [1, 3]}, sharex=True
)
fig.subplots_adjust(hspace=0.06)
ax1.plot(times, mids, color='black', linewidth=1.2)
ax1.set_ylabel('Mid-price')
ax1.set_title(title, fontsize=11)
if len(t_b):
ax2.scatter(t_b, p_b, s=sz(v_b), c='#2ca02c', alpha=0.35,
linewidths=0, label='Bid limit orders')
if len(t_a):
ax2.scatter(t_a, p_a, s=sz(v_a), c='#d62728', alpha=0.35,
linewidths=0, label='Ask limit orders')
ax2.plot(times, mids, color='black', linewidth=1.0,
linestyle='--', alpha=0.7, label='Mid-price')
ax2.set_xlabel('Time (periods)')
ax2.set_ylabel('Price')
ax2.legend(loc='upper right', fontsize=9)
plt.show()
1 Models for order arrival¶
Poisson and Hawkes processes¶
The Poisson process with constant intensity is the baseline model for order arrival: inter-arrival times are and arrivals in disjoint windows are independent. The Hawkes self-exciting process adds a self-excitation term — each event temporarily raises the arrival rate by before decaying at rate :
The process is stationary when , with stationary mean .
Below we simulate both models, compare their inter-arrival distributions, and visualise the autocorrelation in the event counts — the key empirical signature of self-excitation.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from collections import deque
FIGURES_DIR = '../markdown/figures'
def simulate_poisson(lam, T, rng):
"""Homogeneous Poisson process on [0, T]."""
n = rng.poisson(lam * T)
return np.sort(rng.uniform(0, T, n))
def simulate_hawkes(mu, phi, beta, T, rng):
"""
Univariate Hawkes process via Ogata's thinning algorithm.
Stationary when phi < beta.
"""
events = []
t = 0.0
lam_t = mu
while t < T:
# Upper bound: current intensity
lam_hat = lam_t
dt = rng.exponential(1.0 / lam_hat)
t += dt
if t > T:
break
# Recompute exact intensity at proposed time
lam_exact = mu + phi * sum(np.exp(-beta * (t - s)) for s in events)
if rng.random() < lam_exact / lam_hat:
events.append(t)
lam_t = lam_exact + phi # intensity jumps at accepted event
else:
lam_t = lam_exact # no event: update upper bound
return np.array(events)
rng = np.random.default_rng(42)
T = 600.0 # 10-minute window
# Poisson with same mean rate as Hawkes stationary mean
mu, phi, beta = 0.5, 1.5, 2.5
lam_stat = mu * beta / (beta - phi) # stationary mean of Hawkes
ev_poisson = simulate_poisson(lam_stat, T, rng)
ev_hawkes = simulate_hawkes(mu, phi, beta, T, rng)
print(f"Stationary mean rate: {lam_stat:.3f} events/s")
print(f"Poisson realisation: {len(ev_poisson)} events (expected {lam_stat*T:.0f})")
print(f"Hawkes realisation: {len(ev_hawkes)} events")fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))
# Inter-arrival distributions
iat_poisson = np.diff(ev_poisson)
iat_hawkes = np.diff(ev_hawkes)
bins = np.linspace(0, np.percentile(np.concatenate([iat_poisson, iat_hawkes]), 95), 50)
axes[0].hist(iat_poisson, bins=bins, alpha=0.6, density=True, color='#1565C0', label='Poisson')
axes[0].hist(iat_hawkes, bins=bins, alpha=0.6, density=True, color='#C62828', label='Hawkes')
x_plot = np.linspace(1e-6, bins[-1], 300)
axes[0].plot(x_plot, lam_stat * np.exp(-lam_stat * x_plot), '--k', lw=1.5, label='$\\mathrm{Exp}(\\bar{\\lambda})$')
axes[0].set_xlabel('Inter-arrival time (s)'); axes[0].set_ylabel('Density')
axes[0].set_title('Inter-arrival distributions'); axes[0].legend(); axes[0].grid(alpha=0.3)
# Cumulative event count
axes[1].step(ev_poisson, np.arange(1, len(ev_poisson)+1), where='post', color='#1565C0', lw=1.5, label='Poisson')
axes[1].step(ev_hawkes, np.arange(1, len(ev_hawkes)+1), where='post', color='#C62828', lw=1.5, label='Hawkes')
axes[1].set_xlabel('Time (s)'); axes[1].set_ylabel('Cumulative events')
axes[1].set_title('Cumulative event counts'); axes[1].legend(); axes[1].grid(alpha=0.3)
# Autocorrelation of 1-second bucket counts
dt_bin = 1.0
n_bins = int(T / dt_bin)
bins_t = np.arange(n_bins + 1) * dt_bin
counts_p = np.histogram(ev_poisson, bins=bins_t)[0].astype(float)
counts_h = np.histogram(ev_hawkes, bins=bins_t)[0].astype(float)
def acf(x, nlags=30):
x = x - x.mean()
var = (x**2).mean()
return [np.mean(x[k:] * x[:len(x)-k]) / var for k in range(nlags + 1)]
lags = np.arange(31)
axes[2].bar(lags - 0.2, acf(counts_p), width=0.35, color='#1565C0', alpha=0.7, label='Poisson')
axes[2].bar(lags + 0.2, acf(counts_h), width=0.35, color='#C62828', alpha=0.7, label='Hawkes')
axes[2].axhline(0, color='k', lw=0.8)
axes[2].set_xlabel('Lag (s)'); axes[2].set_ylabel('Autocorrelation')
axes[2].set_title('1-second bucket count ACF (clustering signature)'); axes[2].legend(); axes[2].grid(alpha=0.3)
plt.tight_layout()
plt.savefig(f'{FIGURES_DIR}/lob_order_arrival.png', dpi=150, bbox_inches='tight')
plt.show()2 LOB features¶
We compute the standard LOB features — spread, order imbalance, and micro-volatility — from the probabilistic generative simulation defined in the Simulation section below.
The key predictive feature is order imbalance at the best quote:
We verify empirically that predicts the sign of the subsequent mid-price move — a cornerstone result of high-frequency microstructure.
# Build a snapshot dataset from the probabilistic generative model (same params as below)
rng_lob = np.random.default_rng(7)
from collections import defaultdict
class _LOB:
"""Minimal price-time LOB for feature extraction."""
def __init__(self, tick=0.5):
self.tick = tick; self.bids = {}; self.asks = {}
self.best_bid = self.best_ask = self.mid = None
def _snap(self, p): return round(round(p/self.tick)*self.tick, 8)
def _refresh(self):
self.bids = {p:v for p,v in self.bids.items() if v > 1e-9}
self.asks = {p:v for p,v in self.asks.items() if v > 1e-9}
self.best_bid = max(self.bids) if self.bids else None
self.best_ask = min(self.asks) if self.asks else None
if self.best_bid and self.best_ask:
self.mid = (self.best_bid + self.best_ask) / 2
def limit(self, side, price, vol):
price = self._snap(price)
if side == 'buy':
while vol > 1e-9 and self.asks and min(self.asks) <= price:
b = min(self.asks); tr = min(self.asks[b], vol); self.asks[b]-=tr; vol-=tr; self._refresh()
if vol > 1e-9: self.bids[price] = self.bids.get(price, 0) + vol
else:
while vol > 1e-9 and self.bids and max(self.bids) >= price:
b = max(self.bids); tr = min(self.bids[b], vol); self.bids[b]-=tr; vol-=tr; self._refresh()
if vol > 1e-9: self.asks[price] = self.asks.get(price, 0) + vol
self._refresh()
def market(self, side, vol):
if side == 'buy':
while vol > 1e-9 and self.asks:
b = min(self.asks); tr = min(self.asks[b], vol); self.asks[b]-=tr; vol-=tr
else:
while vol > 1e-9 and self.bids:
b = max(self.bids); tr = min(self.bids[b], vol); self.bids[b]-=tr; vol-=tr
self._refresh()
P0, tick = 100.0, 0.5
lob2 = _LOB(tick)
for i in range(10):
lob2.bids[round(P0-(i+.5)*tick,8)] = max(1, int(rng_lob.lognormal(2.2, 0.7)))
lob2.asks[round(P0+(i+.5)*tick,8)] = max(1, int(rng_lob.lognormal(2.2, 0.7)))
lob2._refresh()
snaps = []
for t in range(5000):
if lob2.best_bid is None or lob2.best_ask is None: break
for _ in range(int(rng_lob.poisson(2))):
if lob2.best_bid is None or lob2.best_ask is None: break
side = 'buy' if rng_lob.random() < 0.5 else 'sell'
sz = max(1, int(rng_lob.lognormal(2.2, 0.7)))
if rng_lob.random() < 0.25:
opp = sum(lob2.asks.values()) if side=='buy' else sum(lob2.bids.values())
lob2.market(side, min(sz, max(1, int(opp*0.4))))
else:
n_ticks = min(int(rng_lob.geometric(0.45))-1, 20)
price = (lob2.best_bid - n_ticks*tick) if side=='buy' else (lob2.best_ask + n_ticks*tick)
lob2.limit(side, price, sz)
if t % 5 == 0 and lob2.mid is not None and lob2.bids and lob2.asks:
vb = lob2.bids.get(lob2.best_bid, 0)
va = lob2.asks.get(lob2.best_ask, 0)
snaps.append({'t': t, 'mid': lob2.mid,
'spread': lob2.best_ask - lob2.best_bid,
'OI': (vb - va) / (vb + va) if (vb + va) > 0 else 0})
snaps = snaps[:-1] # drop last (no next mid)
OI = np.array([s['OI'] for s in snaps])
mid = np.array([s['mid'] for s in snaps])
spr = np.array([s['spread'] for s in snaps])
dmid = np.diff(mid) # next mid-price change
print(f"Snapshots: {len(snaps)}, mean spread: {spr.mean():.3f}, mean |OI|: {np.abs(OI[:-1]).mean():.3f}")# Accuracy of OI as a directional predictor
sign_actual = np.sign(dmid)
sign_predicted = np.sign(OI[:-1])
mask = sign_actual != 0 # ignore flat ticks
accuracy = (sign_actual[mask] == sign_predicted[mask]).mean()
print(f"OI → next-tick direction accuracy: {accuracy*100:.1f}%")
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
# Left: scatter OI vs next mid-price change
axes[0].scatter(OI[:-1], dmid, s=4, alpha=0.25, color='#1565C0')
# Fit linear trend
c = np.polyfit(OI[:-1], dmid, 1)
x_line = np.array([-1, 1])
axes[0].plot(x_line, np.polyval(c, x_line), 'r-', lw=2,
label=f'slope = {c[0]:.4f}')
axes[0].axhline(0, color='k', lw=0.8); axes[0].axvline(0, color='k', lw=0.8)
axes[0].set_xlabel('Order imbalance $\\mathrm{OI}_t$')
axes[0].set_ylabel('Next mid-price change $\\Delta M_{t+1}$')
axes[0].set_title('OI as a predictor of mid-price direction')
axes[0].legend(); axes[0].grid(alpha=0.3)
# Right: conditional mean change by OI bucket
n_buckets = 10
oi_edges = np.percentile(OI[:-1], np.linspace(0, 100, n_buckets + 1))
bucket_means = []
bucket_centres = []
for i in range(n_buckets):
mask_b = (OI[:-1] >= oi_edges[i]) & (OI[:-1] < oi_edges[i+1])
if mask_b.sum() > 5:
bucket_means.append(dmid[mask_b].mean())
bucket_centres.append((oi_edges[i] + oi_edges[i+1]) / 2)
axes[1].bar(range(n_buckets), bucket_means, color='#2E7D32', alpha=0.75)
axes[1].axhline(0, color='k', lw=1)
axes[1].set_xticks(range(n_buckets))
axes[1].set_xticklabels([f'{c:.2f}' for c in bucket_centres], rotation=45, fontsize=8)
axes[1].set_xlabel('OI bucket (low → high)')
axes[1].set_ylabel('Mean $\\Delta M_{t+1}$')
axes[1].set_title('Conditional mean mid-price change by OI bucket')
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.savefig(f'{FIGURES_DIR}/lob_oi_prediction.png', dpi=150, bbox_inches='tight')
plt.show()3 Fill probability model¶
A sell limit order placed ticks above mid fills at a Poisson rate . The fill probability within window is:
Calibration from simulated data¶
We generate synthetic limit orders with known , observe whether each fills within a capped window (censored survival data), and recover the parameters by maximum likelihood. The log-likelihood for exponential fill times with censoring is:
from scipy.optimize import minimize
def fill_rate(delta, A, k): return A * np.exp(-k * delta)
def fill_prob(delta, tau, A, k): return 1 - np.exp(-fill_rate(delta, A, k) * tau)
# ── Generate synthetic dataset ────────────────────────────────────────────────
A_true, k_true = 2.0, 1.0
tau_max = 5.0 # observation window (minutes)
N = 3000
rng_fp = np.random.default_rng(99)
deltas = rng_fp.uniform(0, 5, N)
fill_times = rng_fp.exponential(1 / fill_rate(deltas, A_true, k_true))
filled = fill_times <= tau_max
obs_times = np.minimum(fill_times, tau_max)
print(f"Fill rate: {filled.mean()*100:.1f}% ({filled.sum()} of {N})")
# ── MLE via negative log-likelihood ──────────────────────────────────────────
def neg_ll(params):
A, k = np.exp(params) # reparameterise to ensure positivity
lam = fill_rate(deltas, A, k)
return -np.sum(filled * np.log(lam + 1e-12) - lam * obs_times)
res = minimize(neg_ll, x0=[np.log(1.0), np.log(1.0)],
method='Nelder-Mead', options={'xatol': 1e-6, 'fatol': 1e-6, 'maxiter': 5000})
A_hat, k_hat = np.exp(res.x)
print(f"True: A = {A_true:.3f}, k = {k_true:.3f}")
print(f"Estimated: A = {A_hat:.3f}, k = {k_hat:.3f}")fig, axes = plt.subplots(1, 2, figsize=(13, 5))
# Left: fill rate curves
d_plot = np.linspace(0, 5, 200)
k_vals = [0.5, 1.0, 2.0]
clrs = ['#1565C0', '#2E7D32', '#C62828']
for k, col in zip(k_vals, clrs):
axes[0].plot(d_plot, fill_rate(d_plot, 1.0, k), lw=2.2, color=col, label=f'$k={k}$')
axes[0].set_xlabel('Placement depth $\\delta$ (ticks)')
axes[0].set_ylabel('Fill rate $\\lambda(\\delta)$ (fills/min)')
axes[0].set_title('Fill rate vs depth ($A = 1$)')
axes[0].legend(); axes[0].grid(alpha=0.3)
# Right: estimated vs true fill probabilities across depths
d_eval = np.linspace(0, 5, 100)
tau_show = 3.0
axes[1].plot(d_eval, fill_prob(d_eval, tau_show, A_true, k_true), 'k-', lw=2.5, label='True')
axes[1].plot(d_eval, fill_prob(d_eval, tau_show, A_hat, k_hat), 'r--', lw=2.0, label='MLE estimate')
axes[1].set_xlabel('Placement depth $\\delta$ (ticks)')
axes[1].set_ylabel(f'$P(\\mathrm{{fill}}\\mid\\delta,\\tau={tau_show:.0f}\\,\\mathrm{{min}})$')
axes[1].set_title(f'Fill probability: true vs MLE ($N={N}$ orders, {filled.sum()} filled)')
axes[1].legend(); axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.savefig(f'{FIGURES_DIR}/lob_fill_probability.png', dpi=150, bbox_inches='tight')
plt.show()4 Short-term price prediction¶
We fit a simple logistic regression of next mid-price direction on order imbalance to quantify the OI predictability demonstrated above, and show the decision boundary.
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import classification_report
# Features: [OI, OI^2] — quadratic to allow non-linear threshold
X = np.column_stack([OI[:-1], OI[:-1]**2])
y = np.sign(dmid)
# Drop flat ticks
mask_nz = y != 0
X_fit, y_fit = X[mask_nz], y[mask_nz]
lr = LogisticRegression(C=10, max_iter=500)
lr.fit(X_fit, y_fit)
print("Logistic regression on OI → next mid-price direction")
print(classification_report(y_fit, lr.predict(X_fit), target_names=['down','up'], digits=3))
# Decision boundary plot
fig, ax = plt.subplots(figsize=(8, 5))
oi_vals = np.linspace(-1, 1, 300)
prob_up = lr.predict_proba(np.column_stack([oi_vals, oi_vals**2]))[:, 1]
ax.plot(oi_vals, prob_up, 'b-', lw=2.5, label='$P(\\Delta M > 0 \\mid \\mathrm{OI})$')
ax.axhline(0.5, color='k', linestyle='--', lw=1.2, label='Decision boundary')
ax.axvline(0, color='gray', linestyle=':', lw=1.0)
ax.set_xlabel('Order imbalance $\\mathrm{OI}_t$')
ax.set_ylabel('Predicted probability of up-tick')
ax.set_title('Logistic regression: OI → next mid-price direction')
ax.set_ylim(0, 1); ax.legend(); ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()5 Simulation of LOBs¶
Probabilistic generative models¶
Simulation with Poisson order arrivals, log-normal order sizes, and geometric price placement relative to the best quote. The matching engine runs price-time priority continuously.
rng = np.random.default_rng(42)
lob = LimitOrderBook(tick_size=0.5)
P0, tick = 100.0, 0.5
for i in range(10):
lob.bids[round(P0 - (i + 0.5) * tick, 8)] = max(1, int(rng.lognormal(2.2, 0.7)))
lob.asks[round(P0 + (i + 0.5) * tick, 8)] = max(1, int(rng.lognormal(2.2, 0.7)))
lob._refresh()
snapshots_prob = []
for t in range(3000):
if lob.best_bid is None or lob.best_ask is None:
break
for _ in range(int(rng.poisson(2))):
if lob.best_bid is None or lob.best_ask is None:
break
side = 'buy' if rng.random() < 0.5 else 'sell'
sz = max(1, int(rng.lognormal(2.2, 0.7)))
if rng.random() < 0.25:
opp = sum(lob.asks.values()) if side == 'buy' else sum(lob.bids.values())
lob.add_market(side, min(sz, max(1, int(opp * 0.4))))
else:
n_ticks = min(int(rng.geometric(0.45)) - 1, 20)
price = (lob.best_bid - n_ticks * tick) if side == 'buy' \
else (lob.best_ask + n_ticks * tick)
lob.add_limit(side, price, sz)
if t % 8 == 7:
if len(lob.bids) > 20:
lob.cancel('buy', min(lob.bids))
if len(lob.asks) > 20:
lob.cancel('sell', max(lob.asks))
if t % 15 == 0 and lob.mid_price is not None:
bids, asks = lob.snapshot(12)
snapshots_prob.append({'t': t, 'bids': bids, 'asks': asks, 'mid': lob.mid_price})
plot_lob_dynamics(
snapshots_prob,
title='Probabilistic Generative Model – LOB Dynamics\n'
'(Poisson arrivals, log-normal sizes, geometric price placement)'
)
print(f'{len(snapshots_prob)} snapshots, '
f'mid-price: {snapshots_prob[0]["mid"]:.2f} → {snapshots_prob[-1]["mid"]:.2f}')
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[2], line 22
20 if rng.random() < 0.25:
21 opp = sum(lob.asks.values()) if side == 'buy' else sum(lob.bids.values())
---> 22 lob.add_market(side, min(sz, max(1, int(opp * 0.4))))
23 else:
24 n_ticks = min(int(rng.geometric(0.45)) - 1, 20)
Cell In[1], line 57, in LimitOrderBook.add_market(self, side, volume)
55 trade = min(self.asks[best], volume)
56 self.asks[best] -= trade
---> 57 volume -= trade
58 filled += trade
59 else:
KeyboardInterrupt: Agent-based models¶
McGroarty et al. (2019) five-agent simulation: market makers, liquidity consumers, momentum traders, mean reversion traders, and noise traders.
rng = np.random.default_rng(0)
lob = LimitOrderBook(tick_size=0.5)
P0, tick = 100.0, 0.5
T, snap_every = 10000, 50
for i in range(10):
lob.bids[round(P0 - (i + 0.5) * tick, 8)] = max(1, int(rng.lognormal(2.0, 0.5)))
lob.asks[round(P0 + (i + 0.5) * tick, 8)] = max(1, int(rng.lognormal(2.0, 0.5)))
lob._refresh()
# ── agent counts and action probabilities (McGroarty et al. Table 1) ──────────
n_mm, n_lc, n_mr, n_mt, n_nt = 2, 2, 5, 5, 12
d_mm, d_lc, d_mr, d_mt, d_nt = 0.10, 0.10, 0.40, 0.40, 0.75
# market maker state
mm_w = 50
mm_bid_p = [None]*n_mm; mm_bid_v = [0]*n_mm
mm_ask_p = [None]*n_mm; mm_ask_v = [0]*n_mm
order_signs = []
# liquidity consumer state
lc_side = ['buy' if rng.random() < 0.5 else 'sell' for _ in range(n_lc)]
lc_rem = [int(rng.uniform(300, 1200)) for _ in range(n_lc)]
# momentum trader state
mt_nr, mt_kappa, mt_scale = 200, 0.001, 4000
price_hist = [P0] * mt_nr
# mean reversion trader state
mr_alpha, mr_k, mr_vmr = 0.02, 1.5, 8
mr_ema = [P0]*n_mr
mr_var = [0.0]*n_mr
# noise trader state
nt_sz_mu, nt_sz_sig = 1.8, 0.65
nt_lam_m, nt_lam_l = 0.25, 0.60 # lam_c = 0.15
nt_lcrs, nt_linspr, nt_lspr = 0.15, 0.20, 0.40 # loffspr = 0.25
nt_beta, nt_xmin = 3.0, 0.5
nt_orders = [[] for _ in range(n_nt)]
snapshots_abm = []
for t in range(T):
if lob.best_bid is None or lob.best_ask is None:
break
mid = lob.mid_price
t_sgn = 0
# ── market makers ────────────────────────────────────────────────────────
for i in range(n_mm):
if rng.random() >= d_mm:
continue
if lob.best_bid is None or lob.best_ask is None:
continue
pred_buy = (np.mean(order_signs[-mm_w:]) > 0) if len(order_signs) >= mm_w \
else (rng.random() < 0.5)
if mm_bid_p[i] is not None and mm_bid_p[i] in lob.bids:
lob.cancel('buy', mm_bid_p[i], mm_bid_v[i])
if mm_ask_p[i] is not None and mm_ask_p[i] in lob.asks:
lob.cancel('sell', mm_ask_p[i], mm_ask_v[i])
if lob.best_bid is None or lob.best_ask is None:
continue
if pred_buy:
sell_v = max(1, int(rng.uniform(5, 25))); buy_v = 2
else:
buy_v = max(1, int(rng.uniform(5, 25))); sell_v = 2
bp = lob.best_bid; ap = lob.best_ask
lob.bids[bp] = lob.bids.get(bp, 0) + buy_v
lob.asks[ap] = lob.asks.get(ap, 0) + sell_v
mm_bid_p[i], mm_bid_v[i] = bp, buy_v
mm_ask_p[i], mm_ask_v[i] = ap, sell_v
lob._refresh()
# ── liquidity consumers ──────────────────────────────────────────────────
for i in range(n_lc):
if lc_rem[i] <= 0 or rng.random() >= d_lc:
continue
if lob.best_bid is None or lob.best_ask is None:
continue
side = lc_side[i]
opp = lob.asks if side == 'buy' else lob.bids
if not opp:
continue
best_opp = min(opp) if side == 'buy' else max(opp)
vol = min(lc_rem[i], opp[best_opp])
if vol > 0:
lob.add_market(side, vol)
lc_rem[i] -= vol
t_sgn = 1 if side == 'buy' else -1
# ── momentum traders ─────────────────────────────────────────────────────
price_hist.append(mid)
if len(price_hist) > mt_nr + 100:
price_hist.pop(0)
for i in range(n_mt):
if rng.random() >= d_mt or len(price_hist) < mt_nr + 1:
continue
roc = (price_hist[-1] - price_hist[-mt_nr - 1]) / price_hist[-mt_nr - 1]
if abs(roc) < mt_kappa:
continue
side = 'buy' if roc > 0 else 'sell'
vol = max(1, int(abs(roc) * mt_scale))
opp = sum(lob.asks.values()) if side == 'buy' else sum(lob.bids.values())
vol = min(vol, max(1, int(opp * 0.25)))
lob.add_market(side, vol)
t_sgn = 1 if side == 'buy' else -1
# ── mean reversion traders ───────────────────────────────────────────────
for i in range(n_mr):
if rng.random() >= d_mr:
continue
if lob.best_bid is None or lob.best_ask is None:
continue
mr_ema[i] += mr_alpha * (mid - mr_ema[i])
dev = mid - mr_ema[i]
mr_var[i] += mr_alpha * (dev * dev - mr_var[i])
sigma = max(np.sqrt(max(mr_var[i], 0.0)), tick * 0.5)
if dev >= mr_k * sigma:
p = lob.best_ask - tick
if p > lob.best_bid:
lob.add_limit('sell', p, mr_vmr)
else:
lob.add_limit('sell', lob.best_ask, mr_vmr)
elif dev <= -mr_k * sigma:
p = lob.best_bid + tick
if p < lob.best_ask:
lob.add_limit('buy', p, mr_vmr)
else:
lob.add_limit('buy', lob.best_bid, mr_vmr)
# ── noise traders ────────────────────────────────────────────────────────
for i in range(n_nt):
if rng.random() >= d_nt:
continue
if lob.best_bid is None or lob.best_ask is None:
continue
side = 'buy' if rng.random() < 0.5 else 'sell'
sz = max(1, int(rng.lognormal(nt_sz_mu, nt_sz_sig)))
r = rng.random()
if r < nt_lam_m:
opp = sum(lob.asks.values()) if side == 'buy' else sum(lob.bids.values())
lob.add_market(side, min(sz, max(1, int(opp * 0.5))))
t_sgn = 1 if side == 'buy' else -1
elif r < nt_lam_m + nt_lam_l:
rl = rng.random()
if rl < nt_lcrs:
price = lob.best_ask if side == 'buy' else lob.best_bid
elif rl < nt_lcrs + nt_linspr:
price = round(rng.uniform(lob.best_bid, lob.best_ask) / tick) * tick
elif rl < nt_lcrs + nt_linspr + nt_lspr:
price = lob.best_bid if side == 'buy' else lob.best_ask
else:
u = max(rng.random(), 1e-9)
offset = nt_xmin * (1.0 - u) ** (-1.0 / (nt_beta - 1))
offset = min(offset, 20 * tick)
price = (lob.best_bid - offset) if side == 'buy' \
else (lob.best_ask + offset)
price = round(price / tick) * tick
lob.add_limit(side, price, sz)
nt_orders[i].append((side, price))
t_sgn = 1 if side == 'buy' else -1
else:
if nt_orders[i]:
s, p = nt_orders[i].pop(0)
lob.cancel(s, p)
if not lob.bids and lob.best_ask is not None:
lob.bids[lob.best_ask - tick] = 5; lob._refresh()
if not lob.asks and lob.best_bid is not None:
lob.asks[lob.best_bid + tick] = 5; lob._refresh()
order_signs.append(t_sgn)
if t % snap_every == 0 and lob.mid_price is not None:
bids, asks = lob.snapshot(12)
snapshots_abm.append({'t': t, 'bids': bids, 'asks': asks, 'mid': lob.mid_price})
plot_lob_dynamics(
snapshots_abm,
title='Agent-Based Model (McGroarty et al. 2019) – LOB Dynamics\n'
'(Market makers · Liquidity consumers · Momentum · Mean reversion · Noise traders)'
)
print(f'{len(snapshots_abm)} snapshots, '
f'mid-price: {snapshots_abm[0]["mid"]:.2f} → {snapshots_abm[-1]["mid"]:.2f}')