This notebook accompanies the Liquidity Modelling chapter. It provides worked simulations for the main liquidity indicators and aggregation methods discussed in the text.
Sections
Activity-based indicators
Amihud illiquidity ratio
Roll estimator
Price dispersion
Aggregating into a liquidity score
Inventory rotation time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
rng = np.random.default_rng(42)
plt.rcParams.update({
'figure.dpi': 120,
'axes.spines.top': False,
'axes.spines.right': False,
'font.size': 11,
})1. Activity-based indicators¶
We simulate a cross-section of bonds with different liquidity levels and compute the canonical activity-based indicators: number of trades per day, percentage of days traded, turnover, and average trade size.
def simulate_bond_activity(n_bonds: int = 50, n_days: int = 30, seed: int = 42) -> pd.DataFrame:
"""Simulate daily trade records for a cross-section of bonds.
Bonds are parameterised by a latent liquidity level l ~ Uniform(0.05, 1).
More liquid bonds trade more often and in smaller average sizes.
"""
rng = np.random.default_rng(seed)
liquidity = rng.uniform(0.05, 1.0, size=n_bonds) # true latent liquidity
outstanding = rng.uniform(200, 5000, size=n_bonds) # outstanding notional (€M)
rows = []
for i, (liq, out) in enumerate(zip(liquidity, outstanding)):
# Expected trades per day driven by liquidity
lambda_trades = liq * 5 + 0.1
for d in range(n_days):
n_trades = rng.poisson(lambda_trades)
if n_trades == 0:
continue
# Average size inversely related to liquidity
avg_size = rng.lognormal(mean=np.log(10 / liq), sigma=0.5, size=n_trades)
rows.append({
'bond': i,
'day': d,
'n_trades': n_trades,
'volume': avg_size.sum(),
'outstanding': out,
'true_liquidity': liq,
})
return pd.DataFrame(rows)
trades = simulate_bond_activity(n_bonds=60, n_days=30)
# Aggregate to bond level
stats = (
trades.groupby('bond')
.agg(
avg_trades_per_day=('n_trades', 'mean'),
pct_days_traded=('day', lambda x: x.nunique() / 30),
avg_trade_size=('volume', lambda x: x.sum() / trades.loc[x.index, 'n_trades'].sum()),
total_volume=('volume', 'sum'),
outstanding=('outstanding', 'first'),
true_liquidity=('true_liquidity', 'first'),
)
)
stats['turnover'] = stats['total_volume'] / (stats['outstanding'] * 30)
print(stats.sort_values('true_liquidity', ascending=False).head(10).to_string(float_format='{:.3f}'.format))fig, axes = plt.subplots(1, 3, figsize=(13, 4))
metrics = [
('avg_trades_per_day', 'Avg trades per day'),
('pct_days_traded', 'Days traded (%)'),
('turnover', 'Turnover'),
]
for ax, (col, label) in zip(axes, metrics):
ax.scatter(stats['true_liquidity'], stats[col], alpha=0.7, s=30, color='steelblue')
ax.set_xlabel('True latent liquidity')
ax.set_ylabel(label)
corr = stats[['true_liquidity', col]].corr().iloc[0, 1]
ax.set_title(f'{label}\n(r = {corr:.2f})')
fig.suptitle('Activity-based indicators vs. true liquidity', y=1.02, fontsize=13)
plt.tight_layout()
plt.show()2. Amihud illiquidity ratio¶
The Amihud ratio is defined as
We simulate trade-by-trade price paths for instruments with different depths and verify that the Amihud ratio correctly orders them.
def simulate_amihud(n_trades: int, depth: float, sigma: float = 0.001,
seed: int = None) -> float:
"""Simulate trades and return the Amihud ratio.
Price impact per unit volume: impact = volume / depth.
Mid-price also has Gaussian noise with std sigma (information shocks).
"""
rng = np.random.default_rng(seed)
# Trade volumes ~ LogNormal
volumes = rng.lognormal(mean=np.log(5), sigma=0.8, size=n_trades)
# Trade direction
directions = rng.choice([-1, 1], size=n_trades)
# Price impact: linear in volume for simplicity
impact = directions * volumes / depth
# Information shocks
info = rng.normal(0, sigma, size=n_trades)
# Observed price changes
price = 100.0
amihud_vals = []
for v, imp, inf_shock in zip(volumes, impact, info):
dp = imp + inf_shock
ret = abs(dp / price)
amihud_vals.append(ret / v)
price += dp
return np.mean(amihud_vals)
depths = np.logspace(1, 4, 40) # depth from 10 to 10,000
amihud_means = []
amihud_stds = []
n_sims = 50
for d in depths:
vals = [simulate_amihud(200, d, seed=i) for i in range(n_sims)]
amihud_means.append(np.mean(vals))
amihud_stds.append(np.std(vals))
amihud_means = np.array(amihud_means)
amihud_stds = np.array(amihud_stds)
fig, ax = plt.subplots(figsize=(7, 4))
ax.loglog(depths, amihud_means, color='steelblue', lw=2, label='Mean Amihud ratio')
ax.fill_between(depths, amihud_means - amihud_stds, amihud_means + amihud_stds,
alpha=0.25, color='steelblue', label='±1 std')
ax.set_xlabel('Market depth (units of volume)')
ax.set_ylabel('Amihud illiquidity ratio')
ax.set_title('Amihud ratio decreases with market depth')
ax.legend()
plt.tight_layout()
plt.show()Effect of trade size distribution¶
In practice, trade sizes follow a heavy-tailed distribution. We compare the Amihud estimator under log-normal and power-law size distributions.
def simulate_amihud_dist(n_trades: int, depth: float, dist: str,
sigma: float = 0.001, seed: int = None) -> float:
rng = np.random.default_rng(seed)
if dist == 'lognormal':
volumes = rng.lognormal(mean=np.log(5), sigma=0.8, size=n_trades)
else: # power law via inverse transform, exponent 2.5, min size 1
u = rng.uniform(size=n_trades)
volumes = (1 - u) ** (-1 / 1.5)
volumes = np.clip(volumes, 1, 500)
directions = rng.choice([-1, 1], size=n_trades)
price = 100.0
vals = []
for v, d_, info in zip(volumes, directions, rng.normal(0, sigma, n_trades)):
dp = d_ * v / depth + info
ret = abs(dp / price)
vals.append(ret / v)
price += dp
return np.mean(vals)
n_sims = 200
depth = 500
lognorm_vals = [simulate_amihud_dist(200, depth, 'lognormal', seed=i) for i in range(n_sims)]
powerlaw_vals = [simulate_amihud_dist(200, depth, 'powerlaw', seed=i) for i in range(n_sims)]
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
for ax, vals, label in zip(axes, [lognorm_vals, powerlaw_vals], ['Log-normal sizes', 'Power-law sizes']):
ax.hist(vals, bins=30, color='steelblue', edgecolor='white', alpha=0.85)
ax.axvline(np.mean(vals), color='crimson', lw=1.8, label=f'Mean = {np.mean(vals):.2e}')
ax.axvline(np.median(vals), color='darkorange', lw=1.8, linestyle='--',
label=f'Median = {np.median(vals):.2e}')
ax.set_title(f'Amihud distribution — {label}')
ax.set_xlabel('Amihud ratio')
ax.set_ylabel('Count')
ax.legend(fontsize=9)
plt.suptitle('Amihud ratio under different trade size distributions (same depth)', fontsize=12, y=1.02)
plt.tight_layout()
plt.show()3. Roll estimator¶
The Roll estimator recovers the half-spread from the serial correlation of price changes:
We simulate a mid-price random walk with bid-ask bounce and verify that the estimator recovers the true half-spread. We then examine how drift distorts it.
def simulate_roll(n: int, half_spread: float, sigma_mid: float = 0.02,
drift: float = 0.0, seed: int = None):
"""Simulate transaction prices under a random-walk mid + bid-ask bounce."""
rng = np.random.default_rng(seed)
mid = np.zeros(n)
mid[0] = 100.0
for t in range(1, n):
mid[t] = mid[t - 1] + drift + rng.normal(0, sigma_mid)
direction = rng.choice([-1, 1], size=n)
prices = mid + half_spread * direction
return prices, mid
def roll_estimate(prices: np.ndarray) -> float:
dp = np.diff(prices)
cov = np.cov(dp[1:], dp[:-1])[0, 1]
if cov >= 0:
return np.nan # estimator undefined
return np.sqrt(-cov)
# --- Varying true spread, zero drift ---
true_spreads = np.linspace(0.005, 0.20, 25)
n_sims = 100
estimated = []
for s in true_spreads:
ests = [roll_estimate(simulate_roll(500, s, seed=i)[0]) for i in range(n_sims)]
estimated.append(np.nanmean(ests))
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
ax = axes[0]
ax.plot(true_spreads, true_spreads, 'k--', lw=1, label='45° line (true value)')
ax.plot(true_spreads, estimated, 'o-', color='steelblue', ms=5, lw=1.5,
label='Roll estimate (mean)')
ax.set_xlabel('True half-spread')
ax.set_ylabel('Estimated half-spread')
ax.set_title('Roll estimator vs. true spread (zero drift)')
ax.legend()
# --- Effect of drift ---
half_spread = 0.05
drifts = np.linspace(0.0, 0.03, 20)
roll_by_drift = []
for dr in drifts:
ests = [roll_estimate(simulate_roll(500, half_spread, drift=dr, seed=i)[0])
for i in range(n_sims)]
roll_by_drift.append(np.nanmean(ests))
ax = axes[1]
ax.axhline(half_spread, color='k', lw=1, linestyle='--', label='True half-spread')
ax.plot(drifts / half_spread, roll_by_drift, 'o-', color='crimson', ms=5, lw=1.5,
label='Roll estimate')
ax.set_xlabel('Drift / half-spread ratio')
ax.set_ylabel('Estimated half-spread')
ax.set_title('Roll estimator degrades with mid-price drift')
ax.legend()
plt.tight_layout()
plt.show()4. Price dispersion¶
Price dispersion measures how much trade prices deviate from a mid-price benchmark. We simulate bond trades for instruments with different effective spreads and compare the price dispersion estimator with the true spread.
def price_dispersion(trade_prices: np.ndarray, mid_prices: np.ndarray,
volumes: np.ndarray) -> float:
"""Compute volume-weighted price dispersion."""
rel_dev = np.abs((trade_prices - mid_prices) / mid_prices)
return 2 * np.sum(rel_dev * volumes) / np.sum(volumes)
def simulate_bond_trades(n: int, half_spread: float, sigma_mid: float = 0.05,
seed: int = None):
rng = np.random.default_rng(seed)
mid = np.cumsum(rng.normal(0, sigma_mid, n)) + 100.0
direction = rng.choice([-1, 1], size=n)
# Add noise around the quoted spread (simulates partial price improvement)
actual_spread = half_spread + rng.normal(0, half_spread * 0.15, n)
trade_prices = mid + actual_spread * direction
volumes = rng.lognormal(2, 1, size=n)
return trade_prices, mid, volumes
spreads = np.linspace(0.01, 0.40, 30)
n_sims = 80
pd_estimates = []
for s in spreads:
vals = []
for i in range(n_sims):
tp, mp, vol = simulate_bond_trades(300, s, seed=i)
vals.append(price_dispersion(tp, mp, vol))
pd_estimates.append(np.mean(vals))
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(spreads, spreads, 'k--', lw=1, label='True round-trip spread')
ax.plot(spreads, pd_estimates, 'o-', color='steelblue', ms=4, lw=1.5,
label='Price dispersion estimate')
ax.set_xlabel('True half-spread')
ax.set_ylabel('Estimated spread')
ax.set_title('Price dispersion vs. true spread')
ax.legend()
plt.tight_layout()
plt.show()5. Aggregating into a liquidity score¶
We simulate a panel of bonds with multiple noisy liquidity indicators and compare three aggregation methods: ranking, PCA, and a simple Lasso regression predicting future turnover.
from sklearn.linear_model import Lasso
from sklearn.pipeline import Pipeline
def simulate_liquidity_panel(n_bonds: int = 100, noise: float = 0.4,
seed: int = 0) -> pd.DataFrame:
"""Generate a cross-section of bonds with noisy liquidity metrics."""
rng = np.random.default_rng(seed)
liq = rng.uniform(0, 1, n_bonds) # true latent liquidity
def noisy(signal, sign=1):
return sign * signal + rng.normal(0, noise, n_bonds)
df = pd.DataFrame({
'true_liq': liq,
'num_trades': noisy(liq * 10, sign=1), # higher = more liquid
'pct_days_traded': noisy(liq, sign=1), # higher = more liquid
'turnover': noisy(liq * 0.5, sign=1), # higher = more liquid
'amihud': noisy(-liq * 2, sign=1), # higher = less liquid
'bid_ask_spread': noisy(-liq, sign=1), # higher = less liquid
'price_dispersion': noisy(-liq, sign=1), # higher = less liquid
'avg_trade_size': noisy(-liq * 5 + 5, sign=1), # higher = less liquid
# forward turnover (target for ML approach)
'future_turnover': liq * 0.5 + rng.normal(0, 0.1, n_bonds),
})
return df
df = simulate_liquidity_panel(n_bonds=150, noise=0.3, seed=7)
# ------ 1. Ranking aggregation ------
# Metrics where higher = less liquid: invert
higher_is_illiquid = ['amihud', 'bid_ask_spread', 'price_dispersion', 'avg_trade_size']
higher_is_liquid = ['num_trades', 'pct_days_traded', 'turnover']
metric_cols = higher_is_liquid + higher_is_illiquid
ranks = pd.DataFrame(index=df.index)
for col in higher_is_liquid:
ranks[col] = df[col].rank(ascending=True)
for col in higher_is_illiquid:
ranks[col] = df[col].rank(ascending=False) # invert: low amihud = high rank
avg_rank = ranks.mean(axis=1)
df['score_rank'] = (avg_rank - avg_rank.min()) / (avg_rank.max() - avg_rank.min())
# ------ 2. PCA aggregation ------
# For PCA, sign-adjust so that higher always means more liquid
X_signed = df[metric_cols].copy()
for col in higher_is_illiquid:
X_signed[col] = -X_signed[col]
scaler = StandardScaler()
X_std = scaler.fit_transform(X_signed)
pca = PCA(n_components=3)
scores_pca = pca.fit_transform(X_std)
# First PC should align positively with liquidity; flip if needed
if np.corrcoef(scores_pca[:, 0], df['true_liq'])[0, 1] < 0:
scores_pca[:, 0] *= -1
df['score_pca'] = scores_pca[:, 0]
# Normalize to [0, 1] for comparison
df['score_pca'] = (df['score_pca'] - df['score_pca'].min()) / \
(df['score_pca'].max() - df['score_pca'].min())
# ------ 3. Lasso (predict future turnover) ------
pipe = Pipeline([
('scaler', StandardScaler()),
('lasso', Lasso(alpha=0.02, max_iter=5000)),
])
pipe.fit(df[metric_cols], df['future_turnover'])
df['score_lasso'] = pipe.predict(df[metric_cols])
df['score_lasso'] = (df['score_lasso'] - df['score_lasso'].min()) / \
(df['score_lasso'].max() - df['score_lasso'].min())
print('Correlations with true latent liquidity:')
for col in ['score_rank', 'score_pca', 'score_lasso']:
r = df[['true_liq', col]].corr().iloc[0, 1]
print(f' {col:20s}: {r:.3f}')
print()
print('PCA explained variance (first 3 components):',
[f'{v:.1%}' for v in pca.explained_variance_ratio_])fig, axes = plt.subplots(1, 3, figsize=(13, 4))
methods = [
('score_rank', 'Ranking aggregation'),
('score_pca', 'PCA (PC1)'),
('score_lasso', 'Lasso regression'),
]
for ax, (col, title) in zip(axes, methods):
ax.scatter(df['true_liq'], df[col], alpha=0.5, s=25, color='steelblue')
r = df[['true_liq', col]].corr().iloc[0, 1]
ax.set_title(f'{title}\n(r = {r:.3f})')
ax.set_xlabel('True latent liquidity')
ax.set_ylabel('Composite score (normalised)')
fig.suptitle('Comparison of liquidity score aggregation methods', fontsize=13, y=1.02)
plt.tight_layout()
plt.show()# PCA loadings
loadings = pd.Series(
pca.components_[0],
index=[c + (' (inv)' if c in higher_is_illiquid else '') for c in metric_cols]
).sort_values()
fig, ax = plt.subplots(figsize=(8, 4))
colors = ['crimson' if v < 0 else 'steelblue' for v in loadings]
ax.barh(loadings.index, loadings.values, color=colors, edgecolor='white')
ax.axvline(0, color='black', lw=0.8)
ax.set_title('PCA first component loadings\n(positive = aligns with liquidity)')
ax.set_xlabel('Loading')
plt.tight_layout()
plt.show()6. Inventory rotation time¶
Within the Avellaneda-Stoikov framework, the rate at which a market maker’s inventory reverts to zero is driven by the net flow of client requests. We simulate the inventory path under optimal quoting and study how the rotation time depends on the key parameters: arrival rate , average size , volatility , and price elasticity .
The optimal half-spread from the asymptotic Avellaneda-Stoikov solution is
and the win probability for a given inventory level adjusts around this spread asymmetrically to push inventory back toward zero.
def simulate_inventory_rotation(
q0: float = 100.0,
A: float = 10.0, # RfQs per day
x_bar: float = 5.0, # average RfQ size
alpha: float = 2.0, # price elasticity (1/price_unit)
sigma: float = 0.01, # daily volatility
p0: float = 0.5, # win prob at mid
gamma: float = 0.1, # risk aversion
n_days: float = 60.0, # simulation horizon (days)
seed: int = 0,
):
"""Simulate inventory under AS optimal quoting via a daily approximation.
The optimal spread creates an asymmetry: the side of the spread facing
the inventory excess is tightened, increasing hit probability on that side.
"""
rng = np.random.default_rng(seed)
# Optimal base half-spread (asymptotic approximation)
delta_star = 1 / alpha + gamma * sigma ** 2 / (2 * alpha ** 2 * A * p0)
q = q0
history = [q]
dt = 1 / 252 # daily time step (trading year)
steps = int(n_days)
for _ in range(steps):
# Inventory skew: dealer tightens the spread on the side that reduces q
inventory_adjustment = gamma * sigma ** 2 * q / (alpha * A * p0)
# sell-side spread (client buy): tighter if q > 0
delta_sell = max(delta_star - inventory_adjustment, 1e-4)
# buy-side spread (client sell): wider if q > 0
delta_buy = delta_star + inventory_adjustment
# Win probabilities on each side
p_win_sell = p0 * np.exp(-alpha * delta_sell)
p_win_buy = p0 * np.exp(-alpha * delta_buy)
# Expected RfQs today
n_rfq = rng.poisson(A)
for _ in range(n_rfq):
side = rng.choice(['buy', 'sell']) # client request side
size = rng.exponential(x_bar)
if side == 'sell':
# client sells → dealer buys → inventory increases
p_win = p_win_buy
if rng.random() < p_win:
q += size
else:
# client buys → dealer sells → inventory decreases
p_win = p_win_sell
if rng.random() < p_win:
q -= size
history.append(q)
return np.array(history)
# --- Single path ---
path = simulate_inventory_rotation(q0=100, A=10, x_bar=5, seed=42)
threshold = 50.0 # half-inventory level
crossings = np.where(np.abs(path) <= threshold)[0]
t_half = crossings[0] if len(crossings) > 0 else len(path)
fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(path, color='steelblue', lw=1.5)
ax.axhline(threshold, color='crimson', lw=1, linestyle='--', label=f'Half-inventory = {threshold}')
ax.axhline(0, color='black', lw=0.8)
if t_half < len(path):
ax.axvline(t_half, color='darkorange', lw=1.5, linestyle=':',
label=f'Half-life ≈ {t_half} days')
ax.set_xlabel('Days')
ax.set_ylabel('Inventory level')
ax.set_title('Inventory path under optimal Avellaneda-Stoikov quoting')
ax.legend()
plt.tight_layout()
plt.show()# --- Sensitivity analysis: how does T_half depend on each parameter? ---
def expected_half_life(params: dict, n_sims: int = 30, q0: float = 100.0,
n_days: int = 120, seed_offset: int = 0) -> float:
half_lives = []
for s in range(n_sims):
path = simulate_inventory_rotation(q0=q0, n_days=n_days, seed=s + seed_offset,
**params)
cross = np.where(np.abs(path) <= q0 / 2)[0]
half_lives.append(cross[0] if len(cross) > 0 else n_days)
return np.mean(half_lives)
base = dict(A=10.0, x_bar=5.0, alpha=2.0, sigma=0.01, p0=0.5, gamma=0.1)
sweep_params = {
'A (RfQ arrival rate)': ('A', np.linspace(2, 30, 12)),
'x̄ (avg trade size)': ('x_bar', np.linspace(1, 20, 12)),
'σ (volatility)': ('sigma', np.linspace(0.002, 0.05, 12)),
'α (price elasticity)': ('alpha', np.linspace(0.5, 6.0, 12)),
}
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
for ax, (label, (key, values)) in zip(axes.flat, sweep_params.items()):
hl = []
for v in values:
params = {**base, key: v}
hl.append(expected_half_life(params, n_sims=20))
ax.plot(values, hl, 'o-', color='steelblue', ms=5, lw=1.5)
ax.set_xlabel(label)
ax.set_ylabel('Expected half-life (days)')
ax.set_title(f'Inventory half-life vs. {label}')
fig.suptitle('Inventory rotation time — parameter sensitivity', fontsize=13, y=1.01)
plt.tight_layout()
plt.show()# --- On-the-run vs. off-the-run bond comparison ---
scenarios = {
'On-the-run sovereign': dict(A=30, x_bar=10, alpha=3.0, sigma=0.005, p0=0.6, gamma=0.05),
'Off-the-run sovereign': dict(A=8, x_bar=6, alpha=2.0, sigma=0.007, p0=0.45, gamma=0.05),
'Liquid corporate': dict(A=5, x_bar=4, alpha=1.5, sigma=0.012, p0=0.4, gamma=0.1),
'Illiquid corporate': dict(A=1, x_bar=3, alpha=0.8, sigma=0.025, p0=0.3, gamma=0.1),
}
n_sims = 50
n_days = 120
q0 = 100.0
fig, ax = plt.subplots(figsize=(10, 5))
colors = ['steelblue', 'cornflowerblue', 'darkorange', 'crimson']
for (name, params), color in zip(scenarios.items(), colors):
paths = [simulate_inventory_rotation(q0=q0, n_days=n_days, seed=s, **params)
for s in range(n_sims)]
mean_path = np.mean(paths, axis=0)
std_path = np.std(paths, axis=0)
t = np.arange(n_days + 1)
ax.plot(t, mean_path / q0, lw=2, color=color, label=name)
ax.fill_between(t, (mean_path - std_path) / q0, (mean_path + std_path) / q0,
alpha=0.15, color=color)
ax.axhline(0.5, color='black', lw=1, linestyle='--', label='Half-inventory level')
ax.set_xlabel('Days')
ax.set_ylabel('Inventory / initial inventory')
ax.set_title('Inventory rotation across bond types (mean ± 1 std, 50 simulations)')
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()