solver/models/

kelly_hjb.rs

/// Rigorous Kelly (log-utility) HJB model for market making.
///
/// **Rigor: Rigorous.** 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::traits::{ControlOutput, Gradients, Model};

#[derive(Clone)]
pub struct KellyHjb {
    pub sigma: f64,
    pub kappa: f64,
    pub a: f64,
    pub q_min: f64,
    pub q_max: f64,
    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, grads: &Gradients<2>) -> (f64, f64) {
        let base = 1.0 / self.kappa;

        let dv_buy = grads.fwd[0]; // v(q+dx) - v(q)
        let dv_sell = grads.bwd[0]; // v(q) - v(q-dx)

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

impl Model<2> for KellyHjb {
    type Process = ();

    fn process(&self) {}

    fn optimize(&self, state: &[f64; 2], grads: &Gradients<2>) -> ControlOutput<2> {
        let q = state[0];
        let x = state[1];

        // --- Inventory (dim 0): fill intensities ---
        let (delta_bid, delta_ask) = self.get_spreads(grads);

        let lambda_bid = self.a * (-self.kappa * delta_bid).exp();
        let lambda_ask = self.a * (-self.kappa * delta_ask).exp();

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

        // GLT hamiltonian for log-utility (gamma -> 0 limit):
        //   H = (lambda_bid + lambda_ask) / k
        let hamiltonian = (lambda_bid + lambda_ask) / self.kappa;

        // Cancel solver's jump terms to avoid double-counting.
        let jump_correction = self.dq * (lambda_bid * grads.fwd[0] - lambda_ask * grads.bwd[0]);

        // --- 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)
        // Upwind discretization with spacing dx:
        let sigma2 = self.sigma.powi(2);
        let diff_coeff = 0.5 * sigma2 * x.powi(2);
        let drift_x = sigma2 * x;

        let dx: f64 = 1.0; // x-grid spacing, aligned with unit fill jump
        let diff_term = diff_coeff / dx.powi(2);
        let drift_term_abs = drift_x.abs() / dx;

        let mut lambda_x_plus = diff_term;
        let mut lambda_x_minus = diff_term;

        if drift_x > 0.0 {
            lambda_x_plus += drift_term_abs;
        } else {
            lambda_x_minus += drift_term_abs;
        }

        // Constant source term from the Ito generator: -0.5 sigma^2
        let const_source = -0.5 * sigma2;

        ControlOutput {
            lambda_plus: [lambda_bid, lambda_x_plus],
            lambda_minus: [lambda_ask, lambda_x_minus],
            flow: hamiltonian + const_source - jump_correction,
        }
    }

    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 constant_discount_rate(&self) -> Option<f64> {
        Some(0.0)
    }

    fn next_step(&self, current_state: &[f64; 2], dt: f64, noise: &[f64; 2]) -> [f64; 2] {
        let q = current_state[0];
        let x = current_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_diffusion_dimension(&self, dim: usize) -> bool {
        dim == 1 // x diffuses, q is jump-controlled
    }

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

    fn fill_rate_base(&self, _state: &[f64; 2]) -> f64 {
        self.a
    }

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

#[cfg(test)]
mod tests {
    use super::super::traits::{Gradients, Model};
    use super::*;

    fn zero_grads() -> Gradients<2> {
        Gradients {
            fwd: [0.0; 2],
            bwd: [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_grads());
        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 = m.optimize(&[0.0, 5.0], &zero_grads());
        assert!(
            (ctrl.lambda_plus[0] - ctrl.lambda_minus[0]).abs() < 1e-12,
            "bid/ask intensities should be equal at q=0"
        );
    }

    #[test]
    fn positive_flow_at_zero_inventory() {
        let m = default_model();
        let ctrl = m.optimize(&[0.0, 5.0], &zero_grads());
        assert!(
            ctrl.flow > 0.0,
            "flow should be positive at q=0: hamiltonian > source"
        );
    }

    #[test]
    fn terminal_gives_log_x_plus_q() {
        let m = default_model();
        let v = m.terminal(&[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 = m.terminal(&[-4.0, 2.0]);
        assert!(v.is_finite());
        assert!(v < 0.0);
    }

    #[test]
    fn dims_tagged_correctly() {
        let m = default_model();
        assert!(!m.is_diffusion_dimension(0));
        assert!(m.is_diffusion_dimension(1));
        assert!(m.is_integer_dimension(0));
        assert!(!m.is_integer_dimension(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);
    }
}