solver/models/

avellaneda.rs

// src/models/avellaneda.rs
use super::normal_cdf;
use super::traits::{ControlOutput, Gradients, Model};

/// The classic Avellaneda-Stoikov market making model.
///
/// Models a market maker optimizing bid/ask quotes to maximize utility of terminal wealth
/// while penalizing inventory risk.
///
/// # Dynamics
/// * Price: $dS_t = \sigma dW_t$
/// * Inventory: $dq_t = dN^{buy}_t - dN^{sell}_t$
/// * Cash: $dX_t = (S_t + \delta_a) dN^{sell}_t - (S_t - \delta_b) dN^{buy}_t$
#[derive(Clone)]
pub struct AvellanedaStoikov {
    /// Risk aversion parameter ($\gamma$). Controls the penalty for holding inventory.
    pub gamma: f64,
    /// Volatility of the mid-price ($\sigma$).
    pub sigma: f64,
    /// Order filling probability parameter ($\kappa$). Intensity decay rate.
    pub kappa: f64,
    /// Base order arrival intensity ($A$).
    pub a: f64,
}

impl AvellanedaStoikov {
    /// Computes the optimal bid/ask spreads based on current value function gradients.
    ///
    /// # formula
    /// $$ \delta^* = \frac{1}{\gamma} \ln(1 + \frac{\gamma}{\kappa}) \pm \frac{\partial V}{\partial q} $$
    pub fn get_spreads(&self, grads: &Gradients<2>) -> (f64, f64) {
        let base_spread = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();

        let dv_dq_buy = grads.fwd[0];
        let dv_dq_sell = grads.bwd[0];

        let mut delta_bid = base_spread - dv_dq_buy;
        let mut delta_ask = base_spread + dv_dq_sell;

        // Safety cap to prevent e^50 explosions
        let min_spread = -10.0;
        delta_bid = delta_bid.max(min_spread);
        delta_ask = delta_ask.max(min_spread);

        (delta_bid, delta_ask)
    }
}

impl Model<2> for AvellanedaStoikov {
    type Process = ();

    fn process(&self) {}

    fn optimize(&self, state: &[f64; 2], grads: &Gradients<2>) -> ControlOutput<2> {
        let q = state[0];
        // state[1] is Price S, but AS model drift/intensity depends only on inventory q.

        let (d_bid, d_ask) = self.get_spreads(grads);

        // Compute intensities based on spreads: Lambda = A * e^{-k * delta}
        let lambda_bid = self.a * (-self.kappa * d_bid).exp();
        let lambda_ask = self.a * (-self.kappa * d_ask).exp();

        // Hamiltonian for Exponential Utility (Guéant-Lehalle-Tapia transformation):
        // The standard HJB for the transformed variable \theta(t,q) involves:
        // H = sup_d { \lambda(d) [1 - e^{-\gamma(d - \Delta \theta)}] } / \gamma
        // With optimal spreads, this simplifies to term below:
        let hamiltonian_val = (lambda_bid + lambda_ask) / (self.gamma + self.kappa);

        // Running penalty term (risk aversion * volatility)
        let risk_penalty = -0.5 * self.gamma * self.sigma.powi(2) * q.powi(2);

        // Term Cancellation Strategy:
        // The analytical control (Gueant-Lehalle-Tapia) implicitly includes the jump cost terms
        // (lambda * (V(q') - V(q))) due to the utility transformation.
        // However, the BSDE solver explicitly adds these drift/jump terms via its driver update.
        // To avoid double counting, we must subtract the explicit jump terms here.
        let drift_correction = lambda_bid * grads.fwd[0] - lambda_ask * grads.bwd[0];

        // Dim 0 is Inventory.

        // lambda_plus[0] = Buy Intensity (q -> q+1)
        // lambda_minus[0] = Sell Intensity (q -> q-1)
        // Dim 1 is Price (No jumps relevant to value function in reduced coords, but kept for interface)
        ControlOutput {
            lambda_plus: [lambda_bid, 0.0],
            lambda_minus: [lambda_ask, 0.0],
            flow: hamiltonian_val + risk_penalty - drift_correction,
        }
    }

    fn terminal(&self, _state: &[f64; 2]) -> f64 {
        0.0
    }

    fn constant_discount_rate(&self) -> Option<f64> {
        Some(0.0)
    }

    fn next_step(&self, current_state: &[f64; 2], dt: f64, noise: &[f64; 2]) -> [f64; 2] {
        let mut next = *current_state;
        // q (index 0) does not diffuse
        // S (index 1) diffuses with volatility sigma
        // noise is N(0, 1)

        let sqrt_dt = dt.sqrt();
        next[1] += self.sigma * sqrt_dt * noise[1];

        next
    }

    fn next_step_controlled(
        &self,
        current_state: &[f64; 2],
        control: &ControlOutput<2>,
        dt: f64,
        noise: &[f64; 2],
    ) -> [f64; 2] {
        let mut next = *current_state;
        // Dim 0 (q): Bernoulli fill simulation.
        let u = normal_cdf(noise[0]);
        let p_bid = (control.lambda_plus[0] * dt).clamp(0.0, 1.0);
        let p_ask = (control.lambda_minus[0] * dt).clamp(0.0, 1.0);
        if u < p_bid {
            next[0] += 1.0;
        } else if u > 1.0 - p_ask {
            next[0] -= 1.0;
        }
        // Dim 1 (S): diffusive price.
        let sqrt_dt = dt.sqrt();
        next[1] += self.sigma * sqrt_dt * noise[1];
        next
    }

    fn is_diffusion_dimension(&self, dim: usize) -> bool {
        dim == 1
    }

    fn is_integer_dimension(&self, dim: usize) -> bool {
        dim == 0 // q is discrete integer; S is continuous
    }

    fn fill_rate_base(&self, _state: &[f64; 2]) -> f64 {
        self.a
    }

    fn fill_rate_decay(&self) -> f64 {
        self.kappa
    }
}

#[cfg(test)]
mod tests {
    use super::super::traits::{Gradients, Model};
    use super::*;

    fn zero_grads() -> Gradients<2> {
        Gradients {
            fwd: [0.0; 2],
            bwd: [0.0; 2],
        }
    }

    fn default_model() -> AvellanedaStoikov {
        AvellanedaStoikov {
            gamma: 0.1,
            sigma: 0.1,
            kappa: 1.5,
            a: 140.0,
        }
    }

    #[test]
    fn zero_gradient_gives_base_spread() {
        let m = default_model();
        let grads = zero_grads();
        let (bid, ask) = m.get_spreads(&grads);
        let base = (1.0 / m.gamma) * (1.0 + m.gamma / m.kappa).ln();
        assert!((bid - base).abs() < 1e-12);
        assert!((ask - base).abs() < 1e-12);
    }

    #[test]
    fn symmetric_at_zero_inventory() {
        let m = default_model();
        let grads = zero_grads();
        let ctrl = m.optimize(&[0.0, 100.0], &grads);
        assert!(
            (ctrl.lambda_plus[0] - ctrl.lambda_minus[0]).abs() < 1e-12,
            "Bid/ask intensities should be equal at q=0"
        );
    }

    #[test]
    fn positive_flow_at_zero_inventory() {
        let m = default_model();
        let grads = zero_grads();
        let ctrl = m.optimize(&[0.0, 100.0], &grads);
        assert!(
            ctrl.flow > 0.0,
            "Flow should be positive at q=0: hamiltonian > risk_penalty"
        );
    }

    #[test]
    fn risk_penalty_increases_with_inventory() {
        let m = default_model();
        let grads = zero_grads();
        let ctrl_0 = m.optimize(&[0.0, 100.0], &grads);
        let ctrl_5 = m.optimize(&[5.0, 100.0], &grads);
        assert!(
            ctrl_5.flow < ctrl_0.flow,
            "Higher inventory should decrease flow due to risk penalty"
        );
    }

    #[test]
    fn terminal_is_zero() {
        let m = default_model();
        assert_eq!(m.terminal(&[0.0, 100.0]), 0.0);
        assert_eq!(m.terminal(&[5.0, 200.0]), 0.0);
    }

    #[test]
    fn next_step_preserves_inventory() {
        let m = default_model();
        let state = [3.0, 100.0];
        let next = m.next_step(&state, 0.01, &[0.5, 1.0]);
        assert_eq!(
            next[0], 3.0,
            "Inventory should not change in diffusion step"
        );
    }

    #[test]
    fn next_step_diffuses_price() {
        let m = default_model();
        let state = [0.0, 100.0];
        let next = m.next_step(&state, 0.01, &[0.0, 1.0]);
        assert!(
            (next[1] - 100.0).abs() > 1e-10,
            "Price should diffuse with non-zero noise"
        );
    }

    #[test]
    fn dim1_is_diffusion() {
        let m = default_model();
        assert!(!m.is_diffusion_dimension(0));
        assert!(m.is_diffusion_dimension(1));
    }
}