solver/models/

avellaneda_hawkes.rs

use super::TerminalCondition;
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 Hawkes Process for order arrival intensity.
///
/// State Space: [Inventory (q), Intensity (lambda)]
/// - q: Discrete inventory level
/// - lambda: Continuous order arrival intensity
///
/// Dynamics:
/// - q: Jumps ±1 upon fills.
/// - lambda: Mean-reverting jump-diffusion (approximated).
///   d(lambda) = beta * (mu - lambda) * dt + alpha * dN_market
///   where N_market is the market order counting process with intensity lambda.
///
/// The solver uses a grid-based approach. We approximate the Hawkes dynamics
/// on the lambda dimension using an upwind finite difference scheme mapped to
/// the solver's jump interface.
#[derive(Clone)]
pub struct AvellanedaHawkes {
    pub gamma: f64,
    pub sigma: f64,
    pub kappa: f64,
    // Hawkes Parameters
    pub alpha: f64, // Jump size of intensity
    pub beta: f64,  // Decay rate of intensity
    pub mu: f64,    // Baseline intensity level (mean reversion target)
    pub lambda_step: f64,
    /// Grid step size for inventory dimension (set from `grid.dx[0]` for FD).
    pub dq: f64,
    /// Terminal condition V(T, q, lambda).  Defaults to Zero.
    pub terminal_condition: TerminalCondition,
    /// Minimum inventory (lower hard boundary). Defaults to -infinity (no constraint).
    pub q_min: f64,
    /// Maximum inventory (upper hard boundary). Defaults to +infinity (no constraint).
    pub q_max: f64,
    /// Optional terminal liquidation half-spread used when terminal_condition is
    /// LiquidationCost. If None, uses the AS base spread.
    pub terminal_liquidation_half_spread: Option<f64>,
}

impl AvellanedaHawkes {
    pub fn new(gamma: f64, sigma: f64, kappa: f64, alpha: f64, beta: f64, mu: f64) -> Self {
        Self {
            gamma,
            sigma,
            kappa,
            alpha,
            beta,
            mu,
            lambda_step: 1.0,
            dq: 1.0,
            terminal_condition: TerminalCondition::Zero,
            q_min: f64::NEG_INFINITY,
            q_max: f64::INFINITY,
            terminal_liquidation_half_spread: None,
        }
    }

    pub fn with_terminal_liquidation_half_spread(mut self, half_spread: f64) -> Self {
        self.terminal_liquidation_half_spread = Some(half_spread.max(0.0));
        self
    }

    pub fn with_inventory_bounds(mut self, q_min: f64, q_max: f64) -> Self {
        self.q_min = q_min;
        self.q_max = q_max;
        self
    }

    pub fn with_terminal_condition(mut self, terminal_condition: TerminalCondition) -> Self {
        self.terminal_condition = terminal_condition;
        self
    }

    pub fn with_lambda_step(mut self, lambda_step: f64) -> Self {
        self.lambda_step = lambda_step.abs().max(1e-8);
        self
    }

    pub fn with_dq(mut self, dq: f64) -> Self {
        self.dq = dq.abs().max(1e-8);
        self
    }

    /// Computes optimal spreads based on current inventory and intensity.
    /// Note: The base intensity 'A' from standard Avellaneda is replaced by the state variable 'lambda'.
    pub fn get_spreads(&self, _lambda: f64, 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;

        let min_spread = -10.0;
        let max_spread = 10.0;
        delta_bid = delta_bid.max(min_spread).min(max_spread);
        delta_ask = delta_ask.max(min_spread).min(max_spread);

        (delta_bid, delta_ask)
    }

    /// The base arrival intensity is state-dependent: returns the current
    /// intensity `state[1]` (clamped above 0) for intensity-to-spread conversion.
    pub fn fill_rate_base(&self, state: &[f64; 2]) -> f64 {
        state[1].max(1e-10)
    }

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

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

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

        let max_fill_mult = 20.0_f64;
        let lambda_bid = lambda * (-self.kappa * delta_bid).exp().min(max_fill_mult);
        let lambda_ask = lambda * (-self.kappa * delta_ask).exp().min(max_fill_mult);

        let lambda_bid = if q >= self.q_max { 0.0 } else { lambda_bid };
        let lambda_ask = if q <= self.q_min { 0.0 } else { lambda_ask };

        MarketMakingControl::new(lambda_bid, lambda_ask)
    }

    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 {
        // Full-value problem: the Hawkes intensity drift transport is already
        // simulated forward, so only the running reward plus the local
        // inventory-risk source belong in the backward driver. The reward is
        // bounded at `1/dt` per fill side to match the forward fill-probability
        // clamp `lambda * dt <= 1`.
        let q = state[0];
        let local = -0.5 * self.gamma * self.sigma.powi(2) * q.powi(2);
        let rate_cap = 1.0 / dt.max(1e-12);
        let reward = (control.bid_intensity.min(rate_cap) + control.ask_intensity.min(rate_cap))
            / (self.gamma + self.kappa);
        reward + local
    }

    fn generator(
        &self,
        _t: f64,
        state: &[f64; 2],
        _control: &Self::Control,
        derivs: &StateDerivatives<2>,
    ) -> f64 {
        let q = state[0];
        let lambda = state[1];
        let net_drift_lambda = self.beta * (self.mu - lambda) + self.alpha * lambda;

        // Risk penalty plus Hawkes intensity transport, upwinded on the
        // intensity dimension.
        -0.5 * self.gamma * self.sigma.powi(2) * q.powi(2) + net_drift_lambda * derivs.grad[1]
    }

    fn terminal(&self, state: &[f64; 2]) -> f64 {
        match self.terminal_condition {
            TerminalCondition::Zero => 0.0,
            TerminalCondition::LiquidationCost => {
                let q = state[0];
                let base_spread = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();
                let half_spread = self.terminal_liquidation_half_spread.unwrap_or(base_spread);
                -q.abs() * half_spread
            }
        }
    }

    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 lambda = state[1].max(0.0);
        let base_spread = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();
        let lambda_bid = lambda * (-self.kappa * base_spread).exp();
        let lambda_ask = lambda * (-self.kappa * base_spread).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;
        }

        let net_drift_lambda = self.beta * (self.mu - lambda) + self.alpha * lambda;
        next[1] += net_drift_lambda * dt;
        next[1] = next[1].max(0.0);
        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 fill intensities, so the
        // Bernoulli fill events use them 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;
        }

        let lambda = state[1].max(0.0);
        let net_drift_lambda = self.beta * (self.mu - lambda) + self.alpha * lambda;
        next[1] += net_drift_lambda * dt;
        next[1] = next[1].max(0.0);
        next
    }

    fn is_diffusion_dimension(&self, _dim: usize) -> bool {
        false
    }

    fn gradient_step(&self, dim: usize) -> f64 {
        if dim == 1 { self.lambda_step } else { 1.0 }
    }
}

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

    fn transport(
        &self,
        _t: f64,
        state: &[f64; 2],
        control: &Self::Control,
        derivs: &StateDerivatives<2>,
    ) -> Transport<2> {
        let q = state[0];
        let lambda = state[1].max(0.0);
        let net_drift = self.beta * (self.mu - lambda) + self.alpha * lambda;
        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);

        let rate = net_drift.abs() / self.lambda_step.abs().max(1e-12);
        let (plus_lam, minus_lam) = if net_drift >= 0.0 {
            (rate, 0.0)
        } else {
            (0.0, rate)
        };

        Transport::new(
            [control.bid_intensity, plus_lam],
            [control.ask_intensity, minus_lam],
            hamiltonian + local - jump_transport,
        )
    }
}

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

    fn zero_derivs() -> StateDerivatives<2> {
        StateDerivatives::new([0.0; 2], [0.0; 2])
    }

    fn default_model() -> AvellanedaHawkes {
        AvellanedaHawkes::new(0.5, 0.5, 1.5, 0.5, 2.0, 1.0)
            .with_lambda_step(0.1)
            .with_dq(1.0)
    }

    #[test]
    fn zero_gradient_gives_base_spread() {
        let m = default_model();
        let derivs = zero_derivs();
        let (bid, ask) = m.get_spreads(1.0, &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 bsde_driver_excludes_intensity_transport() {
        let m = default_model();
        let state = [2.0, 1.0];
        let control = ControlProblem::optimize(&m, 0.0, &state, &zero_derivs());

        let mut derivs = zero_derivs();
        derivs.grad[1] = 1.0;

        let dt = 0.01;
        let q = state[0];
        let local = -0.5 * m.gamma * m.sigma.powi(2) * q.powi(2);
        let rate_cap = 1.0 / dt;
        let reward = (control.bid_intensity.min(rate_cap) + control.ask_intensity.min(rate_cap))
            / (m.gamma + m.kappa);
        let expected = reward + local;

        assert_eq!(
            ControlProblem::bsde_driver(&m, 0.0, &state, &control, &derivs, dt),
            expected
        );
        assert!(
            ControlProblem::bsde_driver(&m, 0.0, &state, &control, &derivs, dt)
                != ControlProblem::driver(&m, 0.0, &state, &control, &derivs)
        );
    }

    #[test]
    fn next_step_controlled_uses_optimal_intensities() {
        let m = default_model();
        let state = [1.0, 1.0];

        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);

        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 intensity_scales_fill_rates() {
        let m = default_model();
        let derivs = zero_derivs();
        let ctrl_low = ControlProblem::optimize(&m, 0.0, &[0.0, 0.5], &derivs);
        let ctrl_high = ControlProblem::optimize(&m, 0.0, &[0.0, 2.0], &derivs);
        let total_low = ctrl_low.bid_intensity + ctrl_low.ask_intensity;
        let total_high = ctrl_high.bid_intensity + ctrl_high.ask_intensity;
        assert!(
            total_high > total_low,
            "Higher intensity should give higher fill rates: low={}, high={}",
            total_low,
            total_high
        );
    }

    #[test]
    fn subcritical_regime_mean_reverts() {
        let m = AvellanedaHawkes::new(0.5, 0.5, 1.5, 0.5, 2.0, 1.0);
        assert!(m.alpha < m.beta, "Subcritical: alpha < beta");

        // Equilibrium intensity = beta * mu / (beta - alpha)
        let eq = m.beta * m.mu / (m.beta - m.alpha);

        let state_high = [0.0, eq + 1.0];
        let next_high = ControlProblem::next_step(&m, 0.0, &state_high, 0.01, &[0.0, 0.0]);
        assert!(
            next_high[1] < state_high[1],
            "Above equilibrium, intensity should decrease"
        );

        let state_low = [0.0, eq - 0.5];
        let next_low = ControlProblem::next_step(&m, 0.0, &state_low, 0.01, &[0.0, 0.0]);
        assert!(
            next_low[1] > state_low[1],
            "Below equilibrium, intensity should increase"
        );
    }

    #[test]
    fn intensity_stays_non_negative() {
        let m = default_model();
        let state = [0.0, 0.001];
        let next = ControlProblem::next_step(&m, 0.0, &state, 0.01, &[0.0, 0.0]);
        assert!(next[1] >= 0.0, "Intensity must stay non-negative");
    }

    #[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::DeterministicDrift
        );
    }
}