solver/models/

bilateral_hawkes_order_flow_imbalance.rs

use super::normal_cdf;
use super::traits::{ControlOutput, Gradients, Model};

// Re-export so existing paths compile.
pub use super::TerminalCondition;

/// Bilateral Hawkes model with price-impact-driven adverse selection.
///
/// # State Space: `[q, lambda_plus, lambda_minus]`
///
/// | dim | symbol       | description                                      |
/// |-----|--------------|--------------------------------------------------|
/// | 0   | q            | Inventory level (integer, jump-controlled)       |
/// | 1   | lambda_plus  | Buy market-order intensity (ODE drift, no noise) |
/// | 2   | lambda_minus | Sell market-order intensity (ODE drift, no noise)|
///
/// # Dynamics
///
/// `d(lambda+) = beta*(mu - lambda+) dt + alpha dN+`
/// `d(lambda-) = beta*(mu - lambda-) dt + alpha dN-`
/// `dS = sigma dW + eta_ofi (dN+ - dN-)`
///
/// Each buy market order (dN+) moves the mid price UP by `eta_ofi`; each sell market
/// order (dN-) moves it DOWN by `eta_ofi`. The market maker accounts for this
/// price impact when setting optimal spreads.
///
/// # Optimal Spreads Under Price Impact
///
/// Deriving the FOC from the CARA HJB with per-fill wealth effect `(q±1)*eta_ofi` gives:
///
/// $$ \delta_{bid}^* = \delta_{base} - \frac{\partial V}{\partial q} + (q+1)\eta_{OFI} $$
/// $$ \delta_{ask}^* = \delta_{base} + \frac{\partial V}{\partial q} - (q-1)\eta_{OFI} $$
///
/// Economic interpretation:
/// - Long inventory (q > 0): ask narrows (selling into a rising price is favourable),
///   bid widens (avoid accumulating more when sell MOs push price down).
/// - Short inventory (q < 0): ask widens (avoid shorting further into a rising price),
///   bid narrows (cover cheaply when sell MOs push price down).
///
/// Setting `eta_ofi = 0` recovers the plain `BilateralHawkes` (BHK) model exactly.
///
/// In the simulation pass `impact_factor = -eta_ofi` to `SimulatedDataSource` so the
/// price process matches the assumed dynamics `dS = sigma dW + eta_ofi (dN+ - dN-)`.
///
/// # Hawkes Approximation
///
/// Same fluid (mean-field) approximation as `BilateralHawkes`:
///
/// $$ d(\lambda_+)/dt \approx \beta(\mu - \lambda_+) + \alpha \lambda_+ $$
/// $$ d(\lambda_-)/dt \approx \beta(\mu - \lambda_-) + \alpha \lambda_- $$
///
#[derive(Clone, Debug)]
pub struct BilateralHawkesOrderFlowImbalance {
    /// Risk-aversion parameter.
    pub gamma: f64,
    /// Price volatility.
    pub sigma: f64,
    /// Fill-rate decay (Avellaneda-Stoikov kappa).
    pub kappa: f64,
    // --- Hawkes parameters (shared for both sides) ---
    /// Jump size added to the excited side's intensity.
    pub alpha: f64,
    /// Mean-reversion speed.
    pub beta: f64,
    /// Baseline intensity (mean-reversion target).
    pub mu: f64,
    // --- Price-impact parameter ---
    /// Price impact per fill event (eta_ofi in dS = sigma*dW + eta_ofi*(dN+ - dN-)).
    /// Higher eta_ofi = stronger q-dependent spread correction.
    /// Setting eta_ofi = 0 recovers the plain BilateralHawkes (BHK) model.
    pub eta_ofi: f64,
    /// Grid step for both lambda+ and lambda- dimensions.
    pub lambda_step: f64,
    /// Grid step for the inventory dimension.
    pub dq: f64,
    /// Minimum inventory (lower hard boundary). Defaults to -infinity.
    pub q_min: f64,
    /// Maximum inventory (upper hard boundary). Defaults to +infinity.
    pub q_max: f64,
    /// Terminal condition for the value function at T. Defaults to Zero.
    pub terminal_condition: TerminalCondition,
    /// 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 BilateralHawkesOrderFlowImbalance {
    pub fn new(
        gamma: f64,
        sigma: f64,
        kappa: f64,
        alpha: f64,
        beta: f64,
        mu: f64,
        eta_ofi: f64,
    ) -> Self {
        Self {
            gamma,
            sigma,
            kappa,
            alpha,
            beta,
            mu,
            eta_ofi,
            lambda_step: 1.0,
            dq: 1.0,
            q_min: f64::NEG_INFINITY,
            q_max: f64::INFINITY,
            terminal_condition: TerminalCondition::Zero,
            terminal_liquidation_half_spread: None,
        }
    }

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

    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_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 bid/ask half-spreads with price-impact-driven inventory correction.
    ///
    /// When eta_ofi > 0 each fill carries an adverse-selection wealth effect `(q±1)*eta_ofi`.
    /// The FOC on the CARA HJB gives:
    ///   delta_bid = base_spread - dV/dq_fwd + (q+1)*eta_ofi
    ///   delta_ask = base_spread + dV/dq_bwd - (q-1)*eta_ofi
    ///
    /// Long inventory (q > 0): ask narrows (sell into rising price), bid widens.
    /// Short inventory (q < 0): ask widens (avoid adding short), bid narrows.
    /// eta_ofi = 0 recovers the plain BHK spread formula identically.
    pub fn get_spreads(&self, state: &[f64; 3], grads: &Gradients<3>) -> (f64, f64) {
        let q = state[0];
        let base_spread = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();
        let delta_bid = (base_spread - grads.fwd[0] + (q + 1.0) * self.eta_ofi).clamp(-10.0, 10.0);
        let delta_ask = (base_spread + grads.bwd[0] - (q - 1.0) * self.eta_ofi).clamp(-10.0, 10.0);
        (delta_bid, delta_ask)
    }

    /// Upwind transport rates for one intensity dimension.
    #[inline]
    fn intensity_upwind(&self, drift: f64) -> (f64, f64) {
        let rate = drift.abs() / self.lambda_step;
        if drift >= 0.0 {
            (rate, 0.0)
        } else {
            (0.0, rate)
        }
    }
}

impl Model<3> for BilateralHawkesOrderFlowImbalance {
    type Process = ();

    fn process(&self) {}

    fn optimize(&self, state: &[f64; 3], grads: &Gradients<3>) -> ControlOutput<3> {
        let q = state[0];
        let lambda_plus = state[1].max(0.0);
        let lambda_minus = state[2].max(0.0);

        // --- Dimension 0: Inventory controls ---
        let (d_bid, d_ask) = self.get_spreads(state, grads);

        let max_mult = 20.0_f64;
        // Bid filled by sell MOs (lambda-); ask filled by buy MOs (lambda+)
        let lambda_bid_fill = lambda_minus * (-self.kappa * d_bid).exp().min(max_mult);
        let lambda_ask_fill = lambda_plus * (-self.kappa * d_ask).exp().min(max_mult);

        // Shut off forbidden quoting side at inventory bounds.
        let lambda_bid_fill = if q >= self.q_max {
            0.0
        } else {
            lambda_bid_fill
        };
        let lambda_ask_fill = if q <= self.q_min {
            0.0
        } else {
            lambda_ask_fill
        };

        let hamiltonian_inv = (lambda_bid_fill + lambda_ask_fill) / (self.gamma + self.kappa);
        let risk_penalty = -0.5 * self.gamma * self.sigma.powi(2) * q.powi(2);

        // Cancel the solver's q-transport contribution so we can inject the Hamiltonian.
        let drift_correction_q =
            self.dq * (lambda_bid_fill * grads.fwd[0] - lambda_ask_fill * grads.bwd[0]);

        // --- Dimensions 1/2: lambda+ and lambda- dynamics (fluid approximation) ---
        let net_drift_plus = self.beta * (self.mu - lambda_plus) + self.alpha * lambda_plus;
        let net_drift_minus = self.beta * (self.mu - lambda_minus) + self.alpha * lambda_minus;
        let (lp_plus, lp_minus) = self.intensity_upwind(net_drift_plus);
        let (lm_plus, lm_minus) = self.intensity_upwind(net_drift_minus);

        ControlOutput {
            lambda_plus: [lambda_bid_fill, lp_plus, lm_plus],
            lambda_minus: [lambda_ask_fill, lp_minus, lm_minus],
            flow: hamiltonian_inv + risk_penalty - drift_correction_q,
        }
    }

    fn terminal(&self, state: &[f64; 3]) -> 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 constant_discount_rate(&self) -> Option<f64> {
        Some(0.0)
    }

    /// Forward step using the fluid (mean-field) drift.
    fn next_step(&self, current_state: &[f64; 3], dt: f64, _noise: &[f64; 3]) -> [f64; 3] {
        let lambda_plus = current_state[1];
        let lambda_minus = current_state[2];
        let drift_plus = self.beta * (self.mu - lambda_plus) + self.alpha * lambda_plus;
        let drift_minus = self.beta * (self.mu - lambda_minus) + self.alpha * lambda_minus;
        [
            current_state[0],
            (lambda_plus + drift_plus * dt).max(0.0),
            (lambda_minus + drift_minus * dt).max(0.0),
        ]
    }

    /// Coupled step: inventory jumps based on current fill probabilities.
    fn next_step_controlled(
        &self,
        current_state: &[f64; 3],
        control: &ControlOutput<3>,
        dt: f64,
        noise: &[f64; 3],
    ) -> [f64; 3] {
        let mut next = self.next_step(current_state, dt, noise);

        // Bernoulli inventory jump using noise[0] mapped through CDF -> uniform
        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; // sell MO hits bid  -> inventory +1
            next[2] += self.alpha; // lambda- self-excites on sell MO fill
        } else if u > 1.0 - p_ask {
            next[0] -= 1.0; // buy  MO hits ask  -> inventory -1
            next[1] += self.alpha; // lambda+ self-excites on buy MO fill
        }

        next[1] = next[1].max(0.0);
        next[2] = next[2].max(0.0);

        next
    }

    /// Both intensity dimensions follow deterministic ODE drift with no Brownian noise.
    fn is_diffusion_dimension(&self, _dim: usize) -> bool {
        false
    }

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

    fn gradient_step(&self, dim: usize) -> f64 {
        match dim {
            1 | 2 => self.lambda_step.abs().max(1e-8),
            _ => 1.0,
        }
    }

    /// The base arrival intensity is state-dependent: returns the average of
    /// the buy/sell intensity states for intensity-to-spread conversion.
    fn fill_rate_base(&self, state: &[f64; 3]) -> f64 {
        ((state[1] + state[2]) * 0.5).max(1e-10)
    }

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

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

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

    fn zero_state() -> [f64; 3] {
        [0.0, 1.0, 1.0]
    }

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

    #[test]
    fn model_creation() {
        let m = default_model();
        assert_eq!(m.gamma, 0.5);
        assert_eq!(m.sigma, 0.2);
        assert_eq!(m.kappa, 1.5);
        assert_eq!(m.alpha, 0.5);
        assert_eq!(m.beta, 2.0);
        assert_eq!(m.mu, 1.0);
        assert_eq!(m.eta_ofi, 0.5);
    }

    #[test]
    fn zero_inventory_gives_symmetric_spreads() {
        // At q=0, zero grads: delta_bid = base + (0+1)*xi, delta_ask = base - (0-1)*xi
        // Both equal base + xi, so bid == ask (symmetric market)
        let m = default_model();
        let state = [0.0, 1.0, 1.0];
        let (bid, ask) = m.get_spreads(&state, &zero_grads());
        let base = (1.0 / m.gamma) * (1.0 + m.gamma / m.kappa).ln();
        assert!(
            (bid - ask).abs() < 1e-12,
            "bid and ask should be equal at q=0"
        );
        assert!((bid - (base + m.eta_ofi)).abs() < 1e-12);
    }

    #[test]
    fn long_inventory_widens_bid_narrows_ask() {
        // q > 0: MM is long. To offload, widen bid (post higher), narrow ask.
        // delta_bid = base + (q+1)*xi > base+xi, delta_ask = base - (q-1)*xi < base+xi
        let m = default_model(); // xi = 0.5
        let q = 5.0;
        let state_long = [q, 1.0, 1.0];
        let state_zero = [0.0, 1.0, 1.0];
        let (bid_long, ask_long) = m.get_spreads(&state_long, &zero_grads());
        let (bid_zero, ask_zero) = m.get_spreads(&state_zero, &zero_grads());
        assert!(bid_long > bid_zero, "long position widens bid");
        assert!(ask_long < ask_zero, "long position narrows ask");
    }

    #[test]
    fn short_inventory_narrows_bid_widens_ask() {
        // q < 0: MM is short. To cover, narrow bid (easier buy), widen ask.
        // delta_bid = base + (q+1)*xi < base+xi, delta_ask = base - (q-1)*xi > base+xi
        let m = default_model(); // xi = 0.5
        let q = -5.0;
        let state_short = [q, 1.0, 1.0];
        let state_zero = [0.0, 1.0, 1.0];
        let (bid_short, ask_short) = m.get_spreads(&state_short, &zero_grads());
        let (bid_zero, ask_zero) = m.get_spreads(&state_zero, &zero_grads());
        assert!(bid_short < bid_zero, "short position narrows bid");
        assert!(ask_short > ask_zero, "short position widens ask");
    }

    #[test]
    fn larger_xi_gives_larger_spread_corrections() {
        // At q=5 (long), higher xi should produce larger deviations from zero-xi baseline.
        // delta_bid = base + (q+1)*xi, delta_ask = base - (q-1)*xi
        let m1 = BilateralHawkesOrderFlowImbalance::new(0.5, 0.2, 1.5, 0.5, 2.0, 1.0, 0.5)
            .with_lambda_step(0.1);
        let m2 = BilateralHawkesOrderFlowImbalance::new(0.5, 0.2, 1.5, 0.5, 2.0, 1.0, 1.0)
            .with_lambda_step(0.1);
        let m0 = BilateralHawkesOrderFlowImbalance::new(0.5, 0.2, 1.5, 0.5, 2.0, 1.0, 0.0)
            .with_lambda_step(0.1);

        let state = [5.0, 1.0, 1.0]; // q=5 (long)
        let grads = zero_grads();

        let (bid0, ask0) = m0.get_spreads(&state, &grads);
        let (bid1, ask1) = m1.get_spreads(&state, &grads);
        let (bid2, ask2) = m2.get_spreads(&state, &grads);

        // bid monotone increasing in xi at q>0
        assert!(bid1 > bid0);
        assert!(bid2 > bid1);
        // ask monotone decreasing in xi at q>0
        assert!(ask1 < ask0);
        assert!(ask2 < ask1);
    }

    #[test]
    fn terminal_condition_zero() {
        let m = default_model().with_terminal_condition(TerminalCondition::Zero);
        let state = [5.0, 1.0, 1.0];
        assert_eq!(m.terminal(&state), 0.0);
    }

    #[test]
    fn terminal_condition_liquidation_cost() {
        let m = default_model()
            .with_terminal_condition(TerminalCondition::LiquidationCost)
            .with_terminal_liquidation_half_spread(0.1);

        let state_q5 = [5.0, 1.0, 1.0];
        let state_q0 = [0.0, 1.0, 1.0];
        let state_qm3 = [-3.0, 1.0, 1.0];

        assert_eq!(m.terminal(&state_q5), -5.0 * 0.1);
        assert_eq!(m.terminal(&state_q0), 0.0);
        assert_eq!(m.terminal(&state_qm3), -3.0 * 0.1);
    }

    #[test]
    fn optimize_output_structure() {
        let m = default_model();
        let state = zero_state();
        let grads = zero_grads();
        let output = m.optimize(&state, &grads);

        // Check that output dimensions are correctly sized
        assert_eq!(output.lambda_plus.len(), 3);
        assert_eq!(output.lambda_minus.len(), 3);

        // All rates should be non-negative
        for rate in &output.lambda_plus {
            assert!(*rate >= 0.0);
        }
        for rate in &output.lambda_minus {
            assert!(*rate >= 0.0);
        }
    }

    #[test]
    fn inventory_bounds_enforcement() {
        let m = default_model().with_inventory_bounds(-2.0, 2.0);
        let state = [2.0, 2.0, 1.0]; // at upper bound
        let grads = Gradients {
            fwd: [0.0, 0.0, 0.0],
            bwd: [-1.0, 0.0, 0.0], // negative gradient would encourage sells
        };

        let output = m.optimize(&state, &grads);
        // At upper bound, can't increase inventory further
        // lambda_plus[0] = lambda_bid_fill should be zero (no more sell MOs accepted)
        assert_eq!(output.lambda_plus[0], 0.0);
    }

    #[test]
    fn hawkes_dynamics_stability() {
        // Test that Hawkes dynamics remain stable (lambda doesn't grow without bound)
        let m = default_model();

        // Stationary intensity should be mu / (1 - alpha/beta) when both equal
        let rho = m.alpha / m.beta;
        let stationary_intensity = m.mu / (1.0 - rho);

        let mut state = [0.0, stationary_intensity, stationary_intensity];

        // Step forward many times; intensity should stay bounded
        for _ in 0..1000 {
            state = m.next_step(&state, 0.01, &[0.0, 0.0, 0.0]);
        }

        // After many steps, intensity should remain close to stationary level
        assert!(state[1] > 0.0 && state[1] < 10.0 * stationary_intensity);
        assert!(state[2] > 0.0 && state[2] < 10.0 * stationary_intensity);
    }

    #[test]
    fn next_step_consistency() {
        // Verify that next_step uses the same fluid dynamics as optimize()
        let m = default_model();
        let state = [0.0, 1.5, 1.0];
        let dt = 0.01;

        let next = m.next_step(&state, dt, &[0.0; 3]);

        // Inventory should stay the same (no fills in next_step)
        assert_eq!(next[0], state[0]);

        // Lambda+ should decay slightly toward mu (since lambda+ > mu * (1 - rho))
        let drift_plus = m.beta * (m.mu - state[1]) + m.alpha * state[1];
        let expected_lambda_plus = state[1] + drift_plus * dt;
        assert!((next[1] - expected_lambda_plus).abs() < 1e-10);
    }

    #[test]
    fn next_step_controlled_fills_inventory() {
        // Test that next_step_controlled can generate inventory jumps
        let m = default_model();
        let state = [0.0, 2.0, 2.0];

        // Create a control with high bid fill rate
        let control = ControlOutput {
            lambda_plus: [1.0, 0.0, 0.0],   // high bid fill
            lambda_minus: [0.01, 0.0, 0.0], // low ask fill
            flow: 0.0,
        };

        let dt = 0.1;

        // Run many steps and count fills
        let mut up_count = 0;
        let mut down_count = 0;
        for i in 0..100 {
            let noise = [(i as f64 - 50.0) / 10.0; 3]; // sweep through normal values
            let next = m.next_step_controlled(&state, &control, dt, &noise);

            if next[0] > state[0] {
                up_count += 1;
            } else if next[0] < state[0] {
                down_count += 1;
            }
        }

        // With high bid fill rate, we should see many upward inventory moves
        assert!(up_count > down_count);
    }

    #[test]
    fn gradient_step_returns_correct_values() {
        let m = default_model();

        assert_eq!(m.gradient_step(0), 1.0);
        assert_eq!(m.gradient_step(1), 0.1);
        assert_eq!(m.gradient_step(2), 0.1);
    }

    #[test]
    fn is_integer_dimension() {
        let m = default_model();

        assert!(m.is_integer_dimension(0)); // inventory is discrete
        assert!(!m.is_integer_dimension(1)); // lambda+ is continuous
        assert!(!m.is_integer_dimension(2)); // lambda- is continuous
    }

    #[test]
    fn is_diffusion_dimension() {
        let m = default_model();

        assert!(!m.is_diffusion_dimension(0));
        assert!(!m.is_diffusion_dimension(1));
        assert!(!m.is_diffusion_dimension(2));
    }

    #[test]
    fn constant_discount_rate() {
        let m = default_model();
        assert_eq!(m.constant_discount_rate(), Some(0.0));
    }

    #[test]
    fn spread_formula_matches_analytical_expression() {
        // delta_bid(q) = base + (q+1)*xi, delta_ask(q) = base - (q-1)*xi
        let m = default_model(); // xi = 0.5
        let base = (1.0 / m.gamma) * (1.0 + m.gamma / m.kappa).ln();
        for &q in &[-5.0_f64, -1.0, 0.0, 1.0, 5.0] {
            let state = [q, 1.0, 1.0];
            let (bid, ask) = m.get_spreads(&state, &zero_grads());
            let expected_bid = base + (q + 1.0) * m.eta_ofi;
            let expected_ask = base - (q - 1.0) * m.eta_ofi;
            assert!((bid - expected_bid).abs() < 1e-12, "bid mismatch at q={q}");
            assert!((ask - expected_ask).abs() < 1e-12, "ask mismatch at q={q}");
        }
    }

    #[test]
    fn compare_with_balanced_no_imbalance() {
        // When lambdas are balanced, should behave like a symmetric market maker
        let m = default_model();
        let state = [0.0, 1.0, 1.0];
        let grads = Gradients {
            fwd: [0.0, 0.0, 0.0],
            bwd: [0.0, 0.0, 0.0],
        };

        let (bid, ask) = m.get_spreads(&state, &grads);
        // At q=0, zero grads, xi=0.5: bid == ask == base + xi (symmetric)
        assert!((bid - ask).abs() < 1e-12);
    }

    #[test]
    fn extreme_imbalance_clamping() {
        // Spreads should be clamped to [-10, 10] range
        let m = BilateralHawkesOrderFlowImbalance::new(0.5, 0.2, 1.5, 0.5, 2.0, 1.0, 100.0)
            .with_lambda_step(0.1);

        let state = [0.0, 10.0, 0.0]; // extreme imbalance
        let grads = Gradients {
            fwd: [100.0, 0.0, 0.0],
            bwd: [100.0, 0.0, 0.0],
        };

        let (bid, ask) = m.get_spreads(&state, &grads);

        // Spreads should be bounded
        assert!(bid >= -10.0 && bid <= 10.0);
        assert!(ask >= -10.0 && ask <= 10.0);
    }

    #[test]
    fn inventory_control_respects_bounds() {
        let m = default_model().with_inventory_bounds(-5.0, 5.0);
        let state_at_max = [5.0, 2.0, 1.0];
        let state_at_min = [-5.0, 1.0, 2.0];
        let grads = zero_grads();

        let output_max = m.optimize(&state_at_max, &grads);
        let output_min = m.optimize(&state_at_min, &grads);

        // At max inventory, can't increase (lambda_plus[0] = lambda_bid_fill should be zero)
        assert_eq!(output_max.lambda_plus[0], 0.0);

        // At min inventory, can't decrease (lambda_minus[0] = lambda_ask_fill should be zero)
        assert_eq!(output_min.lambda_minus[0], 0.0);
    }

    #[test]
    fn xi_zero_recovers_baseline() {
        // With xi=0, get_spreads reduces to pure AS formula at q=0:
        // delta_bid = base + (0+1)*0 = base, delta_ask = base - (0-1)*0 = base
        let m = BilateralHawkesOrderFlowImbalance::new(0.5, 0.2, 1.5, 0.5, 2.0, 1.0, 0.0)
            .with_lambda_step(0.1);
        let state = [0.0, 2.0, 1.0];
        let (bid, ask) = m.get_spreads(&state, &zero_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 model_builder_methods() {
        let m = BilateralHawkesOrderFlowImbalance::new(0.5, 0.2, 1.5, 0.5, 2.0, 1.0, 0.5)
            .with_lambda_step(0.2)
            .with_dq(0.5)
            .with_inventory_bounds(-10.0, 10.0)
            .with_terminal_condition(TerminalCondition::LiquidationCost)
            .with_terminal_liquidation_half_spread(0.05);

        assert_eq!(m.lambda_step, 0.2);
        assert_eq!(m.dq, 0.5);
        assert_eq!(m.q_min, -10.0);
        assert_eq!(m.q_max, 10.0);
        assert_eq!(m.terminal_liquidation_half_spread, Some(0.05));
    }
}