solver/analytical/drift/
exact.rs

1use crate::analytical::traits::AnalyticalSolution;
2use crate::models::traits::ControlOutput;
3use crate::numeric::ode::LinearSpectralSolver;
4
5/// The exact analytical solution for the Avellaneda-Stoikov model with Drift.
6/// Solves the system of ODEs for v_q(t).
7pub struct AvellanedaDriftExact {
8    pub gamma: f64,
9    pub sigma: f64,
10    pub kappa: f64,
11    pub a: f64,
12    pub mu: f64,  // Drift
13    pub phi: f64, // Inventory Penalty
14    pub terminal_time: f64,
15    pub q_max: usize,
16}
17
18impl AvellanedaDriftExact {
19    #[allow(clippy::too_many_arguments)]
20    pub fn new(
21        gamma: f64,
22        sigma: f64,
23        kappa: f64,
24        a: f64,
25        mu: f64,
26        phi: f64,
27        terminal_time: f64,
28        q_max: usize,
29    ) -> Self {
30        Self {
31            gamma,
32            sigma,
33            kappa,
34            a,
35            mu,
36            phi,
37            terminal_time,
38            q_max,
39        }
40    }
41
42    /// Solves for the optimal spreads at time t given inventory q.
43    pub fn exact_spreads(&self, t: f64, q: f64) -> (f64, f64) {
44        let time_remaining = self.terminal_time - t;
45        if time_remaining <= 1e-9 {
46            return (f64::INFINITY, f64::INFINITY);
47        }
48
49        let n = 2 * self.q_max + 1;
50        // alpha = k * (0.5 * gamma * sigma^2 + phi)
51        let alpha = self.kappa * (0.5 * self.gamma * self.sigma.powi(2) + self.phi);
52        let beta = self.kappa * self.mu;
53        let eta = self.a * (1.0 + self.gamma / self.kappa).powf(-(1.0 + self.kappa / self.gamma));
54
55        let mut d = vec![0.0; n];
56        let mut e = vec![0.0; n];
57
58        for i in 0..n {
59            let q_val = (i as i32) - (self.q_max as i32);
60            let q_f = q_val as f64;
61
62            // Corresponds to Section "Price Trend" in `market_making.tex`.
63            // The system is \dot{v} = -\tilde{M} v.
64            // Solution: v(t) = \exp(\tilde{M} (T-t)) v(T).
65            // We construct the matrix A = \tilde{M} directly to compute exp(A * time_remaining).
66
67            // Text Eq \ref{eq:ode_trend} defines \tilde{M}:
68            // Diagonal: -\tilde{\alpha} q^2 + k \mu q
69            // Off-diagonal: \eta
70
71            d[i] = -alpha * q_f.powi(2) + beta * q_f;
72            if i < n - 1 {
73                e[i] = eta;
74            }
75        }
76
77        // Terminal condition v(T) = 1.
78        let y0 = vec![1.0; n];
79
80        // Pass positive time_remaining to compute exp(M_code * tau)
81        let v_t = LinearSpectralSolver::solve_tridiagonal(d, e, &y0, time_remaining);
82
83        let const_term = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();
84        let q_idx = (q.round() as i32 + self.q_max as i32) as usize;
85
86        if q_idx == 0 || q_idx >= n - 1 {
87            return (999.0, 999.0);
88        }
89
90        let bid = (1.0 / self.kappa) * (v_t[q_idx] / v_t[q_idx + 1]).ln() + const_term;
91        let ask = (1.0 / self.kappa) * (v_t[q_idx] / v_t[q_idx - 1]).ln() + const_term;
92
93        (bid, ask)
94    }
95}
96
97impl AnalyticalSolution<2> for AvellanedaDriftExact {
98    fn value_function(&self, _t: f64, _state: &[f64; 2]) -> f64 {
99        0.0 // Not implemented
100    }
101
102    fn optimal_controls(&self, t: f64, state: &[f64; 2]) -> ControlOutput<2> {
103        let q = state[0];
104        let (d_bid, d_ask) = self.exact_spreads(t, q);
105
106        let lambda_bid = self.a * (-self.kappa * d_bid).exp();
107        let lambda_ask = self.a * (-self.kappa * d_ask).exp();
108
109        let flow = lambda_bid * d_bid + lambda_ask * d_ask;
110
111        ControlOutput {
112            lambda_plus: [lambda_bid, 0.0],
113            lambda_minus: [lambda_ask, 0.0],
114            flow,
115        }
116    }
117}
solver/analytical/drift/exact.rs

solver/analytical/drift/
exact.rs
