solver/models/
bilateral_hawkes.rs

1use super::normal_cdf;
2use super::traits::{ControlOutput, Gradients, Model};
3
4// Re-export so existing `hjb_solver::models::bilateral_hawkes::TerminalCondition` paths compile.
5pub use super::TerminalCondition;
6
7/// Avellaneda-Stoikov model with a bilateral (2-D) Hawkes process.
8///
9/// # State Space  `[q, lambda_plus, lambda_minus]`
10///
11/// | dim | symbol   | description                                      |
12/// |-----|----------|--------------------------------------------------|
13/// | 0   | q        | Inventory level (integer, jump-controlled)       |
14/// | 1   | lambda+  | Buy market-order intensity (ODE drift, no noise) |
15/// | 2   | lambda-  | Sell market-order intensity (ODE drift, no noise)|
16///
17/// # Dynamics
18///
19/// `d(lambda+) = beta*(mu - lambda+) dt + alpha dN+`
20/// `d(lambda-) = beta*(mu - lambda-) dt + alpha dN-`
21///
22/// `N+` counts buy market orders (rate `lambda+`), `N-` counts sell market orders (rate `lambda-`).
23/// Each side self-excites independently — there is no cross-excitation.
24///
25/// # Market-making mapping
26///
27/// | event              | driven by | inventory effect |
28/// |--------------------|-----------|-----------------|
29/// | Sell MO hits bid   | lambda-   | q -> q + 1      |
30/// | Buy  MO hits ask   | lambda+   | q -> q - 1      |
31///
32/// Effective fill rates:
33/// * `Lambda_bid = lambda- * exp(-kappa * delta_bid)`  (sell MO hits MM bid)
34/// * `Lambda_ask = lambda+ * exp(-kappa * delta_ask)`  (buy  MO hits MM ask)
35///
36/// # Hawkes approximation used in the FD solver
37///
38/// The lambda dimensions use the fluid (mean-field) approximation:
39///
40///   `d(lambda+)/dt ~ beta*(mu - lambda+) + alpha * lambda+`
41///   `d(lambda-)/dt ~ beta*(mu - lambda-) + alpha * lambda-`
42///
43/// This is identical to the treatment in `AvellanedaHawkes` and ensures the
44/// stationary mean of each intensity sits at `mu / (1 - alpha/beta)` in the
45/// value-function PDE.  Using discrete fill-triggered jumps instead requires
46/// `lambda_step << alpha` for accuracy; at the grid resolution affordable for
47/// 4-D tables that condition is violated, causing the PDE to see a much lower
48/// effective mean intensity than the actual simulation (`mu` vs `mu/(1-rho)`),
49/// which corrupts the value function and the extracted spreads.
50#[derive(Clone, Debug)]
51pub struct BilateralHawkes {
52    /// Risk-aversion parameter.
53    pub gamma: f64,
54    /// Price volatility.
55    pub sigma: f64,
56    /// Fill-rate decay (Avellaneda-Stoikov kappa).
57    pub kappa: f64,
58    // --- Hawkes parameters (shared for both sides) ---
59    /// Jump size added to the excited side's intensity.
60    pub alpha: f64,
61    /// Mean-reversion speed.
62    pub beta: f64,
63    /// Baseline intensity (mean-reversion target).
64    pub mu: f64,
65    /// Grid step for both lambda+ and lambda- dimensions.
66    pub lambda_step: f64,
67    /// Grid step for the inventory dimension.
68    pub dq: f64,
69    /// Minimum inventory (lower hard boundary). Defaults to -infinity.
70    pub q_min: f64,
71    /// Maximum inventory (upper hard boundary). Defaults to +infinity.
72    pub q_max: f64,
73    /// Terminal condition for the value function at T. Defaults to Zero.
74    pub terminal_condition: TerminalCondition,
75    /// Optional terminal liquidation half-spread used when terminal_condition is
76    /// LiquidationCost. If None, uses the AS base spread.
77    pub terminal_liquidation_half_spread: Option<f64>,
78}
79
80impl BilateralHawkes {
81    pub fn new(gamma: f64, sigma: f64, kappa: f64, alpha: f64, beta: f64, mu: f64) -> Self {
82        Self {
83            gamma,
84            sigma,
85            kappa,
86            alpha,
87            beta,
88            mu,
89            lambda_step: 1.0,
90            dq: 1.0,
91            q_min: f64::NEG_INFINITY,
92            q_max: f64::INFINITY,
93            terminal_condition: TerminalCondition::Zero,
94            terminal_liquidation_half_spread: None,
95        }
96    }
97
98    pub fn with_terminal_condition(mut self, terminal_condition: TerminalCondition) -> Self {
99        self.terminal_condition = terminal_condition;
100        self
101    }
102
103    pub fn with_terminal_liquidation_half_spread(mut self, half_spread: f64) -> Self {
104        self.terminal_liquidation_half_spread = Some(half_spread.max(0.0));
105        self
106    }
107
108    pub fn with_inventory_bounds(mut self, q_min: f64, q_max: f64) -> Self {
109        self.q_min = q_min;
110        self.q_max = q_max;
111        self
112    }
113
114    pub fn with_lambda_step(mut self, lambda_step: f64) -> Self {
115        self.lambda_step = lambda_step.abs().max(1e-8);
116        self
117    }
118
119    pub fn with_dq(mut self, dq: f64) -> Self {
120        self.dq = dq.abs().max(1e-8);
121        self
122    }
123
124    /// Computes optimal bid/ask half-spreads from value-function gradients.
125    ///
126    /// Bid spread: driven by the inventory gradient in the +q direction (sell MO fills).
127    /// Ask spread: driven by the inventory gradient in the -q direction (buy  MO fills).
128    pub fn get_spreads(&self, grads: &Gradients<3>) -> (f64, f64) {
129        let base_spread = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();
130        let delta_bid = (base_spread - grads.fwd[0]).clamp(-10.0, 10.0);
131        let delta_ask = (base_spread + grads.bwd[0]).clamp(-10.0, 10.0);
132        (delta_bid, delta_ask)
133    }
134
135    /// Mean-reversion drift part only (without Hawkes jump term).
136    /// Retained as a utility; currently unused within this module.
137    #[allow(dead_code)]
138    #[inline]
139    fn intensity_mean_reversion_drift(&self, lambda: f64) -> f64 {
140        self.beta * (self.mu - lambda)
141    }
142
143    /// Upwind transport rates for one intensity dimension.
144    #[inline]
145    fn intensity_upwind(&self, drift: f64) -> (f64, f64) {
146        let rate = drift.abs() / self.lambda_step;
147        if drift >= 0.0 {
148            (rate, 0.0)
149        } else {
150            (0.0, rate)
151        }
152    }
153}
154
155impl Model<3> for BilateralHawkes {
156    type Process = ();
157
158    fn process(&self) {}
159
160    fn optimize(&self, state: &[f64; 3], grads: &Gradients<3>) -> ControlOutput<3> {
161        let q = state[0];
162        let lambda_plus = state[1].max(0.0);
163        let lambda_minus = state[2].max(0.0);
164
165        // --- Dimension 0: Inventory controls ---
166        let (d_bid, d_ask) = self.get_spreads(grads);
167
168        let max_mult = 20.0_f64;
169        // Bid filled by sell MOs (lambda-); ask filled by buy MOs (lambda+)
170        let lambda_bid_fill = lambda_minus * (-self.kappa * d_bid).exp().min(max_mult);
171        let lambda_ask_fill = lambda_plus * (-self.kappa * d_ask).exp().min(max_mult);
172
173        // Shut off forbidden quoting side at inventory bounds.
174        let lambda_bid_fill = if q >= self.q_max {
175            0.0
176        } else {
177            lambda_bid_fill
178        };
179        let lambda_ask_fill = if q <= self.q_min {
180            0.0
181        } else {
182            lambda_ask_fill
183        };
184
185        let hamiltonian_inv = (lambda_bid_fill + lambda_ask_fill) / (self.gamma + self.kappa);
186        let risk_penalty = -0.5 * self.gamma * self.sigma.powi(2) * q.powi(2);
187
188        // Cancel the solver's q-transport contribution so we can inject the Hamiltonian.
189        let drift_correction_q =
190            self.dq * (lambda_bid_fill * grads.fwd[0] - lambda_ask_fill * grads.bwd[0]);
191
192        // --- Dimensions 1/2: lambda+ and lambda- dynamics ---
193        // Fluid (mean-field) approximation: replace discrete fill-triggered jumps with a
194        // continuous drift term alpha * lambda.  This gives the same stationary mean as the
195        // true Hawkes process (lambda* = mu / (1-rho)) and avoids the large discretization
196        // error that arises when the FD grid step (lambda_step) is much larger than alpha.
197        // Both sides remain jointly optimised through the 4-D value function.
198        let net_drift_plus = self.beta * (self.mu - lambda_plus) + self.alpha * lambda_plus;
199        let net_drift_minus = self.beta * (self.mu - lambda_minus) + self.alpha * lambda_minus;
200        let (lp_plus, lp_minus) = self.intensity_upwind(net_drift_plus);
201        let (lm_plus, lm_minus) = self.intensity_upwind(net_drift_minus);
202
203        ControlOutput {
204            lambda_plus: [lambda_bid_fill, lp_plus, lm_plus],
205            lambda_minus: [lambda_ask_fill, lp_minus, lm_minus],
206            flow: hamiltonian_inv + risk_penalty - drift_correction_q,
207        }
208    }
209
210    fn terminal(&self, state: &[f64; 3]) -> f64 {
211        match self.terminal_condition {
212            TerminalCondition::Zero => 0.0,
213            TerminalCondition::LiquidationCost => {
214                let q = state[0];
215                let base_spread = (1.0 / self.gamma) * (1.0 + self.gamma / self.kappa).ln();
216                let half_spread = self.terminal_liquidation_half_spread.unwrap_or(base_spread);
217                -q.abs() * half_spread
218            }
219        }
220    }
221
222    fn constant_discount_rate(&self) -> Option<f64> {
223        Some(0.0)
224    }
225
226    /// Forward step using the fluid (mean-field) drift, consistent with `optimize()`.
227    ///
228    /// The drift is `beta*(mu - lambda) + alpha * lambda` (same approximation used in
229    /// the PDE so that the forward simulation matches the value-function's stationary
230    /// lambda regime at `mu / (1 - alpha/beta)`).
231    fn next_step(&self, current_state: &[f64; 3], dt: f64, _noise: &[f64; 3]) -> [f64; 3] {
232        let lambda_plus = current_state[1];
233        let lambda_minus = current_state[2];
234        let drift_plus = self.beta * (self.mu - lambda_plus) + self.alpha * lambda_plus;
235        let drift_minus = self.beta * (self.mu - lambda_minus) + self.alpha * lambda_minus;
236        [
237            current_state[0],
238            (lambda_plus + drift_plus * dt).max(0.0),
239            (lambda_minus + drift_minus * dt).max(0.0),
240        ]
241    }
242
243    /// Coupled step: inventory jumps based on current fill probabilities.
244    fn next_step_controlled(
245        &self,
246        current_state: &[f64; 3],
247        control: &ControlOutput<3>,
248        dt: f64,
249        noise: &[f64; 3],
250    ) -> [f64; 3] {
251        let mut next = self.next_step(current_state, dt, noise);
252
253        // Bernoulli inventory jump using noise[0] mapped through CDF -> uniform
254        let u = normal_cdf(noise[0]);
255        let p_bid = (control.lambda_plus[0] * dt).clamp(0.0, 1.0);
256        let p_ask = (control.lambda_minus[0] * dt).clamp(0.0, 1.0);
257        if u < p_bid {
258            next[0] += 1.0; // sell MO hits bid  -> inventory +1
259            next[2] += self.alpha; // lambda- self-excites on sell MO fill
260        } else if u > 1.0 - p_ask {
261            next[0] -= 1.0; // buy  MO hits ask  -> inventory -1
262            next[1] += self.alpha; // lambda+ self-excites on buy MO fill
263        }
264
265        next[1] = next[1].max(0.0);
266        next[2] = next[2].max(0.0);
267
268        next
269    }
270
271    /// Both intensity dimensions follow deterministic ODE drift with no Brownian noise.
272    fn is_diffusion_dimension(&self, _dim: usize) -> bool {
273        false // no dimension has Brownian noise in this model
274    }
275
276    fn is_integer_dimension(&self, dim: usize) -> bool {
277        dim == 0 // q is discrete integer; lambda_bid and lambda_ask are continuous
278    }
279
280    fn gradient_step(&self, dim: usize) -> f64 {
281        match dim {
282            1 | 2 => self.lambda_step.abs().max(1e-8),
283            _ => 1.0,
284        }
285    }
286
287    /// The base arrival intensity is state-dependent: returns the average of
288    /// the buy/sell intensity states for intensity-to-spread conversion.
289    fn fill_rate_base(&self, state: &[f64; 3]) -> f64 {
290        ((state[1] + state[2]) * 0.5).max(1e-10)
291    }
292
293    fn fill_rate_decay(&self) -> f64 {
294        self.kappa
295    }
296}
297
298#[cfg(test)]
299mod tests {
300    use super::*;
301
302    fn zero_grads() -> Gradients<3> {
303        Gradients {
304            fwd: [0.0; 3],
305            bwd: [0.0; 3],
306        }
307    }
308
309    fn default_model() -> BilateralHawkes {
310        BilateralHawkes::new(0.5, 0.2, 1.5, 0.5, 2.0, 1.0)
311            .with_lambda_step(0.1)
312            .with_dq(1.0)
313    }
314
315    #[test]
316    fn zero_gradient_gives_symmetric_base_spread() {
317        let m = default_model();
318        let (bid, ask) = m.get_spreads(&zero_grads());
319        let base = (1.0 / m.gamma) * (1.0 + m.gamma / m.kappa).ln();
320        assert!((bid - base).abs() < 1e-12);
321        assert!((ask - base).abs() < 1e-12);
322    }
323
324    #[test]
325    fn bid_driven_by_lambda_minus_ask_by_lambda_plus() {
326        // lambda+ high -> ask fill rate high, lambda- high -> bid fill rate high
327        let m = default_model();
328        let g = zero_grads();
329        // High lambda- -> high bid fill
330        let ctrl_high_minus = m.optimize(&[0.0, 1.0, 5.0], &g);
331        // Low  lambda- -> low  bid fill
332        let ctrl_low_minus = m.optimize(&[0.0, 1.0, 0.5], &g);
333        assert!(
334            ctrl_high_minus.lambda_plus[0] > ctrl_low_minus.lambda_plus[0],
335            "bid fill should scale with lambda-"
336        );
337
338        // High lambda+ -> high ask fill
339        let ctrl_high_plus = m.optimize(&[0.0, 5.0, 1.0], &g);
340        let ctrl_low_plus = m.optimize(&[0.0, 0.5, 1.0], &g);
341        assert!(
342            ctrl_high_plus.lambda_minus[0] > ctrl_low_plus.lambda_minus[0],
343            "ask fill should scale with lambda+"
344        );
345    }
346
347    #[test]
348    fn symmetric_state_gives_equal_bid_ask() {
349        let m = default_model();
350        let g = zero_grads();
351        let ctrl = m.optimize(&[0.0, 1.0, 1.0], &g);
352        assert!(
353            (ctrl.lambda_plus[0] - ctrl.lambda_minus[0]).abs() < 1e-12,
354            "bid fill == ask fill when lambda+ == lambda- and q == 0"
355        );
356    }
357
358    #[test]
359    fn independent_intensity_upwinding() {
360        let m = default_model();
361        let g = zero_grads();
362        // Low lambda+ (below mu=1): drift > 0 -> upward transport only
363        let ctrl = m.optimize(&[0.0, 0.1, 1.0], &g);
364        assert!(ctrl.lambda_plus[1] > 0.0, "lambda+ should drift up");
365        assert!(
366            ctrl.lambda_minus[1] < 1e-12,
367            "no downward lambda+ transport"
368        );
369        // High lambda-: the net fluid drift beta(mu-lambda-) + alpha lambda- is negative,
370        // so transport is purely downward in the lambda- dimension.
371        let ctrl2 = m.optimize(&[0.0, 1.0, 3.0], &g);
372        assert!(ctrl2.lambda_minus[2] > 0.0, "lambda- should drift down");
373        assert!(
374            ctrl2.lambda_plus[2] < 1e-12,
375            "no upward lambda- transport when net drift is negative"
376        );
377    }
378
379    #[test]
380    fn hawkes_jump_terms_raise_upward_lambda_rates() {
381        let m = default_model();
382        let g = zero_grads();
383        let ctrl = m.optimize(&[0.0, 1.0, 1.0], &g);
384
385        // lambda+ dimension upward rate gets Hawkes jump from ask fills.
386        assert!(ctrl.lambda_plus[1] > 0.0);
387        // lambda- dimension upward rate gets Hawkes jump from bid fills.
388        assert!(ctrl.lambda_plus[2] > 0.0);
389    }
390
391    #[test]
392    fn controlled_step_applies_side_specific_self_excitation() {
393        let m = default_model();
394        let ctrl = ControlOutput {
395            lambda_plus: [1.0, 0.0, 0.0],
396            lambda_minus: [0.0, 0.0, 0.0],
397            flow: 0.0,
398        };
399
400        // Force bid fill branch (u < p_bid) by choosing very negative noise.
401        let s0 = [0.0, 1.0, 1.0];
402        let s1 = m.next_step_controlled(&s0, &ctrl, 0.1, &[-10.0, 0.0, 0.0]);
403        assert!(s1[2] >= s0[2], "lambda- should increase on bid fill");
404    }
405
406    #[test]
407    fn forward_step_clamped_to_nonnegative() {
408        let m = default_model();
409        let state = [0.0, 0.0, 0.0];
410        let next = m.next_step(&state, 1.0, &[0.0; 3]);
411        assert!(next[1] >= 0.0);
412        assert!(next[2] >= 0.0);
413    }
414
415    #[test]
416    fn terminal_zero_condition_returns_zero() {
417        let m = default_model().with_terminal_condition(TerminalCondition::Zero);
418        let v = m.terminal(&[3.0, 1.0, 1.0]);
419        assert!(v.abs() < 1e-12);
420    }
421
422    #[test]
423    fn terminal_liquidation_cost_matches_formula() {
424        let m = default_model().with_terminal_condition(TerminalCondition::LiquidationCost);
425        let q: f64 = -4.0;
426        let base = (1.0 / m.gamma) * (1.0 + m.gamma / m.kappa).ln();
427        let expected = -q.abs() * base;
428        let got = m.terminal(&[q, 1.2, 0.8]);
429        assert!((got - expected).abs() < 1e-12);
430    }
431}
solver/models/bilateral_hawkes.rs

solver/models/
bilateral_hawkes.rs
