solver/analytical/avellaneda/
exact.rs

1use crate::analytical::traits::AnalyticalSolution;
2use crate::models::TerminalCondition;
3use crate::models::traits::ControlOutput;
4use crate::numeric::ode::LinearSpectralSolver;
5
6/// The exact analytical solution for the Avellaneda-Stoikov model using the spectral method.
7/// Solves the system of ODEs for theta(t, q).
8pub struct AvellanedaExact {
9    pub gamma: f64,
10    pub sigma: f64,
11    pub kappa: f64,
12    pub a: f64,
13    pub terminal_time: f64,
14    pub q_max: usize,
15    /// Terminal condition for the value function at T.  Defaults to Zero.
16    pub terminal_condition: TerminalCondition,
17}
18
19impl AvellanedaExact {
20    pub fn new(
21        gamma: f64,
22        sigma: f64,
23        kappa: f64,
24        a: f64,
25        terminal_time: f64,
26        q_max: usize,
27    ) -> Self {
28        Self {
29            gamma,
30            sigma,
31            kappa,
32            a,
33            terminal_time,
34            q_max,
35            terminal_condition: TerminalCondition::Zero,
36        }
37    }
38
39    pub fn with_terminal_condition(mut self, terminal_condition: TerminalCondition) -> Self {
40        self.terminal_condition = terminal_condition;
41        self
42    }
43
44    /// Terminal value of v(T, q) under the log-transform V = (1/gamma)*ln(v).
45    ///
46    /// Under Zero: V(T,q)=0 => v(T,q)=1.
47    /// Under LiquidationCost: V(T,q)=-|q|*base_spread => v(T,q)=exp(gamma*(-|q|*base_spread)).
48    fn terminal_v(&self, q: i32) -> f64 {
49        match self.terminal_condition {
50            TerminalCondition::Zero => 1.0,
51            TerminalCondition::LiquidationCost => {
52                let base_spread = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();
53                (self.gamma * (-(q.abs() as f64) * base_spread)).exp()
54            }
55        }
56    }
57
58    /// Compute v(t, .) vector: the solution to the linearised HJB system at time t.
59    /// Returns a vector of length 2*q_max+1, indexed by q+q_max.
60    pub fn compute_v(&self, t: f64) -> Vec<f64> {
61        let time_remaining = self.terminal_time - t;
62        let n = 2 * self.q_max + 1;
63        if time_remaining <= 1e-9 {
64            return (0..n)
65                .map(|i| {
66                    let q = (i as i32) - (self.q_max as i32);
67                    self.terminal_v(q)
68                })
69                .collect();
70        }
71
72        let n = 2 * self.q_max + 1;
73        // The text defines \tilde{\alpha} = (k/2)*gamma*sigma^2 + k*phi.
74        // Assuming phi = 0 for pure Avellaneda:
75        let alpha = (self.kappa / 2.0) * self.gamma * self.sigma.powi(2);
76        let eta = self.a * (1.0 + self.gamma / self.kappa).powf(-(1.0 + self.kappa / self.gamma));
77
78        let mut d = vec![0.0; n];
79        let mut e = vec![0.0; n];
80
81        // We solve v(t) = exp(\tilde{M} * (T-t)) * v(T).
82        // The matrix \tilde{M} has diagonal -\tilde{\alpha} q^2 and off-diagonal \eta.
83
84        for i in 0..n {
85            let q_val = (i as i32) - (self.q_max as i32);
86            d[i] = -alpha * (q_val as f64).powi(2);
87            if i < n - 1 {
88                e[i] = eta;
89            }
90        }
91
92        let y0: Vec<f64> = (0..n)
93            .map(|i| {
94                let q = (i as i32) - (self.q_max as i32);
95                self.terminal_v(q)
96            })
97            .collect();
98        // Propagate forward by time_remaining (T-t)
99        LinearSpectralSolver::solve_tridiagonal(d, e, &y0, time_remaining)
100    }
101
102    /// Solves for the optimal spreads at time t given current inventory q.
103    pub fn exact_spreads(&self, t: f64, q: f64) -> (f64, f64) {
104        let v_t = self.compute_v(t);
105        let const_term = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();
106
107        let n = v_t.len();
108        let q_idx = (q.round() as i32 + self.q_max as i32) as usize;
109
110        if q_idx == 0 || q_idx >= n - 1 {
111            return (999.0, 999.0);
112        }
113
114        let bid = (1.0 / self.kappa) * (v_t[q_idx] / v_t[q_idx + 1]).ln() + const_term;
115        let ask = (1.0 / self.kappa) * (v_t[q_idx] / v_t[q_idx - 1]).ln() + const_term;
116
117        (bid, ask)
118    }
119
120    /// Computes the value function u(t, x, q, S) = -exp(-gamma * (x + qS + theta(t, q)))
121    /// Returns theta(t, q).
122    pub fn value_function_theta(&self, t: f64, q: f64) -> f64 {
123        let v_t = self.compute_v(t);
124        let q_idx = (q.round() as i32 + self.q_max as i32) as usize;
125
126        if q_idx >= v_t.len() {
127            return 0.0;
128        }
129
130        // The transformation is v(t,q) = exp(kappa * theta(t,q)).
131        // Thus, theta(t,q) = (1/kappa) * ln(v(t,q)).
132        (1.0 / self.kappa) * v_t[q_idx].ln()
133    }
134}
135
136impl AnalyticalSolution<2> for AvellanedaExact {
137    fn value_function(&self, t: f64, state: &[f64; 2]) -> f64 {
138        self.value_function_theta(t, state[0])
139    }
140
141    fn optimal_controls(&self, t: f64, state: &[f64; 2]) -> ControlOutput<2> {
142        let q = state[0];
143        let (d_bid, d_ask) = self.exact_spreads(t, q);
144
145        let lambda_bid = self.a * (-self.kappa * d_bid).exp();
146        let lambda_ask = self.a * (-self.kappa * d_ask).exp();
147
148        let flow = lambda_bid * d_bid + lambda_ask * d_ask;
149
150        ControlOutput {
151            lambda_plus: [lambda_bid, 0.0],
152            lambda_minus: [lambda_ask, 0.0],
153            flow,
154        }
155    }
156}
solver/analytical/avellaneda/exact.rs

solver/analytical/avellaneda/
exact.rs
