solver/analytical/market_impact/

exact.rs

use crate::analytical::traits::AnalyticalSolution;
use crate::models::traits::ControlOutput;
use crate::numeric::ode::LinearSpectralSolver;

/// The exact analytical solution for the Avellaneda-Stoikov model with Market Impact.
/// Solves the system of ODEs for v_q(t).
pub struct AvellanedaImpactExact {
    pub gamma: f64,
    pub sigma: f64,
    pub kappa: f64,
    pub a: f64,
    pub xi: f64,  // Market Impact
    pub phi: f64, // Inventory Penalty
    pub terminal_time: f64,
    pub q_max: usize,
}

impl AvellanedaImpactExact {
    #[allow(clippy::too_many_arguments)]
    pub fn new(
        gamma: f64,
        sigma: f64,
        kappa: f64,
        a: f64,
        xi: f64,
        phi: f64,
        terminal_time: f64,
        q_max: usize,
    ) -> Self {
        Self {
            gamma,
            sigma,
            kappa,
            a,
            xi,
            phi,
            terminal_time,
            q_max,
        }
    }

    /// Helper to compute z(t) and alpha_trans parameter.
    ///
    /// Corresponds to Section "Market Impact" in `market_making.tex`.
    /// The text defines the system $\dot{\mathbf{v}}(t) = -\tilde{M}_\xi \mathbf{v}(t)$ (Eq \ref{eq:ode_impact}).
    /// To solve this using a symmetric spectral solver, we apply the transformation:
    /// $v_q(t) = \exp(\alpha_{trans} q^2) z_q(t)$ where $\alpha_{trans} = 0.5 k \xi$.
    /// This symmetrizes the transition matrix.
    fn compute_z_and_alpha(&self, t: f64) -> (Vec<f64>, f64) {
        let n = 2 * self.q_max + 1;
        let time_remaining = self.terminal_time - t;
        let alpha_trans = 0.5 * self.kappa * self.xi; // Transformation parameter

        if time_remaining <= 1e-9 {
            // Terminal condition from text: v(T, q) = \exp(-0.5 * k * xi * q^2).
            // Transformation: z(T, q) = \exp(-\alpha_{trans} q^2) v(T, q)
            // = \exp(-0.5 * k * xi * q^2) * \exp(-0.5 * k * xi * q^2)
            // = \exp(-k * xi * q^2).

            let mut z_t = vec![0.0; n];
            for (i, item) in z_t.iter_mut().enumerate() {
                let q_val = (i as i32) - (self.q_max as i32);
                let q = q_val as f64;
                *item = (-self.kappa * self.xi * q.powi(2)).exp();
            }
            return (z_t, alpha_trans);
        }

        // alpha_risk = k * (0.5 * gamma * sigma^2 + phi) (corresponds to \tilde{\alpha} in text)
        let alpha_risk = self.kappa * (0.5 * self.gamma * self.sigma.powi(2) + self.phi);
        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];

        // Symmetrization of the Matrix \tilde{M}_\xi.
        // The text's matrix has off-diagonals \eta e^{k\xi(q-1)} and \eta e^{-k\xi(q+1)}.
        // The transformation $v_q = e^{\alpha_{trans} q^2} z_q$ yields a symmetric system for z
        // with constant off-diagonal: \eta e^{-0.5 * k * xi}.
        let off_diag_val = eta * (-alpha_trans).exp();

        // Note: The text defines \dot{v} = -\tilde{M} v.
        // Our solver computes \exp(A \tau).
        // For the equation \dot{v} = -\tilde{M} v, the solution is v(t) = \exp(\tilde{M}(T-t)) v(T).
        // So we construct A = \tilde{M}_{sym} directly.

        for i in 0..n {
            let q_val = (i as i32) - (self.q_max as i32);
            let q = q_val as f64;
            // Matrix Diagonal from text: -\tilde{\alpha} q^2
            d[i] = -alpha_risk * q.powi(2);

            if i < n - 1 {
                e[i] = off_diag_val;
            }
        }

        // Initial Condition z(T)
        // z(T) = exp(-k * xi * q^2)
        let mut z0 = vec![0.0; n];
        for (i, item) in z0.iter_mut().enumerate() {
            let q_val = (i as i32) - (self.q_max as i32);
            let q = q_val as f64;
            *item = (-self.kappa * self.xi * q.powi(2)).exp();
        }

        // Solve for z(tau) = exp(M_sym * tau) * z0
        let z_t = LinearSpectralSolver::solve_tridiagonal(d, e, &z0, time_remaining);

        (z_t, alpha_trans)
    }

    /// Solves for the optimal spreads at time t given inventory q.
    pub fn exact_spreads(&self, t: f64, q: f64) -> (f64, f64) {
        let (z_t, alpha_trans) = self.compute_z_and_alpha(t);
        let n = z_t.len();

        // Check for validity
        for val in &z_t {
            if val.is_nan() {
                return (f64::NAN, f64::NAN);
            }
        }

        let const_term = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();
        let q_idx = (q.round() as i32 + self.q_max as i32) as usize;

        if q_idx == 0 || q_idx >= n - 1 {
            return (999.0, 999.0);
        }

        // Clamp values to be positive to avoid NaN from log(<=0) due to numerical noise/underflow
        let z_q = z_t[q_idx].max(1e-50);
        let z_qp1 = z_t[q_idx + 1].max(1e-50);
        let z_qm1 = z_t[q_idx - 1].max(1e-50);

        // ln(v_q / v_{q+1}) = ln(z_q / z_{q+1}) + alpha_trans * (q^2 - (q+1)^2)
        // = ln(z_q) - ln(z_{q+1}) - alpha_trans * (2q + 1)
        let ln_ratio_bid = z_q.ln() - z_qp1.ln() - alpha_trans * (2.0 * q + 1.0);

        // ln(v_q / v_{q-1}) = ln(z_q / z_{q-1}) + alpha_trans * (q^2 - (q-1)^2)
        // = ln(z_q) - ln(z_{q-1}) + alpha_trans * (2q - 1)
        let ln_ratio_ask = z_q.ln() - z_qm1.ln() + alpha_trans * (2.0 * q - 1.0);

        let bid = (1.0 / self.kappa) * ln_ratio_bid + (q + 1.0) * self.xi + const_term;
        let ask = (1.0 / self.kappa) * ln_ratio_ask - (q - 1.0) * self.xi + const_term;

        (bid, ask)
    }
}

impl AnalyticalSolution<2> for AvellanedaImpactExact {
    fn value_function(&self, t: f64, state: &[f64; 2]) -> f64 {
        let q = state[0];
        let (z_t, alpha_trans) = self.compute_z_and_alpha(t);
        let q_idx = (q.round() as i32 + self.q_max as i32) as usize;

        if q_idx >= z_t.len() {
            return 0.0;
        }

        // theta(t,q) = (1/kappa) * ln(v(t,q))
        // ln(v) = alpha_trans * q^2 + ln(z)
        (1.0 / self.kappa) * (alpha_trans * q.powi(2) + z_t[q_idx].ln())
    }

    fn optimal_controls(&self, t: f64, state: &[f64; 2]) -> ControlOutput<2> {
        let q = state[0];
        let (d_bid, d_ask) = self.exact_spreads(t, q);

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

        let flow = lambda_bid * d_bid + lambda_ask * d_ask;

        ControlOutput {
            lambda_plus: [lambda_bid, 0.0],
            lambda_minus: [lambda_ask, 0.0],
            flow,
        }
    }
}