"""Optimal execution models.
Implements the Almgren-Chriss (2001) framework for optimal trade
execution, balancing expected market-impact cost against execution risk.
"""
from __future__ import annotations
import numpy as np
from numpy.typing import NDArray
from wraquant.core._coerce import coerce_array
[docs]
def almgren_chriss(
total_qty: float,
sigma: float,
eta: float,
gamma: float,
lambda_risk: float,
n_periods: int,
) -> NDArray[np.floating]:
"""Almgren-Chriss optimal execution trajectory.
Use Almgren-Chriss when you need the optimal rate at which to
liquidate a large position, balancing the cost of trading quickly
(market impact) against the risk of trading slowly (price
uncertainty). This is the foundational model of optimal execution.
Minimises a mean-variance objective:
E[cost] + lambda * Var[cost]
where cost has a temporary component (eta * trade_size^2) and a
permanent component (gamma * trade_size * remaining_position).
The solution is an exponentially decaying trajectory:
x_j = X * sinh(kappa * (T - j)) / sinh(kappa * T)
where kappa = sqrt(lambda * sigma^2 / eta).
Parameters:
total_qty: Total shares to execute.
sigma: Price volatility per period (e.g., daily std of returns
times price).
eta: Temporary impact coefficient (cost per share^2 traded).
gamma: Permanent impact coefficient (permanent price shift per
share traded).
lambda_risk: Risk-aversion parameter. ``lambda_risk = 0``
yields the risk-neutral (minimum-cost) linear trajectory.
Higher values front-load execution (trade faster to reduce
risk).
n_periods: Number of execution periods.
Returns:
Optimal holdings trajectory of length ``n_periods + 1``,
starting at *total_qty* and ending at 0. Take differences
to get the quantity to trade in each period.
Example:
>>> trajectory = almgren_chriss(10_000, sigma=0.02, eta=0.001,
... gamma=0.0001, lambda_risk=1e-4,
... n_periods=20)
>>> trajectory[0]
10000.0
>>> trajectory[-1]
0.0
>>> len(trajectory)
21
References:
- Almgren & Chriss (2001), "Optimal Execution of Portfolio
Transactions"
See Also:
optimal_execution_cost: Compute cost/risk for a given trajectory.
execution_frontier: Sweep lambda to trace cost-risk frontier.
"""
if n_periods <= 0:
raise ValueError("n_periods must be positive")
tau = 1.0 # normalised period length
# Almgren-Chriss kappa
kappa_sq = lambda_risk * sigma**2 / eta if eta > 0 else 0.0
kappa = np.sqrt(kappa_sq) if kappa_sq > 0 else 0.0
trajectory = np.zeros(n_periods + 1, dtype=np.float64)
trajectory[0] = total_qty
if kappa < 1e-12:
# Risk-neutral: linear liquidation
for j in range(1, n_periods + 1):
trajectory[j] = total_qty * (1.0 - j / n_periods)
else:
denom = np.sinh(kappa * n_periods * tau)
if abs(denom) < 1e-15:
# Fallback to linear
for j in range(1, n_periods + 1):
trajectory[j] = total_qty * (1.0 - j / n_periods)
else:
for j in range(1, n_periods + 1):
trajectory[j] = (
total_qty * np.sinh(kappa * (n_periods - j) * tau) / denom
)
# Enforce terminal condition
trajectory[-1] = 0.0
return trajectory
[docs]
def optimal_execution_cost(
trajectory: NDArray[np.floating],
sigma: float,
eta: float,
gamma: float,
) -> dict[str, float]:
"""Compute expected cost and variance of a given execution trajectory.
Use this to evaluate execution strategies by computing their
expected cost and execution risk. Compare different trajectories
(e.g., aggressive vs. passive) to find the right trade-off.
Parameters:
trajectory: Holdings path of length ``n_periods + 1``, starting
at the initial position and ending at 0.
sigma: Price volatility per period.
eta: Temporary impact coefficient.
gamma: Permanent impact coefficient.
Returns:
Dictionary with:
- ``'expected_cost'``: Expected total execution cost (temporary
+ permanent impact). Lower is better.
- ``'variance'``: Variance of execution cost due to price
uncertainty while holding the position.
- ``'std_dev'``: Standard deviation of execution cost. This
is the "execution risk."
Example:
>>> import numpy as np
>>> traj = np.linspace(10_000, 0, 21) # linear liquidation
>>> metrics = optimal_execution_cost(traj, sigma=0.02, eta=0.001, gamma=0.0001)
>>> metrics['expected_cost'] > 0
True
>>> metrics['std_dev'] > 0
True
See Also:
almgren_chriss: Compute the optimal trajectory.
execution_frontier: Trace the cost-risk frontier.
"""
trajectory = coerce_array(trajectory, "trajectory")
n = len(trajectory) - 1
trades = -np.diff(trajectory) # quantities sold each period (positive)
# Temporary impact cost: eta * sum(n_j^2)
temp_cost = eta * float(np.sum(trades**2))
# Permanent impact cost: gamma * sum(n_j * x_j) where x_j is remaining
perm_cost = gamma * float(np.sum(trades * trajectory[:-1]))
expected_cost = temp_cost + perm_cost
# Variance of execution: sigma^2 * sum(x_j^2)
variance = sigma**2 * float(np.sum(trajectory[1:] ** 2))
std_dev = np.sqrt(variance)
return {
"expected_cost": expected_cost,
"variance": variance,
"std_dev": std_dev,
}
[docs]
def execution_frontier(
total_qty: float,
sigma: float,
eta: float,
gamma: float,
n_points: int = 20,
n_periods: int = 20,
) -> dict[str, NDArray[np.floating]]:
"""Efficient frontier of (expected cost, risk) pairs.
Sweeps over a range of risk-aversion parameters to trace out the
trade-off between expected execution cost and execution risk.
Parameters:
total_qty: Total shares to execute.
sigma: Price volatility per period.
eta: Temporary impact coefficient.
gamma: Permanent impact coefficient.
n_points: Number of points on the frontier.
n_periods: Number of execution periods for each trajectory.
Returns:
Dictionary with ``'lambda_values'``, ``'expected_cost'``, and
``'std_dev'`` arrays, each of length *n_points*.
"""
# Sweep lambda from near-zero (risk-neutral) to aggressive risk aversion
lambda_values = np.logspace(-4, 2, n_points)
costs = np.zeros(n_points, dtype=np.float64)
stds = np.zeros(n_points, dtype=np.float64)
for i, lam in enumerate(lambda_values):
traj = almgren_chriss(total_qty, sigma, eta, gamma, lam, n_periods)
metrics = optimal_execution_cost(traj, sigma, eta, gamma)
costs[i] = metrics["expected_cost"]
stds[i] = metrics["std_dev"]
return {
"lambda_values": lambda_values,
"expected_cost": costs,
"std_dev": stds,
}
[docs]
def bertsimas_lo(
total_shares: float,
n_periods: int,
volatility: float,
impact_coeff: float,
risk_aversion: float = 0.0,
) -> dict[str, object]:
"""Bertsimas-Lo (1998) optimal execution with discrete trading.
Use this for a discrete-time optimal execution model that minimises
expected cost plus a risk penalty. Unlike Almgren-Chriss which
models continuous trading, Bertsimas-Lo explicitly works in discrete
periods and provides a closed-form solution under a linear
temporary impact model.
The optimal strategy in each period trades a constant fraction of
the remaining position. With zero risk aversion, this simplifies
to the naive strategy of trading ``total_shares / n_periods`` per
period (linear liquidation).
Temporary impact model:
cost_t = impact_coeff * n_t^2
where n_t is the number of shares traded in period t.
Total expected cost:
E[cost] = impact_coeff * sum(n_t^2)
Parameters:
total_shares: Total shares to liquidate (positive).
n_periods: Number of discrete trading periods.
volatility: Per-period price volatility (standard deviation).
impact_coeff: Temporary impact coefficient. The cost of
trading n shares in one period is ``impact_coeff * n^2``.
risk_aversion: Risk-aversion parameter (default 0.0 for
risk-neutral). Higher values front-load execution to
reduce variance.
Returns:
Dictionary containing:
- **trajectory** (*NDArray*) -- Holdings path of length
``n_periods + 1``, starting at *total_shares* and ending at 0.
- **trades** (*NDArray*) -- Shares traded per period (length
``n_periods``).
- **expected_cost** (*float*) -- Expected total execution cost.
- **cost_variance** (*float*) -- Variance of execution cost
from holding risk.
Example:
>>> result = bertsimas_lo(10_000, n_periods=20, volatility=0.02,
... impact_coeff=0.001, risk_aversion=0.0)
>>> result['trajectory'][0]
10000.0
>>> result['trajectory'][-1]
0.0
>>> len(result['trades'])
20
References:
- Bertsimas & Lo (1998), "Optimal Control of Execution Costs"
*Journal of Financial Markets*, 1(1), 1-50.
See Also:
almgren_chriss: Continuous-time optimal execution.
optimal_execution_cost: Evaluate a given trajectory.
"""
if n_periods <= 0:
raise ValueError("n_periods must be positive")
# Compute optimal fraction to trade each period
# With risk aversion, the Bertsimas-Lo optimal strategy trades
# more aggressively early. The fraction depends on kappa.
if impact_coeff > 0 and risk_aversion > 0:
kappa_sq = risk_aversion * volatility**2 / impact_coeff
kappa = np.sqrt(kappa_sq) if kappa_sq > 0 else 0.0
else:
kappa = 0.0
trajectory = np.zeros(n_periods + 1, dtype=np.float64)
trajectory[0] = total_shares
if kappa < 1e-12:
# Risk-neutral: linear liquidation (trade equal amounts each period)
for j in range(1, n_periods + 1):
trajectory[j] = total_shares * (1.0 - j / n_periods)
else:
# Risk-averse: exponentially decaying trajectory
denom = np.sinh(kappa * n_periods)
if abs(denom) < 1e-15:
for j in range(1, n_periods + 1):
trajectory[j] = total_shares * (1.0 - j / n_periods)
else:
for j in range(1, n_periods + 1):
trajectory[j] = (
total_shares * np.sinh(kappa * (n_periods - j)) / denom
)
trajectory[-1] = 0.0
trades = -np.diff(trajectory) # positive = shares sold
# Expected cost = impact_coeff * sum(n_t^2)
expected_cost = float(impact_coeff * np.sum(trades**2))
# Variance from holding risk: sigma^2 * sum(x_t^2) for remaining positions
cost_variance = float(volatility**2 * np.sum(trajectory[1:] ** 2))
return {
"trajectory": trajectory,
"trades": trades,
"expected_cost": expected_cost,
"cost_variance": cost_variance,
}
[docs]
def estimate_impact_params(
trade_prices: "pd.Series",
volumes: "pd.Series",
directions: "pd.Series",
) -> dict[str, float]:
"""Estimate temporary and permanent impact coefficients from trade data.
Bridges ``wraquant.microstructure.liquidity`` into the Almgren-Chriss
framework by calibrating impact parameters from observed trade data.
The permanent impact coefficient (gamma) is derived from the
microstructure ``price_impact`` function, and the temporary impact
coefficient (eta) is estimated from the effective spread.
Parameters:
trade_prices: Executed trade prices.
volumes: Volume for each trade.
directions: Trade direction (+1 buy, -1 sell).
Returns:
Dictionary with ``eta`` (temporary impact coefficient) and
``gamma`` (permanent impact coefficient) suitable for use with
:func:`almgren_chriss` and :func:`optimal_execution_cost`.
Example:
>>> import pandas as pd, numpy as np
>>> prices = pd.Series([100.0, 100.05, 100.08, 99.95])
>>> vols = pd.Series([1000, 2000, 1500, 1800])
>>> dirs = pd.Series([1, 1, -1, -1])
>>> params = estimate_impact_params(prices, vols, dirs)
>>> 'eta' in params and 'gamma' in params
True
See Also:
almgren_chriss: Uses the estimated parameters.
wraquant.microstructure.liquidity.price_impact: Source of impact data.
"""
import pandas as pd
from wraquant.microstructure.liquidity import price_impact
impact = price_impact(trade_prices, volumes, directions)
clean_impact = impact.dropna()
# Permanent impact: mean absolute impact per unit volume
gamma = float(clean_impact.abs().mean()) if len(clean_impact) > 0 else 1e-4
# Temporary impact: estimated from spread (price changes proportional
# to trade size squared)
vol_arr = coerce_array(volumes, "volumes")
price_changes = trade_prices.diff().dropna().abs()
if len(price_changes) > 0 and len(vol_arr) > 1:
# eta ~ mean(|dp|) / mean(volume^2)
mean_dp = float(price_changes.mean())
mean_v2 = float(np.mean(vol_arr[1:] ** 2))
eta = mean_dp / mean_v2 if mean_v2 > 0 else 1e-4
else:
eta = 1e-4
return {
"eta": eta,
"gamma": gamma,
}