solver/models/

avellaneda_impact.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::pde::PdeProblem;

/// Avellaneda-Stoikov Model with Market Impact
///
/// Models a market maker whose trades permanently impact the mid-price.
/// The inventory penalty is adjusted implicitly by the cost of moving the price against oneself.
///
/// # Dynamics
/// * $S_t \to S_t + \xi$ (Buy) or $S_t - \xi$ (Sell)
/// * Impact is linear and permanent.
#[derive(Clone)]
pub struct AvellanedaImpact {
    pub gamma: f64,
    pub sigma: f64,
    pub kappa: f64,
    pub a: f64,
    /// Permanent market impact parameter ($\xi$).
    /// A trade of size 1 moves the price by $\xi$.
    pub xi: f64,
}

impl AvellanedaImpact {
    /// Computes optimal spreads, adjusted for market impact.
    pub fn get_spreads(&self, derivs: &StateDerivatives<2>, q: f64) -> (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 + (q + 1.0) * self.xi;
        let mut delta_ask = base_spread + dv_dq_sell - (q - 1.0) * self.xi;

        let min_spread = -5.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 AvellanedaImpact {
    type Control = MarketMakingControl;

    fn optimize(&self, _t: f64, state: &[f64; 2], derivs: &StateDerivatives<2>) -> Self::Control {
        let q = state[0];
        let base_spread = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();
        let delta_bid = (base_spread - derivs.fwd[0] + (q + 1.0) * self.xi).max(-5.0);
        let delta_ask = (base_spread + derivs.bwd[0] - (q - 1.0) * self.xi).max(-5.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 {
        (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 {
        self.driver(t, state, control, derivs)
    }

    fn generator(
        &self,
        _t: f64,
        state: &[f64; 2],
        _control: &Self::Control,
        _derivs: &StateDerivatives<2>,
    ) -> f64 {
        let q = state[0];
        -0.5 * self.gamma * self.sigma.powi(2) * q.powi(2)
    }

    fn terminal(&self, state: &[f64; 2]) -> f64 {
        let q = state[0];
        -0.5 * self.xi * q.powi(2)
    }

    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;

        let q = state[0];
        let u = normal_cdf(noise[0]);
        let base_spread = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();
        // Forward proxy for the optimal spreads without the value gradient:
        // the permanent-impact shift is retained so the simulated fill
        // intensities reflect the impact-adjusted quotes.
        let delta_bid = base_spread + (q + 1.0) * self.xi;
        let delta_ask = base_spread - (q - 1.0) * self.xi;
        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;
        }

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

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

        // The optimal control already carries the impact-adjusted fill
        // intensities, so the Bernoulli fill events use them directly instead
        // of a frozen proxy quote.
        let u = normal_cdf(noise[0]);
        let p_bid = (control.bid_intensity * dt).clamp(0.0, 1.0);
        let p_ask = (control.ask_intensity * 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;
        }

        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 AvellanedaImpact {
    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> {
        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,
        )
    }
}

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

    fn model() -> AvellanedaImpact {
        AvellanedaImpact {
            gamma: 0.5,
            sigma: 0.5,
            kappa: 1.5,
            a: 140.0,
            xi: 0.01,
        }
    }

    #[test]
    fn controlled_step_uses_optimal_intensities() {
        let m = model();
        let state = [1.0, 100.0];

        // Zero intensity: no fill can occur, so inventory is unchanged and only
        // the price diffuses.
        let ctrl = MarketMakingControl::new(0.0, 0.0);
        let next = ControlProblem::next_step_controlled(&m, 0.0, &state, &ctrl, 0.01, &[0.0, 0.0]);
        assert_eq!(next[0], 1.0);

        // Saturating intensity with a zero cdf draw forces a bid fill
        // (inventory +1).
        let ctrl = MarketMakingControl::new(1000.0, 0.0);
        let next = ControlProblem::next_step_controlled(&m, 0.0, &state, &ctrl, 0.01, &[0.0, 0.0]);
        assert_eq!(next[0], 2.0);
    }

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