End-to-End Example

This page walks through the full pipeline used in the numerical experiments: solving optimal bid/ask spreads via finite-difference policy iteration on a Hawkes-driven market with order flow imbalance, then backtesting against Gaussian price paths with stochastic matching.

Overview

Three crates, three stages:

solver: solve the HJB PDE on a 3D grid over [q, lambda_plus, lambda_minus] using implicit Euler policy iteration. Extract bid and ask half-spreads from the value-function ratios and store them as 4D lookup tables (adding tau as the fourth axis).
market_model: provide the price process (GBM).
engine: drive a SimulatedDataSource (GBM + bilateral Hawkes matcher) through an Engine, run the strategy, and collect backtest metrics including terminal liquidation.

The model is BilateralHawkesOrderFlowImbalance (OFI). It extends the standard bilateral Hawkes model with a price-impact parameter xi such that each buy market order (filling the ask) raises the mid-price by xi and each sell market order (filling the bid) lowers it by xi. The market maker's optimal spreads account for this adverse selection.

Stage 1: Solve the HJB

The state space is [q, lambda_plus, lambda_minus]. The grid spans inventory [-10, 10] with 21 points (step dq = 1.0), and both intensity axes cover [0, 8] with 11 points. The terminal condition uses LiquidationCost so the solver knows the value of residual inventory at expiry.

#![allow(unused)]
fn main() {
use solver::core::grid::Grid;
use solver::models::bilateral_hawkes_order_flow_imbalance::{
    BilateralHawkesOrderFlowImbalance, TerminalCondition,
};
use solver::numeric::finite_difference::solver::{PolicyIterationSolver, Scheme};

let model = BilateralHawkesOrderFlowImbalance::new(
    0.5,   // gamma  - risk aversion
    0.02,  // sigma  - price volatility
    1.5,   // kappa  - fill-rate decay (AS kappa)
    7.2,   // alpha  - Hawkes self-excitation (rho = alpha / beta = 0.9)
    8.0,   // beta   - Hawkes mean-reversion speed
    2.0,   // mu     - baseline intensity
    0.05,  // xi     - price impact per fill (eta_ofi)
)
.with_terminal_condition(TerminalCondition::LiquidationCost)
.with_terminal_liquidation_half_spread(0.1)
.with_inventory_bounds(-10.0, 10.0)
.with_dq(1.0)
.with_lambda_step(50.0);   // 50 → grid points align with mu=2.0 at index 0

let grid = Grid::<3>::new(
    [21, 11, 11],
    [-10.0, 0.0, 0.0],
    [10.0, 500.0, 500.0],
);

let dt = 1.0 / 200.0;      // 200 time steps
let n_steps = 200;
let solver = PolicyIterationSolver::new(dt).with_scheme(Scheme::Implicit);

// solve_all_steps returns V at tau = 0, tau = dt, ..., tau = T.
let v_history = solver.solve_all_steps(&grid, &model, n_steps);
}

After solving, build the 4D lookup tables from the value function:

#![allow(unused)]
fn main() {
use solver::lookup::{grid::linspace, LookupGrid, LookupTable};

fn build_spread_tables(
    grid: &Grid<3>, model: &BilateralHawkesOrderFlowImbalance,
    v_history: &[Vec<f64>], dt: f64,
) -> (LookupTable, LookupTable) {
    let n_q  = grid.tensor_info.shape[0];
    let n_l  = grid.tensor_info.shape[1];
    let n_t  = v_history.len() - 1;
    let dq   = grid.dx[0];
    let dl   = grid.dx[1];
    let base_spread = (1.0 / model.gamma) * (1.0 + model.gamma / model.kappa).ln();

    let q_axis   = linspace(grid.min[0], grid.max[0], n_q);
    let lp_axis  = linspace(grid.min[1], grid.max[1], n_l);
    let lm_axis  = linspace(grid.min[2], grid.max[2], n_l);
    let tau_axis: Vec<f64> = (0..=n_t).map(|k| k as f64 * dt).collect();

    let axes  = vec![q_axis, lp_axis.clone(), lm_axis, tau_axis];
    let names = vec!["q".into(), "lambda_plus".into(), "lambda_minus".into(), "tau".into()];

    let compute = |coords: &[f64], sign: f64| -> f64 {
        let qi = ((coords[0] - grid.min[0]) / dq).round().clamp(0.0, (n_q-1) as f64) as usize;
        let lpi = ((coords[1] - grid.min[1]) / dl).round().clamp(0.0, (n_l-1) as f64) as usize;
        let lmi = ((coords[2] - grid.min[2]) / dl).round().clamp(0.0, (n_l-1) as f64) as usize;
        let ti  = ((coords[3] / dt).round() as usize).min(n_t);

        let idx = grid.tensor_info.linear_index(&[qi, lpi, lmi]);
        let v   = &v_history[ti];
        let v_q = v[idx].max(1e-300);

        let nb_idx = if sign > 0.0 {
            if qi > 0 { idx - 1 } else { idx }
        } else {
            if qi + 1 < n_q { idx + 1 } else { idx }
        };
        ((1.0 / model.kappa) * (v_q / v[nb_idx].max(1e-300)).ln() + base_spread).max(0.001)
    };

    let bid = LookupGrid::new(axes.clone(), names.clone())
        .generate(|c| compute(c, 1.0));   // V(q) / V(q-1) for ask side
    let ask = LookupGrid::new(axes, names)
        .generate(|c| compute(c, -1.0));  // V(q) / V(q+1) for bid side

    (bid, ask)
}
}

Stage 2: Wire up the market

A GBM drives the mid-price. The matcher uses bilateral Hawkes with the same alpha/beta as the model for consistency, and the data source applies price impact so impact_factor = -xi.

#![allow(unused)]
fn main() {
use std::sync::Arc;
use market_model::process::gbm::GeometricBrownianMotion;
use engine::data_source::SimulatedDataSource;
use engine::matcher::StochasticMatcher;

let sigma = 0.02;
let dt    = 1.0 / 1000.0;   // 1 ms ticks
let n_sim = 2000;           // 2 seconds of trading

let gbm = GeometricBrownianMotion::new(0.0, sigma, 100.0);
let source = SimulatedDataSource::new(
    gbm, 100.0, dt, 0.05, n_sim, -0.05,  // impact_factor = -xi
);

// k = distance decay, a = base fill rate.
let matcher = StochasticMatcher::new(dt, 50.0, 100.0)
    .with_bilateral_hawkes(7.2, 8.0);    // alpha, beta
}

Stage 3: Backtest

Wrap the lookup tables in the OFI strategy and run through the engine. Use run_backtest_with_liquidation to apply the same terminal liquidation cost that was baked into the FDM value function.

#![allow(unused)]
fn main() {
use engine::engine::Engine;
use engine::backtest::run_backtest_with_liquidation;
use engine::strategies::avellaneda_stoikov_bilateral_hawkes_order_flow_imbalance::{
    AvellanedaStoikovBilateralHawkesOrderFlowImbalanceParameters,
    AvellanedaStoikovBilateralHawkesOrderFlowImbalanceStrategy,
};

let params = AvellanedaStoikovBilateralHawkesOrderFlowImbalanceParameters {
    t_horizon: 1.0,
    base_intensity_buy: 2.0,
    base_intensity_sell: 2.0,
};

let strategy = AvellanedaStoikovBilateralHawkesOrderFlowImbalanceStrategy::new(
    params,
    Arc::new(bid_table),
    Arc::new(ask_table),
);

let mut engine = Engine::new(
    matcher, strategy, source,
    10_000.0,  // initial cash
    0.0,       // transaction cost rate (zero for clean spread PnL)
).with_ground_truth(false);   // opaque filter: strategy only sees BBO + portfolio

let result = run_backtest_with_liquidation(
    &mut engine,
    dt,           // dt in years
    0.1,          // terminal liquidation half-spread
);

println!("Sharpe:        {:.2}", result.sharpe_ratio);
println!("Sortino:       {:.2}", result.sortino_ratio);
println!("PnL spread:    {:.4}", result.pnl_spread);
println!("PnL dir:       {:.4}", result.pnl_dir);
println!("Trades:        {}", result.total_trades);
println!("Mean hold:     {:.2} steps", result.mean_hold_time);
println!("Term liq cost: {:.4}", result.terminal_liquidation_cost);
}

What the pipeline demonstrates