solver/models/

kelly_hjb.rs

/// Rigorous Kelly (log-utility) HJB model for market making.
///
/// Solves the true log-utility HJB from first principles without the CARA
/// separation ansatz. Uses the homogeneity reduction
/// $V(t,S,q,X) = \log S + v(t, X/S, q)$ to reduce the state to $(q, x)$.
///
/// # HJB (reduced form)
///
/// $$\begin{aligned}
/// 0 = \partial_t v &+ \tfrac{1}{2}\sigma^2(x^2 v_{xx} + 2x v_x - 1) \\
///   &+ \sup_{\delta^a} A e^{-k\delta^a}[v(x+1, q-1) - v(x, q)] \\
///   &+ \sup_{\delta^b} A e^{-k\delta^b}[v(x-1, q+1) - v(x, q)]
/// \end{aligned}$$
///
/// Terminal: $v(T, x, q) = \log(x + q)$.
///
/// # Optimal Spreads (first-order, $\delta/S$ neglected)
///
/// $$\delta^{a*} = \frac{1}{k} + v(x,q) - v(x+1, q-1)$$
/// $$\delta^{b*} = \frac{1}{k} + v(x,q) - v(x-1, q+1)$$
///
/// Compared to CARA: no risk-aversion parameter $\gamma$; the base spread
/// $1/k$ is the $\gamma \to 0$ limit of the AS formula
/// $\frac{1}{\gamma}\ln(1+\gamma/k)$. Kelly spreads are typically wider
/// than CARA for the same fill parameters because log-utility penalizes
/// wealth drawdowns more severely.
use super::control::{ControlProblem, StateDerivatives};
use super::market_making::MarketMakingControl;

#[derive(Clone)]
pub struct KellyHjb {
    /// Price volatility ($\sigma$).
    pub sigma: f64,
    /// Fill-rate decay ($k$).
    pub kappa: f64,
    /// Base order arrival intensity ($A$).
    pub a: f64,
    /// Minimum inventory (lower hard boundary).
    pub q_min: f64,
    /// Maximum inventory (upper hard boundary).
    pub q_max: f64,
    /// Grid step for the inventory dimension.
    pub dq: f64,
}

impl KellyHjb {
    pub fn new(sigma: f64, kappa: f64, a: f64) -> Self {
        Self {
            sigma,
            kappa,
            a,
            q_min: -10.0,
            q_max: 10.0,
            dq: 1.0,
        }
    }

    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
    }

    /// Compute optimal bid/ask spreads from value-function gradients.
    ///
    /// `fwd[0] = dv/dq = v(q+1) - v(q)`, `bwd[0] = v(q) - v(q-1)` for dq=1.
    /// The x-dimension gradient enters through the finite-difference
    /// stencil on the grid and is handled by the solver's operator.
    pub fn get_spreads(&self, derivs: &StateDerivatives<2>) -> (f64, f64) {
        let base = 1.0 / self.kappa;

        let dv_buy = derivs.fwd[0];
        let dv_sell = derivs.bwd[0];

        let delta_bid = base - dv_buy;
        let delta_ask = base + dv_sell;

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

    /// 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 KellyHjb {
    type Control = MarketMakingControl;

    fn optimize(&self, _t: f64, state: &[f64; 2], derivs: &StateDerivatives<2>) -> Self::Control {
        let q = state[0];
        let (delta_bid, delta_ask) = self.get_spreads(derivs);

        let lambda_bid = if q >= self.q_max {
            0.0
        } else {
            self.a * (-self.kappa * delta_bid).exp()
        };
        let lambda_ask = if q <= self.q_min {
            0.0
        } else {
            self.a * (-self.kappa * delta_ask).exp()
        };

        MarketMakingControl::new(lambda_bid, lambda_ask)
    }

    fn running_reward(&self, _t: f64, _state: &[f64; 2], control: &Self::Control) -> f64 {
        // GLT Hamiltonian for log-utility (gamma -> 0 limit): (lambda_b + lambda_a) / k.
        (control.bid_intensity + control.ask_intensity) / 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 {
        // x dimension (dim 1): diffusion generator.
        // From V = log S + v(x,q), the Ito generator of V gives:
        //   L v = 0.5 sigma^2 (x^2 v_xx + 2x v_x - 1)
        let x = state[1];
        let sigma2 = self.sigma.powi(2);

        0.5 * sigma2 * x.powi(2) * derivs.hessian[1] + sigma2 * x * derivs.grad[1] - 0.5 * sigma2
    }

    fn terminal(&self, state: &[f64; 2]) -> f64 {
        let x = state[1].max(1e-12);
        let q = state[0];
        (x + q.max(-x + 1e-12)).ln()
    }

    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 q = state[0];
        let x = state[1].max(1e-8);
        let sqrt_dt = dt.sqrt();

        // d(x) = d(X/S). With dS = sigma S dW and dX = 0 between fills:
        //   dx = -x * sigma * dW + x * sigma^2 * dt   (Ito)
        let sigma2 = self.sigma.powi(2);
        let dx_drift = x * sigma2 * dt;
        let dx_diff = -x * self.sigma * sqrt_dt * noise[1];
        let x_next = (x + dx_drift + dx_diff).max(1e-8);

        [q, x_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.dq.abs().max(1e-8)
        } else {
            1.0
        }
    }
}

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

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

    fn default_model() -> KellyHjb {
        KellyHjb {
            sigma: 0.3,
            kappa: 1.5,
            a: 140.0,
            q_min: -10.0,
            q_max: 10.0,
            dq: 1.0,
        }
    }

    #[test]
    fn zero_gradient_gives_base_spread() {
        let m = default_model();
        let (bid, ask) = m.get_spreads(&zero_derivs());
        let base = 1.0 / m.kappa;
        assert!((bid - base).abs() < 1e-12);
        assert!((ask - base).abs() < 1e-12);
    }

    #[test]
    fn base_spread_matches_risk_neutral_as() {
        let m = default_model();
        let base = 1.0 / m.kappa;
        let as_base = {
            let gamma = 1e-9;
            (1.0 + gamma / m.kappa).ln() / gamma
        };
        assert!((base - as_base).abs() < 1e-6);
    }

    #[test]
    fn symmetric_at_zero_inventory() {
        let m = default_model();
        let ctrl = ControlProblem::optimize(&m, 0.0, &[0.0, 5.0], &zero_derivs());
        assert!(
            (ctrl.bid_intensity - ctrl.ask_intensity).abs() < 1e-12,
            "bid/ask intensities should be equal at q=0"
        );
    }

    #[test]
    fn positive_reward_at_zero_inventory() {
        let m = default_model();
        let ctrl = ControlProblem::optimize(&m, 0.0, &[0.0, 5.0], &zero_derivs());
        let reward = ControlProblem::running_reward(&m, 0.0, &[0.0, 5.0], &ctrl);
        assert!(
            reward > 0.0,
            "reward should be positive at q=0: hamiltonian > source"
        );
    }

    #[test]
    fn terminal_gives_log_x_plus_q() {
        let m = default_model();
        let v = ControlProblem::terminal(&m, &[2.0, 10.0]);
        assert!((v - 12.0_f64.ln()).abs() < 1e-10);
    }

    #[test]
    fn terminal_handles_negative_total_clamp() {
        let m = default_model();
        // x=2, q=-4 -> x+q = -2. Clamped: q -> -x+eps, x+q -> eps.
        let v = ControlProblem::terminal(&m, &[-4.0, 2.0]);
        assert!(v.is_finite());
        assert!(v < 0.0);
    }

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

    #[test]
    fn fill_rate_methods_return_expected() {
        let m = default_model();
        assert!((m.fill_rate_base(&[0.0, 5.0]) - 140.0).abs() < 1e-12);
        assert!((m.fill_rate_decay() - 1.5).abs() < 1e-12);
    }
}