solver/models/

avellaneda.rs

// src/models/avellaneda.rs
use super::control::{ControlProblem, StateDerivatives};
use super::market_making::MarketMakingControl;
use super::normal_cdf;
use crate::numeric::finite_difference::discretization::{DimensionKind, Transport};
use crate::numeric::finite_difference::elliptic::EllipticControlProblem;
use crate::numeric::finite_difference::pde::PdeProblem;

/// 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, derivs: &StateDerivatives<2>) -> (f64, f64) {
        let base_spread = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();

        let dv_dq_buy = derivs.fwd[0];
        let dv_dq_sell = derivs.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)
    }

    /// The base arrival intensity used for intensity-to-spread conversion.
    pub fn fill_rate_base(&self, _state: &[f64; 2]) -> f64 {
        self.a
    }

    /// The fill-rate decay parameter used for intensity-to-spread conversion.
    pub fn fill_rate_decay(&self) -> f64 {
        self.kappa
    }
}

impl ControlProblem<2> for AvellanedaStoikov {
    type Control = MarketMakingControl;

    fn optimize(&self, _t: f64, _state: &[f64; 2], derivs: &StateDerivatives<2>) -> Self::Control {
        let base_spread = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();
        let delta_bid = (base_spread - derivs.fwd[0]).max(-10.0);
        let delta_ask = (base_spread + derivs.bwd[0]).max(-10.0);

        MarketMakingControl::new(
            self.a * (-self.kappa * delta_bid).exp(),
            self.a * (-self.kappa * delta_ask).exp(),
        )
    }

    fn running_reward(&self, _t: f64, _state: &[f64; 2], control: &Self::Control) -> f64 {
        // Optimized Hamiltonian H* = (lambda_bid + lambda_ask) / (gamma + kappa).
        (control.bid_intensity + control.ask_intensity) / (self.gamma + self.kappa)
    }

    fn bsde_driver(
        &self,
        t: f64,
        state: &[f64; 2],
        control: &Self::Control,
        derivs: &StateDerivatives<2>,
        _dt: f64,
    ) -> f64 {
        // The reduced theta value satisfies d_t theta + H* - 0.5 gamma sigma^2 q^2 = 0.
        // The inventory-risk term is a local source, not forward transport, so
        // the BSDE backward step must include the full reduced driver.
        self.driver(t, state, control, derivs)
    }

    fn generator(
        &self,
        _t: f64,
        state: &[f64; 2],
        _control: &Self::Control,
        _derivs: &StateDerivatives<2>,
    ) -> f64 {
        // Price-diffusion risk penalty: -(1/2) gamma sigma^2 q^2.
        -0.5 * self.gamma * self.sigma.powi(2) * state[0].powi(2)
    }

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

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

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

    fn next_step(&self, _t: f64, state: &[f64; 2], dt: f64, noise: &[f64; 2]) -> [f64; 2] {
        let mut next = *state;

        // Inventory q jumps by +/-1 based on the optimal fill probabilities.
        let u = normal_cdf(noise[0]);
        let base_spread = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();
        let delta_bid = base_spread;
        let delta_ask = base_spread;
        let lambda_bid = self.a * (-self.kappa * delta_bid).exp();
        let lambda_ask = self.a * (-self.kappa * delta_ask).exp();
        let p_bid = (lambda_bid * dt).clamp(0.0, 1.0);
        let p_ask = (lambda_ask * 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;
        }

        // Price diffuses.
        next[1] += self.sigma * dt.sqrt() * noise[1];
        next
    }

    fn is_reduced_value(&self) -> bool {
        true
    }

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

    fn gradient_step(&self, _dim: usize) -> f64 {
        1.0
    }
}

impl PdeProblem<2> for AvellanedaStoikov {
    fn dimension_kind(&self, dim: usize) -> DimensionKind {
        match dim {
            0 => DimensionKind::DiscreteJump,
            _ => DimensionKind::Diffusion,
        }
    }

    fn transport(
        &self,
        _t: f64,
        state: &[f64; 2],
        control: &Self::Control,
        derivs: &StateDerivatives<2>,
    ) -> Transport<2> {
        // The solver adds inventory transport for the two optimal fill rates.
        // The residual source is therefore the full optimized driver minus the
        // transport contribution that the operator will produce.
        let q = state[0];
        let hamiltonian =
            (control.bid_intensity + control.ask_intensity) / (self.gamma + self.kappa);
        let jump_transport =
            control.bid_intensity * derivs.fwd[0] - control.ask_intensity * derivs.bwd[0];
        let local = -0.5 * self.gamma * self.sigma.powi(2) * q.powi(2);
        Transport::new(
            [control.bid_intensity, 0.0],
            [control.ask_intensity, 0.0],
            hamiltonian + local - jump_transport,
        )
    }
}

impl EllipticControlProblem<2> for AvellanedaStoikov {
    fn dimension_kind(&self, dim: usize) -> DimensionKind {
        <Self as PdeProblem<2>>::dimension_kind(self, dim)
    }

    fn transport(
        &self,
        state: &[f64; 2],
        control: &Self::Control,
        derivs: &StateDerivatives<2>,
    ) -> Transport<2> {
        <Self as PdeProblem<2>>::transport(self, 0.0, state, control, derivs)
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    fn zero_derivs() -> StateDerivatives<2> {
        StateDerivatives::new([0.0; 2], [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 derivs = zero_derivs();
        let (bid, ask) = m.get_spreads(&derivs);
        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 derivs = zero_derivs();
        let ctrl = ControlProblem::optimize(&m, 0.0, &[0.0, 100.0], &derivs);
        assert!(
            (ctrl.bid_intensity - ctrl.ask_intensity).abs() < 1e-12,
            "Bid/ask intensities should be equal at q=0"
        );
    }

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

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

    #[test]
    fn dimension_kinds_match_state_space() {
        let m = default_model();
        assert_eq!(
            PdeProblem::dimension_kind(&m, 0),
            DimensionKind::DiscreteJump
        );
        assert_eq!(PdeProblem::dimension_kind(&m, 1), DimensionKind::Diffusion);
    }

    #[test]
    fn transport_source_reconstructs_full_driver() {
        let m = default_model();
        let state = [3.0, 100.0];
        let derivs =
            StateDerivatives::with_directional([0.1, 0.0], [0.0, 0.0], [0.25, 0.0], [0.15, 0.0]);

        let control = ControlProblem::optimize(&m, 0.0, &state, &derivs);
        let transport = PdeProblem::transport(&m, 0.0, &state, &control, &derivs);

        // The discrete operator applies
        //   plus * fwd_delta + minus * bwd_delta + source
        // to reconstruct the full optimized driver. With a unit jump size the
        // forward/bwd directional derivatives are already the discrete deltas.
        let operator_driver = transport.source
            + transport.plus[0] * derivs.fwd[0]
            + transport.minus[0] * (-derivs.bwd[0]);

        let expected = ControlProblem::driver(&m, 0.0, &state, &control, &derivs);
        assert!(
            (operator_driver - expected).abs() < 1e-12,
            "transport decomposition {operator_driver} does not match driver {expected}"
        );
    }
}