solver/models/
avellaneda_impact.rs

1use super::traits::{ControlOutput, Gradients, Model};
2
3/// Avellaneda-Stoikov Model with Market Impact
4///
5/// Models a market maker whose trades permanently impact the mid-price.
6/// The inventory penalty is adjusted implicitly by the cost of moving the price against oneself.
7///
8/// # Dynamics
9/// * $S_t \to S_t + \xi$ (Buy) or $S_t - \xi$ (Sell)
10/// * Impact is linear and permanent.
11#[derive(Clone)]
12pub struct AvellanedaImpact {
13    pub gamma: f64,
14    pub sigma: f64,
15    pub kappa: f64,
16    pub a: f64,
17    /// Permanent market impact parameter ($\xi$).
18    /// A trade of size 1 moves the price by $\xi$.
19    pub xi: f64,
20}
21
22impl AvellanedaImpact {
23    /// Computes optimal spreads, adjusted for market impact.
24    pub fn get_spreads(&self, grads: &Gradients<2>, q: f64) -> (f64, f64) {
25        let base_spread = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();
26
27        let dv_dq_buy = grads.fwd[0];
28        let dv_dq_sell = grads.bwd[0];
29
30        let mut delta_bid = base_spread - dv_dq_buy + (q + 1.0) * self.xi;
31        let mut delta_ask = base_spread + dv_dq_sell - (q - 1.0) * self.xi;
32
33        let min_spread = -5.0;
34        delta_bid = delta_bid.max(min_spread);
35        delta_ask = delta_ask.max(min_spread);
36
37        (delta_bid, delta_ask)
38    }
39}
40
41impl Model<2> for AvellanedaImpact {
42    type Process = ();
43
44    fn process(&self) {}
45
46    fn optimize(&self, state: &[f64; 2], grads: &Gradients<2>) -> ControlOutput<2> {
47        let q = state[0];
48
49        let (d_bid, d_ask) = self.get_spreads(grads, q);
50
51        let lambda_bid = self.a * (-self.kappa * d_bid).exp();
52        let lambda_ask = self.a * (-self.kappa * d_ask).exp();
53
54        // Use the CARA HJB Hamiltonian value: lambda / (gamma + kappa)
55        // This implicitly accounts for the non-linear utility of the spread and impact.
56        let hamiltonian_val = (lambda_bid + lambda_ask) / (self.gamma + self.kappa);
57
58        // Correction for the linear jump term added by the solver
59        // Solver adds: lambda * (theta_new - theta_old)
60        // We want: H(theta)
61        // So flow = H - linear_jump
62        let drift_correction = lambda_bid * grads.fwd[0] - lambda_ask * grads.bwd[0];
63
64        let risk_penalty = -0.5 * self.gamma * self.sigma.powi(2) * q.powi(2);
65
66        ControlOutput {
67            lambda_plus: [lambda_bid, 0.0],
68            lambda_minus: [lambda_ask, 0.0],
69            flow: hamiltonian_val + risk_penalty - drift_correction,
70        }
71    }
72
73    fn terminal(&self, state: &[f64; 2]) -> f64 {
74        // Terminal penalty for liquidation with impact
75        // Usually T condition is just liquidation cost.
76        // v(T, q) = qS - alpha q^2
77        // Here we assume standard liquidation or 0 if we just compare spreads far from T.
78        // The exact solution uses v(T) = exp(-0.5 k xi q^2).
79        // In Mean-Variance, this corresponds to -0.5 * xi * q^2?
80        // Let's add it.
81        let q = state[0];
82        -0.5 * self.xi * q.powi(2)
83    }
84
85    fn constant_discount_rate(&self) -> Option<f64> {
86        Some(0.0)
87    }
88
89    fn next_step(&self, current_state: &[f64; 2], dt: f64, noise: &[f64; 2]) -> [f64; 2] {
90        let mut next = *current_state;
91        let sqrt_dt = dt.sqrt();
92        // q (0) constant
93        // S (1) diffuses
94        next[1] += self.sigma * sqrt_dt * noise[1];
95        next
96    }
97
98    fn fill_rate_base(&self, _state: &[f64; 2]) -> f64 {
99        self.a
100    }
101
102    fn fill_rate_decay(&self) -> f64 {
103        self.kappa
104    }
105}
solver/models/avellaneda_impact.rs

solver/models/
avellaneda_impact.rs
