market-solve

A general-purpose framework for stochastic market simulation, optimal strategy computation, and multi-agent backtesting.

Goal

market-solve is a tool for anyone who needs to simulate financial markets, evaluate trading strategies, or compute optimal policies under uncertainty. The framework treats processes, models, strategies, and exchange mechanics as replaceable components that can be mixed and matched at compile time with zero runtime overhead.

What you can do with it

• Simulate price paths in batch: thousands of parallel Monte Carlo trajectories across 8 stochastic process families, from simple GBM to Hawkes-driven order flow.
• Compute optimal strategies via HJB: solve Hamilton-Jacobi-Bellman equations using finite-difference policy iteration or BSDE-based least-squares Monte Carlo, with dimension-agnostic implementations that compose multi-factor models (inventory + volatility + intensity).
• Backtest against a realistic exchange: run strategies inside a multi-agent engine with configurable order matching (simple, stochastic Poisson-based), exchange mechanics, observation filtering, and portfolio tracking.
• Train RL agents on vectorized environments: run hundreds of independent simulation episodes in parallel via rayon, with reward functions including PnL and differential Sharpe ratio.
• Generate synthetic order book data: produce BBO quotes from a composable price process. Full L2/L3 order book generation is planned.
• Call from Python: the engine exposes strategies and backtest infrastructure via PyO3 bindings.

Workspace structure

The framework is split into three crates with clear responsibility boundaries:

Dependency graph

graph BT
    market_model
    solver --> market_model
    engine --> market_model
    engine --> solver
    py_bindings --> engine

solver depends on market_model for process simulation and SimulationRunner. engine depends on market_model for processes and on solver for HJB-derived optimal strategies.

Quick start

git clone https://github.com/marci/market-solve
cd market-solve
cargo build --release
cargo test --workspace
cargo clippy --workspace -- -D warnings

Generate API documentation:

cargo doc --workspace --no-deps --open

Build the user guide (this book):

cargo install mdbook
cd docs && mdbook build

Example: simulation + backtest

#![allow(unused)]
fn main() {
use market_model::process::gbm::GeometricBrownianMotion;
use market_model::runner::SimulationRunner;
use market_model::types::SimulationConfig;

// 1. Simulate 10,000 GBM paths
let gbm = GeometricBrownianMotion::new(0.05, 0.2, 100.0);
let config = SimulationConfig {
    n_paths: 10_000,
    n_steps: 252,
    dt: 1.0 / 252.0,
    seed: 42,
};
let runner = SimulationRunner::new(gbm, config);
let result = runner.run();

// 2. Access results
let terminal_price = result.terminal(0).unwrap();  // path 0 final price
let path_prices = result.path(1).unwrap();          // all entries for path 1
}

For a full backtest with an order-book-aware strategy, the engine crate wraps the simulated data into an exchange with matching and portfolio tracking. See the engine page for a complete example.

Architecture

The central abstraction: Simulatable

Everything starts from a single trait in market_model:

#![allow(unused)]
fn main() {
pub trait Simulatable: Clone + Send + Sync {
    type Output: Clone + Send + Sync;
    fn dim(&self) -> usize;
    fn step(&mut self, dt: f64, dw: &[f64], rng: &mut impl Rng) -> Self::Output;
    fn current(&self) -> Self::Output;
}
}

dim() tells the runner how many Brownian drivers are needed. step() advances the model by one time step. current() snapshots the current state.

This trait is the only thing SimulationRunner knows about. All complexity above it -- stochastic volatility, jump components, order flow, limit order books -- composes behind implementations of this single interface.

Composition vs. inheritance

Models do not form a class hierarchy. They compose:

graph TD
    S["Simulatable trait"] --> G["GBM: price"]
    S --> H["Heston: price, variance"]
    S --> JD["JumpDiffusion: price"]
    S --> B["Bates: price, variance"]
    S --> W["Hawkes: intensity"]

Each leaf is its own impl Simulatable. The runner is monomorphized at compile time, so there is zero runtime branching on which model is active. When you run a GBM simulation, no Heston or L3 order book code is compiled or linked.

Three-layer architecture

graph TD
    subgraph L1["Layer 1: market_model (Simulation)"]
        SIM["Simulatable trait"]
        SR["SimulationRunner"]
        PS["PriceStrategy trait"]
    end

    subgraph L2["Layer 2: solver (Optimization)"]
        FD["Finite Difference Policy Iteration"]
        BSDE["BSDE LSMC"]
        AN["Analytical solutions"]
    end

    subgraph L3["Layer 3: engine (Backtesting)"]
        ENG["Engine"]
        EX["Exchange"]
        MT["Matcher"]
        ST["Strategy trait"]
        VE["VecEnv"]
    end

    SIM --> FD
    SIM --> BSDE
    SR --> VE
    PS --> ST
    FD --> ST
    BSDE --> ST

Layer 1 generates data. Layer 2 computes optimal controls. Layer 3 runs backtests and trains RL agents. Each layer depends only on the layer below, never above.

Zero-cost principle

This is enforced by the crate dependency graph: solver imports only market_model, and engine imports both. Cargo's dead code elimination strips everything not reachable from your binary.

Implementation Matrix

This page catalogues every configurable option across the three crates and shows which combinations are compatible. Use it to determine what model, solver, and engine setup fits your use case.

Models

The solver crate provides 9 pre-built optimal control models. Each implements Model<N> for a specific state dimension and set of dynamics.

Terminal conditions

Models that do not expose with_terminal_condition have hardcoded terminal values (usually zero, or a model-specific expression like -0.5*xi*q^2 for AvellanedaImpact). To apply uniform terminal liquidation at the backtest level, use run_backtest_with_liquidation from the engine.

Price behaviours captured

Solvers

Both numerical solvers accept any Model<N> implementation.

Engine strategies

All strategies implement Strategy and consume Observation (filtered by ObservationFilter) to emit OrderRequests.

Strategy to model mapping

Matchers

Two matcher implementations determine how limit orders get filled.

StochasticMatcher Hawkes modes

Data sources

Backtesting

BacktestResult metrics

Return, annualised return, volatility, Sharpe, Sortino, max drawdown, total trades, final equity, mean/max/min inventory, adverse selection (bps), realised edge (bps), inventory variance, PnL spread, PnL directional, fill buy/sell counts, mean hold time, terminal liquidation cost.

Observation filtering

VecEnv

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:

1. 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).
2. market_model: provide the price process (GBM).
3. 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

• Model agnosticism: the same PolicyIterationSolver solves 1D Avellaneda-Stoikov and 3D Hawkes OFI with no code changes. The solver only calls Model<N> trait methods.
• Terminal liquidation: the FDM value function is solved with TerminalCondition::LiquidationCost, so the optimal spreads account for the cost of flattening residual inventory at expiry. The engine applies the same half-spread cost via run_backtest_with_liquidation, keeping the strategy and the evaluation consistent.
• Separation of concerns: solver produces lookup tables (offline), market_model produces prices, engine connects them. The strategy never touches the solver; the solver never touches the exchange.
• Performance: the FDM solve runs once and the tables are cached. Each engine tick performs only cheap 4D linear interpolation.

market_model

The simulation layer. Provides stochastic processes, the SimulationRunner, model-agnostic strategies, and synthetic order book generation.

What is currently available

8 stochastic processes, all implementing Simulatable:

3 strategy types implementing PriceStrategy:

Simulation runner with:

• Pre-generated random normals (no RNG in the hot loop)
• SOA (structure-of-arrays) output layout for cache efficiency
• rayon parallel path execution
• Deterministic seeding for reproducibility
• Criterion benchmarks for GBM and Heston at multiple path counts

Order book stubs (order_book/ module): placeholder traits for OrderFlowProcess and OrderBook. L1 BBO generation is available through the engine's SimulatedDataSource. Full L3 generation with Hawkes order flow is planned.

What is not yet available

• L2 or L3 order book depth simulation (L1 BBO generation only)
• Order flow processes (Poisson, Hawkes bilateral)
• GPU-accelerated path simulation
• Multi-asset correlation matrix support in the runner

Processes

All processes in market_model implement the Simulatable trait, which provides a uniform interface for the simulation runner.

Available processes

Geometric Brownian Motion

$$dS_t = \mu S_t dt + \sigma S_t dW_t$$

The simplest model. Log-normal returns, constant volatility. Closed-form expected value and variance.

Ornstein-Uhlenbeck

$$dX_t = \theta(\mu - X_t) dt + \sigma dW_t$$

Mean-reverting process with constant diffusion. Used for interest rates, spreads, or volatility.

Cox-Ingersoll-Ross

$$dv_t = \kappa(\theta - v_t) dt + \sigma \sqrt{v_t} dW_t$$

Non-negative mean-reverting process. Used for variance and intensity modeling. Full truncation enforces the non-negativity constraint.

Heston

$$dS_t = \mu S_t dt + \sqrt{v_t} S_t dW_t^S$$ $$dv_t = \kappa(\theta - v_t) dt + \sigma \sqrt{v_t} dW_t^v$$

Two-factor model with correlated Brownian motions. Captures the volatility smile. The Cholesky decomposition of the correlation matrix is applied in each step.

Jump Diffusion

$$dS_t = \mu S_t dt + \sigma S_t dW_t + dJ_t$$ $$dJ_t = (e^Z - 1) dN_t, \quad Z \sim \mathcal{N}(\mu_J, \sigma_J^2)$$

GBM with compound Poisson jumps and log-normal jump sizes. Models sudden large moves that a pure diffusion misses.

Bates

Heston with Merton-style log-normal jumps on the price process. Reduces to Heston when lambda = 0.

Hawkes

$$\lambda(t) = \mu + \int_0^t \phi(t-s) dN_s$$

Self-exciting intensity process. Two kernel options:

• Exponential: $\phi(t) = \alpha e^{-\beta t}$ -- fast, one decay rate.
• Power law: $\phi(t) = \frac{\alpha}{(c+t)^p}$ -- heavy-tailed decay.

Simulated via Ogata's thinning algorithm. Has dim() = 0 (no Brownian driver) because all randomness comes from the acceptance-rejection step.

Rough Ornstein-Uhlenbeck

Fractional OU with Hurst parameter $H$. The kernel $t^{H-1/2}$ is approximated by a sum of exponentials (Abi Jaber 2019). The number of factors controls the accuracy-speed tradeoff.

Adding a new process

1. Create market_model/src/process/your_model.rs.
2. Implement Simulatable with your state type, dim(), step(), and current().
3. Add pub mod your_model and a re-export to process/mod.rs.
4. Add a test comparing Monte Carlo moments against analytical expectations.
5. Add a criterion benchmark against the GBM baseline.

Strategies

Strategies in market_model operate on price paths only. They produce a Signal at each step given the price history so far. They have zero knowledge of order books, portfolio state, or exchange mechanics.

PriceStrategy trait

#![allow(unused)]
fn main() {
pub trait PriceStrategy: Send + Sync {
    fn signal(&self, history: &[f64]) -> Signal;
}
}

Strategies are stateless: given the same history slice, they always produce the same Signal. Any state (warm-up windows, rolling averages) is derived from the history on each call.

Available strategies

MaCrossover

Long when the fast moving average is above the slow MA, short when below.

#![allow(unused)]
fn main() {
let s = MaCrossover::new(5, 20);
}

RsiStrategy

Long when RSI drops below the oversold threshold, short when above overbought.

#![allow(unused)]
fn main() {
let s = RsiStrategy::new(14, 30.0, 70.0);
}

BuyAndHold / AlwaysFlat

Baseline strategies: always long, and never trade.

Using strategies with the runner

PriceStrategy is separate from Simulatable. To evaluate a strategy on simulated paths, iterate the runner output and call signal at each step:

#![allow(unused)]
fn main() {
let runner = SimulationRunner::new(gbm, config);
let result = runner.run();

for path_idx in 0..result.n_paths {
    let mut history = Vec::new();
    for entry in 0..=result.n_steps {
        let price = result.get(path_idx, entry).unwrap();
        history.push(*price);
        let sig = strategy.signal(&history);
        // accumulate PnL, record signals, etc.
    }
}
}

Adapter for engine

Strategies defined here can be used inside the full exchange engine via an adapter. The adapter converts Signal::Long to a buy market order, Signal::Short to a sell, and Signal::Flat to no action. The adapter lives in the engine crate, not here.

Order Book

The order_book module in market_model provides infrastructure for generating synthetic limit order book data.

Current status

The module exists with traits and placeholder implementations, but full L1/L2/L3 generation is not yet implemented. The planned architecture is described below.

Planned architecture

Three composable layers:

graph LR
    P[PriceProcess] --> OBS[OrderBookSimulator]
    OF[OrderFlowProcess] --> OBS
    OB[OrderBook depth] --> OBS
    OBS --> Runner[SimulationRunner]

OrderBookSimulator implements Simulatable with Output = BookSnapshot, which contains the full book state at each simulation step.

PriceProcess

The true mid-price driver. Any existing Simulatable process (GBM, Heston, etc.) serves as the price process.

OrderFlowProcess

Generates limit orders, market orders, and cancellations each step. Planned implementations:

• Poisson: constant arrival rates.
• Bilateral Hawkes: self-exciting intensities for bid and ask sides.
• NoFlow: no arrivals, just BBO at mid +/- half-spread.

OrderBook depth

Controls the fidelity of the generated book:

What is currently available

• OrderFlowProcess trait and NoFlow implementation.
• OrderBook struct stub.
• The module structure is in place for future implementation.

Simulation Runner

The SimulationRunner is the high-performance batch simulation engine in market_model.

Design

#![allow(unused)]
fn main() {
pub struct SimulationRunner<S: Simulatable> {
    model: S,
    config: SimulationConfig,
}
}

The runner is generic over any Simulatable type. It is monomorphized at compile time, producing machine code specialized to your model's step() implementation.

How it works

1. Pre-generate all random normals: n_paths * n_steps * dim values from a StdRng with fixed seed. No RNG calls in the hot loop.
2. Clone model per path: each parallel path gets its own model instance, seeded RNG, and a slice of the pre-generated random data.
3. Step in parallel: rayon distributes path batches across CPU cores.
4. Store in SOA layout: outputs are interleaved by (path, entry), so sequential access strides contiguously through memory.

Output layout

Each path produces n_steps + 1 entries: the initial state plus one per step. To access:

#![allow(unused)]
fn main() {
let result = runner.run();

// Entry e of path p:
let value = result.get(p, e);  // 0 = initial, n_steps = final

// All entries for path p:
let path_data = result.path(p);

// Only the step outputs for path p (no initial):
let steps_data = result.path_steps(p);

// Final value for path p:
let terminal = result.terminal(p);
}

Configuration

#![allow(unused)]
fn main() {
pub struct SimulationConfig {
    pub n_paths: usize,   // number of independent paths
    pub n_steps: usize,   // steps per path
    pub dt: f64,          // step size in years (e.g. 1/252 for daily)
    pub seed: u64,        // reproducibility seed
}
}

Performance characteristics

• GBM: approximately 5 floating-point operations per step. 100k paths of 252 steps completes in under a second on a modern multi-core CPU.
• Heston: approximately 30 flops per step (2D state, Cholesky, exp). Roughly 6x slower than GBM.
• Jump processes: additional Poisson sampling per step adds ~50% overhead.
• Hawkes: Ogata thinning over the interval; cost depends on the kernel and event intensity.
• Rough OU: cost scales with the number of approximating factors.

Benchmark with: cargo bench -p market_model

solver

The optimization layer. Solves Hamilton-Jacobi-Bellman equations for optimal control problems, with a focus on market making.

What is currently available

Numerical solvers

Finite Difference Policy Iteration (numeric/):

• Discretizes the HJB PDE on a grid over state variables (inventory, variance, intensity).
• Uses policy iteration: alternate between solving the linear system for the value function and optimizing the control.
• Custom sparse linear algebra (linalg/): CSR matrices, SOR solver, eigenvalue decomposition.
• Configurable grid resolution, boundary conditions, and convergence tolerance.

BSDE Least-Squares Monte Carlo (numeric/bsde/):

• Backward stochastic differential equation approach.
• Forward simulation of state trajectories via SimulationRunner.
• Backward regression of value function using polynomial basis functions.
• Dimension-agnostic: works for N-dimensional models via const generics (tested up to N=3 with Heston-Hawkes).
• Basis function scaling: 6 features for N=2, 10 for N=3, 15 for N=4 (degree-2 polynomials).

Analytical solutions

Avellaneda-Stoikov (analytical/):

• Closed-form bid/ask quotes for the canonical market making model.
• Used as a baseline to validate numerical methods.

Models (models/)

Pre-built optimal control models:

Lookup tables (lookup/)

Caches pre-computed grid evaluations for fast interpolation. Used by the engine to evaluate HJB-derived strategies at sub-grid resolution in real-time during backtesting.

Linear algebra (linalg/)

Custom sparse matrix and solver implementations:

• CSR (Compressed Sparse Row) matrix format.
• SOR (Successive Over-Relaxation) iterative solver.
• Eigenvalue decomposition.
• Intel MKL-backed dense linear algebra via lapack.

What is not yet available

• GPU-accelerated linear solvers.
• ADI (Alternating Direction Implicit) splitting for higher-dimensional finite difference grids.
• Additional model types (OU-based, rough volatility, multi-asset).
• Direct integration of solved value functions into the engine as Strategy implementations (infrastructure not yet wired).

Finite Difference Solver

Policy iteration solver for HJB equations discretized on a grid.

Approach

The HJB PDE is discretized on a grid over the state variables. The solver then iterates:

1. Policy evaluation: for the current control policy (bid/ask quotes), solve the linear system for the value function at each grid point. Uses SOR for the inner solve.
2. Policy improvement: at each grid point, find the control that maximizes the HJB residual using the current value function.
3. Repeat until the policy converges.

Example

#![allow(unused)]
fn main() {
use solver::core::grid::Grid;
use solver::models::avellaneda::AvellanedaModel;
use solver::numeric::finite_difference::solver::PolicyIterationSolver;

// 1D grid over inventory q = [-10, 10], 21 points
let grid = Grid::<1>::new([21], [-10.0], [10.0]);

// Avellaneda-Stoikov model with constant intensity
let model = AvellanedaModel::new(
    0.01,   // gamma: risk aversion
    0.2,    // sigma: volatility
    100.0,  // k: order arrival intensity
    1.0,    // T: terminal time
);

let solver = PolicyIterationSolver::new(0.01); // dt = 0.01
let v = solver.solve(&grid, &model, 100); // 100 time steps
// v[i] = value function at grid point i at t=0
}

Configuration

Grid parameters, boundary conditions, convergence tolerance, and SOR relaxation parameters are configurable. See solver/src/numeric/.

When to use

Finite difference is preferred when:

• The state space is low-dimensional (1-3 dimensions).
• You need high accuracy on a fixed grid.
• The dynamics have simple boundary behavior.

It becomes infeasible above 3-4 dimensions due to the curse of dimensionality (grid points grow exponentially).

BSDE Solver

Backward Stochastic Differential Equation solver using least-squares Monte Carlo for high-dimensional HJB problems.

Approach

1. Forward simulation: simulate N state trajectories from t=0 to T using SimulationRunner.
2. Backward regression: starting from the known terminal condition at T, step backward. At each time step, regress the continuation value onto polynomial basis functions of the state variables.
3. Control optimization: at each regression point, evaluate the control (bid/ask intensities) that maximizes the Hamiltonian.

Example

#![allow(unused)]
fn main() {
use solver::numeric::bsde::BsdeSolver;
use solver::models::heston_hawkes::HestonHawkesModel;
use solver::numeric::basis::PolynomialBasis;

// Heston-Hawkes: N=3 model (inventory, variance, intensity)
let model = HestonHawkesModel::new(/* ... */);
let basis = PolynomialBasis::new(2, PolynomialType::Power);
let solver = BsdeSolver::new(10_000, 100, // paths, time steps
    solver::numeric::bsde::InitializationMode::Local(5.0));
let solution = solver.solve(&model, [0.0, 0.04, 0.5], 1.0, &basis);
// solution contains bid/ask intensities and value function
}

Dimension-agnostic design

The solver uses Rust const generics and is fully generic over the state dimension N. The same code works for:

The polynomial basis function count grows as O(N^2) for degree-2 polynomials:

Configuration

• Number of paths and time steps.
• Basis function type (polynomial, scaled).
• Regularization parameter for the regression.
• Convergence tolerance.

When to use

BSDE is preferred when:

• The state dimension is 3 or higher.
• Grid-based methods would be intractable.
• You can tolerate Monte Carlo noise in the solution.
• You need to scale to models with more factors in the future.

Analytical Solutions

Closed-form and semi-analytical solutions for benchmark validation.

Avellaneda-Stoikov

The canonical market making model with constant order arrival intensity. The optimal bid and ask quotes have closed-form expressions as functions of inventory, time remaining, and model parameters.

The reservation price (indifference price) is:

$$r(s, q, t) = s - q \gamma \sigma^2 (T - t)$$

The optimal half-spread is:

$$\delta(t) = \frac{\gamma \sigma^2 (T-t)}{2} + \frac{1}{\gamma} \ln\left(1 + \frac{\gamma}{k}\right)$$

Use as baseline

The analytical solution serves two purposes in the framework:

1. Validation: compare finite-difference and BSDE numerical solutions against the known analytical result to verify correctness.
2. Strategy in the engine: the solution is wrapped as a Strategy in the engine crate and can be backtested against simulated markets. Implemented as AvellanedaStoikovStrategy with optional online volatility estimation.

Limitations

The analytical solution assumes:

• Constant volatility.
• Constant order arrival intensity.
• Linear inventory penalty.

For models with stochastic volatility, self-exciting intensities, or non-linear penalties, use the numerical solvers (finite difference or BSDE) instead.

N-Dimensional Models

The solver crate supports models with arbitrary state dimension through Rust const generics.

Model<N> trait

#![allow(unused)]
fn main() {
pub trait Model<const N: usize> {
    fn optimize(&self, ...) -> ControlOutput<N>;
    fn terminal(&self, ...) -> f64;
    fn next_step(&self, ...) -> Array1<f64>;
    fn is_diffusion_dimension(&self, dim: usize) -> bool;
    fn gradient_step(&self, dim: usize) -> f64;
    fn transform_noise(&self, ...) -> Array1<f64>;
}
}

Conventions

• Dimension 0 is always inventory for market-making models.
• solve_with_spreads and NumericalResult extract lambda_plus[0] and lambda_minus[0] -- convenience methods for the trading dimension only.
• For full N-dimensional control access, use BsdeSolution<N>::evaluate_control().

Adding a new model

1. Create the model struct implementing Model<N>.
2. Register in models/mod.rs.
3. For N above 3, increase num_paths to 30k or more and consider regularization of 1e-4 for stability.

Proven: N=3 Heston-Hawkes

The HestonHawkes model implements Model<3> with state [q, v, lambda]:

• dim 0 (inventory): discrete jumps, controlled via bid/ask intensities.
• dim 1 (variance): CIR diffusion, is_diffusion_dimension(1) = true.
• dim 2 (intensity): deterministic mean-reverting drift (Hawkes).

Run the example: cargo run --release --example heston_hawkes_n3

engine

The backtesting layer. Runs strategies inside a simulated exchange with configurable matching, observation filtering, and multi-agent support.

What is currently available

Core engine

The Engine orchestrates the simulation loop:

1. Advance time (data source produces new MarketState).
2. Build observations per agent (via ObservationFilter).
3. Strategies emit OrderRequests.
4. Exchange resolves matches and updates portfolios.
5. Data source receives aggregate impact.

Supports single-agent and multi-agent configurations with per-agent cash, transaction costs, and observation filters.

Data sources

Simulated (SimulatedDataSource):

• Wraps any Simulatable process from market_model.
• Generates BBO quotes from the process state.
• Supports pre-calculation for zero-impact scenarios (fast path).
• On-the-fly calculation with market impact feedback.

Historical (ParquetDataSource):

• Replays recorded market data from Parquet files.
• Used to backtest strategies on real exchange data.

Exchange (Exchange)

• L1 (BBO-only) exchange with infinite liquidity at top-of-book.
• Per-agent portfolio isolation via HashMap<AgentId, Portfolio>.
• Per-agent CancelAll (cancels only that agent's resting orders).
• Per-agent fill history and last-step fill snapshots.
• Transaction cost deduction per fill.

Matchers

Simple (SimpleMatcher):

• Deterministic immediate fill at BBO.

Stochastic (StochasticMatcher):

• Three-pronged fill evaluation per the Avellaneda-Stoikov model:
  1. Aggressive crossing (limit price crosses BBO).
  2. Deterministic mid-price sweep (mid passed through the limit).
  3. Probabilistic Poisson arrival at distance from mid.

Strategies

All strategy implementations for the engine:

Backtesting

run_backtest produces per-path metrics:

• Total return, annualized return, volatility.
• Sharpe and Sortino ratios.
• Max drawdown, final equity.
• Mean, max, and min inventory.
• Adverse selection (bps) and realized edge (bps).
• PnL decomposition (spread component, directional component).
• Fill counts, mean hold time.
• Terminal liquidation cost.

BacktestSummary aggregates statistics across multiple Monte Carlo runs.

Vectorized environment (VecEnv)

For RL training:

• Runs hundreds of independent Heston + AS engine instances in parallel.
• Configurable reward: PnL or differential Sharpe ratio.
• Returns states, rewards, dones, equities, and positions.
• Fixed configuration via VecEnvConfig with defaults.

Python bindings (optional)

When compiled with --features python, the engine exposes strategies and backtesting via PyO3. See engine/pyproject.toml.

What is not yet available

• L2 or L3 exchange (aggregated price levels, FIFO order-by-order matching).
• Queue priority mechanics.
• Latency profiles per agent (all agents see state simultaneously).
• Transaction fees / maker-taker fee structure.
• Multi-asset / cross-asset data sources.
• Hawkes-process data source in engine (only GBM, Heston, Bates).
• Exogenous signal feeds for agents.
• Endogenous market impact from aggregate fills.
• HJB-solver-driven strategies (optimal policies not yet wired as Strategy).

Engine Strategies

The engine crate defines its own Strategy trait for order-book-aware trading:

#![allow(unused)]
fn main() {
pub trait Strategy {
    fn on_tick(&mut self, obs: &Observation, requests: &mut Vec<OrderRequest>);
}
}

This is a different trait from market_model::PriceStrategy. It operates on Observation (BBO, portfolio, volatility, drift, parameters) and emits OrderRequests (place, cancel, cancel-all) into the exchange.

Available strategies

Avellaneda-Stoikov family

The core market-making strategy with several variants:

• AvellanedaStoikovStrategy: the canonical AS quotes. Optional online volatility estimation (EWMA rolling window) and real-vol passthrough. Optional permanent price impact adjustment.

• AvellanedaStoikovHestonStrategy: AS quotes with Heston stochastic volatility dynamics. The strategy observes the current variance from the Observation filter.

• AvellanedaStoikovHawkesStrategy: AS quotes with Hawkes self-exciting intensity dynamics.

• AvellanedaStoikovBilateralHawkesStrategy: separate Hawkes intensities for bid and ask sides.

• AvellanedaStoikovBilateralHawkesOrderFlowImbalanceStrategy: extends the bilateral Hawkes variant with an order flow imbalance signal for additional adverse selection protection.

• AvellanedaStoikovExactStrategy: uses the exact (non-approximate) AS solution for the optimal spread.

Baseline strategies

• ConstantSymmetricStrategy: posts a fixed spread around the mid price. Good for isolating the effect of spread width.

• ZeroIntelligenceStrategy and RandomStrategy: random order placement. Used as noise-trader agents in multi-agent simulations.

• ExternalStrategy: loads strategy parameters from a configuration file. Allows swapping strategies without recompiling.

Observation filtering

Each agent sees a filtered view of the market via ObservationFilter:

• Transparent: all ground truth visible (volatility, drift, parameters).
• Opaque: only BBO and portfolio visible.
• Partial: whitelist of visible parameters.

This lets you test strategies under realistic information constraints: a market maker that cannot observe true volatility must estimate it from price data.

Matching

The matcher determines how limit orders get filled. Two implementations are available.

Simple matcher

SimpleMatcher: deterministic immediate fill at BBO. If the limit price crosses the BBO (e.g. a buy priced above the ask), the order fills immediately at the BBO price.

Fast and predictable. Use for:

• Low-frequency strategy evaluation.
• PnL attribution.
• When fill probability is not the focus of the test.

Stochastic matcher

StochasticMatcher: three-pronged fill model per the Avellaneda-Stoikov framework. On each tick, every resting order is evaluated:

1. Aggressive crossing

The order's limit price is on the wrong side of the current BBO. Condition (buy): limit >= best_ask. Condition (sell): limit <= best_bid. Fill price: BBO price. Immediate fill with no randomness.

2. Deterministic mid-price sweep

The limit was placed passively, but the mid-price path from $S_t$ to $S_{t+1}$ moved through the limit level. Because the price crossed the order, every market-order arriving during the interval would have filled it. Fill probability: 1.0 (guaranteed). Fill price: limit price.

3. Probabilistic Poisson arrival

The mid has not swept through the limit. Random incoming market orders may still hit the resting order during $[t, t+1]$. The Poisson intensity is computed using the distance from the post-move mid $S_{t+1}$:

$$\lambda = a \cdot e^{-k \cdot \delta}, \quad \delta = |limit - mid_{t+1}|$$

$$p_{fill} = 1 - e^{-\lambda \cdot dt}$$

Fill price: limit price.

The Poisson arrival rate decays exponentially with distance from mid, modeling the fact that a resting order far from the mid is much less likely to be hit by a random market order. The parameters a and k are configurable.

Exchange

The Exchange manages order books, matching, and portfolio state for all agents.

Current implementation: L1

The exchange currently implements level-1 (BBO only) mechanics:

• Tracks only best_bid and best_ask.
• Infinite liquidity at BBO (no depth depletion).
• No queue position or FIFO priority.

This is the fastest model and sufficient for:

• Testing low-frequency strategies.
• PnL attribution.
• Spread capture evaluation.

Per-agent portfolio isolation

Each agent has a separate Portfolio:

#![allow(unused)]
fn main() {
pub struct Portfolio {
    pub cash: f64,
    pub position: f64,
    pub realized_pnl: f64,
}
}

The exchange enforces isolation: agent A's orders cannot see or affect agent B's portfolio. CancelAll cancels only the requesting agent's resting orders. Fill history is tracked per agent.

Planned: L2 and L3

Future implementations will add depth-aware matching:

• L2: aggregated volume per price tick. Market orders walk the book, consuming volume level by level. No queue priority.
• L3: FIFO queues per tick. Strict price-time priority. Per-order identity and agent attribution. Higher computational cost but necessary for testing HFT strategies where queue position matters.

Transaction costs

Transaction costs are deducted per fill as a fraction of the fill value. The cost rate is configurable per agent.

Backtesting

The engine provides a Monte Carlo backtest runner that evaluates a strategy across many independent market paths.

Single-path metrics

run_backtest drives Engine::step for a single path and computes:

Multi-path summary

BacktestSummary aggregates these metrics across multiple Monte Carlo runs, reporting mean and standard deviation for each metric.

Usage

#![allow(unused)]
fn main() {
use engine::backtest::run_backtest;
use engine::data_source::SimulatedDataSource;
use engine::matcher::SimpleMatcher;
use engine::strategies::AvellanedaStoikovStrategy;

let gbm = GeometricBrownianMotion::new(0.0, 0.2, 100.0);
let source = SimulatedDataSource::new(gbm, 100.0, dt, spread, n_steps, 0.0);
let strategy = AvellanedaStoikovStrategy::new(params);

let result = run_backtest(
    &SimpleMatcher::new(),
    &strategy,
    source,
    initial_cash,
    transaction_cost_rate,
);
println!("Sharpe: {:.2}, Sortino: {:.2}", result.sharpe_ratio, result.sortino_ratio);
}

VecEnv for RL training

VecEnv runs multiple engine instances in parallel for reinforcement learning. Each slot is an independent simulation episode. The environment returns states (mid-price, position, time, spread, volatility, intensities), rewards, dones, and metadata.

Reward functions:

• PnL: change in total equity (cash + position * mid).
• DiffSharpe: differential Sharpe ratio, encouraging smooth risk-adjusted returns rather than absolute PnL.

Vectorized Environment

VecEnv runs many independent simulation episodes in parallel for reinforcement learning.

Configuration

VecEnvConfig specifies the full environment:

State space

Each environment step produces a state vector of typically 6-8 elements:

• Mid-price.
• Current position (inventory).
• Time remaining.
• Spread.
• Volatility.
• Intensity (if Hawkes enabled).

Usage

#![allow(unused)]
fn main() {
use engine::vec_env::{VecEnv, VecEnvConfig, RewardType};

let config = VecEnvConfig {
    reward_type: RewardType::DiffSharpe,
    max_steps: 1000,
    ..Default::default()
};

let mut env = VecEnv::new(config, 64); // 64 parallel environments

// Reset to start new episodes.
let states = env.reset();

// Step with actions.
let actions = vec![0.5; 64]; // symmetric spread for all envs
let (next_states, rewards, dones) = env.step(&actions);
}