quant/src/marketmaking.rs
Quoting: the solvable model, and the two costs it does and does not capture.
//! Quoting: the solvable model, and the two costs it does and does not capture.//!//! Three separate things live here, matching the market making chapter's three//! claims.//!//! The Avellaneda-Stoikov problem, whose Hamilton-Jacobi-Bellman equation is//! nonlinear and becomes *linear* under an exponential substitution. The chapter//! makes that the same phenomenon as the solvable models chapter's Riccati//! linearisation, so both routes are computed here and compared: the nonlinear//! equation solved directly, and the linear system solved after substitution.//!//! The winner's curse in a request-for-quote auction, which is adverse selection//! with no informed trader anywhere -- winning means having been the most//! aggressive, and that is a biased place to be.//!//! And the optimal partial hedge, which is where the chapter hands over to risk//! management: the fraction left unhedged has a closed form, and it is not zero. use crate::pathwise::Rng; /// The Avellaneda-Stoikov market making problem.////// A mid price diffuses as `dS = sigma dW`. The maker quotes a bid at `delta_b`/// below the mid and an ask at `delta_a` above it, and fills arrive as Poisson/// processes whose intensity falls off with the distance,////// ```text/// lambda(delta) = A exp(-k delta)./// ```////// Terminal utility is exponential, `-exp(-gamma (cash + inventory * S))`, and/// inventory is capped at `limit` in each direction, which is what makes the/// state space finite.pub struct Quoting { pub sigma: f64, pub gamma: f64, /// Arrival intensity at zero spread. pub intensity: f64, /// How fast arrivals fall off with distance: `k` above. pub decay: f64, pub horizon: f64, pub limit: i32,} impl Quoting { fn inventories(&self) -> usize { (2 * self.limit + 1) as usize } fn index(&self, q: i32) -> usize { (q + self.limit) as usize } /// The linear system the exponential substitution produces. /// /// Writing the value function as /// /// ```text /// u(t,x,q,s) = -exp(-gamma (x + q s)) * v(t,q)^(-gamma/k), /// ``` /// /// the Hamilton-Jacobi-Bellman equation becomes /// /// ```text /// -dv/dt(q) = -alpha q^2 v(q) + eta (v(q-1) + v(q+1)), /// alpha = k gamma sigma^2 / 2, /// eta = A (1 + gamma/k)^(-(1 + k/gamma)), /// ``` /// /// which is linear in `v` --- a tridiagonal system indexed by inventory, with /// no maximisation left in it. Solved here by stepping backwards from /// `v(T,q) = 1`. /// /// Returned as `v[q]` at time zero, indexed by [`Self::index`]. pub fn linear_solution(&self, steps: usize) -> Vec<f64> { let n = self.inventories(); let alpha = self.decay * self.gamma * self.sigma * self.sigma / 2.0; let eta = self.intensity * (1.0 + self.gamma / self.decay).powf(-(1.0 + self.decay / self.gamma)); let mut v = vec![1.0; n]; let dt = self.horizon / steps as f64; for _ in 0..steps { let previous = v.clone(); for i in 0..n { let q = i as i32 - self.limit; let neighbours = (if i > 0 { previous[i - 1] } else { 0.0 }) + (if i + 1 < n { previous[i + 1] } else { 0.0 }); // Backwards in time, so the sign of dt flips. v[i] = previous[i] + dt * (-alpha * (q * q) as f64 * previous[i] + eta * neighbours); } } v } /// The optimal quote distances, from the linear solution. /// /// ```text /// delta_ask(q) = (1/k) ln( v(q) / v(q-1) ) + (1/gamma) ln(1 + gamma/k), /// delta_bid(q) = (1/k) ln( v(q) / v(q+1) ) + (1/gamma) ln(1 + gamma/k). /// ``` /// /// The second term is common to both and is the spread a maker would charge /// with no inventory concern at all. The first is the skew: it is what makes /// a long position quote a tighter ask than bid. /// Returns `(bid distance, ask distance)`, both measured from the mid as /// positive numbers. The order is bid first, and it is worth saying so /// explicitly: a caller that reads it the other way round gets a model whose /// long book shows a worse offer, which is economically backwards and looks /// plausible enough on a chart to survive review. pub fn quotes(&self, q: i32, steps: usize) -> (f64, f64) { let v = self.linear_solution(steps); let base = (1.0 + self.gamma / self.decay).ln() / self.gamma; let i = self.index(q); // Being hit on the bid takes the inventory up, so the bid is priced // against `v` one step higher; being taken on the offer takes it down. let bid = if q < self.limit { (v[i] / v[i + 1]).ln() / self.decay + base } else { f64::INFINITY }; let ask = if q > -self.limit { (v[i] / v[i - 1]).ln() / self.decay + base } else { f64::INFINITY }; (bid, ask) } /// The same problem solved without the substitution, by iterating the /// nonlinear Hamilton-Jacobi-Bellman equation on the value function itself. /// /// Kept because the chapter's structural claim is that the substitution /// *linearises* the problem rather than approximating it, and that claim is /// worth checking against the equation it started from. The maximisation over /// the quote distance is done in closed form --- for `lambda = A exp(-k d)` /// the first order condition gives `exp(-gamma(d + jump)) = k/(k+gamma)`, and /// the value at that optimiser carries a factor `gamma/(k+gamma)` --- so /// what is being compared is the two routes, not two optimisers. /// /// Works with `h(t,q)` defined by `u = -exp(-gamma(x + q s)) exp(gamma h(t,q))`, /// in which the equation reads /// /// ```text /// dh/dt(q) = (gamma sigma^2 q^2 / 2) /// - (1/gamma) sum over sides of lambda(d*) / k , /// ``` /// /// with `d*` the optimiser on each side. pub fn nonlinear_solution(&self, steps: usize) -> Vec<f64> { let n = self.inventories(); let mut h = vec![0.0; n]; let dt = self.horizon / steps as f64; for _ in 0..steps { let previous = h.clone(); for i in 0..n { let q = i as i32 - self.limit; let mut gain = 0.0; // Selling one unit: inventory falls to q - 1. if q > -self.limit { let jump = previous[i - 1] - previous[i]; let d = 1.0 / self.gamma * (1.0 + self.gamma / self.decay).ln() - jump; gain += self.intensity * (-self.decay * d).exp() * self.gamma / (self.decay + self.gamma); } // Buying one unit: inventory rises to q + 1. if q < self.limit { let jump = previous[i + 1] - previous[i]; let d = 1.0 / self.gamma * (1.0 + self.gamma / self.decay).ln() - jump; gain += self.intensity * (-self.decay * d).exp() * self.gamma / (self.decay + self.gamma); } let inventory_cost = 0.5 * self.gamma * self.sigma * self.sigma * (q * q) as f64; h[i] = previous[i] + dt * (-inventory_cost + gain / self.gamma); } } h } /// Quote distances read off the nonlinear solution, for comparison. pub fn quotes_nonlinear(&self, q: i32, steps: usize) -> (f64, f64) { let h = self.nonlinear_solution(steps); let base = (1.0 + self.gamma / self.decay).ln() / self.gamma; let i = self.index(q); let ask = if q > -self.limit { base - (h[i - 1] - h[i]) } else { f64::INFINITY }; let bid = if q < self.limit { base - (h[i + 1] - h[i]) } else { f64::INFINITY }; (bid, ask) } /// Avellaneda and Stoikov's asymptotic total spread, /// /// ```text /// gamma sigma^2 (T - t) + (2/gamma) ln(1 + gamma/k), /// ``` /// /// which is the quoted form and holds when the inventory limit is far away. pub fn asymptotic_spread(&self, time_left: f64) -> f64 { self.gamma * self.sigma * self.sigma * time_left + 2.0 / self.gamma * (1.0 + self.gamma / self.decay).ln() } /// The reservation price offset, `-q gamma sigma^2 (T-t)`: how far below the /// mid a maker holding `q` values the asset. pub fn reservation_offset(&self, q: i32, time_left: f64) -> f64 { -(q as f64) * self.gamma * self.sigma * self.sigma * time_left }} /// The expected cost of winning a request-for-quote auction against `dealers`/// competitors, when every dealer's valuation error is independent and normal/// with standard deviation `error`.////// This is adverse selection with nobody informed. The dealer who wins is the one/// whose error was most aggressive, so conditioning on having won selects a/// biased draw --- and the bias grows with the number of competitors, roughly as/// `error * sqrt(2 ln n)`.////// Returns the mean signed error of the winner, which is negative: the winner has/// underpriced.pub fn winners_curse(dealers: usize, error: f64, trials: usize, seed: u64) -> f64 { let mut rng = Rng::new(seed); let mut total = 0.0; for _ in 0..trials { let mut best = f64::INFINITY; for _ in 0..dealers { best = best.min(error * rng.next_normal()); } total += best; } total / trials as f64} /// The optimal fraction of a position to hedge, when hedging costs and not/// hedging carries variance.////// Hedging a fraction `h` of a unit exposure costs `cost * h` and leaves residual/// variance `variance * (1-h)^2`, penalised at `risk_aversion / 2`. Minimising////// ```text/// cost * h + (risk_aversion / 2) * variance * (1 - h)^2/// ```////// gives `1 - h = cost / (risk_aversion * variance)`, so////// ```text/// h* = max(0, 1 - cost / (risk_aversion * variance))./// ```////// The unhedged fraction is the transaction cost divided by the risk penalty. It/// is zero only when hedging is free, which is the market making chapter's point:/// a partial hedge is a decision with an optimum, not a failure to finish the/// job, and the residual is a position that has to be carried and capitalised.pub fn optimal_hedge_fraction(cost: f64, risk_aversion: f64, variance: f64) -> f64 { let unhedged = cost / (risk_aversion * variance); (1.0 - unhedged).max(0.0)} #[cfg(test)]mod tests { use super::*; fn desk() -> Quoting { Quoting { sigma: 0.30, gamma: 0.1, intensity: 140.0, decay: 1.5, horizon: 1.0, limit: 20, } } const STEPS: usize = 4000; #[test] fn the_substitution_linearises_rather_than_approximates() { // The chapter's structural claim, and the reason this model belongs in a // book about solvability. The Hamilton-Jacobi-Bellman equation has a // maximisation inside it and is nonlinear; the exponential substitution // removes both. // // Testing that against a fixed tolerance would be the wrong test. The two // routes are explicit Euler schemes on different variables -- one on v, // one on log v -- so they disagree at first order in the step whatever // else is true. What distinguishes a linearisation from an approximation // is that the disagreement goes to *zero*: refine the grid and an // approximation leaves a floor while a change of variables does not. // // So the test is the convergence order. The gap should halve each time // the step does, and extrapolate to nothing. let m = desk(); for q in [-8, -3, 0, 2, 7] { let gap = |steps: usize| { let (bid_linear, ask_linear) = m.quotes(q, steps); let (bid_direct, ask_direct) = m.quotes_nonlinear(q, steps); (bid_linear - bid_direct).abs().max((ask_linear - ask_direct).abs()) }; let (coarse, fine, finer) = (gap(4000), gap(8000), gap(16000)); assert!( (coarse / fine - 2.0).abs() < 0.05, "q={q}: halving the step changed the gap by {:.3}, not 2", coarse / fine ); assert!( (fine / finer - 2.0).abs() < 0.05, "q={q}: and again by {:.3}", fine / finer ); // Richardson extrapolation of a first order scheme: 2 * fine - coarse // removes the leading error, and what is left is not a difference of // equations. let extrapolated = (2.0 * finer - fine).abs(); assert!( extrapolated < 1e-6, "q={q}: extrapolated gap {extrapolated:.3e} should vanish" ); } } #[test] fn inventory_skews_the_quotes_and_does_not_widen_them() { // What the model actually says, and the part practitioners take from it. // A long position lowers both quotes rather than widening the pair: the // maker wants to sell, so it shows a better ask and a worse bid, and the // total spread is almost unchanged. let m = desk(); let (bid_flat, ask_flat) = m.quotes(0, STEPS); let (bid_long, ask_long) = m.quotes(6, STEPS); assert!(ask_long < ask_flat, "a long book should show a tighter ask"); assert!(bid_long > bid_flat, "and a worse bid"); let widening = (bid_long + ask_long) - (bid_flat + ask_flat); assert!( widening.abs() < 0.1 * (bid_flat + ask_flat), "the total spread moved by {widening:.5}, which is more than a skew" ); // And the skew is symmetric in the sign of the inventory. let (bid_short, ask_short) = m.quotes(-6, STEPS); assert!((ask_short - bid_long).abs() < 1e-9, "the problem is symmetric in q"); assert!((bid_short - ask_long).abs() < 1e-9); } #[test] fn the_asymptotic_spread_is_the_flat_book_answer() { // Avellaneda and Stoikov's quoted formula is an approximation valid away // from the inventory limits. Checked against the solved system at zero // inventory, where it should be closest. let m = desk(); let (bid, ask) = m.quotes(0, STEPS); let solved = bid + ask; let asymptotic = m.asymptotic_spread(m.horizon); assert!( (solved - asymptotic).abs() < 0.25 * asymptotic, "solved {solved:.5} against asymptotic {asymptotic:.5}" ); // The floor is the term that survives as the horizon shrinks: with no // inventory risk left there is still a spread, set by how price sensitive // the flow is. let floor = 2.0 / m.gamma * (1.0 + m.gamma / m.decay).ln(); assert!(m.asymptotic_spread(0.0) > 0.0); assert!((m.asymptotic_spread(0.0) - floor).abs() < 1e-12); } #[test] fn the_winners_curse_grows_with_the_number_of_dealers() { // Adverse selection with nobody informed. The measured bias should track // the expected minimum of n normals, which grows like sqrt(2 ln n). let error = 0.5; let mut previous = 0.0; for dealers in [2usize, 3, 5, 10] { let bias = winners_curse(dealers, error, 400_000, 20260808); assert!(bias < 0.0, "the winner has underpriced"); assert!(bias < previous, "more competition should cost more"); previous = bias; } // Two dealers is the one case with a clean closed form: the expected // minimum of two standard normals is -1/sqrt(pi). let two = winners_curse(2, 1.0, 400_000, 20260808); let exact = -1.0 / std::f64::consts::PI.sqrt(); assert!((two - exact).abs() < 0.01, "two dealers gave {two:.4} against {exact:.4}"); // And the growth is slow, which is the practical point: going from two // dealers to ten roughly doubles the curse rather than multiplying it by // five. let ten = winners_curse(10, 1.0, 400_000, 20260808); let ratio = ten / two; assert!((1.5..3.0).contains(&ratio), "ten against two is {ratio:.2}"); } #[test] fn the_unhedged_fraction_is_the_cost_over_the_risk_penalty() { // The handover to risk management. Hedging is never complete unless it // is free, and what is left over has a closed form. let (variance, aversion) = (0.04, 2.0); // Free hedging: hedge everything. assert!((optimal_hedge_fraction(0.0, aversion, variance) - 1.0).abs() < 1e-12); // A cost leaves a residual, and the residual is exactly the cost divided // by the risk penalty. for cost in [0.001, 0.01, 0.04] { let h = optimal_hedge_fraction(cost, aversion, variance); let residual = 1.0 - h; assert!( (residual - cost / (aversion * variance)).abs() < 1e-12, "cost={cost}: residual {residual}" ); assert!(h < 1.0, "a costly hedge is never complete"); } // And past a point it is not worth hedging at all: an illiquid risk with // a wide bid-ask is warehoused, which is the case the risk chapter has to // deal with. assert_eq!(optimal_hedge_fraction(0.2, aversion, variance), 0.0); // Verified against a direct minimisation of the objective, which knows // nothing about the formula. let objective = |h: f64| 0.01 * h + 0.5 * aversion * variance * (1.0 - h) * (1.0 - h); let mut best = (f64::INFINITY, 0.0); for i in 0..=100_000 { let h = i as f64 / 100_000.0; let value = objective(h); if value < best.0 { best = (value, h); } } assert!((best.1 - optimal_hedge_fraction(0.01, aversion, variance)).abs() < 1e-4); }}