This notebook accompanies the Execution Tactics chapter. It implements:
Single-venue optimal tactic — HJB backward induction and aggressiveness matrix.
Smart order routing (SOR) — multi-venue benefit via aggregated fill rates.
Reinforcement learning — tabular Q-learning for adaptive execution.
Note: The fill probability model is developed in full in the
lob_modelsnotebook (Section 3). Here we use it directly as a building block.
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
np.random.seed(42)
FIGURES_DIR = '../markdown/figures'
plt.rcParams.update({'font.size': 11, 'axes.titlesize': 12})1 Single-venue optimal tactic¶
1.1 Hamilton–Jacobi–Bellman equation under CARA utility¶
The tactic maximises CARA utility of total proceeds :
where is the risk aversion coefficient (the aggressiveness control) and is the calibrated terminal liquidation cost per unit (a market property, not a free parameter).
Working with the certainty equivalent and disutility , the Bellman equation yields the optimal placement depth:
The risk-comfort depth replaces of the risk-neutral case, and is strictly decreasing in : higher risk aversion → shallower placement → faster execution. As it recovers .
Substituting back gives the reduced ODE in remaining time :
The structure is identical to the risk-neutral ODE (same backward induction); only the driving constant (which reduces to as ) changes.
Analytical solution for (with , since ):
1.2 Backward induction¶
We integrate the reduced ODE forward in from the terminal condition .
def compute_hjb_cara(A_list, k_list, gamma, b, T, Q_max=20, N_tau=600):
"""
CARA-utility HJB backward induction.
Parameters
----------
gamma : float
Risk aversion coefficient — controls aggressiveness (tactic parameter).
b : float
Terminal liquidation cost per unit — calibrated from market data, not a free parameter.
Returns
-------
H_arr : (N_tau+1, Q_max+1) certainty-equivalent value function; index 0 = terminal
delta_arr : (N_tau+1, Q_max+1, n_venues) optimal depths
tau_grid : (N_tau+1,)
"""
dtau = T / N_tau
H = np.array([-b * q for q in range(Q_max + 1)], dtype=float)
H_arr = np.zeros((N_tau + 1, Q_max + 1))
delta_arr = np.zeros((N_tau + 1, Q_max + 1, len(A_list)))
H_arr[0] = H.copy()
# Per-venue CARA constants: C_gamma = A/(k+gamma) * (k/(k+gamma))^(k/gamma)
C_gamma = [A / (k + gamma) * (k / (k + gamma)) ** (k / gamma)
for A, k in zip(A_list, k_list)]
# Risk-comfort depth per venue: (1/gamma) * log(1 + gamma/k)
comfort = [(1.0 / gamma) * np.log(1.0 + gamma / k) for k in k_list]
for step in range(N_tau):
H_new = np.zeros(Q_max + 1)
H_new[0] = 0.0
for q in range(1, Q_max + 1):
DH = np.clip(H[q] - H[q - 1], -30.0, 5.0)
rhs = sum(Cg * np.exp(-k * DH) for Cg, k in zip(C_gamma, k_list))
H_new[q] = H[q] + dtau * rhs
for v, (k_v, cd) in enumerate(zip(k_list, comfort)):
delta_arr[step + 1, q, v] = cd + DH
H = H_new.copy()
H_arr[step + 1] = H.copy()
return H_arr, delta_arr, np.linspace(0, T, N_tau + 1)
# ── Aggressiveness matrices for three risk-aversion levels ────────────────
A_mat, k_mat, b_mat, T_mat, Q_max = 1.0, 1.0, 1.0, 10.0, 20
gamma_vals = [0.2, 2.0, 10.0]
titles = [
r'Low risk aversion ($\gamma = 0.2$)',
r'Moderate ($\gamma = 2$)',
r'High risk aversion ($\gamma = 10$)',
]
fig, axes = plt.subplots(1, 3, figsize=(16, 5.5))
for ax, gamma, title in zip(axes, gamma_vals, titles):
H_arr, d_arr, _ = compute_hjb_cara(
[A_mat], [k_mat], gamma, b_mat, T_mat, Q_max=Q_max, N_tau=600)
D_plot = np.flip(d_arr[1:, 1:, 0].T, axis=1)
vmin, vmax = -0.5, min(float(D_plot.max()), 7.0)
im = ax.imshow(D_plot, aspect='auto', origin='lower',
extent=[0, 1, 1/Q_max, 1], cmap='RdYlGn_r', vmin=vmin, vmax=vmax)
t_r = np.linspace(0, 1, D_plot.shape[1])
q_r = np.linspace(1/Q_max, 1, Q_max)
Tg, Qg = np.meshgrid(t_r, q_r)
try:
ax.contour(Tg, Qg, D_plot, levels=[0.0], colors='black', linewidths=1.8, linestyles='--')
except Exception:
pass
ax.set_xlabel('Time elapsed $t/T$')
if ax is axes[0]:
ax.set_ylabel('Remaining fraction $q/Q$')
ax.set_title(title)
plt.colorbar(im, ax=ax, label='$\\delta^*$ (ticks)')
fig.suptitle(
r'CARA-utility aggressiveness matrices ($b = 1$ tick, calibrated)'
'\n(dashed = market-order boundary $\\delta^* = 0$)',
fontsize=12, y=1.02
)
plt.tight_layout()
fig.savefig(f'{FIGURES_DIR}/tact_aggressiveness_matrix.png', dpi=150, bbox_inches='tight')
plt.show()
print('Risk-comfort depths:')
for gamma in gamma_vals:
cd = (1.0 / gamma) * np.log(1.0 + gamma / k_mat)
print(f' gamma={gamma:4.1f}: risk-comfort depth = {cd:.3f} (vs 1/k = {1/k_mat:.3f})')
print('\nAggressiveness matrix at start (tau=T, gamma=2.0):')
_, d_arr2, _ = compute_hjb_cara([A_mat], [k_mat], 2.0, b_mat, T_mat, Q_max=Q_max, N_tau=600)
for q in [1, 5, 10, 15, 20]:
print(f' q={q:2d}: delta* = {d_arr2[-1, q, 0]:.3f}')
Risk-comfort depths:
gamma= 0.2: risk-comfort depth = 0.912 (vs 1/k = 1.000)
gamma= 2.0: risk-comfort depth = 0.549 (vs 1/k = 1.000)
gamma=10.0: risk-comfort depth = 0.240 (vs 1/k = 1.000)
Aggressiveness matrix at start (tau=T, gamma=2.0):
q= 1: delta* = 1.379
q= 5: delta* = -0.090
q=10: delta* = -0.428
q=15: delta* = -0.450
q=20: delta* = -0.451
1.3 Simulation of the single-venue tactic¶
We simulate executions by reading from the look-up table at each step, drawing Poisson inter-fill times, and applying the calibrated terminal liquidation cost to any unfilled units. The risk aversion controls how aggressively the tactic places orders throughout the window; is fixed from market calibration.
def simulate_tactic(H_arr, delta_arr_v0, A, k, b, Q, T, n_sims=5000, seed=0):
"""
Single-venue CARA tactic simulation.
delta* < 0 triggers a market order (zero price improvement).
Residual at T costs b per unit (calibrated market cost).
"""
N_tau = H_arr.shape[0] - 1
rng = np.random.default_rng(seed)
proceeds = np.zeros(n_sims)
for s in range(n_sims):
q, tau, total = Q, T, 0.0
while q > 0 and tau > 1e-9:
tau_idx = int(np.clip(tau / T * N_tau, 0, N_tau))
q_idx = min(q, delta_arr_v0.shape[1] - 1)
d = float(delta_arr_v0[tau_idx, q_idx, 0])
if d < 0:
q -= 1 # market order at mid, zero price improvement
continue
lam = A * np.exp(-k * d)
wait = rng.exponential(1.0 / lam) if lam > 1e-12 else np.inf
if wait < tau:
tau -= wait; total += d; q -= 1
else:
tau = 0
if q > 0:
total -= b * q # b is the calibrated terminal cost, not a tactic parameter
proceeds[s] = total
return proceeds
# Compare three risk-aversion levels; b=1.0 is the calibrated market parameter
A_sim, k_sim, b_mkt, T_sim = 1.0, 1.0, 1.0, 10.0
fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))
for ax, (gamma, Q_test) in zip(axes, [(0.2, 10), (2.0, 10), (10.0, 10)]):
H_sim, d_sim, _ = compute_hjb_cara([A_sim], [k_sim], gamma, b_mkt, T_sim,
Q_max=20, N_tau=600)
proc = simulate_tactic(H_sim, d_sim, A_sim, k_sim, b_mkt, Q_test, T_sim, n_sims=5000)
ax.hist(proc, bins=40, color='#1565C0', alpha=0.7, density=True)
ax.axvline(proc.mean(), color='#C62828', lw=2, linestyle='--',
label=f'mean = {proc.mean():.2f}')
ax.set_title(f'$\\gamma = {gamma}$, $Q = {Q_test}$')
ax.set_xlabel('Total proceeds (ticks)')
ax.set_ylabel('Density')
ax.legend(); ax.grid(True, alpha=0.3)
fig.suptitle('Proceeds distribution — single-venue CARA tactic ($b = 1$, calibrated)', fontsize=12)
plt.tight_layout()
plt.show()
2 Smart Order Routing (SOR)¶
When orders can be posted simultaneously on venues, the total fill rate is
Each venue contributes independently, so the combined HJB simply sums per-venue contributions:
and the optimal depth for venue is (same for all venues).
Key SOR insight¶
Adding a second venue increases , reducing expected waiting time for the next fill. This lowers the probability of hitting the terminal deadline with residual inventory, improving execution quality even when the second venue is illiquid.
def simulate_sor_fixed(A_list, k_list, delta_fixed, b, Q, T, n_sims=10_000, seed=0):
"""
Post at a fixed depth delta_fixed on every venue simultaneously.
Total fill rate = sum of individual venue rates.
Terminal residual at T costs b per unit.
"""
lams = [A * np.exp(-k * delta_fixed) for A, k in zip(A_list, k_list)]
total_lam = sum(lams)
rng = np.random.default_rng(seed)
proceeds = np.zeros(n_sims)
for s in range(n_sims):
q, tau, total = Q, T, 0.0
while q > 0 and tau > 1e-9:
wait = rng.exponential(1.0 / total_lam) if total_lam > 1e-12 else np.inf
if wait < tau:
tau -= wait; total += delta_fixed; q -= 1
else:
tau = 0
if q > 0:
total -= b * q
proceeds[s] = total
return proceeds
# Parameters
A1, k1 = 2.0, 0.5 # liquid venue
A2, k2 = 0.5, 1.5 # illiquid venue
b_sor, Q_sor, T_sor, delta_sor = 2.0, 15, 10.0, 1.0
proc_v1 = simulate_sor_fixed([A1], [k1], delta_sor, b_sor, Q_sor, T_sor, seed=1)
proc_v2 = simulate_sor_fixed([A2], [k2], delta_sor, b_sor, Q_sor, T_sor, seed=2)
proc_sor = simulate_sor_fixed([A1, A2], [k1, k2], delta_sor, b_sor, Q_sor, T_sor, seed=3)
print('Fixed-depth strategy (delta=1 tick), Q=15, T=10 min')
print(f' Venue 1 only (A={A1}, k={k1}): mean={proc_v1.mean():.1f}, std={proc_v1.std():.1f}')
print(f' Venue 2 only (A={A2}, k={k2}): mean={proc_v2.mean():.1f}, std={proc_v2.std():.1f}')
lam_tot = A1*np.exp(-k1*delta_sor) + A2*np.exp(-k2*delta_sor)
print(f' Two-venue SOR (lambda_tot={lam_tot:.2f}/min): mean={proc_sor.mean():.1f}, std={proc_sor.std():.1f}')Fixed-depth strategy (delta=1 tick), Q=15, T=10 min
Venue 1 only (A=2.0, k=0.5): mean=5.2, std=8.4
Venue 2 only (A=0.5, k=1.5): mean=-26.6, std=3.2
Two-venue SOR (lambda_tot=1.32/min): mean=7.4, std=7.8
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
# Left: fill rate curves
d_plt = np.linspace(0, 5, 200)
axes[0].plot(d_plt, A1 * np.exp(-k1 * d_plt), lw=2.5, color='#1565C0',
label=f'Venue 1 (liquid, $A={A1}$, $k={k1}$)')
axes[0].plot(d_plt, A2 * np.exp(-k2 * d_plt), lw=2.5, color='#C62828',
label=f'Venue 2 (illiquid, $A={A2}$, $k={k2}$)')
axes[0].axvline(delta_sor, color='gray', lw=1.5, linestyle=':',
label=f'$\\delta = {delta_sor}$ (strategy)')
axes[0].set_xlabel('Depth $\\delta$ (ticks)')
axes[0].set_ylabel('Fill rate (fills/min)')
axes[0].set_title('Venue fill rates')
axes[0].legend(); axes[0].grid(True, alpha=0.3)
# Right: P&L distributions
lo = min(proc_v1.min(), proc_v2.min(), proc_sor.min()) - 1
hi = max(proc_v1.max(), proc_v2.max(), proc_sor.max()) + 1
bins = np.linspace(lo, hi, 60)
for proc, lbl, col in [
(proc_v1, 'Venue 1 only', '#1565C0'),
(proc_v2, 'Venue 2 only', '#C62828'),
(proc_sor, 'Two-venue SOR', '#2E7D32'),
]:
axes[1].hist(proc, bins=bins, alpha=0.40, color=col, label=lbl, density=True)
axes[1].axvline(proc.mean(), color=col, lw=2.0, linestyle='--')
axes[1].set_xlabel(f'Total proceeds (ticks, $Q={Q_sor}$, $\\delta={delta_sor}$)')
axes[1].set_ylabel('Density')
axes[1].set_title('Execution P&L distributions (10,000 simulations)')
axes[1].legend(); axes[1].grid(True, alpha=0.3)
plt.tight_layout()
fig.savefig(f'{FIGURES_DIR}/tact_sor.png', dpi=150, bbox_inches='tight')
plt.show()3 Reinforcement learning for execution tactics¶
Execution tactics can be framed as a Markov Decision Process:
| MDP element | Execution interpretation |
|---|---|
| State | Elapsed time, remaining inventory |
| Action | Placement depth (or market order) |
| Reward | Price improvement per fill; at terminal time |
| Transition | Poisson fill event; time step |
Q-learning maintains a table and updates on each transition:
The -greedy policy starts fully exploratory and decays toward the greedy policy, recovering the HJB-optimal solution as episode count grows.
# MDP parameters — b_mkt is the calibrated terminal cost; gamma_rl controls RL aggressiveness
# The RL reward uses the linear (risk-neutral) objective, approximating gamma -> 0
A_rl, k_rl, b_mkt_rl, T_rl, Q_rl = 1.0, 1.0, 1.0, 10.0, 10
DT_RL = T_rl / 20 # 0.5-min time steps per episode tick
ACTIONS = np.array([-1, 0, 1, 2, 3, 4]) # -1 = market order (proceeds = 0)
N_T_BINS = 10
N_Q_BINS = Q_rl + 1
def rl_step(tau, q, action_idx, rng):
d = int(ACTIONS[action_idx])
if d < 0 or tau <= DT_RL + 1e-9:
new_q, new_tau = q - 1, max(tau - DT_RL, 0.0)
done = new_tau <= 0 or new_q == 0
reward = -b_mkt_rl * new_q if (done and new_q > 0) else 0.0
return new_tau, new_q, reward, done
lam = A_rl * np.exp(-k_rl * d)
filled = rng.random() < (1 - np.exp(-lam * DT_RL))
new_q = (q - 1) if filled else q
new_tau = tau - DT_RL
done = new_tau <= 0 or new_q == 0
reward = (float(d) if filled else 0.0) + (-b_mkt_rl * new_q if (done and new_q > 0) else 0.0)
return new_tau, new_q, reward, done
# Q-learning
Q_table = np.zeros((N_T_BINS, N_Q_BINS, len(ACTIONS)))
alpha, gamma_disc = 0.12, 0.99 # gamma_disc is the RL discount factor (not risk aversion)
eps_start, eps_end, N_EP = 1.0, 0.05, 6000
eps_decay = (eps_start - eps_end) / N_EP
rng_rl = np.random.default_rng(42)
ep_rewards = []
for ep in range(N_EP):
eps = max(eps_end, eps_start - ep * eps_decay)
tau, q, total_r = T_rl, Q_rl, 0.0
for _ in range(120):
if q == 0 or tau <= 0: break
tb = min(int(tau / T_rl * N_T_BINS), N_T_BINS - 1)
qb = min(q, N_Q_BINS - 1)
a = (rng_rl.integers(len(ACTIONS)) if rng_rl.random() < eps
else int(np.argmax(Q_table[tb, qb])))
ntau, nq, r, done = rl_step(tau, q, a, rng_rl)
total_r += r
ntb = min(int(ntau / T_rl * N_T_BINS), N_T_BINS - 1)
nqb = min(nq, N_Q_BINS - 1)
target = r if done else r + gamma_disc * np.max(Q_table[ntb, nqb])
Q_table[tb, qb, a] += alpha * (target - Q_table[tb, qb, a])
tau, q = ntau, nq
if done: break
ep_rewards.append(total_r)
print(f'Training done. Final 500-ep average reward: {np.mean(ep_rewards[-500:]):.2f}')
print(f'(RL uses linear reward — risk-neutral limit; b={b_mkt_rl} is calibrated market cost)')
Training done. Final 500-ep average reward: 1.50
(RL uses linear reward — risk-neutral limit; b=1.0 is calibrated market cost)
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
# Left: learning curve
window = 300
smooth = np.convolve(ep_rewards, np.ones(window)/window, mode='valid')
ep_ax = np.arange(len(smooth)) + window // 2
axes[0].scatter(np.arange(N_EP)[::5], np.array(ep_rewards)[::5],
s=3, alpha=0.15, color='#90CAF9')
axes[0].plot(ep_ax, smooth, lw=2, color='#1565C0', label=f'{window}-ep rolling avg')
axes[0].axhline(np.mean(ep_rewards[-500:]), color='#C62828', lw=1.5,
linestyle='--', label=f'Converged ≈ {np.mean(ep_rewards[-500:]):.1f}')
axes[0].set_xlabel('Episode')
axes[0].set_ylabel('Total episodic reward')
axes[0].set_title('Q-learning training curve (risk-neutral reward)')
axes[0].legend(); axes[0].grid(True, alpha=0.3)
# Right: greedy policy — best action mapped to delta value
best_a = np.argmax(Q_table, axis=2)
best_d = ACTIONS[best_a]
plot_d = np.flip(best_d[1:, 1:].T, axis=1)
im = axes[1].imshow(plot_d, aspect='auto', origin='lower',
extent=[0, 1, 1/Q_rl, 1],
cmap='RdYlGn_r', vmin=-1.5, vmax=4.5)
plt.colorbar(im, ax=axes[1], label='Greedy action $\\delta^*$ (ticks)')
axes[1].set_xlabel('Time elapsed $t/T$')
axes[1].set_ylabel('Remaining fraction $q/Q$')
axes[1].set_title(f'Learned greedy policy ($b={b_mkt_rl}$ calibrated, $A=1$, $k=1$)')
plt.tight_layout()
fig.savefig(f'{FIGURES_DIR}/tact_rl_convergence.png', dpi=150, bbox_inches='tight')
plt.show()
print('\nGreedy RL policy at start (t=0) vs CARA-HJB (gamma=2.0):')
tb0 = N_T_BINS - 1
_, d_hjb, _ = compute_hjb_cara([A_rl], [k_rl], 2.0, b_mkt_rl, T_rl, Q_max=Q_rl, N_tau=600)
for q_test in [1, 3, 5, 10]:
qb = min(q_test, N_Q_BINS - 1)
a_rl = int(np.argmax(Q_table[tb0, qb]))
d_cara = d_hjb[-1, q_test, 0]
print(f' q={q_test:2d}: delta_RL={ACTIONS[a_rl]:2d} delta_CARA={d_cara:.2f}')
Greedy RL policy at start (t=0) vs CARA-HJB (gamma=2.0):
q= 1: delta_RL=-1 delta_CARA=1.38
q= 3: delta_RL=-1 delta_CARA=0.34
q= 5: delta_RL=-1 delta_CARA=-0.09
q=10: delta_RL=-1 delta_CARA=-0.43
Exercises¶
Exercise 1 (Analytical ODE for , CARA utility): For a single venue, the certainty-equivalent ODE for reduces to where .
(a) Verify that solves this ODE with initial condition .
(b) Show that as , and the solution reduces to the risk-neutral form .
(c) Verify numerically against compute_hjb_cara for , , , .
Exercise 2 (Risk-comfort depth): For parameters , , compute the risk-comfort depth for and verify the limiting cases and . At what does the risk-comfort depth fall below half of the risk-neutral benchmark ?
Exercise 3 (Multi-venue independence): In the CARA SOR model with venues, show that the optimal depth for venue is , where is venue-independent. How does the cross-venue structure change compared with the risk-neutral case?
Exercise 4 (Dark pool routing): Extend simulate_tactic to include a dark pool (fills at
mid-price, zero price improvement, with fill probability per time step) and show the
P&L improvement for a large order units.
Exercise 5 (DQN extension): Replace the Q-table in the RL section with a two-layer neural network using PyTorch. Train on the same MDP and compare the greedy policy to the CARA-HJB solution for .