solver/models/

avellaneda_lot_size.rs

//! Avellaneda-Stoikov market making with a configurable lot size.
//!
//! A validating instance for [`DimensionKind::Jump`] and [`JumpKernel`]. The
//! inventory is an integer lattice with spacing `lot_size`, so a single fill
//! moves inventory by `+lot_size` (bid fill) or `-lot_size` (ask fill) rather
//! than one unit. With `lot_size = 1` the model reduces exactly to the base
//! Avellaneda-Stoikov reduced value problem, providing a closed-form oracle for
//! the general-jump machinery.
//!
//! # Dynamics
//!
//! The inventory `q` is measured in units of the lot, so the jump amplitude is
//! `lot_size` state units per fill. The reduced value obeys
//!
//! ```text
//! d_t theta + H*(q) - 0.5 gamma sigma^2 q^2 = 0
//! ```
//!
//! where `q` ranges over the integer lattice `lot_size * Z` and the jump
//! operator is `lambda_bid (theta(q+lot) - theta(q)) +
//! lambda_ask (theta(q-lot) - theta(q))`.

use crate::models::control::{ControlProblem, StateDerivatives};
use crate::models::market_making::MarketMakingControl;
use crate::models::normal_cdf;
use crate::numeric::finite_difference::discretization::{
    DimensionKind, JumpKernel, JumpTransition, Transport,
};
use crate::numeric::finite_difference::pde::PdeProblem;

/// Avellaneda-Stoikov with lot size `lot_size` on the inventory dimension.
#[derive(Clone, Copy, Debug)]
pub struct AvellanedaLotSize {
    /// Risk aversion.
    pub gamma: f64,
    /// Mid-price volatility.
    pub sigma: f64,
    /// Fill-intensity decay.
    pub kappa: f64,
    /// Base order arrival intensity.
    pub a: f64,
    /// Number of inventory units moved per fill (positive integer).
    pub lot_size: i32,
}

impl AvellanedaLotSize {
    /// Creates the model.
    ///
    /// # Panics
    ///
    /// Panics if `lot_size` is not positive.
    ///
    /// # Examples
    ///
    /// ```
    /// use solver::models::avellaneda_lot_size::AvellanedaLotSize;
    /// let model = AvellanedaLotSize::new(0.5, 0.5, 1.5, 140.0, 2);
    /// assert_eq!(model.lot_size, 2);
    /// ```
    pub fn new(gamma: f64, sigma: f64, kappa: f64, a: f64, lot_size: i32) -> Self {
        assert!(lot_size > 0, "lot size must be positive");
        Self {
            gamma,
            sigma,
            kappa,
            a,
            lot_size,
        }
    }

    fn base_spread(&self) -> f64 {
        (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln()
    }
}

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

    fn optimize(&self, _t: f64, _state: &[f64; 2], derivs: &StateDerivatives<2>) -> Self::Control {
        let base = self.base_spread();
        let delta_bid = (base - derivs.fwd[0]).max(-10.0);
        let delta_ask = (base + 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 {
        (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 {
        -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;
        let u = normal_cdf(noise[0]);
        let base = self.base_spread();
        let lambda_bid = self.a * (-self.kappa * base).exp();
        let lambda_ask = self.a * (-self.kappa * base).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] += self.lot_size as f64;
        } else if u > 1.0 - p_ask {
            next[0] -= self.lot_size as f64;
        }
        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;
        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] += self.lot_size as f64;
        } else if u > 1.0 - p_ask {
            next[0] -= self.lot_size as f64;
        }
        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 {
        if dim == 0 { self.lot_size as f64 } else { 1.0 }
    }
}

impl PdeProblem<2> for AvellanedaLotSize {
    fn dimension_kind(&self, dim: usize) -> DimensionKind {
        match dim {
            0 => DimensionKind::Jump,
            _ => 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,
        )
    }

    fn jump_kernel(
        &self,
        _t: f64,
        _state: &[f64; 2],
        control: &Self::Control,
        _derivs: &StateDerivatives<2>,
        _transport: &Transport<2>,
    ) -> JumpKernel<2> {
        let mut kernel = JumpKernel::empty();
        if control.bid_intensity > 0.0 {
            kernel.dims[0].push(JumpTransition::new(self.lot_size, control.bid_intensity));
        }
        if control.ask_intensity > 0.0 {
            kernel.dims[0].push(JumpTransition::new(-self.lot_size, control.ask_intensity));
        }
        kernel
    }
}

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

    #[test]
    fn lot_size_one_reduces_to_unit_jump_kernel() {
        let m = AvellanedaLotSize::new(0.5, 0.5, 1.5, 140.0, 1);
        let derivs = StateDerivatives::new([0.0; 2], [0.0; 2]);
        let state = [0.0, 0.0];
        let control = ControlProblem::optimize(&m, 0.0, &state, &derivs);
        let kernel =
            PdeProblem::jump_kernel(&m, 0.0, &state, &control, &derivs, &Transport::zero());
        assert_eq!(kernel.dims[0][0].amplitude, 1);
        assert_eq!(kernel.dims[0][1].amplitude, -1);
    }

    #[test]
    fn lot_size_two_jumps_by_two() {
        let m = AvellanedaLotSize::new(0.5, 0.5, 1.5, 140.0, 2);
        let derivs = StateDerivatives::new([0.0; 2], [0.0; 2]);
        let state = [0.0, 0.0];
        let control = ControlProblem::optimize(&m, 0.0, &state, &derivs);
        let kernel =
            PdeProblem::jump_kernel(&m, 0.0, &state, &control, &derivs, &Transport::zero());
        assert_eq!(kernel.dims[0][0].amplitude, 2);
        assert_eq!(kernel.dims[0][1].amplitude, -2);
    }

    #[test]
    fn next_step_controlled_jumps_by_lot_size() {
        let m = AvellanedaLotSize::new(0.5, 0.5, 1.5, 140.0, 3);
        let state = [0.0, 0.0];
        let control = MarketMakingControl::new(1000.0, 0.0);
        let next =
            ControlProblem::next_step_controlled(&m, 0.0, &state, &control, 0.01, &[0.0, 0.0]);
        assert_eq!(next[0], 3.0);
    }
}