Execution Algorithms

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

FIGURES_DIR = '../markdown/figures'

np.random.seed(42)

def make_price(t, p0, drift=0.18, noise=0.06, seed=None):
    """Random walk price trajectory."""
    rng = np.random.RandomState(seed)
    n = len(t)
    dt = t[1] - t[0]
    increments = rng.randn(n) * noise * np.sqrt(dt) + drift * dt
    p = np.zeros(n); p[0] = p0
    for i in range(1, n):
        p[i] = p[i - 1] + increments[i]
    return p

t = np.linspace(0, 1, 300)
mid = make_price(t, 0.42, drift=0.28, noise=0.07, seed=10)
ask = mid + 0.07
bid = mid - 0.07

fig, ax = plt.subplots(figsize=(10, 5.5))
ax.set_xlim(-0.1, 1.1)
ax.set_ylim(-0.05, 1.15)
ax.axis('off')

ax.plot(t, ask, color='steelblue', lw=1.4, label='best ask')
ax.plot(t, mid, color='steelblue', lw=1.4, ls='--', label='mid-price')
ax.plot(t, bid, color='steelblue', lw=1.4, ls=':', label='best bid')

# vertical markers for decision, arrival, execution times
td, ta, te = 0.0, 0.38, 0.82
for x, lbl in [(td, 'decision time'), (ta, 'arrival time'), (te, 'execution time')]:
    ax.axvline(x, color='gray', lw=0.8, ls=':')
    ax.text(x, -0.04, lbl, ha='center', va='top', fontsize=8.5, color='#444')

# key price points
pd_val = mid[0]
pa_val = mid[int(ta * 299)]
pe_val = ask[int(te * 299)] + 0.08   # market impact lifts price above ask

ax.plot(td, pd_val, 'o', color='#2ca02c', ms=7, zorder=5)
ax.plot(ta, pa_val, 'o', color='#2ca02c', ms=7, zorder=5)
ax.plot(te, pe_val, 'o', color='#d62728', ms=7, zorder=5)
ax.text(td - 0.015, pd_val, r'$p_d$', ha='right', va='center', fontsize=9)
ax.text(ta - 0.015, pa_val, r'$p_a$', ha='right', va='center', fontsize=9)
ax.text(te + 0.015, pe_val, r'$p_{exec}$', ha='left', va='center', fontsize=9)

def double_arrow(ax, x, y1, y2, label, side='right', color='#444'):
    cx = x + (0.045 if side == 'right' else -0.045)
    ax.annotate('', xy=(cx, y1), xytext=(cx, y2),
                arrowprops=dict(arrowstyle='<->', color=color, lw=1.2))
    ax.text(cx + (0.015 if side == 'right' else -0.015), (y1 + y2) / 2,
            label, ha='left' if side == 'right' else 'right',
            va='center', fontsize=8.5, color=color)

double_arrow(ax, -0.07, pd_val, pa_val, 'delay cost',    color='#8c6d31')
double_arrow(ax, ta + 0.04, pa_val, ask[int(ta * 299)], 'spread', color='#1f77b4')
double_arrow(ax, te - 0.04, ask[int(te * 299)], pe_val, 'market\nimpact', color='#d62728')
double_arrow(ax, -0.085, pd_val, mid[int(te * 299)], 'price risk', color='#9467bd')

ax.legend(loc='lower right', fontsize=9, framealpha=0.8)
ax.set_title('Transaction cost components along the execution timeline (buy order)', fontsize=11)
fig.tight_layout()
fig.savefig(f'{FIGURES_DIR}/exec_cost_timeline.png', dpi=150, bbox_inches='tight')
plt.close(fig)
print('Saved exec_cost_timeline.png')

fig, ax = plt.subplots(figsize=(7, 5.5))

# hyperbola-like frontier: high market impact at low timing risk, vice versa
x = np.linspace(0.15, 2.0, 300)
a, b = 0.30, 0.05
y = a / (x - b) + 0.1

ax.plot(x, y, color='steelblue', lw=2.5)
ax.fill_between(x, y, 2.5, alpha=0.10, color='gray')
ax.text(1.0, 1.8, 'Infeasible\nregion', fontsize=9, color='gray', ha='center')

# key points A, B, C on frontier, D off frontier
pts = {'A': 0.22, 'B': 0.60, 'C': 1.60}
for name, xp in pts.items():
    yp = a / (xp - b) + 0.1
    ax.plot(xp, yp, '*', color='#d62728', ms=14, zorder=5)
    ax.text(xp - 0.07, yp + 0.07, name, fontsize=11, fontweight='bold', color='#444')

ax.plot(1.0, 1.2, '*', color='#d62728', ms=14, zorder=5)
ax.text(1.07, 1.25, 'D — "Irrational"\nstrategy', fontsize=8.5, color='#444')

ax.text(0.30, 1.8, 'A: High Impact,\nLow Risk', fontsize=8, color='#444')
ax.text(0.68, 0.75, 'B: Medium Impact,\nMedium Risk', fontsize=8, color='#444')
ax.text(1.20, 0.35, 'C: Low Impact,\nHigh Risk', fontsize=8, color='#444')

ax.set_xlabel('Timing Risk (variance of execution cost)', fontsize=10)
ax.set_ylabel('Market Impact (expected execution cost)', fontsize=10)
ax.set_title('Efficient Trading Frontier\n(Almgren & Chriss, 2000)', fontsize=11)
ax.set_xlim(0.0, 2.1); ax.set_ylim(0.0, 2.5)
ax.set_xticks([]); ax.set_yticks([])
ax.spines[['top', 'right']].set_visible(False)
fig.tight_layout()
fig.savefig(f'{FIGURES_DIR}/exec_efficient_frontier.png', dpi=150, bbox_inches='tight')
plt.close(fig)
print('Saved exec_efficient_frontier.png')

# Square-root market impact model vs linear model
q_over_v = np.logspace(-5, -2, 300)   # relative order size Q/V
alpha = 1.0

# normalise both models to the same impact at the midpoint
q_mid = np.sqrt(q_over_v[0] * q_over_v[-1])
impact_sqrt   = alpha * np.sqrt(q_over_v)                       # δ = 0.5
impact_linear = alpha * q_over_v / np.sqrt(q_mid)               # δ = 1.0, calibrated

fig, ax = plt.subplots(figsize=(8, 5))
ax.loglog(q_over_v, impact_sqrt,   lw=2,   label=r'Square-root ($\delta=0.5$)')
ax.loglog(q_over_v, impact_linear, lw=2, ls='--', label=r'Linear ($\delta=1.0$)')

ax.set_xlabel(r'Relative order size $Q/V$', fontsize=11)
ax.set_ylabel(r'Impact / volatility', fontsize=11)
ax.set_title('Square-root market impact model\n(Grinold & Kahn, 2000)', fontsize=11)
ax.legend(fontsize=10)
ax.grid(True, which='both', alpha=0.3)
fig.tight_layout()
plt.show()

def pretrade_market_impact(Q, V, sigma, alpha=1.0, spread_bps=5):
    """Estimate pre-trade market impact in basis points using the square-root model.

    Parameters
    ----------
    Q     : order size (shares or notional)
    V     : average daily volume (same units as Q)
    sigma : average daily volatility (as a fraction, e.g. 0.015 for 1.5%)
    alpha : model constant (calibrated, default 1.0)
    spread_bps : bid-ask spread in basis points

    Returns
    -------
    impact_bps : expected market impact in basis points
    """
    spread_cost = spread_bps / 2           # half-spread cost (crossing the spread once)
    impact_bps = alpha * sigma * np.sqrt(Q / V) * 1e4  # convert fraction to bps
    return spread_cost + impact_bps

# Example: order = 5% of ADV, daily vol = 1.5%, spread = 5 bps
Q, V, sigma, spread_bps = 100_000, 2_000_000, 0.015, 5
impact = pretrade_market_impact(Q, V, sigma, spread_bps=spread_bps)
print(f'Order size:      {Q:,} ({100*Q/V:.1f}% of ADV)')
print(f'Estimated impact: {impact:.1f} bps')

# Sensitivity: impact vs order size as % of ADV
fractions = np.array([0.01, 0.02, 0.05, 0.10, 0.20, 0.30, 0.50])
impacts = [pretrade_market_impact(f * V, V, sigma, spread_bps=spread_bps)
           for f in fractions]
for f, imp in zip(fractions, impacts):
    print(f'  Q/V = {f:.0%}  →  {imp:.1f} bps')

# Synthetic intraday volume profiles: typical day (U-shape) and event day
bins = np.arange(9, 17.5, 0.5)    # 30-min bins from 9:00 to 17:00
n = len(bins)
t_norm = np.linspace(0, 1, n)      # normalised time 0→1

# Typical day: U-shape — high at open/close, low at mid-day
vol_typical = 0.5 * np.exp(-6 * t_norm) + 0.5 * np.exp(-6 * (1 - t_norm)) + 0.15
vol_typical += np.abs(np.random.RandomState(0).randn(n)) * 0.02
vol_typical /= vol_typical.sum()          # normalise to relative fractions

# Event day (e.g. option expiry): spike at mid-morning (11:30 fix)
vol_event = vol_typical.copy()
spike_idx = np.argmin(np.abs(bins - 11.5))
vol_event[spike_idx] += 0.15
vol_event /= vol_event.sum()

bin_labels = [f'{int(h)}:{int((h%1)*60):02d}' for h in bins]

fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=False)
for ax, vol, title in zip(
        axes,
        [vol_typical, vol_event],
        ['Typical day (U-shaped volume profile)', 'Option expiration day']):
    ax.bar(range(n), vol * 100, color='steelblue', alpha=0.8)
    ax.set_xticks(range(0, n, 2))
    ax.set_xticklabels(bin_labels[::2], rotation=45, ha='right', fontsize=8)
    ax.set_xlabel('Time of Day')
    ax.set_ylabel('Volume (%)')
    ax.set_title(title)
    ax.set_ylim(0, max(vol_event.max() * 100 * 1.15, vol_typical.max() * 100 * 1.15))

fig.suptitle('Intraday volume profiles — VWAP static curve', fontsize=12)
fig.tight_layout()
plt.show()

def simulate_execution(Q, n_bins, sigma_per_bin, k_impact, vol_curve,
                        strategy='twap', pov_rate=0.10, seed=0):
    """Simulate a buy execution and compute benchmark metrics.

    Parameters
    ----------
    Q           : total order size
    n_bins      : number of time bins
    sigma_per_bin : price volatility per bin
    k_impact    : linear market impact coefficient (price per unit of q)
    vol_curve   : relative volume curve (sums to 1), length n_bins
    strategy    : 'twap', 'vwap', or 'pov'
    pov_rate    : target participation rate for 'pov' strategy

    Returns
    -------
    dict with execution metrics
    """
    rng = np.random.RandomState(seed)
    p0 = 100.0                            # arrival price

    # market volumes per bin (Poisson-distributed around expected)
    avg_market_vol = 1_000_000            # total daily market volume
    mkt_vols = rng.poisson(vol_curve * avg_market_vol).astype(float)

    # price random walk
    price = np.zeros(n_bins + 1)
    price[0] = p0
    for i in range(1, n_bins + 1):
        price[i] = price[i - 1] * (1 + rng.normal(0, sigma_per_bin))

    # determine slice sizes by strategy
    if strategy == 'twap':
        slices = np.full(n_bins, Q / n_bins)
    elif strategy == 'vwap':
        slices = Q * vol_curve           # proportional to volume curve
    elif strategy == 'pov':
        slices = pov_rate * mkt_vols     # track market volume
        slices = np.minimum(slices, Q)   # can't execute more than order
    else:
        raise ValueError(f'Unknown strategy: {strategy}')

    # clip and track remaining quantity
    executed = np.zeros(n_bins)
    remaining = Q
    exec_prices = np.zeros(n_bins)
    for i in range(n_bins):
        q_i = min(slices[i], remaining)
        impact_i = k_impact * q_i
        exec_prices[i] = price[i] + impact_i    # linear temporary impact
        executed[i] = q_i
        remaining -= q_i
        if remaining <= 0:
            break

    # aggregate execution metrics
    total_exec = executed.sum()
    p_avg = (executed * exec_prices).sum() / total_exec if total_exec > 0 else np.nan

    # IS benchmark: arrival price p0
    opportunity_cost = remaining * (price[-1] - p0)
    IS_cost = total_exec * (p_avg - p0) + opportunity_cost
    IS_pnl_bps = (p0 - p_avg) / p0 * 1e4        # positive = better than arrival

    # TWAP benchmark
    twap = price[:-1].mean()                     # simple average of bin open prices
    TWAP_cost = total_exec * (p_avg - twap)
    TWAP_pnl_bps = (twap - p_avg) / twap * 1e4

    # VWAP benchmark (volume-weighted market price)
    vwap = (mkt_vols * price[:-1]).sum() / mkt_vols.sum()
    VWAP_cost = total_exec * (p_avg - vwap)
    VWAP_pnl_bps = (vwap - p_avg) / vwap * 1e4

    # PoV benchmark
    pov_actual = (executed / mkt_vols).mean()

    return {
        'strategy': strategy, 'p_avg': p_avg, 'fill_rate': total_exec / Q,
        'IS_cost': IS_cost, 'IS_pnl_bps': IS_pnl_bps,
        'TWAP_cost': TWAP_cost, 'TWAP_pnl_bps': TWAP_pnl_bps,
        'VWAP_cost': VWAP_cost, 'VWAP_pnl_bps': VWAP_pnl_bps,
        'pov_actual': pov_actual, 'twap': twap, 'vwap': vwap,
    }

# Run all three strategies on the same simulated day
params = dict(
    Q=50_000, n_bins=len(vol_typical), sigma_per_bin=0.002,
    k_impact=1e-6, vol_curve=vol_typical, seed=42
)

results = {s: simulate_execution(strategy=s, **params)
           for s in ['twap', 'vwap', 'pov']}

print(f"{'Strategy':<8}  {'Avg price':>10}  {'Fill rate':>9}  "
      f"{'IS (bps)':>9}  {'TWAP (bps)':>11}  {'VWAP (bps)':>11}")
print('-' * 68)
for s, r in results.items():
    print(f"{r['strategy'].upper():<8}  {r['p_avg']:>10.4f}  {r['fill_rate']:>9.1%}  "
          f"{r['IS_pnl_bps']:>9.2f}  {r['TWAP_pnl_bps']:>11.2f}  {r['VWAP_pnl_bps']:>11.2f}")

# Visualise execution schedule and price path for TWAP vs VWAP
rng = np.random.RandomState(42)
avg_market_vol = 1_000_000
mkt_vols = rng.poisson(vol_typical * avg_market_vol).astype(float)

n_b = len(vol_typical)
price = np.zeros(n_b + 1); price[0] = 100.0
for i in range(1, n_b + 1):
    price[i] = price[i - 1] * (1 + rng.normal(0, 0.002))

Q = 50_000
twap_slices = np.full(n_b, Q / n_b)
vwap_slices = Q * vol_typical

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 7), sharex=True)

bins_x = np.arange(n_b)
ax1.bar(bins_x - 0.2, twap_slices, width=0.35, label='TWAP', alpha=0.8, color='steelblue')
ax1.bar(bins_x + 0.2, vwap_slices, width=0.35, label='VWAP', alpha=0.8, color='darkorange')
ax1.bar(bins_x, mkt_vols / 20, width=0.05, label='Market vol (scaled)', color='gray', alpha=0.5)
ax1.set_ylabel('Child order size (shares)')
ax1.set_title('Execution schedule: TWAP vs VWAP')
ax1.legend()

ax2.plot(np.arange(n_b + 1), price, lw=2, color='black', label='mid-price')
ax2.axhline(price[0], ls='--', color='green', alpha=0.7, label=f'Arrival price p₀ = {price[0]:.2f}')
twap_bench = price[:n_b].mean()
ax2.axhline(twap_bench, ls=':', color='steelblue', alpha=0.9, label=f'TWAP = {twap_bench:.3f}')
vwap_bench = (mkt_vols * price[:n_b]).sum() / mkt_vols.sum()
ax2.axhline(vwap_bench, ls=':', color='darkorange', alpha=0.9, label=f'VWAP = {vwap_bench:.3f}')
ax2.set_ylabel('Price')
ax2.set_xlabel('Time bin')
ax2.legend(fontsize=9)

fig.tight_layout()
plt.show()

def pretrade_tca(Q, V, sigma_daily, n_bins, spread_bps=5, alpha=1.0,
                 risk_aversion=0.0):
    """Pre-trade TCA: estimate expected cost and timing risk for TWAP vs IS.

    Under linear temporary impact p_exec = p + k*q_i, the expected IS cost
    of TWAP (equal slices q_i = Q/n) is:
      E[IS] = k * Q^2 / n + spread_cost
    Timing risk (variance) scales as:
      Var[IS] ≈ sigma^2 * (Q^2 / n) * T   (T normalised to 1)
    """
    k = alpha * sigma_daily * np.sqrt(1 / V)   # linear impact coefficient
    q_slice = Q / n_bins                       # TWAP slice size
    sigma_bin = sigma_daily / np.sqrt(252 * n_bins)  # per-bin volatility

    # expected impact cost (deterministic under linear impact)
    impact_cost = k * Q * q_slice * n_bins      # = k * Q^2
    spread_cost = spread_bps / 2 / 1e4 * Q     # total spread cost
    expected_cost = impact_cost + spread_cost

    # timing risk (std dev of price movement on remaining inventory)
    # sum over bins: Var ~ sigma_bin^2 * sum_{i=0}^{n-1} ((n-i)/n * Q)^2
    remaining_fracs = np.array([(n_bins - i) / n_bins for i in range(n_bins)])
    timing_var = sigma_bin**2 * np.sum((remaining_fracs * Q)**2)
    timing_risk_bps = np.sqrt(timing_var) / 100 * 1e4   # convert to bps

    return {
        'expected_cost_bps': expected_cost / (Q * 100) * 1e4,
        'timing_risk_bps': timing_risk_bps,
        'impact_cost_bps': impact_cost / (Q * 100) * 1e4,
    }

# Pre-trade TCA for different execution horizons (number of bins)
Q_ex = 100_000
V_ex = 2_000_000
sigma_ex = 0.015

print(f"Pre-trade TCA: Q={Q_ex:,} ({100*Q_ex/V_ex:.0f}% ADV), σ={sigma_ex:.1%}")
print(f"{'Horizon (bins)':<16} {'E[cost] (bps)':>14} {'Timing risk (bps)':>18} {'Impact (bps)':>13}")
print('-' * 65)
for n in [4, 8, 16, 32, 64]:
    tca = pretrade_tca(Q_ex, V_ex, sigma_ex, n)
    print(f"{n:<16} {tca['expected_cost_bps']:>14.2f} "
          f"{tca['timing_risk_bps']:>18.2f} {tca['impact_cost_bps']:>13.2f}")

# Visualise the cost-risk frontier for different horizons
horizons = np.arange(2, 65)
costs, risks = [], []
for n in horizons:
    tca = pretrade_tca(Q_ex, V_ex, sigma_ex, n)
    costs.append(tca['expected_cost_bps'])
    risks.append(tca['timing_risk_bps'])

fig, ax = plt.subplots(figsize=(8, 5))
sc = ax.scatter(risks, costs, c=horizons, cmap='plasma', s=30, zorder=3)
ax.plot(risks, costs, '-', color='gray', alpha=0.4, lw=1)
plt.colorbar(sc, ax=ax, label='Number of time bins (execution horizon)')

# annotate a few horizons
for n in [4, 16, 64]:
    tca = pretrade_tca(Q_ex, V_ex, sigma_ex, n)
    ax.annotate(f'n={n}', (tca['timing_risk_bps'], tca['expected_cost_bps']),
                textcoords='offset points', xytext=(6, 4), fontsize=9)

ax.set_xlabel('Timing Risk — std dev of execution cost (bps)', fontsize=10)
ax.set_ylabel('Expected execution cost (bps)', fontsize=10)
ax.set_title('Pre-trade efficient frontier:\n'
             'cost vs risk as execution horizon varies', fontsize=11)
ax.grid(True, alpha=0.3)
fig.tight_layout()
plt.show()

def ac_trajectory(kappa, T, X=1.0, n=300):
    """Almgren-Chriss optimal inventory trajectory x(t)."""
    t = np.linspace(0, T, n)
    if kappa < 1e-8:
        x = X * (1.0 - t / T)
    else:
        x = X * np.sinh(kappa * (T - t)) / np.sinh(kappa * T)
    return t, x


def ac_trading_rate(kappa, T, X=1.0, n=300):
    """Almgren-Chriss optimal trading rate v(t) = -dx/dt."""
    t = np.linspace(0, T, n)
    if kappa < 1e-8:
        v = np.full(n, X / T)
    else:
        v = kappa * X * np.cosh(kappa * (T - t)) / np.sinh(kappa * T)
    return t, v


T, X = 1.0, 1.0
kappa_T_vals = [0.0, 0.5, 1.0, 2.0, 4.0]
colors_ac    = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd']

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

for kT, col in zip(kappa_T_vals, colors_ac):
    kappa = kT / T
    t, x = ac_trajectory(kappa, T, X)
    t, v = ac_trading_rate(kappa, T, X)
    lbl = r'TWAP ($\kappa T=0$)' if kT == 0.0 else rf'$\kappa T = {kT}$'
    ls  = '--' if kT == 0.0 else '-'
    ax1.plot(t, x, color=col, ls=ls, lw=2, label=lbl)
    ax2.plot(t, v, color=col, ls=ls, lw=2, label=lbl)

for ax, ylabel, title in [
    (ax1, r'Remaining inventory $x_t/X$',    'Optimal inventory trajectory'),
    (ax2, r'Trading rate $v_t \cdot T / X$', 'Optimal trading rate'),
]:
    ax.set_xlabel(r'Normalised time $t/T$', fontsize=11)
    ax.set_ylabel(ylabel, fontsize=11)
    ax.set_title(title, fontsize=11)
    ax.legend(fontsize=9)
    ax.set_xlim(0, 1)
    ax.grid(True, alpha=0.3)
ax1.set_ylim(0, 1)

fig.tight_layout()
fig.savefig(f'{FIGURES_DIR}/ac_trajectories.png', dpi=150, bbox_inches='tight')
plt.close(fig)
print('Saved ac_trajectories.png')

def ac_cost(kappa, sigma, eta, gamma, X, T):
    """Expected IS and Var[IS] under the Almgren-Chriss optimal strategy."""
    perm = 0.5 * gamma * X**2
    if kappa < 1e-8:
        return perm + eta * X**2 / T,  sigma**2 * X**2 * T / 3.0
    coth_half = 1.0 / np.tanh(kappa * T / 2.0)
    E   = perm + 0.5 * eta * kappa * X**2 * coth_half
    Var = (sigma**2 * X**2) / (2.0 * kappa) * coth_half
    return E, Var


sigma_f, eta_f, gamma_f, X_f, T_f = 1.0, 1.0, 0.0, 1.0, 1.0

kappa_range = np.logspace(-3, 1.5, 500)
E_f, V_f    = zip(*[ac_cost(k, sigma_f, eta_f, gamma_f, X_f, T_f)
                     for k in kappa_range])
E_f, std_f  = np.array(E_f), np.sqrt(np.array(V_f))

kappa_pts = [0.001, 0.5, 1.0, 2.0, 5.0]
pt_labels = [
    r'$\kappa T\approx0$ (TWAP)',
    r'$\kappa T=0.5$', r'$\kappa T=1$',
    r'$\kappa T=2$',   r'$\kappa T=5$ (aggressive)',
]

fig, ax = plt.subplots(figsize=(9, 6))
sc = ax.scatter(std_f, E_f, c=np.log10(kappa_range),
                cmap='plasma', s=8, zorder=3)
plt.colorbar(sc, ax=ax, label=r'$\log_{10}(\kappa T)$')
ax.plot(std_f, E_f, '-', color='gray', alpha=0.25, lw=1)

for km, lbl in zip(kappa_pts, pt_labels):
    e, v = ac_cost(km, sigma_f, eta_f, gamma_f, X_f, T_f)
    ax.plot(np.sqrt(v), e, 'o', color='black', ms=8, zorder=5)
    ax.annotate(lbl, (np.sqrt(v), e),
                textcoords='offset points', xytext=(8, 2), fontsize=8.5)

ax.set_xlabel(r'$\sqrt{\operatorname{Var}[IS^*]}$ — timing risk', fontsize=11)
ax.set_ylabel(r'$\mathbb{E}[IS^*]$ — expected impact cost', fontsize=11)
ax.set_title('Almgren–Chriss efficient frontier', fontsize=11)
ax.grid(True, alpha=0.3)
fig.tight_layout()
fig.savefig(f'{FIGURES_DIR}/ac_frontier.png', dpi=150, bbox_inches='tight')
plt.close(fig)
print('Saved ac_frontier.png')

def ac_discrete_schedule(X, T, N, kappa):
    """Discrete Almgren-Chriss schedule: inventory x_k and child orders n_k."""
    tau = T / N
    t   = np.arange(N + 1) * tau
    if kappa < 1e-8:
        x = X * (1.0 - t / T)
    elif kappa * T > 100:           # sinh overflows; use exponential asymptote
        x = X * np.exp(-kappa * t)
        x[-1] = 0.0
    else:
        x = X * np.sinh(kappa * (T - t)) / np.sinh(kappa * T)
    n = x[:-1] - x[1:]
    return t, x, n


X_ex, T_ex, N_ex = 1_000_000, 1.0, 12
sigma_is, eta_is  = 0.015, 5e-9
lambda_vals       = [0, 1e4, 5e4, 2e5]
kappa_vals_ex     = [np.sqrt(lam * sigma_is**2 / eta_is) if lam > 0 else 0.0
                     for lam in lambda_vals]
labels_ex         = [r'TWAP ($\lambda=0$)', r'$\lambda=10^4$',
                     r'$\lambda=5\times10^4$', r'$\lambda=2\times10^5$']
colors_ex         = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728']
bin_labels        = [f"{9+i//2}:{'00' if i%2==0 else '30'}" for i in range(N_ex)]

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(11, 8), sharex=True)
width = 0.18
for idx, (kappa, lbl, col) in enumerate(zip(kappa_vals_ex, labels_ex, colors_ex)):
    t_d, x_d, n_d = ac_discrete_schedule(X_ex, T_ex, N_ex, kappa)
    ax1.bar(np.arange(N_ex) + (idx - 1.5) * width, n_d / 1e3,
            width=width, label=lbl, alpha=0.85, color=col)
    ax2.plot(np.arange(N_ex + 1) / N_ex, x_d / 1e6,
             '-o', ms=4, color=col, label=lbl, lw=1.8)

ax1.set_ylabel('Child order size (k shares)', fontsize=10)
ax1.set_title(f'Almgren–Chriss IS schedule — {X_ex:,} shares, {N_ex} bins', fontsize=11)
ax1.legend(fontsize=9)
ax1.set_xticks(range(N_ex))
ax1.set_xticklabels(bin_labels, rotation=30, ha='right', fontsize=8)
ax2.set_xlabel('Normalised time', fontsize=10)
ax2.set_ylabel('Remaining inventory (M shares)', fontsize=10)
ax2.set_title('Inventory trajectories', fontsize=11)
ax2.legend(fontsize=9); ax2.grid(True, alpha=0.3); ax2.set_xlim(0, 1)
fig.tight_layout()
plt.show()

_, _, n_ag = ac_discrete_schedule(X_ex, T_ex, N_ex, kappa_vals_ex[-1])
cumul = np.cumsum(n_ag)
print(f'\nSchedule (kT={kappa_vals_ex[-1]*T_ex:.1f}, most aggressive):')
print(f"{'Bin':>4}  {'Time':>6}  {'n_k':>14}  {'Cumul%':>7}")
for i, (ni, ci) in enumerate(zip(n_ag, cumul)):
    print(f'{i+1:>4}  {bin_labels[i]:>6}  {ni:>14,.0f}  {100*ci/X_ex:>6.1f}%')

def simulate_vwap_day(Pi, vol_true, vol_pred, sigma_bin, seed=42):
    """Simulate static and dynamic VWAP on one day.

    Dynamic strategy: q_n = Pi*(V_n + E_n[v_n])/E_n[V_M] - Q_n
    Static strategy : q_n = Pi * vol_pred[n] / sum(vol_pred)
    """
    rng         = np.random.RandomState(seed)
    n_bins_loc  = len(vol_true)
    avg_mkt     = 1_000_000
    realised    = rng.poisson(vol_true * avg_mkt).astype(float)
    V_M_true    = realised.sum()
    ret         = rng.normal(0, sigma_bin, n_bins_loc)
    prices      = np.cumprod(np.concatenate([[1.0], 1.0 + ret]))
    mid_loc     = prices[:-1]
    vwap_mkt    = (realised * mid_loc).sum() / V_M_true

    results = {}
    for strategy in ('static', 'dynamic'):
        exec_qty = np.zeros(n_bins_loc)
        Q_n = V_n = 0.0
        for n in range(n_bins_loc):
            remaining = Pi - Q_n
            if strategy == 'static':
                q_n = Pi * vol_pred[n] / vol_pred.sum()
            else:
                E_v_n = vol_pred[n] * avg_mkt
                E_V_M = V_n + vol_pred[n:].sum() * avg_mkt
                q_n   = max(Pi * (V_n + E_v_n) / max(E_V_M, 1) - Q_n, 0.0)
            q_n = min(q_n, remaining)
            if n == n_bins_loc - 1:
                q_n = remaining          # dump any residual
            exec_qty[n] = q_n
            Q_n += q_n;  V_n += realised[n]
        p_exec = (exec_qty * mid_loc).sum() / Pi
        results[strategy] = {
            'exec_qty':    exec_qty,
            'slippage_bps': (p_exec - vwap_mkt) / vwap_mkt * 1e4,
        }
    results.update(mid=mid_loc, realised_vol=realised, vwap_mkt=vwap_mkt)
    return results


n_bins_v, Pi_v, sigma_bv = len(vol_typical), 100_000, 0.0015

for vol_true, day_name in [(vol_typical, 'typical'), (vol_event, 'event')]:
    res = simulate_vwap_day(Pi_v, vol_true, vol_typical, sigma_bv, seed=7)
    print(f"{day_name.capitalize()} day  static: {res['static']['slippage_bps']:+.2f} bps"
          f"  |  dynamic: {res['dynamic']['slippage_bps']:+.2f} bps")

n_sims  = 2_000
mc_rng  = np.random.RandomState(0)
slippage_mc = {'static': [], 'dynamic': []}

for _ in range(n_sims):
    noise    = np.exp(mc_rng.randn(n_bins_v) * 0.3 - 0.3**2 / 2)
    vol_true = vol_typical * noise
    vol_true /= vol_true.sum()
    seed_s   = mc_rng.randint(0, 100_000)
    res = simulate_vwap_day(Pi_v, vol_true, vol_typical, sigma_bv, seed=seed_s)
    slippage_mc['static'].append(res['static']['slippage_bps'])
    slippage_mc['dynamic'].append(res['dynamic']['slippage_bps'])

fig, axes = plt.subplots(1, 2, figsize=(12, 5))
palette   = {'static': 'steelblue', 'dynamic': 'darkorange'}
for ax, strat in zip(axes, ('static', 'dynamic')):
    bps  = np.array(slippage_mc[strat])
    rmse = np.sqrt(np.mean(bps**2))
    ax.hist(bps, bins=60, color=palette[strat], alpha=0.85, density=True)
    ax.axvline(bps.mean(), color='black', lw=1.5, ls='--',
               label=f'Mean = {bps.mean():.2f} bps')
    ax.axvline(np.percentile(bps,  5), color='red', lw=1.2, ls=':'
               , label=f'5th pct = {np.percentile(bps, 5):.2f} bps')
    ax.axvline(np.percentile(bps, 95), color='red', lw=1.2, ls=':'
               , label=f'95th pct = {np.percentile(bps, 95):.2f} bps')
    ax.set_title(f'{strat.capitalize()} VWAP — RMSE = {rmse:.2f} bps', fontsize=11)
    ax.set_xlabel('VWAP slippage (bps)', fontsize=10)
    ax.set_ylabel('Density', fontsize=10)
    ax.legend(fontsize=9)

fig.suptitle(f'Static vs Dynamic VWAP slippage ({n_sims:,} Monte Carlo paths)',
             fontsize=12)
fig.tight_layout()
plt.show()

s_rmse = np.sqrt(np.mean(np.array(slippage_mc['static'])**2))
d_rmse = np.sqrt(np.mean(np.array(slippage_mc['dynamic'])**2))
print(f'RMSE improvement: {100*(s_rmse-d_rmse)/s_rmse:.1f}%  ({s_rmse:.2f} -> {d_rmse:.2f} bps)')