solver/models/

stationary_avellaneda.rs

//! Infinite-horizon Avellaneda-Stoikov market making in inventory space.
//!
//! This is a validating instance of the ergodic stationary path
//! ([`crate::numeric::ergodic`]). It is deliberately not an
//! [`crate::numeric::finite_difference::elliptic::EllipticControlProblem`]:
//! the undiscounted market-making HJB is translation invariant and has no
//! unique solution of `sup_u { f + L^u V } = 0`. Its exact
//! `T -> infinity` limit is the principal eigenvector of the GLT operator.
//!
//! # Problem and conditions
//!
//! State: inventory `q in Z`. Control: bid/ask half-spreads
//! `(delta_b, delta_a)`. Fill intensities are
//!
//! ```text
//! lambda_b = a exp(-kappa delta_b)
//! lambda_a = a exp(-kappa delta_a)
//! ```
//!
//! The finite-horizon CARA value reduces to `theta(t, q)`; see
//! `docs/src/reference/soc_exact.md`.
//!
//! # Exact solution
//!
//! Let `v_q = exp(kappa theta_q)`. The GLT substitution gives the linear
//! system
//!
//! ```text
//! v_q_t = alpha q^2 v_q - eta (v_{q-1} + v_{q+1})
//! alpha = (kappa / 2) gamma sigma^2
//! eta   = a (1 + gamma/kappa)^-(1 + kappa/gamma)
//! ```
//!
//! As `T -> infinity` the solution concentrates on the principal (Perron)
//! eigenvector of the negated operator, assembled here as
//!
//! ```text
//! A_qq = -alpha q^2,   A_{q,q+1} = A_{q,q-1} = +eta
//! ```
//!
//! The eigenvector is normalized by `theta(0) = 0`:
//!
//! ```text
//! theta(q) = (1/kappa) ln(v_q / v_0)
//! ```
//!
//! Optimal stationary spreads are recovered as
//!
//! ```text
//! delta_b = base + theta(q) - theta(q+1)
//! delta_a = base + theta(q) - theta(q-1)
//! base    = (1/gamma) ln(1 + gamma/kappa)
//! ```
//!
//! See [`crate::analytical::avellaneda::approximations::gueant::AvellanedaGueant`]
//! for the reference implementation used in tests.
use crate::core::grid::Grid;
use crate::linalg::csr::CsrMatrix;
use crate::numeric::ergodic::StationaryEigenProblem;

/// Infinite-horizon Avellaneda-Stoikov market making model in inventory space.
///
/// See the module-level documentation for the formulation and exact solution.
#[derive(Clone, Debug)]
pub struct StationaryAvellaneda {
    /// Risk aversion for inventory holding.
    pub gamma: f64,
    /// Order filling intensity decay (`kappa`).
    pub kappa: f64,
    /// Base order arrival intensity (`A`).
    pub a: f64,
    /// Mid-price volatility.
    pub sigma: f64,
}

impl StationaryAvellaneda {
    /// Creates a stationary Avellaneda model.
    ///
    /// # Examples
    ///
    /// ```
    /// use solver::models::stationary_avellaneda::StationaryAvellaneda;
    /// let model = StationaryAvellaneda::new(0.1, 1.5, 140.0, 0.1);
    /// assert_eq!(model.gamma, 0.1);
    /// ```
    pub fn new(gamma: f64, kappa: f64, a: f64, sigma: f64) -> Self {
        Self {
            gamma,
            kappa,
            a,
            sigma,
        }
    }

    /// The exact base half-spread at zero inventory.
    pub fn base_spread(&self) -> f64 {
        (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln()
    }

    /// The `alpha` diagonal coefficient of the GLT operator.
    pub fn alpha(&self) -> f64 {
        (self.kappa / 2.0) * self.gamma * self.sigma.powi(2)
    }

    /// The `eta` off-diagonal coefficient of the GLT operator.
    pub fn eta(&self) -> f64 {
        self.a * (1.0 + self.gamma / self.kappa).powf(-(1.0 + self.kappa / self.gamma))
    }

    /// Builds the symmetric GLT operator with the Perron sign convention.
    fn glt_operator(&self, grid: &Grid<1>) -> CsrMatrix {
        let n = grid.total_size();
        let mut mat = CsrMatrix::new(n, n * 3);
        let q0 = grid.min[0].round() as i32;
        let alpha = self.alpha();
        let eta = self.eta();

        for i in 0..n {
            let q = q0 + i as i32;
            if i > 0 {
                mat.add_entry(i - 1, eta);
            }
            mat.add_entry(i, -alpha * (q as f64).powi(2));
            if i + 1 < n {
                mat.add_entry(i + 1, eta);
            }
            mat.finish_row();
        }

        mat
    }
}

impl StationaryEigenProblem<1> for StationaryAvellaneda {
    fn operator(&self, grid: &Grid<1>) -> CsrMatrix {
        self.glt_operator(grid)
    }
}

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

    fn default_model() -> StationaryAvellaneda {
        StationaryAvellaneda::new(0.1, 1.5, 140.0, 0.1)
    }

    #[test]
    fn base_spread_matches_reference() {
        let model = default_model();
        let base = (1.0 / model.gamma) * (1.0 + model.gamma / model.kappa).ln();
        assert!((model.base_spread() - base).abs() < 1e-12);
    }

    #[test]
    fn glt_coefficients_match_reference() {
        let model = default_model();
        let alpha = (model.kappa / 2.0) * model.gamma * model.sigma.powi(2);
        let eta =
            model.a * (1.0 + model.gamma / model.kappa).powf(-(1.0 + model.kappa / model.gamma));
        assert!((model.alpha() - alpha).abs() < 1e-12);
        assert!((model.eta() - eta).abs() < 1e-12);
    }

    #[test]
    fn operator_is_symmetric_tridiagonal() {
        let model = default_model();
        let grid = Grid::<1>::new([9], [-4.0], [4.0]);
        let mat = model.operator(&grid);

        assert_eq!(mat.size, 9);
        for i in 0..mat.size {
            let start = mat.row_ptr[i];
            let end = mat.row_ptr[i + 1];
            assert!(end - start <= 3, "row {i} has too many entries");
            for idx in start..end {
                let col = mat.col_indices[idx];
                assert!(
                    (col as i32 - i as i32).abs() <= 1,
                    "row {i} has non-tridiagonal column {col}"
                );
            }
        }
    }
}