solver/analytical/avellaneda/approximations/

gueant.rs

use crate::analytical::traits::AnalyticalSolution;
use crate::linalg::eigen::eigen_symmetric_tridiagonal;
use crate::models::market_making::MarketMakingControl;

/// The asymptotic approximation (infinite horizon) for the Avellaneda-Stoikov model.
/// Based on the Gueant-Lehalle-Tapia closed-form solution.
///
/// Reference: Guéant, O., Lehalle, C. A., & Tapia, J. (2013). Dealing with the inventory risk: a solution to the market making problem.
pub struct AvellanedaGueant {
    pub gamma: f64,
    pub sigma: f64,
    pub kappa: f64,
    pub a: f64,
    /// Inventory truncation used when evaluating the value reference.
    pub q_max: usize,
}

impl AvellanedaGueant {
    pub fn new(gamma: f64, sigma: f64, kappa: f64, a: f64) -> Self {
        Self {
            gamma,
            sigma,
            kappa,
            a,
            q_max: 50,
        }
    }

    /// Sets the inventory truncation used for the value reference.
    pub fn with_q_max(mut self, q_max: usize) -> Self {
        assert!(q_max > 0, "q_max must be positive");
        self.q_max = q_max;
        self
    }

    /// Guéant-Lehalle-Tapia closed-form approximation for optimal spreads.
    ///
    /// This approximation assumes infinite time horizon ($T \to \infty$) and small spreads,
    /// providing a stationary strategy that depends on inventory $q$ but not time $t$.
    ///
    /// The spreads are symmetric around the mid-price, shifted by inventory:
    /// $\delta^\pm(q) \approx \text{RiskNeutral} \pm \text{InventorySkew} \cdot q$
    ///
    /// Reference: Guéant, O., Lehalle, C. A., & Tapia, J. (2013). Dealing with the inventory risk.
    pub fn approximate_spreads(&self, q: f64) -> (f64, f64) {
        let risk_factor = (self.sigma.powi(2) * self.gamma / (2.0 * self.kappa * self.a)).sqrt();
        let const_term = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();
        let correction =
            (1.0 + self.gamma / self.kappa).powf(0.5 * (1.0 + self.kappa / self.gamma));

        let delta_bid = const_term + (2.0 * q + 1.0) / 2.0 * risk_factor * correction;
        let delta_ask = const_term - (2.0 * q - 1.0) / 2.0 * risk_factor * correction;

        (delta_bid, delta_ask)
    }

    /// Stationary value function `theta(q)` for the reduced Avellaneda model.
    ///
    /// The exact stationary solution of the GLT linearized HJB is the principal
    /// eigenvector `v(q)` of the symmetric tridiagonal operator with diagonal
    /// `-alpha q^2` and off-diagonal `eta`, where
    ///
    /// ```text
    /// alpha = (kappa / 2) gamma sigma^2
    /// eta   = a (1 + gamma/kappa)^{-(1 + kappa/gamma)}
    /// ```
    ///
    /// The reduced solver works directly in `theta` space, with
    /// `v(q) = exp(kappa theta(q))`. The value is normalized so `theta(0) = 0`,
    /// which removes the one-parameter translation invariance of the stationary
    /// HJB. Only differences `theta(q) - theta(q')` are invariant under the
    /// normalization and therefore meaningful.
    ///
    /// # Panics
    ///
    /// Panics if `q` is outside `[-q_max, q_max]`.
    pub fn stationary_theta(&self, q_max: usize, q: i32) -> f64 {
        let n = 2 * q_max + 1;
        assert!(
            q >= -(q_max as i32) && q <= q_max as i32,
            "inventory q must lie within [-q_max, q_max]"
        );
        let idx = (q + q_max as i32) as usize;
        assert!(idx < n, "index {idx} out of range for n={n}");

        let alpha = (self.kappa / 2.0) * self.gamma * self.sigma.powi(2);
        let eta = self.a * (1.0 + self.gamma / self.kappa).powf(-(1.0 + self.kappa / self.gamma));

        let mut d = vec![0.0; n];
        let mut e = vec![0.0; n - 1];
        for i in 0..n {
            let q_val = (i as i32) - (q_max as i32);
            d[i] = -alpha * (q_val as f64).powi(2);
            if i < n - 1 {
                e[i] = eta;
            }
        }

        let (_evals, evecs) = eigen_symmetric_tridiagonal(d, e);
        // dstev returns eigenvalues ascending; the principal eigenvector is the
        // one with the largest eigenvalue.
        let k = n - 1;
        let v_q = evecs[idx + k * n];
        let v_0 = evecs[q_max + k * n];

        (1.0 / self.kappa) * (v_q / v_0).ln()
    }

    /// Stationary bid/ask spreads implied by the principal-eigenvector value.
    ///
    /// The bid spread is `base + (theta(q) - theta(q+1)) / 2` adjusted for the
    /// half-inventory convention, and the ask spread is the mirror around the
    /// base spread. This matches the `approximate_spreads` small-gamma limit.
    pub fn stationary_spreads(&self, q_max: usize, q: i32) -> (f64, f64) {
        let base = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();
        let theta_q = self.stationary_theta(q_max, q);

        // Use one-sided differences at the inventory boundary, matching the
        // solver's boundary derivative convention.
        let theta_up = if q < q_max as i32 {
            self.stationary_theta(q_max, q + 1)
        } else {
            theta_q
        };
        let theta_dn = if q > -(q_max as i32) {
            self.stationary_theta(q_max, q - 1)
        } else {
            theta_q
        };

        let bid = base + (theta_q - theta_up);
        let ask = base + (theta_q - theta_dn);
        (bid, ask)
    }
}

impl AnalyticalSolution<2> for AvellanedaGueant {
    fn value_function(&self, _t: f64, state: &[f64; 2]) -> f64 {
        self.stationary_theta(self.q_max, state[0].round() as i32)
    }

    fn optimal_controls(&self, _t: f64, state: &[f64; 2]) -> MarketMakingControl {
        let q = state[0];
        let (d_bid, d_ask) = self.approximate_spreads(q);

        let lambda_bid = self.a * (-self.kappa * d_bid).exp();
        let lambda_ask = self.a * (-self.kappa * d_ask).exp();

        MarketMakingControl::new(lambda_bid, lambda_ask)
    }
}

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

    fn default_gueant() -> AvellanedaGueant {
        AvellanedaGueant::new(0.1, 0.1, 1.5, 140.0).with_q_max(50)
    }

    #[test]
    fn stationary_theta_is_zero_at_zero_inventory() {
        let gueant = default_gueant();
        let theta = gueant.stationary_theta(50, 0);
        assert!((theta - 0.0).abs() < 1e-12);
    }

    #[test]
    fn stationary_theta_decreases_with_inventory_magnitude() {
        let gueant = default_gueant();
        let theta0 = gueant.stationary_theta(50, 0);
        let theta1 = gueant.stationary_theta(50, 1);
        let theta2 = gueant.stationary_theta(50, 2);
        assert!(theta1 < theta0);
        assert!(theta2 < theta1);
    }

    #[test]
    fn stationary_spreads_match_small_gamma_approximation() {
        let gueant = default_gueant();
        let q = 3;
        let (bid_approx, ask_approx) = gueant.approximate_spreads(q as f64);
        let (bid, ask) = gueant.stationary_spreads(50, q);

        // The GLT approximation and the principal-eigenvector shape agree in the
        // small-gamma (small-spread) limit. A loose tolerance is intentional.
        assert!((bid - bid_approx).abs() < 1e-2, "bid {bid} vs {bid_approx}");
        assert!((ask - ask_approx).abs() < 1e-2, "ask {ask} vs {ask_approx}");
    }
}