Optimal Execution — The Almgren-Chriss Framework

Companion notebook to the Optimal Execution article (Optimal Execution, Module 1). Uses Python stdlib (math) and numpy only — no external dependencies beyond Pyodide's bundled packages.

§1 — Imports and named constants
§2 — Optimal trajectory for a range of risk-aversion parameters; expected cost and variance
§3 — Efficient frontier: sweep λ, compute (expected cost, std) pairs
§4 — TWAP limit: verify λ→0 convergence analytically and numerically
§5 — Validation summary: pass/fail for all article claims

Press Shift+Enter or click ▶ Run on each cell in order. §2–§5 depend on §1.

Conventions used throughout: Time is measured in trading days. Volatility σ is daily (not annualised). Inventory X is in shares. Impact parameters η (temporary) and γ (permanent) have units of $ per share per (million shares / day) and $ per share respectively. All costs are in dollars.

§1 — Imports and Model Parameters

All model parameters are declared as named constants. No magic numbers appear in subsequent cells.

Parameter interpretation: ETA (temporary impact) measures the price discount per unit trading rate; GAMMA (permanent impact) measures the permanent price shift per share traded; SIGMA is daily price volatility; LAMBDA_VALS are the risk-aversion parameters whose units are 1/($ per share).

import math
import numpy as np

# ── Problem parameters ────────────────────────────────────────────────────────
X_SHARES    = 1_000_000      # total shares to liquidate
T_DAYS      = 5              # liquidation horizon (trading days)
N_INTERVALS = 5              # number of trading intervals
TAU         = T_DAYS / N_INTERVALS  # interval length (days), = 1 here

# ── Market parameters ─────────────────────────────────────────────────────────
SIGMA       = 0.02           # daily price volatility (decimal, not %)
S0          = 50.0           # initial mid-price ($)

# ── Impact parameters ─────────────────────────────────────────────────────────
# ETA: temporary impact coefficient
#   Units: $ per share per (shares/day). For 1M shares/day rate, impact = ETA * 1e6.
#   ETA = 0.1 means selling at 1M shares/day moves execution price by 0.1 * (1M/1) = $0.1M
#   Per-share: 0.1 * rate. Here we express rate as shares so ETA has units $/share^2 * day.
ETA         = 0.1e-6         # temporary impact: $ cost per share per (share/day) of rate
GAMMA_PERM  = 0.0            # permanent impact coefficient (set to 0 to isolate temporary)

# ── Risk-aversion parameters to examine ──────────────────────────────────────
LAMBDA_VALS = [0.0, 1e-6, 1e-5, 5e-5, 1e-4]

print("Almgren-Chriss Model — Parameters")
print("=" * 50)
print(f"  Shares to liquidate X:   {X_SHARES:>12,}")
print(f"  Horizon T:               {T_DAYS:>12} days")
print(f"  Intervals N:             {N_INTERVALS:>12}")
print(f"  Interval length tau:     {TAU:>12.2f} days")
print(f"  Daily volatility sigma:  {SIGMA:>12.4f}  ({SIGMA*100:.2f}%)")
print(f"  Initial price S0:        ${S0:>11.2f}")
print(f"  Temporary impact eta:    {ETA:>12.2e}  $/share/(share/day)")
print(f"  Permanent impact gamma:  {GAMMA_PERM:>12.4f}  $/share")
print(f"  Lambda values:           {LAMBDA_VALS}")
print()
print("Ready. Run §2 to compute optimal trajectories.")

§2 — Optimal Trajectories and Cost Decomposition

For each value of the risk-aversion parameter λ, we compute:

The urgency parameter κ = √(λσ²/η)
The optimal inventory path x(t) = X ⋅ sinh(κ(T−t)) / sinh(κT)
The trade sequence n_k = x(t_{k-1}) − x(t_k)
Expected temporary impact cost E[temp] = η ∑ n_k² / τ
Variance of IS = σ² τ ∑ x_{k-1}²

The λ = 0 case requires the limiting formula x(t) = X(T−t)/T (TWAP), since sinh(κT) → 0 as κ → 0 and direct evaluation is undefined.

# ── Helper: compute optimal Almgren-Chriss trajectory ────────────────────────
def ac_trajectory(lam, x_total, t_days, n_intervals, tau, sigma, eta):
    """
    Returns (times, inventory, trades, expected_cost, variance) for the
    Almgren-Chriss optimal schedule.

    Parameters
    ----------
    lam       : risk-aversion parameter (>= 0)
    x_total   : total shares to liquidate
    t_days    : liquidation horizon (days)
    n_intervals: number of intervals
    tau       : interval length (days)
    sigma     : daily volatility
    eta       : temporary impact coefficient

    Returns
    -------
    times     : array of shape (n_intervals+1,), time points t_k = k*tau
    inventory : array of shape (n_intervals+1,), inventory x_k
    trades    : array of shape (n_intervals,), shares sold n_k in each interval
    exp_cost  : expected temporary impact cost (permanent cost excluded)
    variance  : variance of implementation shortfall
    kappa     : urgency parameter
    """
    times = np.array([k * tau for k in range(n_intervals + 1)])

    if lam == 0.0:
        # TWAP limit: linear inventory decay
        kappa = 0.0
        inventory = x_total * (t_days - times) / t_days
    else:
        kappa = math.sqrt(lam * sigma**2 / eta)
        sinh_kT = math.sinh(kappa * t_days)
        if abs(sinh_kT) < 1e-12:
            # Numerical fallback to TWAP for very small kappa*T
            kappa = 0.0
            inventory = x_total * (t_days - times) / t_days
        else:
            inventory = x_total * np.sinh(kappa * (t_days - times)) / sinh_kT

    # Enforce exact boundary conditions
    inventory[0]  = x_total
    inventory[-1] = 0.0

    # Trade sizes: n_k = x_{k-1} - x_k
    trades = inventory[:-1] - inventory[1:]

    # Expected temporary impact cost: eta * sum(n_k^2 / tau)
    exp_cost = eta * float(np.sum(trades**2)) / tau

    # Variance of IS: sigma^2 * tau * sum(x_{k-1}^2)
    variance = sigma**2 * tau * float(np.sum(inventory[:-1]**2))

    return times, inventory, trades, exp_cost, variance, kappa


# ── Compute for each lambda ───────────────────────────────────────────────────
results = {}
for lam in LAMBDA_VALS:
    t, inv, tr, ec, var, kap = ac_trajectory(
        lam, X_SHARES, T_DAYS, N_INTERVALS, TAU, SIGMA, ETA
    )
    results[lam] = {
        "times": t, "inventory": inv, "trades": tr,
        "exp_cost": ec, "variance": var, "kappa": kap
    }

# ── Print trajectories ────────────────────────────────────────────────────────
print("Optimal Trajectories — Inventory Schedule (shares)")
print("=" * 75)

# Header
header_parts = ["  k  |  t_k  "]
for lam in LAMBDA_VALS:
    lam_label = "TWAP" if lam == 0.0 else ("lam={:.0e}".format(lam))
    header_parts.append(lam_label)
header = " | ".join(header_parts)
print(header)
print("  " + "-" * 72)

for k in range(N_INTERVALS + 1):
    row_parts = ["  {:>2} | {:>5.1f}  ".format(k, k * TAU)]
    for lam in LAMBDA_VALS:
        inv_k = results[lam]["inventory"][k]
        row_parts.append("{:>12,.0f}".format(inv_k))
    print(" | ".join(row_parts))

print()

# ── Print trade sequences ─────────────────────────────────────────────────────
print("Trade Sequence n_k (shares sold each interval)")
print("=" * 75)
header2_parts = ["  k  |  t_k  "]
for lam in LAMBDA_VALS:
    lam_label = "TWAP" if lam == 0.0 else ("lam={:.0e}".format(lam))
    header2_parts.append(lam_label)
print(" | ".join(header2_parts))
print("  " + "-" * 72)

for k in range(N_INTERVALS):
    row_parts = ["  {:>2} | {:>5.1f}  ".format(k + 1, (k + 1) * TAU)]
    for lam in LAMBDA_VALS:
        n_k = results[lam]["trades"][k]
        row_parts.append("{:>12,.0f}".format(n_k))
    print(" | ".join(row_parts))

print()

# ── Print cost summary ────────────────────────────────────────────────────────
print("Cost and Variance Summary")
print("=" * 75)
col_header = "  {:>14} | {:>14} | {:>14} | {:>14}".format(
    "lambda", "kappa (1/day)", "E[temp cost]($)", "std(IS) ($)"
)
print(col_header)
print("  " + "-" * 62)
for lam in LAMBDA_VALS:
    r = results[lam]
    lam_str = "0 (TWAP)" if lam == 0.0 else "{:.1e}".format(lam)
    print("  {:>14} | {:>14.6f} | {:>14,.2f} | {:>14,.2f}".format(
        lam_str, r["kappa"], r["exp_cost"], math.sqrt(r["variance"])
    ))

§3 — Efficient Frontier

We sweep λ over a logarithmic grid from near-zero to large values, computing the optimal (expected cost, std of IS) pair for each λ. This traces the efficient frontier in execution space.

The frontier should be a curve from the TWAP point (low cost, high variance) to the instant-liquidation point (high cost, zero variance). Any execution algorithm that falls above and to the right of this frontier is suboptimal.

Desk interpretation: Each point on the frontier is an optimal execution strategy parameterised by risk tolerance. The client's risk-aversion parameter λ maps directly to a position on the frontier and an associated algorithm parameter set (e.g. participation rate, urgency level).

# ── Efficient frontier: sweep lambda ─────────────────────────────────────────
N_FRONTIER_PTS = 40   # number of lambda values on the frontier grid

# Build lambda grid: 0 (TWAP) + logarithmic range up to fast execution
lam_grid = [0.0] + list(np.logspace(-7, -3, N_FRONTIER_PTS - 1))

frontier_pts = []
for lam in lam_grid:
    _, _, _, ec, var, kap = ac_trajectory(
        lam, X_SHARES, T_DAYS, N_INTERVALS, TAU, SIGMA, ETA
    )
    frontier_pts.append({
        "lambda": lam,
        "kappa":  kap,
        "exp_cost": ec,
        "std_is":   math.sqrt(var),
        "mv_obj":   ec + lam * var if lam > 0 else ec
    })

# ── Print selected frontier points ────────────────────────────────────────────
print("Efficient Frontier — Selected Points")
print("=" * 78)
print("  {:>10} | {:>12} | {:>16} | {:>16} | {:>10}".format(
    "lambda", "kappa(1/day)", "E[temp cost] ($)", "std(IS) ($)", "M-V obj"
))
print("  " + "-" * 74)

# Print every 5th point plus the endpoints
selected_indices = [0] + list(range(4, N_FRONTIER_PTS, 5)) + [N_FRONTIER_PTS - 1]
seen = set()
for idx in selected_indices:
    if idx in seen or idx >= len(frontier_pts):
        continue
    seen.add(idx)
    fp = frontier_pts[idx]
    lam_str = "0 (TWAP)" if fp["lambda"] == 0.0 else "{:.2e}".format(fp["lambda"])
    print("  {:>10} | {:>12.6f} | {:>16,.2f} | {:>16,.2f} | {:>10,.2f}".format(
        lam_str, fp["kappa"], fp["exp_cost"], fp["std_is"], fp["mv_obj"]
    ))

print()

# ── Verify frontier properties ────────────────────────────────────────────────
print("Frontier Structural Checks")
print("=" * 50)

# Check 1: expected cost is increasing in lambda (faster execution = more cost)
costs  = [fp["exp_cost"] for fp in frontier_pts]
stds   = [fp["std_is"]   for fp in frontier_pts]

cost_increasing = all(costs[i] <= costs[i+1] for i in range(len(costs)-1))
std_decreasing  = all(stds[i]  >= stds[i+1]  for i in range(len(stds)-1))

print("  Expected cost non-decreasing in lambda?  {}".format(
    "YES" if cost_increasing else "NO (check model)"
))
print("  Std(IS) non-increasing in lambda?        {}".format(
    "YES" if std_decreasing else "NO (check model)"
))

# Print TWAP vs fastest-execution comparison
twap_pt  = frontier_pts[0]
fast_pt  = frontier_pts[-1]
print()
print("  TWAP point:  E[cost]={:>12,.0f}  std={:>12,.0f}".format(
    twap_pt["exp_cost"], twap_pt["std_is"]
))
print("  Fast point:  E[cost]={:>12,.0f}  std={:>12,.0f}".format(
    fast_pt["exp_cost"], fast_pt["std_is"]
))
print()
print("  Cost increase from TWAP to fast: {:,.0f} ({:.1f}x)".format(
    fast_pt["exp_cost"] - twap_pt["exp_cost"],
    fast_pt["exp_cost"] / twap_pt["exp_cost"] if twap_pt["exp_cost"] > 0 else float("inf")
))
print("  Std reduction from TWAP to fast: {:,.0f} ({:.1f}% reduction)".format(
    twap_pt["std_is"] - fast_pt["std_is"],
    100.0 * (twap_pt["std_is"] - fast_pt["std_is"]) / twap_pt["std_is"]
    if twap_pt["std_is"] > 0 else 0.0
))

# Store for validation
frontier_cost_increasing_stored = cost_increasing
frontier_std_decreasing_stored  = std_decreasing
twap_cost_stored  = twap_pt["exp_cost"]
twap_std_stored   = twap_pt["std_is"]
fast_cost_stored  = fast_pt["exp_cost"]
fast_std_stored   = fast_pt["std_is"]

§4 — TWAP Limit: λ → 0 Convergence

The article claims that the optimal trajectory converges to the TWAP (linear) schedule as λ → 0. We verify this in two ways:

Analytically: Taylor expansion sinh(κu) ≈ κu for small κu gives x(t) → X(T−t)/T.
Numerically: compute the trajectory for λ = 10⁻⁸ and compare to the exact TWAP formula.

The Taylor expansion gives the first-order correction term, which quantifies how much the optimal schedule deviates from TWAP for small but non-zero λ. This correction always causes front-loading (sell more early) for any λ > 0.

# ── TWAP trajectory (exact linear formula) ───────────────────────────────────
t_grid = np.array([k * TAU for k in range(N_INTERVALS + 1)])
twap_inventory = X_SHARES * (T_DAYS - t_grid) / T_DAYS

# ── Near-TWAP trajectory: lambda = 1e-8 ──────────────────────────────────────
LAM_SMALL = 1e-8
_, near_twap_inv, _, _, _, kap_small = ac_trajectory(
    LAM_SMALL, X_SHARES, T_DAYS, N_INTERVALS, TAU, SIGMA, ETA
)

# ── First-order Taylor correction ─────────────────────────────────────────────
# x_opt(t) = X * sinh(kappa*(T-t)) / sinh(kappa*T)
# For small kappa: sinh(kappa*u) = kappa*u + (kappa*u)^3/6 + ...
# So x_opt(t) = X * [(T-t) + (kappa^2/6)*(T-t)^3 + ...] / [T + (kappa^2/6)*T^3 + ...]
# = X*(T-t)/T * [1 + (kappa^2/6)*((T-t)^2 - T^2) + O(kappa^4)]
# First correction: x_opt(t) approx X*(T-t)/T - X*(kappa^2/6)*((T-t) - (T-t)^3/T^2)
def twap_correction(t, x_total, t_horizon, kap):
    """First-order Taylor correction to TWAP for small kappa."""
    u = t_horizon - t
    twap_term = x_total * u / t_horizon
    # Correction: -X * (kappa^2/6) * (u - u^3/T^2)
    correction = -x_total * (kap**2 / 6.0) * (u - u**3 / t_horizon**2)
    return twap_term + correction

taylor_correction = np.array([
    twap_correction(t, X_SHARES, T_DAYS, kap_small) for t in t_grid
])

# ── Print comparison ──────────────────────────────────────────────────────────
print("TWAP Limit Verification")
print("=" * 75)
print("kappa for lambda={:.1e}: {:.8f} (1/day)".format(LAM_SMALL, kap_small))
print("kappa*T = {:.6f}  (should be << 1 for near-TWAP)".format(kap_small * T_DAYS))
print()
print("  {:>4} | {:>14} | {:>14} | {:>14} | {:>10}".format(
    "k", "TWAP (linear)", "AC (near-TWAP)", "Taylor approx", "|AC - TWAP|"
))
print("  " + "-" * 65)
for k in range(N_INTERVALS + 1):
    diff = abs(near_twap_inv[k] - twap_inventory[k])
    print("  {:>4} | {:>14,.2f} | {:>14,.2f} | {:>14,.2f} | {:>10.4f}".format(
        k, twap_inventory[k], near_twap_inv[k], taylor_correction[k], diff
    ))

print()

# ── Verify correction direction: should be front-loaded relative to TWAP ──────
# At t=0 both are X; check t=tau (first interval): AC should be LESS than TWAP
# (has sold MORE in first interval, so inventory is lower)
ac_first = near_twap_inv[1]
twap_first = twap_inventory[1]
front_loaded = ac_first < twap_first  # AC inventory lower => sold more early

print("Front-loading verification (lambda={:.0e}):".format(LAM_SMALL))
print("  Inventory after interval 1 (t={}):".format(TAU))
print("    TWAP:  {:>12,.0f} shares".format(twap_first))
print("    AC:    {:>12,.0f} shares".format(ac_first))
print("  AC sells more early (front-loaded)?  {}".format(
    "YES" if front_loaded else "NO"
))
print()

# ── Closeness check ───────────────────────────────────────────────────────────
max_abs_diff = float(np.max(np.abs(near_twap_inv - twap_inventory)))
rel_diff = max_abs_diff / X_SHARES
print("Max absolute deviation |AC - TWAP|: {:,.2f} shares".format(max_abs_diff))
print("Relative deviation:                 {:.6f}%".format(rel_diff * 100))
converges_ok = rel_diff < 1e-3   # less than 0.1% relative deviation

print()
print("Convergence to TWAP for lambda={:.0e}?  {}".format(
    LAM_SMALL, "YES (within 0.1% relative)" if converges_ok else "NO"
))

# ── Store for validation ──────────────────────────────────────────────────────
twap_convergence_stored   = converges_ok
front_loaded_stored       = front_loaded
max_rel_diff_stored       = rel_diff

§5 — Validation Summary

Final cell: verify all key claims from the article numerically. Each check is a precise assertion with a stated tolerance.

print("=" * 65)
print("VALIDATION SUMMARY — Almgren-Chriss Framework")
print("=" * 65)
print()

checks = []

# ── 1. TWAP schedule has uniform trades ───────────────────────────────────────
twap_trades = results[0.0]["trades"]
twap_uniform = float(np.std(twap_trades)) < 1.0  # all equal up to float rounding
checks.append(("TWAP (lambda=0) has uniform trade sizes", twap_uniform))
print("1. TWAP trade std = {:.4f} (should be < 1)  {}".format(
    float(np.std(twap_trades)), "✓" if twap_uniform else "✗"
))

# ── 2. Boundary conditions satisfied for all lambda ──────────────────────────
bc_ok = all(
    abs(results[lam]["inventory"][0] - X_SHARES) < 1.0 and
    abs(results[lam]["inventory"][-1]) < 1.0
    for lam in LAMBDA_VALS
)
checks.append(("Boundary conditions x(0)=X and x(T)=0 satisfied for all lambda", bc_ok))
print("2. Boundary conditions x(0)=X, x(T)=0 for all lambda?  {}".format(
    "✓" if bc_ok else "✗"
))

# ── 3. Trades are non-negative (no buying back) ───────────────────────────────
trades_nonneg = all(
    float(np.min(results[lam]["trades"])) >= -1.0   # allow tiny float rounding
    for lam in LAMBDA_VALS
)
checks.append(("All trades n_k >= 0 (no short-selling during liquidation)", trades_nonneg))
print("3. All trades n_k >= 0?  {}".format("✓" if trades_nonneg else "✗"))

# ── 4. Inventory is monotonically decreasing ─────────────────────────────────
inventory_decreasing = all(
    all(results[lam]["inventory"][k] >= results[lam]["inventory"][k+1] - 1.0
        for k in range(N_INTERVALS))
    for lam in LAMBDA_VALS
)
checks.append(("Inventory x_k is non-increasing for all lambda", inventory_decreasing))
print("4. Inventory non-increasing for all lambda?  {}".format(
    "✓" if inventory_decreasing else "✗"
))

# ── 5. Efficient frontier: cost non-decreasing in lambda ─────────────────────
checks.append(("Efficient frontier: expected cost non-decreasing in lambda", frontier_cost_increasing_stored))
print("5. Expected cost non-decreasing in lambda (frontier)?  {}".format(
    "✓" if frontier_cost_increasing_stored else "✗"
))

# ── 6. Efficient frontier: std non-increasing in lambda ──────────────────────
checks.append(("Efficient frontier: std(IS) non-increasing in lambda", frontier_std_decreasing_stored))
print("6. std(IS) non-increasing in lambda (frontier)?  {}".format(
    "✓" if frontier_std_decreasing_stored else "✗"
))

# ── 7. TWAP convergence for small lambda ─────────────────────────────────────
checks.append(("lambda->0 trajectory converges to TWAP (within 0.1%)", twap_convergence_stored))
print("7. Convergence to TWAP for lambda=1e-8?  {}  (max rel diff = {:.4e})".format(
    "✓" if twap_convergence_stored else "✗", max_rel_diff_stored
))

# ── 8. Front-loading: higher lambda => sell more early ───────────────────────
checks.append(("Positive lambda causes front-loading relative to TWAP", front_loaded_stored))
print("8. Front-loading for lambda > 0?  {}".format(
    "✓" if front_loaded_stored else "✗"
))

# ── 9. Kappa scaling: kappa increases with lambda ────────────────────────────
kappa_vals = [results[lam]["kappa"] for lam in LAMBDA_VALS if lam > 0]
kappa_increasing = all(kappa_vals[i] <= kappa_vals[i+1] for i in range(len(kappa_vals)-1))
checks.append(("Urgency kappa increases with lambda (kappa = sqrt(lambda*sigma^2/eta))", kappa_increasing))
print("9. Kappa increasing with lambda?  {}".format(
    "✓" if kappa_increasing else "✗"
))

# ── 10. TWAP minimises expected cost over all uniform-block schedules ─────────
# For fixed X and N, TWAP gives n_k = X/N for all k, which minimises sum(n_k^2)
# by Cauchy-Schwarz. Check that TWAP cost <= any other schedule's expected cost.
# Compare to the most aggressive schedule (lambda=1e-4).
twap_cost = results[0.0]["exp_cost"]
aggr_cost = results[1e-4]["exp_cost"]
twap_cheaper = twap_cost <= aggr_cost
checks.append(("TWAP has lower expected cost than aggressive schedule (lambda=1e-4)", twap_cheaper))
print("10. TWAP cost <= aggressive schedule cost?  {}  TWAP={:,.0f}  aggr={:,.0f}".format(
    "✓" if twap_cheaper else "✗", twap_cost, aggr_cost
))

print()
passed = sum(1 for _, ok in checks if ok)
print("=" * 65)
print("Results: {}/{} claims verified  {}".format(
    passed, len(checks), "✓" if passed == len(checks) else "— check failures above"
))
if passed == len(checks):
    print("All validations passed — article claims confirmed numerically.")