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 serve --open

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 ControlProblem<N> for a specific state dimension and set of dynamics; models solved on a diffusive grid also implement the FD-specific PdeProblem<N> transport contract. All models are solved rigorously — they discretize the exact HJB PDE without approximation. See Rigor Reference for the full mapping of HJB to method.

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 ControlProblem<N> implementation. The grid-based PolicyIterationSolver additionally requires PdeProblem<N> for the implicit, Crank-Nicolson, and Strang-ADI schemes; explicit Euler only needs ControlProblem.

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_grid_pde_history returns the terminal condition at index 0 and the
// value at tau = 0, tau = dt, ..., tau = T as the remaining snapshots.
let v_history = solver.solve_grid_pde_history(&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 ControlProblem and PdeProblem 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.

Validation

Validation of solver results is recorded in three pages:

• Classical methods - the finite-difference (FD) and least-squares Monte Carlo (BSDE) paths in the Rust solver crate, validated inside cargo test.
• Neural network methods - the deep BSDE, deep HJB (DGM), and neural operator methods in the neural_solver package (JAX), validated by the Python test suite and against committed Rust reference values.
• BSDE precision - the one-time high-precision BSDE verification against exact solutions, targeting \(\leq 0.1\%\) relative error.

Canonical error metric

All numerical-versus-reference comparisons use the canonical error metric in solver::validation. The primary quantity is the normalized relative error

\[ \varepsilon = \frac{|x - x^{\text{ref}}|}{\max(|x^{\text{ref}}|, f)} \]

where f is the small absolute floor (DEFAULT_ABS_FLOOR = 1e-12) that prevents division by zero when the reference is at or near zero. A percentage is \(100\,\varepsilon\). The helpers relative_error, percent_error, assert_within, and assert_close are shared across every solver test and benchmark so that a reported tolerance means the same thing regardless of the magnitude of the quantity under test.

Classical Methods

This page records how the classical solver results are validated: the finite-difference (FD) and least-squares Monte Carlo (BSDE) paths in the Rust solver crate. For every benchmark problem it lists the reference used, the numerical method checked against it, and the tolerance at which the check currently passes. It complements the narrative solver pages (PDE Solving Methods and Elliptic PDE Solving Methods) by compressing the same information into machine-checkable tables. The neural methods are validated separately on the neural network methods page.

Validation policy

The solver has two kinds of validation problems:

• Exact-solution problems have a closed-form or semi-analytical reference. They are validated numerical-versus-exact. Both the finite-difference (FD) and least-squares Monte Carlo (BSDE) paths are checked against the same reference where the method is applicable.
• No-exact-solution problems (Hawkes, Heston, and higher-dimensional market-making models) have no reduced closed form. They are validated numerical-versus-numerical: FD against BSDE, plus degeneracy tests that reduce the model to a model with an exact solution.

Tolerances are recorded as obtained, not as goals. They reflect the assertions currently present in the tests. Every entry states whether the check is deterministic (FD, closed-form consistency) or stochastic (BSDE).

BSDE precision policy

BSDE is Monte Carlo and does not converge to high precision cheaply. Its standard-suite assertions therefore bound coarse, smoke-level error, not a tight tolerance. High-precision BSDE runs (targeting relative error at or below \(0.1\%\)) are a one-time exercise, not part of cargo test. They are run on demand and their results are recorded in the BSDE table below; they are not re-run on every build.

The standard BSDE checks divide into:

• Tier 1 (default suite): fast smoke tests and zero-time-step consistency. These assert well-formedness, symmetry, and agreement at the terminal boundary, not tight finite-horizon accuracy.
• Tier 2 (slow-tests): finite-horizon BSDE checks with bounded, generous tolerances driven by Monte Carlo noise.
• Tier 3 (one-time exercise): high-precision BSDE runs targeting \(\leq 0.1\%\) relative error. Run explicitly, never in the default suite.

Reading the tolerance column

• A percentage is a relative error bound.
• An absolute figure is a spread or value error bound.
• "Machine" means the check is at or near floating-point round-off.
• "Finite/well-formed" means the assertion is qualitative (finite, positive, symmetric), not a numerical tolerance.

Canonical error metric

All numerical-versus-reference comparisons use the canonical error metric defined on the Validation page.

Exact-solution benchmarks: FD path

Exact-solution benchmarks: BSDE path

BSDE is Monte Carlo. The standard assertions are smoke-level; the high-precision target is a one-time exercise whose obtained results are documented in BSDE Precision.

The BSDE high-precision target is not asserted in the default suite. The one-time high-precision verification and its obtained precision are recorded in BSDE Precision.

No-exact-solution benchmarks: numerical vs numerical

These models have no reduced closed form. FD and BSDE are checked against each other, and degeneracy limits reduce each model to an exact-solution reference. The expanded models (Heston, AvellanedaHawkes, HestonHawkes, BilateralHawkes, BilateralHawkesOrderFlowImbalance) are full-value problems: their bsde_driver is running_reward + local_source (factor transport excluded, reward rate bounded at 1/dt per fill side) and their next_step_controlled simulates inventory fills under the optimal control. Because their inventory dynamics depend on the control, finite-horizon BSDE checks run in coupled (Picard) mode; the decoupled mode is correct only at the terminal boundary.

Consistency and convergence checks

Beyond exact-vs-numerical and numerical-vs-numerical comparisons, the following checks pin down the analytical layer and solver convergence.

Test suite tiers

Neural Network Methods

This page records how the neural methods are validated. The neural_solver package (JAX) implements mesh-free deep BSDE, deep HJB (DGM), and neural operator methods. Unlike the classical solver, its validation is not part of cargo test: it is exercised by the Python test suite in neural_solver/tests/ and its measured results are documented in the neural empirical results. This page is the single-place summary; the per-architecture numbers are stated once in the empirical pages and are not duplicated here beyond the headline figures.

The neural methods are validated in two ways:

• Against exact references - the same closed forms the classical solver uses (Black-Scholes, Merton, LQ, jump-LQ), and for the jump path a Riccati closed form that exercises the learned jump integrand directly.
• Against the classical solver's committed values - reference agreement with no training: the JAX closed forms reproduce the Rust solver values to float32 precision.

Cross-architecture cost and error (deep BSDE vs DGM vs operator, including the no-exact-solution targets) are in comparison; which architecture fits which exact-solution problem is in coverage. The reference schema the agreement check consumes is documented in usage.

BSDE Precision

This page documents the one-time high-precision verification of the least-squares Monte Carlo (BSDE) solver against exact solutions. It is a companion to the Classical Methods page, which records the tolerances asserted by the automated test suite. The standard suite bounds coarse, smoke-level BSDE error; the target here is the tighter \(\leq 0.1\%\) relative error that the BSDE solver should reach when it is correct.

Scope

The verification covers every exact-solution benchmark the BSDE path can solve:

• Merton log-utility portfolio (one dimension).
• Finite-horizon linear-quadratic regulator (one dimension).
• Avellaneda-Stoikov market making and its drift and impact variants.

For each problem the exact closed form is stated once in Exact Solutions in Stochastic Optimal Control. This page records only the concrete verification parameter set, the resulting reference value, and the measured error.

Metric

All errors are reported with the canonical metric defined on the Validation page.

Merton portfolio

Reference closed form in Exact Solutions in Stochastic Optimal Control. With \(r = 0\), \(\mu = 0.1\), \(\sigma = 0.3\), \(T = 0.5\), \(x_0 = 10\) the exact value is \(V(0, 10) = 2.330363\).

A fixed-seed replicate study (solver/examples/bsde_noise_vs_bias) at dt = 0.005, 100_000 paths, five replicates, Hermite(4), reg = 1e-3 gives a mean relative error of 0.04% with a standard error of 0.01%. The mean is within three standard errors of zero, so the residual is Monte Carlo noise: there is no detectable deterministic bias.

Finite-horizon linear-quadratic regulator

Reference Riccati closed form in Exact Solutions in Stochastic Optimal Control. With \(A = -0.5\), \(B = 1\), \(C = 0.1\), \(Q = 1\), \(Q_T = 1\), \(R = 1\), \(T = 0.5\), \(x_0 = 0.5\), the Riccati value used by LqRegulator::exact_value is approximately -0.18666439.

A fixed-seed replicate study (solver/examples/bsde_noise_vs_bias) at dt = 0.005, 100_000 paths, five replicates, Hermite(4), reg = 1e-3 gives a mean relative error of 0.47% with a standard error of 0.05%. The mean is roughly nine standard errors from zero, so this is a deterministic bias, not noise. It is Euler-Maruyama discretization error: halving the time step to dt = 0.0025 lowers the mean to 0.30%. The \(\leq 0.1\%\) target is therefore time-step limited; reaching it for the LQ value requires dt below about 0.001 at this path count.

Avellaneda-Stoikov market making

Reference matrix-exponential closed form in Exact Solutions in Stochastic Optimal Control. With \(\gamma = 0.5\), \(k = 1.5\) the base total spread is \(2\delta_0 = 1.15073\), and the exact total spread at \(q = 1\) is approximately 1.18164 (the inventory skew contributes 0.03091).

A single-replicate run at dt = 0.005, 20_000 paths gives a total spread within 0.02% of the exact 1.18164; the inventory skew is recovered.

Drift and impact variants

Reference spectral closed forms in Exact Solutions in Stochastic Optimal Control. The degeneracy limits are \(\mu = 0\) (drift) and \(\xi = 0\) (impact), both of which reduce to the base Avellaneda-Stoikov exact spread. The reduced forward pass holds the inventory (jump) dimension fixed; the reason the reduced value must not simulate inventory jumps is in BSDE Solver.

The drift total spread recovers the exact skew to within Monte Carlo noise (0.02% mean, 0.008% SEM at 100_000 paths). The impact total spread matches the exact solution to within Monte Carlo noise (0.0000% at the same settings).

Conclusion

A fixed-seed replicate study (solver/examples/bsde_noise_vs_bias, dt = 0.005, 100_000 paths, five replicates) gives:

The <= 0.1% target is met for Merton, Avellaneda-Stoikov, drift, and impact. The LQ value is time-step limited (needs dt <= ~0.001), tracked as I18.

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. The mathematical definitions of these processes — their driving equations, parameters, and moments — are stated once in Stochastic Processes and are not restated here.

Available processes

Each process implements dim(), step(), and current(); the Hawkes process has dim() = 0 because its randomness comes entirely from the thinning (Ogata) acceptance-rejection step rather than a Brownian driver.

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 the analytical expectations in Stochastic Processes.
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: it solves Hamilton-Jacobi-Bellman (HJB) equations for optimal-control problems, with a focus on market making.

This page indexes the implementations against their canonical mathematics. Per the ownership rule, exact solutions and mathematical background live in the mathematical reference, not here: solver pages document implementation (contracts, configuration, usage) and link back to the reference for the mathematics they realize. The single exception is the neural theory section, which carries its own method theory.

Implementation-to-solution index

Each exact solution is stated once in the reference; the table maps the code that realizes it to that statement. The stochastic-optimal-control closed forms are on the SOC page; the pure PDE benchmarks (parabolic and elliptic) are on the PDE page. The SOC models without a closed form (Heston, Hawkes, and their combinations) are stated on SOC Models Without Exact Solutions.

Analytical solutions

Optimal-control models (ControlProblem<N>)

Numerical solvers

• Finite Difference — policy iteration on a grid, validated against the exact solutions above under grid refinement.
• PDE Solving Methods — the transport/source contract and dimension kinds.
• Elliptic PDE Solving Methods — the stationary path; its benchmark exact solutions are on the PDE page.
• BSDE Solver — least-squares Monte Carlo.
• N-Dimensional Models — the ControlProblem<N> contract.
• Economic Invariants — qualitative sign/ordering/symmetry checks.

Analytical solutions

The closed-form Avellaneda-Stoikov quotes are the reservation price and half-spread

\[r(s,q,t) = s - q\gamma\sigma^2(T-t), \qquad \delta(t) = \frac{\gamma\sigma^2(T-t)}{2} + \frac{1}{\gamma}\ln\Big(1 + \frac{\gamma}{k}\Big),\]

the small-\(\tau\) (near-maturity) limit of the exact solution in the SOC page. They serve two purposes: validating the numerical solvers, and driving the engine AvellanedaStoikovStrategy. The exact (matrix-exponential) form is AvellanedaExact.

Neural methods

The neural_solver package implements GPU-native BSDE and HJB methods in JAX and is validated against the analytical and numerical references above. See Neural Methods.

Finite Difference Solver

Policy iteration solver for HJB equations discretized on a grid.

The mathematical contract for the grid-based methods is documented in PDE Solving Methods. This page covers the solver's role and configuration; the FD-specific transport contract and dimension kinds live on that page.

Approach

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

1. Policy improvement: at each grid node, find the control that maximizes the HJB residual using the current value function.
2. Policy evaluation: solve the resulting linear system for the value function. Implicit, Crank-Nicolson, and Strang-ADI paths use the FD transport operator; explicit Euler uses the scalar driver.
3. Repeat backward in time from the terminal condition to t = 0.

Solvers

• PolicyIterationSolver::solve_grid_control solves a generic [ControlProblem] with an explicit Euler step.
• The FD-specific PdeProblem contract supplies the transport operator for the implicit, Crank-Nicolson, and Strang-ADI integrators.

See solver/src/numeric/finite_difference/ for the implementation.

Configuration

Grid bounds and resolution are configured through Grid<N>. Time stepping, SOR tolerance, and scheme selection are configured on PolicyIterationSolver. 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).

PDE Solving Methods

Ground truth for the grid-based HJB solvers. This page fixes the mathematical contract that the finite-difference path implements.

Scope and relationship to ControlProblem

Two contracts are kept separate:

• ControlProblem<N> is the generic stochastic-control contract. It supplies the running reward, the scalar infinitesimal generator, the terminal condition, and the control optimizer. It is the contract consumed by the BSDE regression solver and by an explicit-Euler grid step.

• The FD-specific contract (PdeProblem<N>) is consumed only by grid PDE solvers. It supplies the transport operator, not a scalar generator. This separation exists because a scalar driver is sufficient for an explicit step and for the BSDE backward pass, but it is not sufficient to assemble the implicit operator \((I - dt\, L) V_{new} = V_{old} + dt\, s\).

The generic contract is the one a user implements when they only need BSDE or an explicit FD smoke test. A model that must be solved accurately on a diffusive grid also implements the FD-specific contract.

HJB equation

For state \(x \in \mathbb{R}^n\), control \(u \in U\), running reward \(f\), and terminal reward \(g\), the value function satisfies

\[ 0 = \partial_t V + \sup_{u \in U} \lbrace f(t,x,u) + \mathcal{L}^u V(t,x) \rbrace, \qquad V(T,x) = g(x). \]

The generator splits into a transport part and a purely local source:

\[ \mathcal{L}^u V = \mathcal{T}^u V + s^u, \]

where T collects the terms that couple a grid node to its neighbours and s collects the terms evaluated at the node itself.

Dimension kinds

A reliable FD scheme must know how to discretize each coordinate. The framework distinguishes three kinds:

The kind is declarative, not inferred from whether a model happens to use fwd or bwd. The generic trait exposes only is_diffusion_dimension and gradient_step; the FD contract adds a dimension_kind(dim) method so the solver selects the stencil without model-specific branching.

Finite-difference stencil

Let \(h_i\) be the grid spacing in dimension \(i\) and \(V_{i,+}\) / \(V_{i,-}\) the value at the forward and backward neighbours. The transport operator is the sum of one-dimensional contributions

\[ \mathcal{T} V = \sum_i \big[ a_i^+ (V_{i,+} - V) + a_i^- (V_{i,-} - V) \big]. \]

Diffusion dimension

Write the diffusion coefficient as \(D_i = 0.5 \sigma_i^2\) and the drift as \(b_i\). Central second differences and upwind first differences give

\[ a_i^+ = \frac{D_i}{h_i^2} + \frac{\max(b_i, 0)}{h_i}, \qquad a_i^- = \frac{D_i}{h_i^2} + \frac{\max(-b_i, 0)}{h_i}. \]

The central term recovers

\[ \frac{D_i}{h_i^2}(V_{i,+} - 2V + V_{i,-}) = D_i\, \partial_{x_i x_i} V + O(h_i^2). \]

Deterministic drift dimension

With no diffusion, \(D_i = 0\) and only the upwind drift terms remain:

\[ a_i^+ = \frac{\max(b_i, 0)}{h_i}, \qquad a_i^- = \frac{\max(-b_i, 0)}{h_i}. \]

This is the stable upwind approximation of \(b_i \partial_{x_i} V\).

Correlated diffusion (cross-derivative terms)

When two coordinates are driven by correlated Brownian shocks, the generator contains mixed second derivatives weighted by the off-diagonal covariance. For a diffusion matrix \(\sigma(x)\) with instantaneous covariance \(D = \sigma\sigma^{\top}\), the diffusion term is

\[ \tfrac{1}{2}\,\mathrm{tr}\!\big(D\,\mathrm{Hess}\, V\big) = \tfrac{1}{2}\sum_{i,j} D_{ij}\,\partial_{x_i x_j} V. \]

The diagonal \(i = j\) recovers the single-coordinate diffusion terms already handled by the per-dimension stencils. The off-diagonal \(i \neq j\) entries are the correlated-diffusion cross terms this section adds.

The mixed second derivative is approximated on the grid by the centered four-corner stencil

\[ \partial_{x_i x_j} V(x) \approx \frac{V_{++} - V_{+-} - V_{-+} + V_{--}}{4 h_i h_j}, \]

where the four evaluations offset the state by \(\pm h_i e_i\) and \(\pm h_j e_j\). This is second-order accurate in each spacing and is only defined at interior nodes where both coordinates have both neighbours.

The cross term is consumed through the full symmetric Hessian carried by StateDerivatives::hessian_full. A model whose generator sums \(\tfrac12 \sum_{i,j} D_{ij}\, \partial_{x_i x_j} V\) reads that field directly; the explicit-Euler ControlProblem path therefore handles correlated noise without any ad hoc drift correction. See solver::models::lq_regulator::LqRegulator for the validating model with a full covariance matrix.

Discrete jump dimension

For an integer state with unit jumps, the forward and backward transition intensities are used directly:

\[ a_i^+ = \lambda_i^{+}, \qquad a_i^- = \lambda_i^{-}, \]

with \(h_i = 1\). The contribution \(\lambda_i^+(V(q+1)-V(q)) + \lambda_i^-(V(q-1)-V(q))\) is the exact infinitesimal generator of the jump process, not a Taylor approximation.

General jump dimension

A Jump dimension relaxes the unit-amplitude assumption. The jump kernel of dimension \(i\) is a finite set of transitions \(\lbrace(\Delta_j, \lambda_j)\rbrace\), where \(\Delta_j\) is an integer amplitude and \(\lambda_j\) the arrival intensity at which that transition fires. The transport contribution is the sum over the kernel

\[ \mathcal{T}_i V = \sum_j \lambda_j \, \big( V(x + \Delta_j e_i) - V(x) \big), \]

with the amplitude mapped directly to a grid-node offset. The state is interpreted on an integer lattice, so one amplitude unit equals one grid node: \(V(x + \Delta_j e_i)\) is the value at the node \(\Delta_j\) positions along coordinate \(i\). This is the exact infinitesimal generator of the competing-Poisson jump process, not a finite-difference approximation.

The DiscreteJump kind is the specialization whose kernel is exactly \(\lbrace(+1, \lambda_i^+), (-1, \lambda_i^-)\rbrace\). The FD contract exposes the kernel through jump_kernel; its default folds Transport::plus/minus into that unit pair, so a DiscreteJump model does not override it. A Jump model overrides jump_kernel to return its full kernel, and the explicit-Euler path sums the kernel rather than using the single forward/backward pair.

The general kernel is consumed by the explicit-Euler path. The implicit, Crank-Nicolson, and Strang-ADI integrators still assemble a tridiagonal operator from plus/minus and therefore remain limited to unit jumps; a general kernel is tridiagonal-compatible only when every transition is a \(\pm 1\) neighbour.

Transport contract

The FD-specific contract decomposes the controlled generator into a transport term and a local residual source:

\[ (T^u V, s^u). \]

T is assembled from the per-dimension plus and minus coefficients. s is the scalar right-hand-side term chosen so that, at the current value, the transport plus the source reproduces the full optimized driver:

\[ T^{u} V + s^{u} = \sup_u \lbrace f + \mathcal{L}^u V \rbrace. \]

This is the standard linearized operator form. The transport goes into the implicit operator and the source goes into the right-hand side. The source is a scalar at each node but may depend on the current derivatives when a model chooses to move part of a nonlinear or already-optimized term into the implicit operator for stability.

This contract is deliberately not the legacy market-making triple \((\lambda_+, \lambda_-, flow)\). The coefficients have explicit meaning: transport coefficients define T; the source is the residual that makes the scheme consistent. The legacy flow mixed a Hamiltonian value with a cancellation correction and hid the operator semantics.

Time integrators

At each backward step the linear problem is

\[ (I - \theta\, \Delta t\, \mathcal{T}^{u^*}) V^n = (I + (1-\theta)\Delta t\, \mathcal{T}^{u^*}) V^{n+1} + \Delta t\, s^{u^*}, \]

with theta = 0 explicit Euler, theta = 1 implicit Euler, and theta = 1/2 Crank-Nicolson. The discount term enters the diagonal as an additional \(+ dt \, r\).

• Explicit Euler is simple but CFL-limited: diffusive models require dt small enough to keep the operator non-negative.
• Implicit Euler is unconditionally stable and is the fallback for stiff diffusive models.
• Crank-Nicolson is second-order in time but can oscillate on non-smooth terminal data. Rannacher smoothing applies two implicit steps at startup.
• Strang ADI is the two-dimensional splitting used when one dimension is stiff (for example Heston variance) and the other is transport-like (for example inventory). It runs three tridiagonal Thomas sweeps per step.

The FD-specific contract is required for the implicit, Crank-Nicolson, and ADI paths. The scalar ControlProblem::driver is sufficient only for explicit Euler. Hyperbolic equations such as (u_{tt} = c^2 u_{xx}) require a second-order-in-time integrator and are out of scope.

Boundary convention

The time-dependent transport solver treats a missing neighbour as a zero-flux boundary: it drops the forward/backward coefficient at the edge. Exact solutions on unbounded or non-reflecting domains are therefore compared only in the interior. First-class Dirichlet, Neumann, and Robin conditions are available on the stationary path; see Elliptic PDE Solving Methods.

Market-making Hamiltonian transport cancellation

In the CARA market-making HJB the jump supremum is already maximized by the optimal spread. The optimized Hamiltonian

\[ H^* = \frac{\lambda^{a*} + \lambda^{b*}}{\gamma + \kappa} \]

contains the jump contribution

\[ \lambda^{a*}(V(q-1)-V(q)) + \lambda^{b*}(V(q+1)-V(q)). \]

The implemented transport/source decomposition keeps the optimal intensities as inventory transport and subtracts the same jump contribution from the local source:

• inventory dimension: \(plus = \lambda^{b*}\), \(minus = \lambda^{a*}\),
• local source: \(H^* + risk\_penalty - jump\_transport\),

where

\[ jump\_transport = \lambda^{b*} (V(q+1)-V(q)) + \lambda^{a*} (V(q-1)-V(q)). \]

This is algebraically equivalent to folding the jump term entirely into the source and setting the inventory transport to zero. The implemented form is chosen because it keeps the jump coupling in the implicit operator, where it improves stability. The subtraction is mandatory: without it the jump contribution is double-counted.

The legacy Model<N> code emulated the folded form with a drift_correction that subtracted the inventory transport out of flow. That cancellation is unnecessary under the clean transport/source split and must not be carried into the new contract.

Configuring discrete dimensions

Grid dimensions are typed by kind through the FD-specific contract. The generic trait still exposes is_diffusion_dimension and gradient_step for mesh-free stencils; PdeProblem adds a declarative dimension_kind(dim) method so that:

• DiscreteJump dimensions use transition rates and unit grid spacing,
• Jump dimensions use an arbitrary jump kernel, whose amplitudes map directly to grid-node offsets and are summed by the explicit scheme,
• Diffusion dimensions use diffusion plus upwind drift,
• DeterministicDrift dimensions use upwind drift only.

This removes the need for model-specific fwd/bwd bookkeeping inside the solver and is what enables the stable implicit and ADI integrators for Heston, Hawkes, and American-put models.

Elliptic PDE Solving Methods

Ground truth for the stationary finite-difference path. Elliptic problems have no time variable, so the boundary conditions are the problem: the solution is the direct solve of

\[ -\mathcal{T} u + \rho(x) u = f \]

rather than a backward march from terminal data. Here T is the transport operator built from Transport<N> and \(\rho\) is the reaction coefficient.

Relationship to the time-dependent contract

The parabolic and elliptic paths share the spatial discretization in finite_difference/discretization.rs. Both consume DimensionKind to select a stencil and Transport<N> to supply the per-dimension plus/minus coefficients and local source. They differ only in the outer solve strategy:

• PdeProblem<N> marches backward in time.
• EllipticProblem<N> assembles the operator once and solves it directly.

This keeps the stencils in one place while making the stationary strategy explicit. The stationary work lives inside finite_difference, not as a sibling module.

Elliptic problem contract

EllipticProblem<N> supplies:

• dimension_kind(dim): the stencil kind for coordinate dim.
• transport(state, derivs): the plus/minus transport coefficients. The returned source is ignored by the elliptic assembler; forcing belongs in rhs.
• rhs(state): the forcing f.
• reaction(state): a local zero-order coefficient added to the diagonal.
• boundary_conditions(): per-dimension lower and upper conditions.

The operator is assembled as

\[ -\sum_i \big[ a_i^{+} (u_{i,+} - u) + a_i^{-} (u_{i,-} - u) \big] + \rho(x)\, u = f(x), \]

where \(\rho\) is the reaction term. The negative sign follows the standard elliptic convention: for diffusion coefficient \(D\), \(T u = D u''\) and the problem reads \(-D u'' + \rho u = f\). The same DimensionKind stencil table as in solver_pde.md is used for \(a_i^{+}\) and \(a_i^{-}\).

Elliptic control problem contract

For infinite-horizon stochastic control, the solver consumes EllipticControlProblem<N>, which extends crate::models::control::ControlProblem<N> with the same transport /discretization data. Its stationary HJB equation is

\[ 0 = \sup_u \lbrace f(x,u) + \mathcal{L}^u V - r V \rbrace. \]

After optimizing over u, the equation is assembled in the same form as the pure elliptic problem:

\[ -\mathcal{T}^{u^*} V + r V = s^{u^*}, \]

where \(T V + s = \sup_u \lbrace f + \mathcal{L}^u V \rbrace\). The transport method returns plus, minus, and source = s; the discount \(r\) becomes the reaction diagonal.

StationarySolver::solve_control performs policy iteration:

1. Compute derivatives of the current value field.
2. Optimize the control at every node.
3. Assemble the operator and right-hand side from the optimized transport.
4. Solve the linear system.
5. Repeat until the value function stops changing.

The pure EllipticProblem path uses rhs for the forcing and ignores Transport::source; the control path uses Transport::source as the right-hand side.

Boundary conditions

BoundaryCondition is defined in discretization.rs:

• Dirichlet(value): \(u = \text{value}\).
• Neumann(value): \(du/dx = \text{value}\), with the derivative taken along the increasing coordinate direction \(x\), not an outward normal.
• Robin { value, alpha, beta }: \(\alpha u + \beta\, du/dx = \text{value}\).

BoundaryConditions<N> bundles the lower and upper conditions for every dimension. It defaults to Neumann(0) on every side.

Boundary condition implementation status

Solution strategies

StationarySolver assembles a sparse CSR operator and solves it:

• 1D problems use the tridiagonal LAPACK solve.
• N-dimensional problems use successive over-relaxation.

Dirichlet conditions replace the boundary row with the identity. Neumann and Robin conditions use a ghost point. The missing neighbour coefficient is folded into the existing interior neighbour and produces a right-hand-side (and diagonal, for Robin) correction. The derivative in Neumann and Robin is taken along the increasing coordinate direction \(x\).

Reference problems

The implementation is validated against closed forms in solver/tests/elliptic_benchmarks.rs (pure elliptic) and solver/tests/elliptic_control.rs (the control path). The exact solutions are stated once in Exact Solutions to PDEs; the problems are:

• 1D Laplace (Dirichlet), Poisson, and reaction-diffusion;
• 2D Poisson with separable forcing;
• 1D Laplace with Neumann and Robin lower boundaries.

These six problems cover constant, polynomial, separable multi-dimensional, reaction-dominated, Neumann, and Robin boundary behaviour, and verify that the assembler solves the bare elliptic operator \(-\mathcal{T} u + \rho u = f\) rather than an identity-shifted operator.

The control path is validated against the stationary (algebraic Riccati) LQ solution in Exact Solutions in Stochastic Optimal Control, with \(a = -1\), \(b = 1\), \(q = r = 1\), \(\rho = 0.1\), and \(c = 0\). Dirichlet boundary conditions match the exact value at the grid endpoints, so the test isolates the interior operator and the policy-iteration loop.

Economic Invariants

This page is the source of truth for qualitative, economics-driven tests of the solver. It complements exact-solution tests and convergence/consistency tests. Exact-solution tests verify a known closed form; economic-invariant tests verify that the numerical result has the sign, ordering, symmetry, or limit that the model's mathematics predicts.

The exact model formulations and derivations live in Exact Solutions in Stochastic Optimal Control. This page only concerns qualitative properties of the resulting values and controls.

Why these tests exist

Exact-solution tests cannot catch every error:

• A model may solve a well-posed but wrong equation.
• A solver may converge to the wrong sign, monotonicity, or boundary behavior in a regime with no closed form.
• A control optimizer may return finite but economically nonsensical quotes.
• A stochastic solver may produce unstable sign changes hidden by averaging.

Economic invariants are high-signal, low-cost guards against those failures. They are deterministic where possible; stochastic solvers use fixed seeds or repeated runs with a statistical tolerance.

Invariant taxonomy

Every economic test is assigned one or more of these labels.

R1. Symmetry and parity

The value and controls respect the symmetries of the model.

R2. Inventory skew and no-crossing

Quotes create a price incentive to reduce inventory. For the convention used throughout this codebase, a bid fill increases inventory and an ask fill decreases inventory.

R3. Parameter monotonicity and directional sensitivity

A parameter change moves the value or control in the direction predicted by the economics.

R4. Degenerate and limit reduction

A model parameter set to a limit reduces the model to a simpler reference.

R5. Terminal and boundary behavior

The solution satisfies known terminal or boundary facts.

R6. Well-posedness and stability

Results are finite, sign-correct, and repeatable.

R7. Cross-solver consistency

Independent numerical methods agree where their approximations overlap.

Current test layout

The current economic-style tests are spread across the existing model-specific folders. The structure is kept as-is.

Per-model invariants and their mathematics

Each invariant below includes the economic argument that justifies the assertion. The argument follows the model code, not a generic reference.

Avellaneda-Stoikov

The value function separates as

\[ V(t, q, S) = -\exp(-\gamma (X + q S + \theta(t, q))) \]

The optimal half-spreads are

\[ \delta_b = \text{base} - fwd[0] \] \[ \delta_a = \text{base} + bwd[0] \] \[ \text{base} = (1/\gamma) \ln(1 + \gamma/\kappa) \]

where \(fwd[0] = V(q+1) - V(q)\) and \(bwd[0] = V(q) - V(q-1)\).

\(bid(0) = ask(0)\)

At \(q = 0\) the model is symmetric in the jump directions. The inventory gradient term vanishes at zero inventory because \(\theta(t, q)\) is even in \(q\). Hence \(fwd[0] = bwd[0]\) and both spreads equal the base spread.

\(bid(q) = ask(-q)\) and \(ask(q) = bid(-q)\)

The AS model has no drift and no permanent impact, so \(\theta(t, q) = \theta(t, -q)\). Substituting \(q \to -q\) swaps the forward and backward differences, which swaps the two spread formulas.

Zero volatility gives the base spread

When \(\sigma = 0\), the inventory risk term \(-0.5 \gamma \sigma^2 q^2\) vanishes, so \(\theta(t, q)\) is flat and both directional derivatives are zero. The spread therefore equals base.

Positive zero-inventory spread

The base spread is positive because \(\gamma > 0\), \(\kappa > 0\), and \(\ln(1 + \gamma/\kappa) > 0\).

Avellaneda-Stoikov with drift and impact

The reduction tests check that the extensions reproduce the base model in the limit:

• AvellanedaDrift(mu = 0) has no drift term in the diagonal of the GLT matrix.
• AvellanedaImpact(xi = 0) has no asymmetric off-diagonal.

Both therefore coincide with AvellanedaStoikov. The tests compare values and spreads at \(q = 0\).

Hawkes (unilateral)

The unilateral Hawkes intensity evolves as

\[ d\lambda = \beta (\mu - \lambda) \, dt + \alpha \lambda \, dt \]

in the deterministic PDE approximation.

\(alpha = 0\) reduces to Avellaneda-Stoikov

With \(\alpha = 0\) and \(\lambda = \mu\), the intensity stays at \(\mu\) and the model is exactly AvellanedaStoikov(a = mu). The tests compare FD spreads against AvellanedaExact.

Positive value at zero inventory

At \(q = 0\) and \(\lambda = \mu\) the market maker earns the base spread on incoming orders, so \(V > 0\).

Inventory symmetry

For a fixed \(\lambda = \mu\), the unilateral model has no inventory drift or asymmetry; \(V(q) = V(-q)\).

Control intensity scales with state intensity

The fill intensities are proportional to the state intensity lambda:

\[ \lambda_b = \lambda \exp(-\kappa \delta_b) \] \[ \lambda_a = \lambda \exp(-\kappa \delta_a) \]

The test verifies that the total control intensity reacts to the state intensity, with finite and non-negative bounds.

Bilateral Hawkes

The value is a function of inventory \(q\) and two intensities \((\lambda_+, \lambda_-)\).

Positive value at zero inventory

Same argument as unilateral Hawkes.

Symmetry in inventory

At \(\lambda_+ = \lambda_-\), the model is symmetric in the buy/sell directions, so \(V(q, \lambda_+, \lambda_-) = V(-q, \lambda_+, \lambda_-)\).

Side sensitivity

A higher sell order-flow intensity \(\lambda_-\) increases the probability of a bid fill (inventory increases and the market maker earns the bid spread), raising \(V\) at \(q = 0\). A higher buy order-flow intensity \(\lambda_+\) has the mirror effect. The tests assert strict monotonicity in each intensity at \(q = 0\).

Terminal conditions

The terminal value under TerminalCondition::Zero is identically zero. Under TerminalCondition::LiquidationCost, the terminal value is

\[ V(T, q) = -|q| \cdot \text{base\_spread} \]

where \(base\_spread = (1/\gamma) \ln(1 + \gamma/\kappa)\). The tests check both forms at zero time steps, so the solver returns the terminal data exactly.

Heston

The Heston market-making model has state (q, v), where v is the variance process. The value function solves a PDE with the generator

\[ \mu_v = v_{\kappa} (v_{\theta} - v) - \gamma \cdot \text{price\_scale} \cdot \rho \cdot v_{\xi} \cdot v \cdot q \] \[ \sigma2_v = v_{\xi}^2 v \]

and the inventory risk penalty

\[ -0.5 \gamma \cdot \text{price\_scale}^2 v q^2 \]

Higher instantaneous variance widens spreads

At \(q = 0\) the spread is symmetric, and the inventory risk penalty is proportional to \(v q^2\). Increasing \(v\) deepens the value function's inventory curvature, so both fwd[0] and bwd[0] grow in magnitude and both spreads widen.

Higher long-run variance widens spreads

The drift \(v_\kappa (v_\theta - v)\) pulls future variance toward \(v_\theta\). Higher \(v_\theta\) means higher expected future variance and therefore higher expected inventory risk, which widens the spreads.

Faster mean reversion tightens spreads when \(v > v_{\theta}\)

If the current variance is above its long-run mean, a larger \(v_\kappa\) pulls it back down faster, so expected future risk is lower. The spreads tighten.

Higher vol-of-vol widens spreads

\(v_\xi\) controls the variance-of-variance. Increasing \(v_\xi\) increases the diffusion term \(0.5 v_\xi^2 v V_{vv}\) in the variance generator, which increases the value function's curvature and therefore the inventory gradient. The test uses a low liquidity \(a = 10\) so the effect is visible.

rho has no effect at terminal

At \(T = 0\) the value function and all of its gradients are zero, so the spreads equal the base spread. The \(\rho\) term enters only through the variance drift \(-\gamma \rho \xi v q\) and therefore has no effect at terminal.

No-crossing inventory skew

At a positive inventory \(q\), the market maker quotes a lower ask and a higher bid to shed inventory; at a negative inventory the skew reverses. Therefore the bid spread is increasing in \(q\) and the ask spread is decreasing in \(q\), with no crossing. This is checked for a range of \(\rho\).

BSDE stability

BSDE Monte Carlo results must be finite, positive, and have bounded coefficient of variation across repeated runs. The rho sweep checks all three.

American put

The American put value satisfies

\[ \min(-V_t - L V + r V, \, V - (K - S)_+) = 0 \]

Deep out-of-the-money

For \(S \gg K\), the payoff \((K - S)_+\) is zero and immediate exercise is never optimal, so \(V = 0\).

Deep in-the-money

For \(S \ll K\), immediate exercise is optimal and \(V = K - S\).

Merton

The Merton problem has log-utility value

\[ V(t, x) = \ln x + \left[ r + 0.5 (\mu - r)^2 / \sigma^2 \right] (T - t) \]

and constant optimal fraction

\[ u^* = (\mu - r) / \sigma^2 \]

The policy is independent of wealth and time; the tests check the closed form and that the optimizer reproduces it.

Separability

The value is ln x + g(t); the tests check the closed form directly.

Linear-quadratic regulator

The finite-horizon LQ value is

\[ V(t, x) = -x' P(t) x - q(t) \]

and the optimal control is

\[ u^*(t, x) = -R^{-1} B' P(t) x \]

The running cost is positive definite, so the value is negative and quadratic in x.

Linear negative feedback

The optimal control is linear in x with a stabilizing negative gain.

Stationary finite value

For the discounted stationary LQ problem with rho > 0, the value \(V(x) = -P x^2 - c^2 P / \rho\) is finite and negative.

Algebraic Riccati root

The Riccati equation residual is checked directly.

Stationary Avellaneda

The stationary value is the principal eigenvector of the GLT operator \(A = -\alpha q^2\) with off-diagonal \(+\eta\).

Positive Perron eigenvector

By Perron-Frobenius theory the principal eigenvector is positive.

\(theta(0) = 0\) normalization

The gauge is fixed by normalizing \(\theta(0) = 0\).

theta decreases away from zero inventory

The diagonal \(-\alpha q^2\) penalizes large inventory, so the principal eigenvector is largest near \(q = 0\) and decreases as \(|q|\) grows. Since \(\theta = (1/\kappa) \ln v\), \(\theta\) is decreasing in \(|q|\).

Implied spreads match the GLT reference

The spreads are computed from the eigenvector and must agree with AvellanedaGueant::stationary_spreads.

Expanding the suite

To add economic coverage for a model or parameter:

1. Add a row to the per-model table above with the invariant, the expected direction, and the one-sentence economic argument.
2. Add a test in the existing model-specific folder.
3. For deterministic FD tests, assert a strict ordering or equality.
4. For stochastic BSDE tests, average enough repeats or use a fixed seed and assert the invariant with a statistical tolerance.
5. For reductions, compare the reduced model and reference model on a representative state, not only at q = 0.
6. Update the table to mark the invariant as covered.

An invariant test must fail when the sign flips. A loose "close enough" assertion on a qualitative sign is not an economic invariant.

BSDE Solver

The BSDE solver approximates the solution of a stochastic optimal control problem by Monte Carlo regression rather than by discretizing the HJB equation on a grid. This page states the mathematical problem the solver approximates, the least-squares projection it applies, and the reduced form used for the market-making models; implementation notes follow the mathematics.

Problem setting

The controlled state \(X_t \in \mathbb{R}^N\) follows

\[dX_t = b(t, X_t, u_t)\,dt + \sigma(t, X_t)\,dW_t,\]

where \(W_t\) is a standard Brownian motion and \(u_t\) is the control. The objective is

\[V(t, x) = \sup_{u} \mathbb{E}\!\left[g(X_T) + \int_t^T f(s, X_s, u_s)\,ds \;\middle|\; X_t = x\right],\]

with running reward \(f\) and terminal value \(g\). The value function is the solution \(Y_t = V(t, X_t)\) of the backward stochastic differential equation

\[-dY_t = f(t, X_t, u^*_t)\,dt - Z_t\,dW_t, \qquad Y_T = g(X_T),\]

where \(u^*\) is the optimal control and \(Z_t\) is the integrand dual to the forward noise [@pardoux1990adapted, @el1997backward]. The control appears in both the forward drift \(b\) and the running reward \(f\); the process is decoupled when \(b\) does not depend on \(u^*\), and coupled otherwise.

The backward driver is the running reward \(f\) alone. The infinitesimal generator of the controlled forward process,

\[\mathcal{L}^u V = b\cdot\nabla V + \tfrac12\operatorname{tr}\!\big(\sigma\sigma^\top D^2 V\big),\]

is already accounted for by simulating the forward SDE under the control; adding it to the backward step would double-count the drift and diffusion.

Recovering the Hessian from the regression fit

The value is reconstructed from the regression fit, so its first and second derivatives are finite differences of the fitted surface, not of the Monte Carlo samples. For a coordinate \(i\) with step \(h_i\), the diagonal second derivative uses the standard three-point stencil

\[ \partial_{x_i x_i} V(x) \approx \frac{V(x + h_i e_i) - 2V(x) + V(x - h_i e_i)}{h_i^2}. \]

For correlated diffusions the mixed second derivative between coordinates \(i\) and \(j\) uses the centered four-corner stencil

\[ \partial_{x_i x_j} V(x) \approx \frac{V_{++} - V_{+-} - V_{-+} + V_{--}}{4 h_i h_j}, \]

where the four offsets are \(\pm h_i e_i\) combined with \(\pm h_j e_j\). This is the same stencil used by the finite-difference path, so the two solvers agree on the meaning of a mixed derivative. The full symmetric Hessian is carried in StateDerivatives::hessian_full and consumed by models whose noise is correlated across coordinates.

Least-squares projection

Simulate \(M\) forward paths \(\lbrace X_n^{(m)} \rbrace\) on a time grid \(0 = t_0 < \dots < t_N = T\) with step \(\Delta t\). Set \(Y_N^{(m)} = g\big(X_N^{(m)}\big)\) and, for \(n = N-1, \dots, 0\), form the target

\[T_n^{(m)} = Y_{n+1}^{(m)} + f\big(t_n, X_n^{(m)}, u^*_n\big)\,\Delta t.\]

The continuation value is projected onto a finite basis \(\lbrace \psi_k \rbrace_{k=1}^{K}\) by ordinary least squares,

\[c^n = \arg\min_{c \in \mathbb{R}^K} \sum_{m=1}^{M}\Big(T_n^{(m)} - \textstyle\sum_{k=1}^{K} c_k\,\psi_k(X_n^{(m)})\Big)^2,\]

and the value is updated to the regression fit \(Y_n^{(m)} = \sum_k c_k^n\,\psi_k(X_n^{(m)})\). The estimator is biased for finite \(M\) and \(K\); its error splits into a projection error from the truncated basis and a Monte Carlo error that scales as \(M^{-1/2}\) [@gobet2005empirical].

Reduced form for market making

The CARA market-making value admits the ansatz

\[V(t, S, q, X) = -\exp\!\big(-\gamma (X + qS + \theta(t, q))\big),\]

so the value is determined by the reduced function \(\theta(t, q)\), which depends on inventory alone and solves a scalar-in-state HJB in \(\theta\). The full derivation and the reduced HJB are in Exact Solutions in Stochastic Optimal Control; only the pieces the BSDE solver needs are restated here.

The optimal half-spreads are the first-order conditions of the reduced HJB,

\[ \begin{aligned} \delta^{b*}(t, q) &= \delta_0 + \theta(t, q) - \theta(t, q+1), \\ \delta^{a*}(t, q) &= \delta_0 + \theta(t, q) - \theta(t, q-1), \end{aligned} \]

with base spread \(\delta_0 = \tfrac{1}{\gamma}\ln\!\big(1 + \tfrac{\gamma}{k}\big)\). The corresponding fill intensities are \(\lambda^{b*} = A e^{-k\delta^{b*}}\) and \(\lambda^{a*} = A e^{-k\delta^{a*}}\).

The reduced value satisfies the backward equation

\[\partial_t\theta + H^*(t, q) - \tfrac{1}{2}\gamma\sigma^2 q^2 = 0, \qquad \theta(T, q) = \theta_T(q),\]

where \(H^*\) is the optimized Hamiltonian of the two fill events and the inventory-risk term \(-\tfrac12\gamma\sigma^2 q^2\) is a local source, not a transport term. The BSDE solver therefore regresses the reduced quantity \(\theta\) and advances each path with the full reduced driver

\[f_\theta(t, q) = H^*(t, q) - \tfrac{1}{2}\gamma\sigma^2 q^2,\]

rather than with \(H^*\) alone; omitting the source term removes the inventory skew. The terminal \(\theta_T(q)\) is model-specific: it vanishes for the base model and equals \(-\tfrac{1}{2}\xi q^2\) for permanent impact, where \(\xi\) is the price shift per fill.

A crucial consistency condition for the reduced value is that the forward pass must not simulate the inventory jumps. Equation \(\partial_t\theta + H^* - \tfrac12\gamma\sigma^2 q^2 = 0\) is a pure backward (deterministic-in-\(q\)) PDE: the jump operator has been eliminated by the CARA ansatz and is contained entirely in \(H^*\). If the forward simulation instead advances \(q\) by its fill events and the regression then fits \(\theta(t+\Delta t, q_{t+\Delta t})\) onto \(q_t\), it fits the jump-convolved continuation rather than \(\theta(t+\Delta t, q_t)\). The result is a deterministic, \(q\)-dependent bias whose size is controlled by the jump probabilities times the inventory curvature of \(\theta\). The solver therefore holds the jump dimensions fixed after each forward step whenever the problem declares itself reduced; the diffusion dimensions may still diffuse because \(\theta\) is independent of them.

Because the reduced value depends only on inventory, the regression basis should in principle span inventory alone. The spread is a second difference of \(\theta\),

\[\delta^{b*}+\delta^{a*} = 2\delta_0 - \big(\theta(q+1) - 2\theta(q) + \theta(q-1)\big),\]

so the spread error equals the negative of the error in the discrete inventory curvature of the fitted value.

This second-difference structure is the source of the solver's principal limitation. Let \(\hat\theta\) be the regression fit and write \(\hat\theta = \theta + \varepsilon\) with estimation error \(\varepsilon\). The curvature error is

\[\hat\theta(q+1) - 2\hat\theta(q) + \hat\theta(q-1) - \big[\theta(q+1) - 2\theta(q) + \theta(q-1)\big] = \varepsilon(q+1) - 2\varepsilon(q) + \varepsilon(q-1),\]

a second difference of \(\varepsilon\). When the three values of \(\varepsilon\) fluctuate independently with scale \(\sigma_\varepsilon\), the curvature error has scale \(\sqrt{6}\,\sigma_\varepsilon\), larger than \(\sigma_\varepsilon\) by a constant factor. The relative spread error is then

\[\varepsilon_{\text{spread}} = \frac{\varepsilon(q+1) - 2\varepsilon(q) + \varepsilon(q-1)}{2\delta_0 - \Delta^2\theta},\]

where \(\Delta^2\theta\) denotes the exact discrete curvature. The numerator is exactly the absolute spread error, and it is \(\sqrt{6}\) times the value noise scale: differentiating a value fit with a second-difference operator amplifies the Monte Carlo noise in \(\varepsilon\).

The projection is unbiased for the value \(\theta\) itself in the \(M \to \infty\), \(K \to \infty\) limit, but differentiating a noisy fit is ill-conditioned in the sense that the second-difference operator is a high-pass filter: it suppresses the smooth, well-identified part of the fit and retains the fluctuating, poorly identified part. This is not a representability limit of the basis — a basis that spans \(\theta\) exactly still leaves a curvature error proportional to the Monte Carlo variance of the regression coefficients — but a variance limit of estimating a derivative from a Monte Carlo regression. It does not vanish by raising the basis degree or by adding paths alone, because the coefficient variance decays only as \(M^{-1/2}\) while the high-pass amplification is order-independent.

Regressing the inventory difference directly

The amplification can be removed by changing the regressed quantity from the value \(\theta\) to its inventory difference. The spread is linear in the difference, so estimating the difference directly converts a second-difference error into a first-difference error.

Define the forward inventory difference

\[\Phi(t, q) = \theta(t, q+1) - \theta(t, q).\]

Then the optimal half-spreads of the reduced value are

\[\delta^{b*}(t, q) = \delta_0 - \Phi(t, q), \qquad \delta^{a*}(t, q) = \delta_0 + \Phi(t, q-1),\]

and the total spread is

\[\delta^{b*}(t, q) + \delta^{a*}(t, q) = 2\delta_0 - \Phi(t, q) + \Phi(t, q-1).\]

If the backward regression estimates \(\Phi\) instead of \(\theta\), then the spread is computed from the estimated \(\hat\Phi\) directly, without numerical differentiation.

The gain is visible by repeating the error analysis. Write \(\hat\Phi = \Phi + \eta\) with estimation error \(\eta\). The total-spread error is

\[\varepsilon_{\text{spread}} = -\eta(q) + \eta(q-1),\]

a difference of two estimation errors, not a second difference. When the two values fluctuate independently with scale \(\sigma_\eta\), the spread error has scale \(\sqrt{2}\,\sigma_\eta\), whereas the indirect estimate has scale \(\sqrt{6}\,\sigma_\varepsilon\). The high-pass filter has been replaced by a first-order difference, so the error is controlled by the first-order estimation noise \(\sigma_\eta\) rather than by the second-order difference \(\sigma_\varepsilon\).

Both forms carry the same denominator — the total spread — so the improvement is entirely in the numerator. The indirect form recovers the spread as a second difference of the value fit, whose scale is \(\sqrt{6}\,\sigma_\varepsilon\); the direct form recovers it as a first-order difference of \(\hat\Phi\), whose scale is \(\sqrt{2}\,\sigma_\eta\). For comparable regression noise the direct form therefore reduces the spread error by a factor of order \(\sqrt{3}\), and it does so by removing the second-difference amplification rather than by increasing the basis or the path count.

This is not a statement about any particular basis. It holds whenever \(\Phi\) is regressed as a first-order object rather than recovered as a second difference of \(\theta\). The cost is that \(\Phi\) must be propagated backward through the reduced equation consistently: the backward step for \(\Phi\) couples neighbouring inventory levels, so the resulting regression is no longer a per-level projection of a single value. See [@gobet2005empirical] for the regression formulation underlying the least-squares step; the direct-difference variant is a control-level reformulation of that step.

This direct-difference formulation is not currently implemented. It addresses variance amplification in the second-difference, which is orthogonal to the jump-convolution bias fixed by the reduced forward pass. The current solver regresses \(\theta\) directly and holds the jump dimensions fixed in the forward pass.

Full-value (expanded) market-making models

The expanded models -- Heston, Hawkes, Heston-Hawkes, bilateral Hawkes, and bilateral Hawkes with order-flow imbalance -- are full-value problems, not reduced ones. They do not collapse the value onto inventory through a CARA ansatz; every state dimension evolves and the value is regressed in the full state space. The reduced-driver argument above therefore does not apply, and a different forward/backward separation is required.

Let \(x = (q, z)\) split the state into the inventory \(q\) and the remaining factors \(z\) (variance, one or two Hawkes intensities). The HJB driver of a full-value market-making problem is

\[f(x, u) + \mathcal{L}^u V = H^*(q, z) + \ell(q, z) + \mathcal{T}V,\]

where \(H^*\) is the optimized fill Hamiltonian, \(\ell\) is the local inventory-risk source

\[\ell(q, z) = -\tfrac{1}{2}\gamma\,s(z)^2 q^2,\]

and \(\mathcal{T}V\) is the transport operator of the factor process (the variance drift and diffusion, and the Hawkes intensity drift). The running reward \(f\) equals \(H^*\); the transport \(\mathcal{T}V\) is carried by the forward simulation.

Because these are full-value problems, the forward pass simulates every dimension under the control, so the backward driver must add only the terms that are not already accounted for by the forward transport. Those are exactly the running reward and the local source:

\[f_{\text{BSDE}}(x, u) = H^*(q, z) + \ell(q, z).\]

Including \(\mathcal{T}V\) in the backward step would double-count the factor drift and diffusion, exactly as in the base problem setting, but now the transport is the factor process rather than the reduced inventory jump. The reduced models override bsde_driver to running_reward + generator only because their generator is the pure local source \(\ell\) with no transport; the expanded models must instead override bsde_driver to running_reward + \ell directly, since their generator mixes \(\ell\) with \(\mathcal{T}V\).

A second consistency condition holds for these models. Their forward fill intensities depend on the optimal control, so the forward pass of the coupled (Picard) solver must use the control recovered from the backward regression, not a frozen proxy quote. The expanded models therefore override next_step_controlled so the coupled forward pass draws fill events from the current optimal intensities. In decoupled mode the proxy next_step is still used, which is the correct choice only when the control is assumed exogenous for the forward process.

The inventory dynamics of a full-value market-making model are control-dependent: the fill intensities \(\lambda^{b*}, \lambda^{a*}\) are functions of the optimal control \(u^*\) itself, not exogenous rates. The forward process therefore reads

\[dq_t = dN^{b}_t - dN^{a}_t,\]

where \(N^b, N^a\) are counting processes whose intensities are \(\lambda^{b*}_t\) and \(\lambda^{a*}_t\). A decoupled forward pass replaces these optimal intensities with a frozen proxy (the base-spread symmetric fill), while the backward regression still recovers \(u^*\). The forward and backward steps then disagree on the jump measure, and the regression fits a jump-convolved continuation that biases the recovered inventory gradient. The coupled (Picard) mode removes this inconsistency by re-simulating the forward inventory under the current optimal control each iteration, so the forward jump measure and the backward driver are both evaluated at \(u^*\). The control-dependent forward process is therefore solved in coupled mode, and the decoupled mode is reserved for models whose forward drift and jump intensity are exogenous in the control.

A final consistency condition is discretization-level. Each forward Euler step represents a Bernoulli fill with probability \(\min(\lambda^*\Delta t, 1)\) per side, which saturates at \(1\) when \(\lambda^* \Delta t > 1\). The backward driver must not count a fill rate the forward step cannot produce, so the running reward is bounded at the forward-representable rate \(1/\Delta t\) per fill side before it enters the backward step:

\[f_{\text{BSDE}}(x, u) = \frac{\min(\lambda^{b*}, 1/\Delta t) + \min(\lambda^{a*}, 1/\Delta t)}{\gamma + \kappa} + \ell(q, z).\]

This bound is applied only inside bsde_driver; the raw intensities returned to the spread inversion are untouched, because inverting the bound would corrupt the recovered half-spread. The bound makes the backward reward consistent with the forward fill-probability clamp without changing the optimal control's meaning.

Arbitrary jump sampling

The forward pass of the BSDE solver draws state increments through the model's next_step (decoupled) or next_step_controlled (coupled) methods; the generic solver does not construct jump kernels itself. For a model with arbitrary-amplitude jumps, the shared helper solver::numeric::finite_difference::discretization::sample_jump_kernel selects a transition from a kernel by mapping standard-normal noise through the normal CDF to a uniform draw and picking the first competing exponential to fire. A general-jump model overrides its forward step to use that helper (or the equivalent Bernoulli per side), so the forward dynamics match the kernel assembled for the finite-difference path.

The worked AvellanedaLotSize model exercises this: its inventory jumps by lot_size state units per fill, and its next_step/next_step_controlled apply that amplitude directly. With lot_size = 1 the model reduces to the base Avellaneda-Stoikov reduced-value problem, providing a closed-form oracle against which the general-jump forward and backward machinery is validated.

Currently implemented

• The generic solver is solver::numeric::bsde::BsdeSolver, configured with with_basis, with_regularization, with_initial_range, and with_seed; solve_full_control returns a BsdeSolution whose evaluate reconstructs the fitted value at any state. The regression is Householder QR (faer) with a ridge floor for rank-deficient designs.
• The reduced market-making models AvellanedaStoikov, AvellanedaDrift, and AvellanedaImpact implement bsde_driver to include the local inventory-risk source, so the backward step uses the full reduced driver \(f_\theta\). They also override is_reduced_value to true, so the forward pass holds the inventory (jump) dimension fixed and regresses \(\theta(t+\Delta t, q_t)\) onto \(q_t\).
• The full-value market-making models Heston, AvellanedaHawkes, HestonHawkes, BilateralHawkes, and BilateralHawkesOrderFlowImbalance override bsde_driver to running_reward + \ell (running reward plus the local inventory-risk source), excluding the factor transport, and bound the reward rate at 1/dt per fill side to match the forward fill-probability clamp. They also override next_step_controlled so the coupled forward pass simulates inventory fills under the current optimal control rather than a frozen proxy quote. Their control-dependent inventory dynamics require coupled (Picard) mode for finite-horizon accuracy; the decoupled mode is correct only at the terminal boundary (zero time steps).
• The polynomial basis is a full multi-dimensional polynomial over all \(N\) state components, not inventory alone; see solver/src/numeric/basis.rs.
• Coupled (Picard) mode is implemented for problems whose forward dynamics depend on the control; the decoupled mode remains the default.

Implementation: solver::numeric::bsde

N-Dimensional Models

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

ControlProblem<N> trait

The generic contract is [solver::models::control::ControlProblem]. It is the solver-agnostic abstraction consumed by both the BSDE regression solver and explicit finite-difference steps:

#![allow(unused)]
fn main() {
pub trait ControlProblem<const N: usize> {
    type Control;
    fn optimize(&self, t, state, derivs) -> Self::Control;
    fn running_reward(&self, t, state, control) -> f64;
    fn generator(&self, t, state, control, derivs) -> f64;
    fn driver(&self, t, state, control, derivs) -> f64;
    fn terminal(&self, state) -> f64;
    fn apply_constraint(&self, state, value) -> f64;
    fn discount_rate(&self, state) -> f64;
    fn constant_discount_rate(&self) -> Option<f64>;
    fn next_step(&self, t, state, dt, noise) -> [f64; N];
    fn is_diffusion_dimension(&self, dim) -> bool;
    fn gradient_step(&self, dim) -> f64;
}
}

StateDerivatives<N> carries first and diagonal second derivatives plus forward/backward directional differences. Continuous controls consume grad/hessian; jump controls consume fwd/bwd. For correlated diffusions it also carries the full symmetric Hessian hessian_full, whose off-diagonal entries are the mixed second derivatives d^2 V / dx_i dx_j.

Correlated diffusions and the Heston rho term

The generator of a correlated diffusion is \(\tfrac12\,\mathrm{tr}(D\, \mathrm{Hess}\, V)\) with covariance \(D = \sigma\sigma^{\top}\), so the off-diagonal entries of hessian_full are weighted by the off-diagonal covariance. The Heston model does not use this general path: its value is reduced to the coordinates \([q, v]\) by the CARA ansatz, so the spot \(S\) is not a grid dimension. The spot-variance correlation therefore appears as a closed-form drift correction

\[ \rho\,\xi\,v\,\partial_{v}\partial_{S}V = -\gamma\,q\,\rho\,\xi\,v\,\partial_{v}V, \]

folded into the variance drift. This is exact for the reduced model, not an ad hoc substitution, and it does not need hessian_full because the correlated coordinate is collapsed. The general cross term is exercised instead by models with two (or more) diffusive grid coordinates, such as [solver::models::lq_regulator::LqRegulator] with a full diffusion matrix.

The FD-specific transport contract for implicit and ADI schemes is documented in PDE Solving Methods. It is a separate trait, not part of ControlProblem.

Dimension kinds

The generic trait currently exposes two per-dimension signals:

• is_diffusion_dimension(dim): the coordinate is driven by Brownian diffusion.
• gradient_step(dim): physical finite-difference step for mesh-free stencils.

This is not sufficient to select a stable FD stencil. The PDE contract adds a declarative dimension_kind(dim) with three values: DiscreteJump, Diffusion, and DeterministicDrift. See PDE Solving Methods for the stencils.

Conventions

• Dimension 0 is always inventory for market-making models.
• The generic result carries the model's associated Control type directly; there is no inventory-specific to_spreads on the generic result.
• Spread conversion for market-making lives in solver::models::market_making::{MarketMakingControl, SpreadResult}.

Proven: N=3 Heston-Hawkes

HestonHawkes implements ControlProblem<3> with state [q, v, lambda]:

• dim 0 (inventory): discrete jump-controlled dimension.
• dim 1 (variance): CIR diffusion, is_diffusion_dimension(1) = true.
• dim 2 (intensity): deterministic mean-reverting drift.

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

Neural Methods

This section covers the neural_solver Python package, the mesh-free, differentiable companion to the Rust solver crate. It implements neural methods for backward stochastic differential equations (BSDE) and Hamilton-Jacobi-Bellman (HJB) equations in JAX.

Neural methods are a separate package, not part of the Rust workspace. They are coupled to the Rust solver only through committed reference values: the Rust solver emits exact values, and the JAX tests assert agreement against them. There is no live call graph between the two codebases.

Theory and empirical results are separate

This folder keeps the mathematical formulation and the measured results in distinct files so the two are not conflated:

The Rust numerical methods (finite difference, least-squares Monte Carlo BSDE, analytical solutions) are documented separately in the sections above; they are the rigorous reference against which the neural methods are validated.

Network architectures

Three architectures are implemented. All are built on a single feedforward block, Mlp (Xavier init, selectable activation tanh (default), softplus, relu).

Deep BSDE

Forward-backward SDE: learns a scalar initial value Y0 plus an integrand network on \((t, x)\).

• Z net (Mlp) — the diffusion integrand.
• \Gamma net (Mlp) — the predictable jump integrand, added by the jump extension.

Applied to Black-Scholes and Merton (diffusion), and to jump-LQ and Merton-jump (jump).

• Random-jump-size sampling — the \Gamma path currently samples fixed-amplitude jumps; the log-normal MertonJumpLognormal is reference-only with no bsde_problem. (planned)

Deep HJB (DGM)

Residual minimization over the value function, parameterized directly.

• value net (Mlp) on \((t, x)\).

Applied to Black-Scholes and jump-LQ.

• Full DGM residual for arbitrary problems — only the Black-Scholes and jump-LQ residuals are implemented; the general HJB residual is not. (planned)
• Adaptive activation — a tanh with a learnable per-unit scale and shift, to sharpen near kinks while staying differentiable. (planned)
• Partition-of-unity / local-basis networks — piecewise-smooth global functions for kink capture, discussed on the jump page. (planned)

Neural operator

Learns the solution operator over a family of instances (amortized), rather than one instance at a time.

• FNO — FourierOperator, a stack of spectral_layer Fourier-multiplier layers (a learned symbol).
• DeepONet — a branch Mlp + trunk Mlp factorization.

Applied to the jump-LQ family amortized over the jump amplitude.

• Exact-reduction match — the FNO/DeepONet primitives and the amortized training loop exist, but the stage validating the learned symbol against the exact jump-generator symbol is still in progress (07-neural-operators). (planned)

The architecture and activation details are on the BSDE/DGM and jump theory pages, and the operator architectures on the operators page.

Other planned architectures

Worked out in the neural-architectures research problem:

• Hamiltonian saddle network — value/control as a min-max of a learned Lagrangian, so the control supremum is exact by construction. (planned)
• Monotone neural operator — input-convex-plus-affine layers carrying the viscosity comparison principle. (planned)
• Measure transformer / Wasserstein attention — permutation-invariant attention over the BSDE particle cloud (dimension-free). (planned)
• Backward-causal operator — time-masked attention over reversed-time slices. (planned)
• Characteristic / Hamiltonian-flow transformer — parameterizes along the Hamiltonian characteristics. (planned)
• Actor-critic / learned policy — a learned control network replacing the analytic controls, for problems with no closed-form policy. (planned)
• Graph / equivariant network — a permutation-equivariant head for multi-asset inventory. (planned)
• Physics-informed neural network (PINN) — strong-form residual with explicit boundary/terminal conditioning. (planned)
• Score-based / diffusion generative model — for sampling the forward measure. (planned)
• Kolmogorov-Arnold network (KAN) — an MLP alternative. (planned)
• RNN / GRU sequence network — the cheap alternative to masked attention for the backward recursion. (planned)

Status

Implementation state is tracked in the neural_solver workstream. In short: the package, forward SDE simulation, closed-form reference models, the DGM diffusion-term helpers, the full deep BSDE (Z process) method (validated against Black-Scholes and Merton), and the jump-aware deep BSDE ($\Gamma$ process, validated against the jump-LQ and Merton jump problems) exist. The remaining items are the "(planned)" entries under each architecture in the Network architectures section above.

Neural Method Theory

This section contains the mathematical formulations and method designs for the neural BSDE and HJB approaches. It is the single exception to the ownership rule: unlike the rest of the solver, this section carries its own mathematical theory rather than delegating it to the mathematical reference. It contains no measurements; the measured results live in the empirical results section.

Neural BSDE and DGM Formulation

This page describes the theory of the neural BSDE and deep HJB (DGM) methods. Measured results are in the empirical results section.

The BSDE formulation

For a forward process \(X_t\) and a terminal condition \(g\), the value process \(Y_t\) and its martingale integrand \(Z_t\) satisfy the backward equation

\[ dY_t = -f(t, X_t, Y_t, Z_t)\, dt + Z_t\, dW_t, \qquad Y_T = g(X_T). \]

A deterministic value function \(V(t, x)\) is recovered when \(Y_t = V(t, X_t)\) and \(Z_t = \sigma(t, X_t)\, \nabla V(t, X_t)\), but the neural method does not require solving for \(V\) on a grid. It simulates \(M\) forward paths and learns the functional relationship directly.

Deep BSDE method

The deep BSDE method (Han, Jentzen, E 2018) parameterizes the unknown initial value \(Y_0\) and the \(Z\) process (and, in a controlled problem, the policy) as neural networks. The forward pass simulates the SDE under the current parameters, accumulates the backward equation, and compares the reconstructed terminal value against \(g(X_T)\). Training minimizes the mean squared terminal matching error.

The key implementation detail for performance is that the whole time march is a single compiled, device-resident program. jax.lax.scan runs the sequential backward pass without a Python loop, and the loss plus optimizer step are fused so the host reads only a scalar.

Network architectures

The Z network is a feedforward Mlp with a selectable activation. The choice is exposed rather than hard-coded because it changes both representational capacity and whether a second derivative exists for a residual-based method.

Adaptive activation functions (a tanh with learnable per-unit scale and shift) and partition-of-unity/local-basis networks are refinements that trade smoothness for kink capture; they are discussed further on the jump processes page and are not yet implemented.

Deep HJB (DGM) method

A closely related method parameterizes \(V(t, x)\) directly and minimizes the HJB or PDE residual

\[ d_t V + \sup_u \big\lbrace f(t,x,u) + \mathcal{L}^u V \big\rbrace \]

using first and second spatial derivatives obtained by automatic differentiation. This is the mesh-free complement to the finite-difference path. The diffusion-term helpers (the spatial Hessian diagonal and the \(\tfrac{1}{2}\operatorname{tr}(\sigma\sigma^\top \nabla^2 V)\) term) exist in neural_solver/dgm.py, and neural_solver/dgm.py::train_dgm trains a value network against the residual.

Two residuals are implemented in neural_solver/loss.py:

• black_scholes_residual for the linear Black-Scholes PDE \(V_t + r s V_s + \tfrac12 \sigma^2 s^2 V_{ss} - r V = 0\);
• jump_lq_residual for the jump-LQ PIDE, including the jump operator \(\lambda[V(x+\xi) - V(x)]\) and the analytic control \(u = \tfrac{b}{2r} V_x\).

Each residual is checked by substituting its closed-form solution (which must annihilate it to float32 tolerance), then the trained value network is compared against the closed form. On CPU the DGM path recovers the Black-Scholes call to ~1.1% and the jump-LQ value to ~1.3% relative error.

Neural Methods for Jump Processes

Stochastic control problems with jumps (Poisson, compound Poisson, or more general Levy processes) cannot be handled by the diffusion-only neural BSDE or DGM formulations as-is. A jump adds a nonlocal term to the backward driver and to the HJB equation. This page records the theory; measured results live in the empirical results section.

Why jumps are structurally different

For a diffusion, the infinitesimal generator is a local differential operator. For a jump-diffusion with jump measure \(\nu\), the generator gains a nonlocal integral:

\[ \mathcal{L} v = \mu \cdot \nabla v + \tfrac{1}{2}\operatorname{tr}(\sigma\sigma^\top \nabla^2 v) + \int \big[ v(x+j) - v(x) - j \cdot \nabla v\, \mathbf{1}_{|j|<1} \big] \, \nu(dj). \]

The integral connects the value at \(x\) to the value at \(x + j\) across a continuum of jump sizes. No amount of network depth or width removes it; it must be computed or sampled explicitly.

Jump architecture is two decisions, not one

The phrase "what architecture handles jumps best" conflates two independent choices. They are documented separately here because they have separate trade-offs.

Value and policy networks: represent the kinks

Jump control problems produce value functions with kinks at exercise, continuation, and constraint boundaries. Smooth activations used by DGM and PINNs (softplus, tanh) force global smoothness and approximate kinks poorly, while non-smooth activations (relu, leaky_relu) cannot supply the second derivative the residual requires. The two workable resolutions are:

• Adaptive activation functions (tanh with learnable scale and shift), which sharpen near kinks while remaining differentiable.
• Partition-of-unity or local-basis networks, which are smooth locally while the global function is only piecewise smooth.

The jump operator: compute it, do not learn it

Once the jump measure \(\nu\) is known, the jump operator is a known, convolution-like integral. It is computed directly rather than learned:

• Finite-activity, fixed jump sizes: a finite sum over the jump grid.
• Infinite-activity Levy measures: Gauss-Laguerre quadrature, or truncation plus quadrature.
• Forward simulation: sample the jump from \(\nu\) inside the SDE; no quadrature is needed (natural for jax.lax.scan).

Deep BSDE with jumps

The martingale representation for a jump BSDE adds a predictable jump term \(\Gamma\) alongside the diffusion integrand \(Z\). The backward equation becomes

\[ dY_t = -f(t, X_t, Y_t, Z_t, \Gamma_t)\, dt + Z_t\, dW_t + \Gamma_t\, d\tilde N_t, \qquad Y_T = g(X_T), \]

where \(\tilde N_t\) is the compensated jump process. The neural method learns \((Y, Z, \Gamma)\) jointly; \(\Gamma\) is the jump analogue of the hedge or control and is exactly what a controlled jump problem needs.

Deep HJB (DGM) with jumps

The DGM method parameterizes \(V(t, x)\) directly and minimizes the residual of the partial integro-differential equation (PIDE), which is the HJB equation with the jump integral added to the generator. The jump integral is computed by quadrature over \(\nu\) at each residual evaluation. For the jump-LQ problem (deterministic amplitude xi) this is implemented as neural_solver/loss.py::jump_lq_residual and validated against the closed form; neural_solver/dgm.py::train_dgm trains the value network against that residual.

Jump taxonomy and which method fits

For stochastic control specifically, the deep BSDE route is preferred: the \(\Gamma\) process is the direct control analogue, and the method scales to high state dimension without a grid. The DGM route is preferred when the full value surface on a low- or medium-dimensional domain is wanted.

Relationship to the Rust solver

The Rust solver crate already has an arbitrary jump-kernel stage (docs/project/plan/workstreams/solver/stages/04-arbitrary-jump-processes.md) that generalizes unit inventory jumps to arbitrary amplitudes and rates. neural_solver reuses the same jump-kernel convention and validates its jump BSDE against the Rust finite-difference and analytical reference values, not only against itself.

The minimal jump benchmark implemented in neural_solver is the jump-LQ regulator (neural_solver/src/neural_solver/models/jump_lq.py): scalar state, a single Poisson source, and a closed-form Riccati/affine value, so the Gamma path is validated against an exact reference rather than a sampled one. Two Merton jump targets are also implemented where the optimal control itself depends on the jump:

• merton_jump.py: a deterministic multiplicative jump, with a closed-form quadratic policy and value, validated both against the Rust reference and by the jump deep BSDE (which recovers the value and the constant Gamma).
• merton_jump_lognormal.py: a log-normal jump size, whose policy is the root of a transcendental first-order condition computed by Gauss-Hermite quadrature and bisection. This is a reference-only model: the deep BSDE loop does not yet sample random jump sizes, so no bsde_problem is provided.

Neural Operators

A neural operator learns a map between function spaces rather than between finite-dimensional vectors. This page states the concept and its two main architectures — DeepONet and the Fourier neural operator — on mathematical grounds, and explains why the latter is the natural parameterization for the repository's jump-diffusion control problems. No measurements appear here; they belong in the empirical results section, and the jump connection is worked out in full in the research note.

From functions to function maps

A plain neural network approximates a map \(f : \mathbb{R}^{d_x} \to \mathbb{R}^{d_y}\) between finite-dimensional spaces. An operator-learning problem replaces both ends by function spaces: given an input function \(a \in \mathcal{A}\) (terminal data, a coefficient field, or a jump measure) and an output function \(u \in \mathcal{U}\) (the value function or the control), the object of interest is the solution operator

\[ G : \mathcal{A} \to \mathcal{U}, \qquad G(a) = u . \]

A neural operator \(G_\theta\) parameterizes this map directly. Two properties distinguish it from solving one instance at a time. First, discretization invariance: the same weights act on the input function regardless of the grid or sensor set on which it is represented, so the approximation is a property of the function, not of one mesh. Second, amortization: a single trained operator evaluates a whole family of instances (every terminal condition, every coefficient, every jump measure) at inference cost, which is the pay-off that justifies the training cost. The general theory is in [@kovachki2023neural].

Universal approximation for operators

The theoretical foundation is a classical result of Chen and Chen: a nonlinear continuous operator \(G\) on a compact set of inputs can be uniformly approximated by a finite expansion of the form

\[ G(a)(y) \;\approx\; \sum_{k=1}^{p} b_k(a)\,\phi_k(y), \]

where the \(b_k\) are continuous functionals of the input function \(a\) and the \(\phi_k\) are functions of the output variable \(y\). The two architectures below are two different ways of realizing this factorization with a neural network: DeepONet learns the \(b_k\) and \(\phi_k\) explicitly as two subnetworks, while the Fourier neural operator realizes the operator as a composition of integral kernels rather than a truncated expansion.

Special case 1: DeepONet (branch/trunk)

DeepONet [@lu2021deeponet] realizes the Chen–Chen factorization with two networks and an inner product. A branch net encodes the input function from its values at a fixed set of sensor points \(x_1,\dots,x_m\) into the coefficients \(b_1,\dots,b_p\); a trunk net encodes the query point \(y\) into the basis functions \(\phi_1,\dots,\phi_p\). The operator is then

\[ G_\theta(a)(y) = \big\langle \mathrm{branch}(a(x_1),\dots,a(x_m)),\; \mathrm{trunk}(y) \big\rangle . \]

The output is evaluated at arbitrary query points \(y\): once the branch latent vector is computed, the trunk can be evaluated at any \(y\) without a grid, which is what makes DeepONet mesh-free on the output side. Its input, however, is a fixed set of sensor values, so the input function is represented pointwise rather than spectrally.

Special case 2: Fourier neural operator (spectral)

The Fourier neural operator [@li2021fourier] takes the alternative route of building the operator as a composition of integral-kernel layers separated by pointwise nonlinearities. Each layer maps

\[ v \;\mapsto\; \sigma\!\big(W v + \mathcal{K} v\big), \qquad (\mathcal{K}v)(x) = \int \kappa(x, y)\, v(y)\, dy, \]

where \(W\) is a pointwise (local) linear map and \(\mathcal{K}\) is a nonlocal integral operator with kernel \(\kappa\). The decisive step is to take a translation-invariant kernel \(\kappa(x,y) = \kappa(x-y)\), so that \(\mathcal{K}\) is a convolution. A convolution is diagonalized by the Fourier transform,

\[ \widehat{\mathcal{K} v}(\xi) = \hat\kappa(\xi)\,\hat v(\xi), \]

so the layer reduces to a spectral multiplier with a learned symbol \(R_\theta\):

\[ v \;\mapsto\; \mathcal{F}^{-1}\!\big(R_\theta \cdot \mathcal{F}(v)\big). \]

This is efficient (the kernel is parameterized by a few Fourier modes) and it bakes in the translation-invariance that convolution expresses. The cost is that both input and output must live on a grid on which the transform is defined, so FNO is grid-bound where DeepONet is mesh-free.

Why the jump generator is a Fourier multiplier

The repository's jump-diffusion control problems are exactly the setting in which the FNO's structure is not an approximation but an identity. For a jump measure \(\nu\), the jump generator is

\[ \mathcal{J}v(x) = \int_{\mathcal{X}} \big(v(x+z) - v(x)\big)\,\nu(dz) = (\nu \ast v)(x) - v(x), \]

a convolution with \(\nu\) minus the identity, whose Fourier symbol is \(\widehat{\mathcal{J}v} = (\hat\nu - 1)\,\hat v\) (see the research note for the derivation and its lattice \(\mathbb{Z}^m\) and group \(\mathbb{Z}_2^m\) counterparts). A single FNO layer with \(R_\theta = \hat\nu - 1\) therefore represents the jump generator exactly. The nonlocal term that makes classical grid methods expensive is the term a spectral operator represents natively; the diffusion term (a local differential operator) and the control supremum (a pointwise nonlinearity) are the parts FNO handles less naturally, and they are the reason a jump-aware operator is a composition rather than one layer.

Choosing between the two

They are complementary: the spectral layer captures the convolution/jump structure, and the branch/trunk head handles mesh-free evaluation and non-function inputs. The stage that turns this into code is 07-neural-operators.

Empirical Results

Measured results for the neural methods. The theory behind these methods is in the theory section.

Architecture-specific results

Each page focuses on finding the optimal network size/shape for each problem within a single architecture. Numbers are single-seed and indicative, not a rigorous hyperparameter optimum.

Cross-architecture results

Problems considered

The no-closed-form cases (Heston, Hawkes, correlated-default L3, unequal-\(k\) market making) are not yet in neural_solver; they are the target of the neural-operator stage and are documented in the neural operator theory.

Empirical Results

Measured results for the neural methods. The theory behind these methods is in the theory section. All numbers below are CPU measurements with JAX 0.11.x, recorded on 2026-08-23; they are indicative of correctness, not of GPU throughput.

Deep BSDE validation

The deep BSDE Z-process method is validated against closed forms. In all runs the network is an Mlp((2, 32, 32, 1)) and training uses Adam with learning rate \(10^{-2}\) and a fixed Monte Carlo dataset.

Black-Scholes European call

Parameters: \(S_0 = 100\), \(K = 100\), \(r = 0.05\), \(\sigma = 0.2\), \(T = 1\). Exact price \(10.4506\). 512 paths, 25 steps, 3000 iterations.

Merton log-utility portfolio

Parameters: \(\mu = 0.08\), \(r = 0.03\), \(\sigma = 0.2\), \(X_0 = 100\), \(T = 1\). The known optimal fraction is \(\pi^* = 1.25\); the exact initial value is \(\log(100) + 0.06125 = 4.66642\). The deep BSDE recovers this within the 10% test tolerance.

Architecture comparison

The activation choice was compared on the Black-Scholes problem above (seed 0, same hyperparameters). Reproduced by neural_solver/benchmarks/architecture_study.py.

softplus wins for this smooth problem, consistent with its role as the standard DGM activation. The ordering is not universal: kinked value functions (which jump and control-boundary problems produce) are expected to favor relu or an adaptive activation, as discussed in the jump theory.

Deep HJB (DGM) validation

Black-Scholes European call

The DGM value network minimizes the Black-Scholes PDE residual plus a terminal-condition term. Parameters: \(S_0 = 100\), \(K = 100\), \(r = 0.05\), \(\sigma = 0.2\), \(T = 1\). Exact price \(10.4506\). 1024 interior and 256 terminal collocation points, Mlp((2, 32, 32, 1)) with softplus, 3000 iterations, learning rate \(10^{-3}\).

Jump linear-quadratic regulator

The jump-LQ problem is the minimal benchmark for the Gamma path: the state and control are scalar, a single Poisson jump source drives the jump term, and the value has a closed-form Riccati/affine solution. Parameters: \(a = -0.5\), \(b = 1\), \(\sigma = 0.2\), \(q = 1\), \(r = 1\), \(q_T = 1\), \(\lambda = 1\), \(\xi = 0.25\), \(x_0 = 0.5\), \(T = 1\). 512 paths, 25 steps, 3000 iterations.

The recovered values are negative as required, and the relative errors are comparable to the diffusion deep BSDE on a smooth problem. The exact value -0.37630 comes from the Riccati solution P = 0.65345, m = 0.21258, n = 0.10664 at tau = 1; the asymmetric jump (xi = 0.25) is what makes the linear coefficient m non-zero and therefore exercises the affine correction in the control.

Integrand and correction validation

Matching the scalar initial value alone cannot detect a network that learns the no-jump -P x^2 part and drops the jump correction, because the affine m x + n term is small relative to the total value. Two stronger checks isolate the jump path directly:

The learned Gamma is strictly negative and non-zero, so the jump correction is not silently absorbed into the value. These are evaluated on the reachable region (near \(x_0\), where the Monte Carlo paths concentrate); off-manifold points degrade as expected because the network never saw them. The checks live in neural_solver/tests/test_bsde_jump.py and use a documented 10% mean relative-error band (the terminal-matching loss does not supervise Z or Gamma directly).

Merton jump portfolio

The Merton log-utility portfolio with a deterministic multiplicative jump is the first jump target where the optimal control itself depends on the jump. Parameters: \(\mu = 0.08\), \(r = 0.03\), \(\sigma = 0.3\), \(\lambda = 1\), \(y = 1.1\), \(X_0 = 100\), \(T = 1\). The optimal fraction rises from the no-jump value \(0.5556\) to \(u^* = 1.5201\) because the upward jump makes the risky asset more attractive. The exact initial value is \(4.74870\).

Because Z and Gamma are constants for this problem (\(Z = \sigma u^*\) and \(\Gamma = \ln(1 + u^*(y-1))\)), the method converges far more tightly than on the jump-LQ problem. The learned Gamma is positive and tracks the exact constant to within 10%, confirming the jump correction is captured rather than absorbed. See neural_solver/tests/test_merton_jump.py.

For the log-normal jump size (\(\ln Y \sim \mathcal{N}(m, \delta^2)\) with \(m = -0.02\), \(\delta = 0.05\)), the optimal fraction is a numerical root, \(u^* = 0.3386\) (below the no-jump \(0.5556\), since the jumps are downward on average), and the exact initial value is \(4.64049\). This is a reference-only model: it is validated against the Rust reference (neural_solver/tests/test_merton_jump_lognormal.py) but is not trained by the deep BSDE loop, which does not yet sample random jump sizes.

Reference agreement (no training)

The JAX closed forms reproduce the Rust solver reference values within 1e-6 (float32) for Merton policy and value, and Black-Scholes positivity. See neural_solver/tests/test_reference_match.py.

Neural Solver Usage

The neural_solver package is a GPU-native, mesh-free companion to the Rust solver crate. It implements neural methods for backward stochastic differential equations (BSDE) and Hamilton-Jacobi-Bellman (HJB) equations in JAX.

This page covers the package itself. The theory is on the theory pages and measured results are in the empirical results section.

The Rust solver crate is the rigorous reference: it solves the same problems with finite-difference policy iteration, least-squares Monte Carlo BSDE, and analytical solutions, all validated against closed forms. neural_solver adds a differentiable, device-resident path that targets problems where a fixed grid becomes impractical (high state dimension).

Why a separate package

neural_solver is a Python package, not a Cargo workspace member. JAX targets Linux with NVIDIA (CUDA) or AMD (ROCm) accelerators, independent of the Rust toolchain. The two codebases are coupled only through committed reference values, not a live call graph: the Rust solver emits exact values, and the JAX tests assert agreement against them.

Design principles

• Device-resident training. The forward simulation, backward driver, loss, and optimizer step all run on-device; the host reads only scalars.
• Differentiate through the problem. Value and control networks are differentiated with JAX autodiff, not finite differences.
• Validated against the Rust solver. Every neural method is checked against exact or committed reference values.

Layout

• sde.py: differentiable forward SDE simulation via jax.lax.scan. The diffusion-only euler_maruyama and the jump-aware euler_maruyama_jump (single compound Poisson source) are implemented.
• bsde.py: deep BSDE solver. train_deep_bsde handles diffusions (full Z-process method); train_deep_bsde_jump adds a learned Gamma process and jump sampling.
• networks.py: tanh/softplus/relu Mlp used by the deep BSDE and DGM methods.
• loss.py: terminal-matching loss.
• dgm.py: deep HJB (DGM) helpers: the correct spatial Hessian diagonal and the diffusion term. The full DGM residual is not yet implemented.
• models/: one file per reference problem (Black-Scholes, Merton, Merton with deterministic jumps, Merton with log-normal jumps, LQ, jump-LQ, Avellaneda-Stoikov). jump_lq.py provides the Poisson-jump LQ regulator with its Riccati closed form and a decoupled BSDE problem for validating the Gamma path. merton_jump.py provides the deterministic-jump deep BSDE target; merton_jump_lognormal.py is reference-only because the deep BSDE loop does not yet sample random jump sizes.

Runnable Python examples live in neural_solver/README.md; they are not embedded here because the mdbook examples must remain valid Rust.

Installation

See neural_solver/README.md for the full installation matrix. In short:

• CPU only: python -m pip install -e .
• NVIDIA (CUDA 12): python -m pip install -e ".[cuda]"
• AMD (ROCm 7.14, Linux): install rocm[libraries,device-gfx1201] and jax_rocm7_plugin/jax_rocm7_pjrt from the AMD wheel index, then jax==0.10.0/jaxlib==0.10.0, then python -m pip install -e ..

JAX-on-ROCm is Linux-only. Native Windows and WSL2 are not currently distributed for JAX, even though the ROCm runtime itself supports Windows for the RX 9070 (gfx1201).

An unattended Ubuntu Desktop 24.04 installer that provisions amdgpu, ROCm 7.14.0, GitHub Desktop, and Zed is provided in neural_solver/autoinstall/; see its README.md.

Reference agreement

Neural methods are validated against exact reference values produced by the Rust solver crate. The first benchmark in the ladder is pure agreement, with no training: the JAX closed forms for Merton, Merton with jumps, LQ, and Avellaneda-Stoikov are compared against the Rust exact values. Black-Scholes is validated against its own JAX closed form; the Rust generator emits no Black-Scholes entry.

Reference schema

The committed fixture neural_solver/tests/data/reference.json has five entries, each produced by solver/examples/reference/emit_reference.rs:

The scalar LQ entry uses the Rust LqRegulator<1,1> with a=-0.5, b=1, c=0.2, q=1, q_terminal=1, r=1, horizon=1. The Avellaneda-Stoikov entry uses AvellanedaExact with gamma=0.5, sigma=0.5, kappa=1.5, a=140, terminal_time=1, q_max=10. The theta values are gauge-relative (see the issue recorded in docs/src/project/issues.md); the spreads are gauge-invariant and are the physical outputs.

Generating reference samples

Reference values are produced by the Rust example solver/examples/reference/emit_reference.rs, which prints a single JSON object to stdout. To regenerate the committed fixture, run from the repository root:

cargo run -p solver --example emit_reference > neural_solver/tests/data/reference.json

Then verify the JAX closed forms still agree:

cd neural_solver && python -m pytest tests/test_reference_match.py

The fixture neural_solver/tests/data/reference.json is checked in so the JAX tests run without requiring a Rust toolchain. Regenerate it whenever the reference problems, parameters, or exact solutions change.

See the neural_solver workstream for the stages and their definitions of done.

Usage

Runnable usage examples live in neural_solver/README.md; they are not embedded here because the mdbook examples must remain valid Rust. The main entry points are:

• neural_solver.sde.euler_maruyama for differentiable forward simulation.
• neural_solver.bsde.train_deep_bsde for the full deep BSDE training loop.
• neural_solver.dgm.hessian_diagonal and neural_solver.dgm.diffusion_term for the DGM diffusion helpers.
• neural_solver.models for closed-form reference values.

See neural_solver/README.md for installation and a runnable example.

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.

Signal-based strategies

• SignalEngineStrategy: generic engine strategy driven by one or more PriceStrategy instances. Wraps market_model signals (MaCrossover, RsiStrategy) and translates them to exchange orders. Supports:
  • Multiple signal sources with AND/OR/Majority voting (SignalCombinator).
  • Configurable position sizing: fixed quantity or fraction of wealth.
  • Convenience constructors: new_ma(fast, slow, qty), new_rsi(window, oversold, overbought, qty). This strategy is a heuristic — no HJB or optimal-control derivation.

Baseline strategies

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

• KellyStrategy: growth-optimal (Kelly) market making. CARA-based approximation — not a rigorous log-utility HJB solution. Estimates drift and volatility online from recent prices via EWMA. Computes target inventory q* = mu/(gamma * sigma^2) and shifts the AS reservation price toward q* instead of zero. Uses CARA spread dynamics with Kelly-inspired inventory bias.

• KellyRigorousStrategy: rigorous Kelly (log-utility) market making. Rigorous solution. Uses precomputed 3D lookup tables [q, x, tau] from solving the true log-utility HJB via finite difference policy iteration. No CARA approximation. x = cash/mid is the wealth ratio. Contrasts with KellyStrategy which uses the CARA ansatz with a Kelly-inspired target shift.

• 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. Supports GBM, Heston, and Bates market processes.

Configuration

VecEnvConfig specifies the full environment. The process field selects the market process via the VecEnvProcess enum:

#![allow(unused)]
fn main() {
pub enum VecEnvProcess {
    GBM {
        mu: f64,
        sigma: f64,
        initial_price: f64,
    },
    Heston {
        mu: f64,
        kappa: f64,
        theta: f64,
        sigma_v: f64,
        rho: f64,
        initial_price: f64,
        initial_variance: f64,
    },
    Bates {
        mu: f64,
        kappa: f64,
        theta: f64,
        sigma_v: f64,
        rho: f64,
        lambda_jump: f64,
        mu_jump: f64,
        sigma_jump: f64,
        initial_price: f64,
        initial_variance: f64,
    },
}
}

Other configuration fields:

State space

Depends on the process type:

• GBM (5D): [mid, inventory, hawkes_buy, hawkes_sell, time]
• Heston (6D): [mid, inventory, variance, hawkes_buy, hawkes_sell, time]
• Bates (6D): [mid, inventory, variance, hawkes_buy, hawkes_sell, time]

Usage

Rust

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

let config = VecEnvConfig {
    process: VecEnvProcess::Heston {
        mu: 0.0,
        kappa: 2.0,
        theta: 0.04,
        sigma_v: 0.3,
        rho: -0.7,
        initial_price: 100.0,
        initial_variance: 0.04,
    },
    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: (64, 2) flattened row-major [bid_dist, ask_dist, ...]
let actions = vec![0.01f32; 128];
let (next_states, rewards, dones) = env.step(&actions);

// State dimension
assert_eq!(env.state_dim(), 6); // Heston: 6D
}

Python

from market_solve import MarketVecEnv

# Heston (default)
env = MarketVecEnv(num_envs=64, max_steps=10000)

# GBM
env_gbm = MarketVecEnv(
    num_envs=64,
    process_type="gbm",
    mu=0.05,
    sigma=0.2,
)

# Bates
env_bates = MarketVecEnv(
    num_envs=64,
    process_type="bates",
    lambda_jump=0.1,
    mu_jump=-0.02,
)

states = env.reset()                # (64, state_dim) float32
states, rewards, dones = env.step(actions)  # actions: (64, 2)

Mathematical Reference

Ownership rule

Exact solutions and mathematical background belong in this section. A derivation, closed form, sign convention, process definition, or numerical method's mathematics is stated once here and referenced from elsewhere; it is never re-derived or restated. The solver, market_model, and engine pages document implementation only — contracts, configuration, and usage — and link back here for the mathematics they realize.

The single exception is the neural theory section, which carries its own method theory (deep BSDE, DGM, jump processes, neural operators) and is not folded into this section.

Canonical pages

The canonical mathematics lives in four self-contained pages, each stating its results once and without implementation detail:

• Stochastic Processes — the driving equations, parameters, and moments of every process (geometric Brownian motion, Ornstein-Uhlenbeck, CIR, Heston, jump diffusion, Bates, Hawkes, rough OU).
• Exact Solutions in Stochastic Optimal Control — the Avellaneda-Stoikov market-making family (CARA separation, reduced HJB, matrix-exponential solution, drift/impact and stationary limits), the Merton portfolio (diffusion, deterministic jumps, log-normal jumps), and the linear-quadratic regulator (correlated and with Poisson jumps).
• SOC Models Without Exact Solutions — SOC models without a closed form: the Heston stochastic-volatility and Hawkes order-flow reduced PDE/ODE systems and the American-put optimal-stopping problem, solved numerically and validated by degeneracy reduction.
• Exact Solutions to PDEs — the parabolic (heat, Black-Scholes, convection-diffusion) and elliptic (Laplace, Poisson, reaction-diffusion) benchmark problems used to validate the numerical solvers.

These pages carry the sign conventions and validation targets that the numerical solvers must reproduce. When a test, doc comment, or strategy needs a canonical formula, it links here rather than restating it.

Supplementary pages:

• Optimal Criteria — asymptotic approximations and the Kelly criterion, which are not exact solutions.
• Numerical Schemes — the finite-difference and BSDE methods whose convergence is validated against the exact solutions above.
• Backtest Framework — the simulation-layer accounting.
• Rigor Reference — the per-model mapping from problem to method and rigor level.
• Symmetries in Stochastic Optimal Control — a catalog of the symmetries and limit-induced symmetries of the SOC problems, cross-referencing the exact-solutions page.
• Limits of Stochastic Optimal Control Models — a catalog of the models obtained as limits of a more canonical model, cross-referencing the exact-solutions page.
• Cointegration — integrated and cointegrated series, the vector error-correction model and its Gaussian rank inference, the exact law of a mean-reverting relation, switching filters, the high-dimensional spectral laws, and the optimal stopping of a spread.
• Cointegration: Proof Obligations — the unproven statements the cointegration programme depends on, each with what must be proven, why, and what would refute it.
• Bibliography — citations.

The implementation-to-solution mapping (which code realizes which exact solution) is in the solver index.

Stochastic Processes

This page states the stochastic processes used throughout the project as mathematical objects: their driving equations, parameters, and the closed-form moments used for validation. Each process is defined once here. Simulation schemes and implementation notes are documented in the solver and market_model sections, not here.

Throughout, \(W_t\) is a standard Brownian motion and \(dW_t\) its increment. See [@shreve2004stochastic, @oksendal2013stochastic] for the stochastic calculus foundations.

Geometric Brownian motion

The log-normal price model with drift \(\mu \in \mathbb{R}\) and volatility \(\sigma > 0\):

\[dS_t = \mu S_t\,dt + \sigma S_t\,dW_t.\]

The strong solution is

\[S_t = S_0 \exp\!\Big(\big(\mu - \tfrac12 \sigma^2\big) t + \sigma W_t\Big),\]

with moments

\[\mathbb{E}[S_t] = S_0 e^{\mu t}, \qquad \operatorname{Var}(S_t) = S_0^2 e^{2\mu t}\big(e^{\sigma^2 t} - 1\big).\]

Ornstein-Uhlenbeck

The mean-reverting Gaussian diffusion with speed \(\theta > 0\), long-run mean \(\mu\), and volatility \(\sigma > 0\):

\[dX_t = \theta(\mu - X_t)\,dt + \sigma\,dW_t.\]

The solution is Gaussian with moments

\[\mathbb{E}[X_t] = X_0 e^{-\theta t} + \mu\big(1 - e^{-\theta t}\big), \qquad \operatorname{Var}(X_t) = \frac{\sigma^2}{2\theta}\big(1 - e^{-2\theta t}\big).\]

As \(t \to \infty\) the distribution converges to \(\mathcal{N}\big(\mu, \tfrac{\sigma^2}{2\theta}\big)\). See [@vasicek1977equilibrium] for the original term-structure application.

Cox-Ingersoll-Ross

The non-negative mean-reverting diffusion used for variance and intensity, with speed \(\kappa > 0\), long-run level \(\theta > 0\), and vol-of-vol \(\sigma > 0\):

\[dv_t = \kappa(\theta - v_t)\,dt + \sigma\sqrt{v_t}\,dW_t.\]

The process stays non-negative when the Feller condition \(2\kappa\theta \ge \sigma^2\) holds. Its moments are

\[\mathbb{E}[v_t] = v_0 e^{-\kappa t} + \theta\big(1 - e^{-\kappa t}\big),\]

\[ \begin{aligned} \operatorname{Var}(v_t) &= v_0\frac{\sigma^2}{\kappa}\big(e^{-\kappa t} - e^{-2\kappa t}\big)\\ &\quad + \frac{\theta\sigma^2}{2\kappa}\big(1 - e^{-\kappa t}\big)^2. \end{aligned} \]

Heston stochastic volatility

A two-factor model coupling a log-normal price \(S_t\) to a CIR variance \(v_t\), with price drift \(\mu\), mean-reversion speed \(\kappa\), long-run variance \(\theta\), vol-of-vol \(\sigma\), and spot-variance correlation \(\rho \in [-1,1]\):

\[ \begin{aligned} 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, \end{aligned} \qquad d\langle W^S, W^v \rangle_t = \rho\,dt. \]

The variance moments are the CIR moments above; the price has \(\mathbb{E}[S_t] = S_0 e^{\mu t}\). The correlation enters only through the quadratic covariation of the two Brownian drivers. See [@heston1993closed] for the original model and [@lord2010comparison] for simulation-scheme comparisons.

Jump diffusion (Merton)

Geometric Brownian motion with compound Poisson jumps. The price follows

\[ \begin{aligned} dS_t &= \mu S_t\,dt + \sigma S_t\,dW_t + dJ_t,\\ dJ_t &= (e^{Z} - 1)\,dN_t, \qquad Z \sim \mathcal{N}(\mu_J, \sigma_J^2), \qquad N_t \sim \operatorname{Poisson}(\lambda t), \end{aligned} \]

where \(N_t\) is a Poisson process with intensity \(\lambda\) and the jump sizes \(e^{Z} - 1\) are log-normal with \(Z\) normal of mean \(\mu_J\) and variance \(\sigma_J^2\), independent of the diffusion.

Bates

Heston stochastic volatility with Merton log-normal jumps on the price:

\[dS_t = \mu S_t\,dt + \sqrt{v_t}\,S_t\,dW_t^S + dJ_t,\]

with the same jump structure as jump diffusion above and the same variance dynamics as Heston. The model reduces to Heston when \(\lambda = 0\).

Hawkes process

A self-exciting point process with conditional intensity

\[\lambda^*(t) = \mu + \sum_{t_i < t} \phi(t - t_i) = \mu + \int_{-\infty}^{t-} \phi(t - s)\,dN_s,\]

where \(\mu > 0\) is the baseline intensity and the kernel \(\phi \ge 0\) encodes self-excitation. Two kernels are used:

• Exponential: \(\phi(t) = \alpha e^{-\beta t}\) with \(\alpha, \beta > 0\), equivalent to the SDE representation \(d\lambda(t) = \beta\big(\mu - \lambda(t)\big)\,dt + \alpha\,dN_t\).
• Power law: \(\phi(t) = \dfrac{\alpha}{(c + t)^p}\) with \(p > 1\).

Stationarity requires the branching ratio to be subcritical,

\[\int_0^\infty \phi(u)\,du < 1,\]

which for the exponential kernel is \(\alpha/\beta < 1\). The stationary mean intensity of the exponential kernel is

\[\lambda^* = \frac{\mu}{1 - \alpha/\beta}.\]

See [@hawkes2018hawkes] for the finance review and [@ogata1988statistical] for the thinning (Ogata) simulation algorithm.

Rough Ornstein-Uhlenbeck

Fractional Ornstein-Uhlenbeck with Hurst parameter \(H \in (0,1)\) and mean-reversion speed \(\theta\):

\[dX_t = -\theta X_t\,dt + \sigma\,dW_t^H,\]

where \(W_t^H\) is a fractional Brownian motion. The driving kernel \(t^{H - 1/2}\) is approximated by a finite sum of exponentials (Abi Jaber 2019), yielding a multi-factor Markovian representation whose accuracy is controlled by the number of factors.

Exact Solutions in Stochastic Optimal Control

This page derives the exact solutions to the stochastic optimal control problems used as validation targets: the Avellaneda-Stoikov market-making family (base, risk-neutral, drift, impact, and stationary), the Merton portfolio, and the linear-quadratic regulator. Each problem is stated with its dynamics, objective, and exact value function and policy. The problems without a closed form (the Heston and Hawkes market-making extensions and the American put) are on SOC Models Without Exact Solutions. The pure (control-free) PDE benchmarks — parabolic and elliptic — are on the PDE page; the underlying processes are on the processes page.

The general framework follows [@pham2009continuous, @fleming2006controlled, @yong1999stochastic].

Market making (Avellaneda-Stoikov)

Notation

Setting

The market maker controls the bid and ask half-spreads \(\delta_t^b, \delta_t^a\) relative to the mid-price \(S_t\):

\[S_t^{\text{bid}} = S_t - \delta_t^b, \qquad S_t^{\text{ask}} = S_t + \delta_t^a.\]

The state is \([S, q, X]\), with mid-price \(S\), inventory \(q\), and cash \(X\), driven by

\[ \begin{aligned} dS_t &= \mu(t, S_t)\,dt + \sigma(t, S_t)\,dW_t, \\ dq_t &= dN_t^b - dN_t^a, \\ dX_t &= (S_t + \delta_t^a)\,dN_t^a - (S_t - \delta_t^b)\,dN_t^b, \end{aligned} \]

where \(N_t^b, N_t^a\) are point processes with intensities \(\lambda^b(\delta^b)\) and \(\lambda^a(\delta^a)\). The standard specification is exponential,

\[\lambda^b(\delta^b) = A e^{-k\delta^b}, \qquad \lambda^a(\delta^a) = A e^{-k\delta^a}.\]

The objective is the Bolza-form expected terminal utility with a running inventory penalty \(\phi\, q^2\):

\[V(t, S, q, X) = \sup_{\delta^a, \delta^b} \mathbb{E}\!\left[ U\!\left(X_T + q_T S_T - \int_t^T \phi\, q_s^2\, ds\right) \;\middle|\; S_t=S, q_t=q, X_t=X\right].\]

CARA separation and reduced PDE

With CARA utility \(U(x) = -e^{-\gamma x}\) the value separates through the ansatz

\[V(t, S, q, X) = -\exp\!\big(-\gamma (X + qS + \theta(t, q))\big),\]

reducing the problem to a PDE in \(\theta(t,q)\) alone. This reduction is a translation symmetry of the state: the dynamics and terminal payoff depend on the mid-price \(S\) and cash \(X\) only through the mark-to-market wealth \(Y = X + qS\), so the value is invariant under the one-parameter shift \((S, X) \mapsto (S + c,\; X - qc)\) for any constant \(c\). The value therefore depends on \(S\) and \(X\) through \(Y\) alone, and the CARA terminal condition forces the exponential factor in \(Y\), leaving only the inventory correction \(\theta(t,q)\). To see this, set \(\mu = 0\), constant \(\sigma\), and exponential intensities. The HJB is [@bellman1952theory, @bellman1966dynamic]

\[ 0 = \partial_t V + \tfrac12 \sigma^2 \partial_S^2 V + \sup_{\delta^a} A e^{-k\delta^a}\big[V(t,S,q-1,X+S+\delta^a) - V\big] + \sup_{\delta^b} A e^{-k\delta^b}\big[V(t,S,q+1,X-S+\delta^b) - V\big]. \]

Substituting the ansatz, the price-diffusion term yields the quadratic inventory penalty

\[ \partial_t V = -\gamma\,(\partial_t\theta)\,V, \qquad \tfrac12\sigma^2 \partial_S^2 V = -\tfrac12\gamma^2\sigma^2 q^2\,V, \]

and each jump term factors across the state shift. For the ask side,

\[ V(t,S,q-1,X+S+\delta^a) = -e^{-\gamma(X + qS + \delta^a + \theta(t,q-1))} = V\, e^{-\gamma(\delta^a + \theta_{q-1} - \theta_q)}. \]

Writing \(\Delta_a = \theta(t,q-1) - \theta(t,q)\) and \(\Delta_b = \theta(t,q+1) - \theta(t,q)\), the \(-\gamma V\) factor cancels and the HJB reduces to an equation in \(\theta\) alone:

\[ \partial_t\theta - \tfrac12\gamma\sigma^2 q^2 + \sup_{\delta^a} \frac{A e^{-k\delta^a}}{\gamma}\big[1 - e^{-\gamma(\delta^a + \Delta_a)}\big] + \sup_{\delta^b} \frac{A e^{-k\delta^b}}{\gamma}\big[1 - e^{-\gamma(\delta^b + \Delta_b)}\big] = 0. \]

Each supremum is a scalar optimization. Differentiating \(f(\delta) = e^{-k\delta}\big[1 - e^{-\gamma(\delta + \Delta)}\big]\) and setting the derivative to zero gives \(k = (k+\gamma)\, e^{-\gamma(\delta^* + \Delta)}\), hence the optimal spreads

\[ \delta^{a*} = \frac{1}{\gamma}\ln\!\Big(1 + \frac{\gamma}{k}\Big) - \Delta_a = \frac{1}{\gamma}\ln\!\Big(1 + \frac{\gamma}{k}\Big) + \theta(t,q) - \theta(t,q-1), \]

\[ \delta^{b*} = \frac{1}{\gamma}\ln\!\Big(1 + \frac{\gamma}{k}\Big) - \Delta_b = \frac{1}{\gamma}\ln\!\Big(1 + \frac{\gamma}{k}\Big) + \theta(t,q) - \theta(t,q+1). \]

The optimal intensities are \(\lambda^{a*} = A e^{-k\delta^{a*}}\) and \(\lambda^{b*} = A e^{-k\delta^{b*}}\). At optimality \(e^{-\gamma(\delta^* + \Delta)} = \frac{k}{k+\gamma}\), so \(1 - e^{-\gamma(\delta^* + \Delta)} = \frac{\gamma}{k+\gamma}\) and the optimized Hamiltonian collapses to

\[H^* = \frac{\lambda^{a*} + \lambda^{b*}}{\gamma + k}.\]

Substituting back yields the reduced PDE for \(\theta\) (with the running penalty \(\phi\) restored):

\[ \partial_t\theta - \big(\tfrac12\gamma\sigma^2 + \phi\big) q^2 + \frac{A}{k+\gamma}\Big(1 + \frac{\gamma}{k}\Big)^{-k/\gamma} \Big[e^{-k(\theta_q - \theta_{q-1})} + e^{-k(\theta_q - \theta_{q+1})}\Big] = 0, \]

with terminal condition \(\theta(T, q) = 0\) (Mayer form). For the base model \(\phi = 0\).

Finite-horizon exact solution (matrix exponential)

Define \(v_q(t) = \exp(k\,\theta(t,q))\). The nonlinear PDE becomes the linear ODE system

\[ \dot{v}_q(t) = \tilde{\alpha}\, q^2 v_q(t) - \eta\big(v_{q-1}(t) + v_{q+1}(t)\big), \]

with constants

\[ \tilde{\alpha} = \tfrac{k}{2}\gamma\sigma^2 + k\phi, \qquad \eta = \frac{kA}{k+\gamma}\Big(1 + \frac{\gamma}{k}\Big)^{-k/\gamma}. \]

For the base model \(\phi = 0\), so \(\tilde{\alpha} = \tfrac{k}{2}\gamma\sigma^2\), and the terminal condition is \(v_q(T) = 1\).

Stacking \(v_q\) over \(q \in [-q_{\max}, q_{\max}]\) gives a linear system \(\dot{\mathbf{v}}(t) = G\,\mathbf{v}(t)\) with tridiagonal generator \(G_{ii} = \tilde{\alpha} q_i^2\), \(G_{i,i\pm1} = -\eta\). The solution is expressed in terms of the negated generator \(M = -G\), with

\[ M_{i,j} = \begin{cases} -\tilde{\alpha}\, q_i^2, & i = j, \\ +\eta, & |i-j| = 1, \\ 0, & \text{otherwise}, \end{cases} \]

as

\[ \mathbf{v}(t) = e^{M\,(T-t)}\,\mathbf{v}(T). \]

This negated form is the sign convention used throughout: the operator has negative diagonal and positive off-diagonal, and the matrix exponential is evaluated directly in \(T - t\).

Optimal spreads

From \(v_q(t)\), the optimal half-spreads are

\[ \begin{aligned} \delta^{a*}(t, q) &= \delta_0 + \frac{1}{k}\ln\!\Big(\frac{v_q(t)}{v_{q-1}(t)}\Big), \\ \delta^{b*}(t, q) &= \delta_0 + \frac{1}{k}\ln\!\Big(\frac{v_q(t)}{v_{q+1}(t)}\Big), \end{aligned} \qquad \delta_0 = \frac{1}{\gamma}\ln\!\Big(1 + \frac{\gamma}{k}\Big), \]

and the fill intensities are \(\lambda^{a*} = A e^{-k\delta^{a*}}\), \(\lambda^{b*} = A e^{-k\delta^{b*}}\). The base spread \(\delta_0\) is the asymptotic value at \(q = 0\) as \(t \to -\infty\); inventory skew enters through the \(\ln(v_q/v_{q\pm1})\) terms. As \(\gamma \to 0\) the base spread tends to \(1/k\), recovering the risk-neutral (linear-utility) solution below.

Value recovery

\[\theta(t,q) = \frac{1}{k}\ln v_q(t), \qquad V(t,S,q,X) = -\exp\!\big(-\gamma(X + qS + \theta(t,q))\big).\]

Terminal conditions

Under the liquidation-cost condition, the terminal value penalizes nonzero inventory by the base spread, matching the cost of crossing the spread to unwind.

Risk-neutral (linear utility) reduction

Under linear (risk-neutral) utility \(U(x) = x\) the objective is the expected terminal mark-to-market wealth,

\[V(t, S, q, X) = \sup_{\delta^a, \delta^b} \mathbb{E}\big[X_T + q_T S_T \;\big|\; S_t=S, q_t=q, X_t=X\big].\]

The value separates as \(V = X + qS + \theta(t, q)\), the same translation symmetry with an affine (translation-invariant) rather than exponential value. Because \(V\) is affine in \(S\), the price-diffusion term \(\tfrac12\sigma^2\,\partial_S^2 V\) vanishes and with it the inventory-risk penalty. The HJB reduces to

\[0 = \partial_t\theta + \sup_{\delta^a} A e^{-k\delta^a}\big(\delta^a + \Delta_a\big) + \sup_{\delta^b} A e^{-k\delta^b}\big(\delta^b + \Delta_b\big),\]

with \(\Delta_a = \theta(t,q-1) - \theta(t,q)\) and \(\Delta_b = \theta(t,q+1) - \theta(t,q)\) as before. The first-order condition of \(f(\delta) = e^{-k\delta}(\delta + \Delta)\) is \(k(\delta^* + \Delta) = 1\), giving the optimal spreads

\[\delta^{a*} = \frac1k + \theta(t,q) - \theta(t,q-1), \qquad \delta^{b*} = \frac1k + \theta(t,q) - \theta(t,q+1),\]

with base spread \(\delta_0 = 1/k\). The log-linearization \(v_q = e^{k\,\theta(t,q)}\) turns the reduced equation into the linear system

\[\dot{v}_q(t) = -A e^{-1}\big(v_{q-1}(t) + v_{q+1}(t)\big),\]

with no diagonal term, or in the negated form

\[M_{i,j} = \begin{cases} +A e^{-1}, & |i-j| = 1, \\ 0, & \text{otherwise}, \end{cases} \qquad \mathbf{v}(t) = e^{M(T-t)}\,\mathbf{v}(T),\]

with \(v_q(T) = 1\) and value recovery \(V = X + qS + \tfrac1k\ln v_q\). This is the \(\gamma \to 0\) limit of the CARA solution: the diagonal \(\tilde\alpha = \tfrac{k}{2}\gamma\sigma^2\) vanishes and the off-diagonal \(\eta = \tfrac{kA}{k+\gamma}(1+\gamma/k)^{-k/\gamma}\) tends to \(A e^{-1}\), while the base spread \(\tfrac1\gamma\ln(1+\gamma/k)\) tends to \(1/k\).

Drift extension

A constant mid-price drift \(dS_t = \mu\,dt + \sigma\,dW_t\) adds a linear \(-k\mu q\) term to the diagonal of the generator, so in the negated form

\[ M_{i,j} = \begin{cases} -\tilde{\alpha}\, q_i^2 + k\mu\, q_i, & i = j, \\ +\eta, & |i-j| = 1. \end{cases} \]

The optimal target inventory is \(q^* = \dfrac{k\mu}{2\tilde{\alpha}}\); the drift breaks the \(q \mapsto -q\) reflection symmetry of the base generator, shifting the target away from zero.

Permanent-impact extension

A permanent price impact \(\xi\) shifts the mid-price by \(\pm\xi\) on each fill, \(dS_t = \sigma\,dW_t + \xi\,dN_t^a - \xi\,dN_t^b\), producing asymmetric off-diagonals and a nontrivial terminal condition:

\[ M_{i,j} = \begin{cases} -\tilde{\alpha}\, i^2, & j = i, \\ +\eta\, e^{k\xi(i-1)}, & j = i-1, \\ +\eta\, e^{-k\xi(i+1)}, & j = i+1, \end{cases} \qquad v_q(T) = \exp\!\Big(-\tfrac{k}{2}\xi\, q^2\Big). \]

The impact \(\xi\) breaks the reflection symmetry \(q \mapsto -q\) of the base generator: the up- and down-couplings \(e^{k\xi(i-1)}\) and \(e^{-k\xi(i+1)}\) differ once \(\xi \neq 0\), reflecting that a fill moves the mid-price in the direction of the trade. The optimal spreads gain impact adjustments

\[ \begin{aligned} \delta^{a*} &= \delta^{a*}_{\text{base}} - \xi(q - 1), \\ \delta^{b*} &= \delta^{b*}_{\text{base}} + \xi(q + 1). \end{aligned} \]

Stationary (infinite-horizon) solution

As \(T \to \infty\) the finite-horizon solution concentrates on the principal (Perron) eigenvector of \(M\). The undiscounted HJB is translation invariant and has no unique solution of \(\sup_u\lbrace f + \mathcal{L}^u V\rbrace = 0\); its exact limit is the eigenpair

\[M\,\mathbf{v} = \lambda_{\max}\,\mathbf{v},\]

where \(M\) is the symmetric tridiagonal operator with diagonal \(-\tilde{\alpha}\, q^2\) and off-diagonal \(+\eta\) defined above.

Gauge fixing

The eigenvector is defined up to a positive scalar, a scale symmetry \(\mathbf{v} \mapsto c\,\mathbf{v}\) of the eigenproblem \(M\mathbf{v} = \lambda_{\max}\mathbf{v}\). The gauge is fixed by normalizing

\[\theta(0) = 0 \quad\Longleftrightarrow\quad \theta(q) = \frac{1}{k}\ln\frac{v_q}{v_0}.\]

Stationary spreads

\[ \begin{aligned} \delta^{b*}(q) &= \delta_0 + \theta(q) - \theta(q+1), \\ \delta^{a*}(q) &= \delta_0 + \theta(q) - \theta(q-1). \end{aligned} \]

Perron-Frobenius guarantees

The operator \(M\) has non-negative off-diagonals and is irreducible on the truncated inventory grid, so by Perron-Frobenius:

• the principal eigenvalue \(\lambda_{\max}\) is real and simple;
• the principal eigenvector \(\mathbf{v}\) is strictly positive;
• \(\theta(q) = \tfrac1k\ln(v_q/v_0)\) is even in \(q\) and decreasing in \(|q|\), since the diagonal \(-\tilde{\alpha}\, q^2\) penalizes large inventory; the evenness is the reflection symmetry \(q \mapsto -q\) of the generator.

Merton portfolio problem

An investor allocates a fraction \(u \in [0,1]\) of wealth \(x\) to a risky asset with drift \(\mu\) and volatility \(\sigma\), the remainder to a risk-free asset at rate \(r\), to maximize expected log terminal wealth.

Dynamics and objective

\[ \begin{aligned} \frac{dS}{S} &= \mu\,dt + \sigma\,dW, \\ dx &= x\big[r + u(\mu - r)\big]\,dt + x u \sigma\,dW, \\ J(u) &= \mathbb{E}\big[\ln x_T\big]. \end{aligned} \]

HJB equation

The value function \(V(t,x) = \sup_u J(u)\) satisfies

\[ 0 = \partial_t V + \sup_u \Big\lbrace x\big(r + u(\mu - r)\big)\partial_x V + \tfrac12 (x u \sigma)^2 \partial_x^2 V \Big\rbrace, \qquad V(T,x) = \ln x. \]

Exact solution

The log-utility ansatz \(V(t,x) = \ln x + B\,(T-t)\) is exact. The first-order condition in \(u\) gives the constant Merton fraction

\[ u^* = \frac{\mu - r}{\sigma^2}, \]

and substitution yields

\[ V(t,x) = \ln x + \left[ r + \frac{(\mu - r)^2}{2\sigma^2} \right](T - t). \]

The policy is independent of wealth and time. The wealth independence is a scale symmetry of log utility: \(x \mapsto c x\) shifts the value by the additive constant \(\ln c\), so the wealth scale drops out of the first-order condition.

Merton portfolio with deterministic jumps

A jump-diffusion extension in which the risky asset has a deterministic multiplicative jump \(y > 0\).

Dynamics

The risky asset and wealth follow

\[ \begin{aligned} \frac{dS}{S} &= \mu\,dt + \sigma\,dW + (y - 1)\,dN_t, \\ dx &= x\big[r + u(\mu - r)\big]\,dt + x u \sigma\,dW + x u (y - 1)\,dN_t, \end{aligned} \]

where \(N_t\) is a Poisson process of intensity \(\lambda\), so a jump scales wealth by \(1 + u(y-1)\). The objective remains \(J(u) = \mathbb{E}[\ln x_T]\).

HJB equation

\[ 0 = \partial_t V + \sup_u \Big\lbrace x\big[r + u(\mu-r)\big] V_x + \tfrac12 x^2 u^2 \sigma^2 V_{xx} + \lambda\big[ V\!\big(x(1 + u(y-1))\big) - V(x) \big] \Big\rbrace. \]

Exact solution

With \(V(t,x) = \ln x + B\,(T-t)\) the jump term contributes the additive constant \(\lambda \ln\big(1 + u(y-1)\big)\). The value coefficient is

\[ B = r + u^*(\mu - r) - \tfrac12 (u^*\sigma)^2 + \lambda\ln\big(1 + u^*(y-1)\big), \]

and the optimal fraction \(u^*\) solves the quadratic first-order condition

\[ 0 = (\mu - r) - \sigma^2 u^* + \lambda\,\frac{y - 1}{1 + u^*(y-1)}. \]

Setting \(\lambda = 0\) or \(y = 1\) recovers the no-jump Merton fraction \(u^* = (\mu - r)/\sigma^2\) and its value.

Merton portfolio with log-normal jumps

When the jump multiplier \(Y\) is log-normal, \(\ln Y \sim \mathcal{N}(m, \delta^2)\), the jump distribution is a continuum and the policy has no closed form.

Dynamics

\[ dx = x\big[r + u(\mu-r)\big]\,dt + x u \sigma\,dW + x u (Y - 1)\,dN_t, \]

with \(Y\) independent of the Poisson process. Log utility requires \(0 \le u \le 1\) so that wealth stays positive for any jump.

Exact solution (semi-closed form)

The ansatz \(V(t,x) = \ln x + B\,(T-t)\) remains exact, with

\[ B = r + u^*(\mu - r) - \tfrac12 (u^*\sigma)^2 + \lambda\,\mathbb{E}\big[\ln(1 + u^*(Y-1))\big], \]

and \(u^*\) solves the transcendental first-order condition

\[ 0 = (\mu - r) - \sigma^2 u^* + \lambda\,\mathbb{E}\!\left[\frac{Y - 1}{1 + u^*(Y-1)}\right]. \]

There is no closed form for \(u^*\). The expectations are evaluated by Gauss-Hermite quadrature on the standard normal \(z = (\ln Y - m)/\delta\), and \(u^*\) is bracketed on \([0,1]\) and refined by bisection. With \(\delta = 0\) the jump degenerates to the deterministic multiplier \(y = e^m\), reducing to the deterministic-jump case above.

Linear-quadratic regulator (correlated)

A finite-horizon LQ regulator with \(n\)-dimensional state and additive noise through a constant diffusion matrix \(C\).

Dynamics and objective

\[ \begin{aligned} dx &= (A x + B u)\,dt + C\,dW, \\ J(u) &= \mathbb{E}\!\left[\int_0^T \big(x^\top Q x + u^\top R u\big)\,dt + x_T^\top Q_T x_T\right], \end{aligned} \]

where the components of \(W\) are independent standard Brownian motions. The problem is to minimize \(J\), equivalently to maximize \(-J\). \(Q, Q_T\) are positive semidefinite and \(R\) is positive definite. The covariance of the noise is \(D = C C^\top\), so a full (non-diagonal) \(C\) yields correlated coordinates with \(d\langle x_i, x_j\rangle = D_{ij}\,dt\).

HJB equation

\[ 0 = \partial_t V + \sup_u \Big\lbrace -\big(x^\top Q x + u^\top R u\big) + (Ax + Bu)^\top \nabla V + \tfrac12\,\operatorname{tr}\!\big(C C^\top \operatorname{Hess} V\big) \Big\rbrace. \]

Exact solution

The value is quadratic,

\[ V(t,x) = -x^\top P(t)\,x - q(t), \qquad u^*(t,x) = -R^{-1} B^\top P(t)\,x, \]

where \(P(t)\) solves the backward Riccati ODE

\[ \dot{P} = -A^\top P - P A - Q + P B R^{-1} B^\top P, \qquad P(T) = Q_T, \]

and \(q(t)\) is the scalar, state-independent correction

\[ \dot{q} = -\operatorname{tr}\big(C C^\top P\big), \qquad q(T) = 0. \]

Because the noise is additive, \(P(t)\) and the feedback gain are independent of \(C\); the correlation enters the value only through \(\operatorname{tr}(C C^\top P)\) in \(q(t)\), and it cancels out of the optimal control. Setting \(C\) block-diagonal recovers the uncorrelated reduction.

Stationary (algebraic Riccati) reduction

In the infinite-horizon limit with discount rate \(\rho > 0\), the value is stationary and the Riccati ODE reduces to the algebraic Riccati equation. For the scalar problem \(dx = (a x + b u)\,dt + c\,dW\) with running cost \(q x^2 + r u^2\), the value is \(V(x) = -P x^2\) and the control is \(u^*(x) = -\tfrac{b}{r} P x\), where \(P\) is the positive root of

\[ \frac{b^2}{r} P^2 + (\rho - 2a) P - q = 0. \]

Linear-quadratic regulator with Poisson jumps

A scalar jump extension providing the minimal jump-diffusion control benchmark. A single Poisson source drives the jump term.

Dynamics and objective

\[ dx = (a x + b u)\,dt + \sigma\,dW + \xi\,dN_t, \]

where \(N_t\) is Poisson with intensity \(\lambda\), and at each jump the state increments by the amplitude \(\xi\), a random variable with law \(\mu\), mean \(\bar{\xi} = \mathbb{E}[\xi]\), and second moment \(\mathbb{E}[\xi^2]\). The amplitude may be deterministic. The running and terminal costs are \(q x^2 + r u^2\) and \(q_T x_T^2\) with \(q, q_T \ge 0\) and \(r > 0\); the problem maximizes the negative cost, so \(V\) below is the maximized value.

HJB equation

\[ 0 = \partial_t V + \sup_u \big\lbrace -(q x^2 + r u^2) + (a x + b u) V_x \big\rbrace + \tfrac12 \sigma^2 V_{xx} + \lambda \int \big[ V(x+\xi) - V(x) \big]\, \mu(d\xi), \]

with terminal \(V(T,x) = -q_T x^2\).

Exact solution

The ansatz \(V(t,x) = -P(t) x^2 - m(t) x - n(t)\) is exact, since the jump integral maps quadratics to quadratics. Matching powers of \(x\) gives three decoupled ODEs. The quadratic coefficient satisfies the no-jump Riccati equation

\[ \dot{P} = -2 a P - q + \frac{b^2}{r} P^2, \qquad P(T) = q_T, \]

so jumps do not enter \(P\). The linear and constant coefficients satisfy

\[ \dot{m} = -\big(a - \tfrac{b^2}{r} P(t)\big) m - 2\lambda P(t)\, \bar{\xi}, \qquad m(T) = 0, \]

\[ \dot{n} = -\sigma^2 P(t) + \frac{b^2}{4r} m(t)^2 - \lambda P(t)\, \mathbb{E}[\xi^2] - \lambda m(t)\, \bar{\xi}, \qquad n(T) = 0. \]

The optimal control keeps the linear feedback form plus an affine jump correction,

\[ u^*(t,x) = -\frac{b}{r} P(t)\, x - \frac{b}{2r} m(t). \]

For a reflection-symmetric jump law (\(\mu\) symmetric about zero, so \(\bar{\xi} = 0\)) the linear term vanishes, \(m \equiv 0\), so the control is the pure linear feedback \(u^* = -\tfrac{b}{r} P(t)\,x\) and the jump enters the value only through the constant

\[ \dot{n} = -\sigma^2 P(t) - \lambda P(t)\, \mathbb{E}[\xi^2], \qquad n(T) = 0. \]

Setting \(\lambda = 0\) recovers the scalar no-jump regulator: \(m \equiv 0\) and \(\dot n = -\sigma^2 P\) reproduce the diffusion-only correction with \(C = \sigma\). The validation targets are the sign and reduction checks: \(P(t) \ge 0\) throughout (since \(q, q_T \ge 0\)), \(V(t,x) \le 0\), the policy reduces to the no-jump feedback law when \(\lambda = 0\) or \(\bar{\xi} = 0\), and the jump contributes only a negative constant shift \(\lambda P\,\mathbb{E}[\xi^2]\) to the symmetric-case value.

SOC Models Without Exact Solutions

These problems have no closed-form value or policy in general. The Heston and Hawkes market-making extensions remove the closed form of a model that has one; the American put is an optimal-stopping problem whose free boundary has no closed form at all. Each HJB, PDE, or obstacle system is stated here as the ground truth the numerical solvers discretize. Each is validated numerical-versus-numerical (finite difference against least-squares Monte Carlo) and by a degeneracy reduction to a model with a closed form, recorded in Limits of Stochastic Optimal Control Models. The numerical methods are on Numerical Schemes; the closed-form problems are on Exact Solutions in Stochastic Optimal Control.

Heston stochastic-volatility extension

With stochastic volatility, the mid-price variance is a CIR process. Writing \(\nu_t\) for the variance, \(\kappa\) for its mean-reversion speed (distinct from the fill-decay \(k\)), \(\theta_\nu\) its long-run level, \(\xi_\nu\) its vol-of-vol, and \(\rho_\nu\) the spot-variance correlation:

\[ dS_t = \sqrt{\nu_t}\,dW_t^S, \qquad d\nu_t = \kappa(\theta_\nu - \nu_t)\,dt + \xi_\nu\sqrt{\nu_t}\,dW_t^\nu, \qquad d\langle W^S, W^\nu\rangle_t = \rho_\nu\,dt. \]

The multiplicative ansatz \(V(t,x,q,S,\nu) = -\exp\big(-\gamma(x + qS)\big)\, w_q(t,\nu)\) leaves a linear parabolic PDE system for \(w_q(\tau, \nu)\) in \(\tau = T - t\):

\[ \begin{aligned} \partial_\tau w_q &= \tfrac12 \xi_\nu^2 \nu\, \partial_\nu^2 w_q + \big[\kappa(\theta_\nu - \nu) + \rho_\nu \xi_\nu \gamma \nu q\big]\,\partial_\nu w_q - \tfrac12 \gamma^2 \nu q^2 w_q \\ &\quad + A\Big(1 + \tfrac{\gamma}{k}\Big)^{-(1 + k/\gamma)} \Big(w_{q-1}\mathbf{1}_{q > -Q} + w_{q+1}\mathbf{1}_{q < Q}\Big), \end{aligned} \]

with initial condition \(w_q(0, \nu) = 1\). The optimal spreads are

\[ \begin{aligned} \delta^{a*}(t,q,\nu) &= \delta_0 + \frac{1}{\gamma}\ln\!\Big(\frac{w_{q-1}(t,\nu)}{w_q(t,\nu)}\Big), \\ \delta^{b*}(t,q,\nu) &= \delta_0 + \frac{1}{\gamma}\ln\!\Big(\frac{w_{q+1}(t,\nu)}{w_q(t,\nu)}\Big), \end{aligned} \qquad \delta_0 = \frac{1}{\gamma}\ln\!\Big(1 + \frac{\gamma}{k}\Big). \]

See [@heston1993closed] for the model and [@lord2010comparison] for the simulation schemes.

Hawkes order-flow extension

Self-exciting arrival intensities couple the fill process to its own history [@hawkes2018hawkes]. The intensity is governed by the process defined on the processes page.

Unilateral Hawkes

Each side has its own self-exciting intensity

\[d\lambda_t^j = \beta\big(\mu - \lambda_t^j\big)\,dt + \alpha\,dN_t^j, \qquad j \in \lbrace a, b\rbrace,\]

and the effective fill rate is \(\lambda_t^j e^{-k\delta^j}\). The fluid (mean-field) approximation replaces the discrete jumps with a continuous drift,

\[\frac{d\lambda_t}{dt} = \beta(\mu - \lambda_t) + \alpha\lambda_t,\]

with stationary mean \(\lambda^* = \dfrac{\mu}{1 - \alpha/\beta}\) and stability condition \(\alpha/\beta < 1\).

Optimal quotes (bilateral)

With state \((q, \lambda^a, \lambda^b)\) the affine ansatz

\[\theta(t, q, \lambda^a, \lambda^b) = C(\tau, q) + h_a(\tau, q)\,\lambda^a + h_b(\tau, q)\,\lambda^b\]

gives coupled ODEs for \(C\), \(h_a\), \(h_b\). The optimal quotes decompose into a baseline inventory-control term plus an excitation adjustment:

\[ \begin{aligned} \delta^{a*} &= \delta_0 + C(\tau, q) - C(\tau, q-1) + h_a(\tau, q)\,\lambda^a - h_a(\tau, q-1)(\lambda^a + \alpha), \\ \delta^{b*} &= \delta_0 + C(\tau, q) - C(\tau, q+1) + h_b(\tau, q)\,\lambda^b - h_b(\tau, q+1)(\lambda^b + \alpha), \end{aligned} \]

where \(\delta_0 = \tfrac1\gamma\ln(1 + \tfrac{\gamma}{k})\).

Order-flow imbalance

The mid-price responds to net order flow with sensitivity \(\eta_{\text{OFI}}\):

\[dS_t = \sigma\,dW_t + \eta_{\text{OFI}}\big(dN_t^b - dN_t^a\big).\]

American put (optimal stopping)

The American put is the optimal-stopping value

\[ V(t, S) = \sup_{\tau \in [t, T]} \mathbb{E}\!\Big[e^{-r(\tau - t)}(K - S_\tau)_+ \;\Big|\; S_t = S\Big], \]

under risk-neutral geometric Brownian motion \(dS/S = r\,dt + \sigma\,dW\). The value satisfies the linear-complementarity (obstacle) formulation

\[ \min\!\Big( -\partial_t V - \mathcal{L} V + r V,\; V - (K - S)_+ \Big) = 0, \]

with generator \(\mathcal{L} V = rS\,\partial_S V + \tfrac12\sigma^2 S^2\,\partial_S^2 V\) and terminal condition \(V(T, S) = (K - S)_+\).

There is no closed-form value for general \(r, \sigma, T\); a reference solution is obtained with a fine-grid finite-difference scheme. In the special case \(r = 0\), early exercise is never strictly optimal and the American value coincides with the Black-Scholes European put [@black1973pricing]:

\[ \begin{aligned} V_{\mathrm{BS}}(S) &= K\,\Phi(-d_2) - S\,\Phi(-d_1), \\ d_1 &= \frac{\ln(S/K) + (r + \tfrac12\sigma^2)T}{\sigma\sqrt{T}}, \\ d_2 &= d_1 - \sigma\sqrt{T}. \end{aligned} \]

The early-exercise boundary \(S^*(t)\), below which immediate exercise is optimal, is obtained numerically rather than in closed form.

Exact Solutions to PDEs

This page states the exact solutions to the pure (control-free) PDEs used to validate the numerical solvers. They isolate the spatial discretization, time integrators, and boundary handling from any control or jump term. The stochastic-optimal-control closed forms — market making, Merton, and linear-quadratic regulator — are on the SOC page.

The problems split into two classes: parabolic problems, which march backward in time from terminal data, and elliptic problems, which have no time variable and are solved directly with boundary conditions.

Parabolic problems

1D heat equation

On \(x \in [0,1]\), with diffusion coefficient \(D\):

\[ \partial_t V + D\, V_{xx} = 0, \qquad V(T, x) = \cos(\pi x). \]

The cosine mode has zero derivative at both endpoints, so the terminal data is compatible with a zero-flux boundary. The exact solution is

\[ V(t, x) = \cos(\pi x)\exp\!\big(-D \pi^2 (T - t)\big). \]

2D separable heat equation

On \((x, y) \in [0,1]^2\) with a single diffusion coefficient \(D\):

\[ \partial_t V + D\,(V_{xx} + V_{yy}) = 0, \qquad V(T, x, y) = \cos(\pi x)\cos(\pi y). \]

Each cosine factor has zero derivative at both endpoints, so the problem is the separable product of two independent 1D heat modes. The exact solution is

\[ V(t, x, y) = \cos(\pi x)\cos(\pi y)\exp\!\big(-2D \pi^2 (T - t)\big). \]

The factor \(2\) in the exponent is the sum of the two identical per-dimension decay rates. Because the two directions share the same coefficient \(D\) and the same terminal mode, the solution is symmetric under the exchange \(x \leftrightarrow y\).

Black-Scholes European call

On \(S \in [0, S_{\max}]\), with risk-free rate \(r\) and volatility \(\sigma\):

\[ \partial_t V + r S\, V_S + \tfrac12 \sigma^2 S^2\, V_{SS} - r V = 0, \qquad V(T, S) = \max(S - K, 0). \]

The exact solution is the Black-Scholes call formula [@black1973pricing]:

\[ V(t, S) = S\,\Phi(d_1) - K e^{-r\tau}\Phi(d_2), \qquad d_1 = \frac{\ln(S/K) + (r + \tfrac12\sigma^2)\tau}{\sigma\sqrt{\tau}}, \qquad d_2 = d_1 - \sigma\sqrt{\tau}, \]

with \(\tau = T - t\). The value is homogeneous of degree one in \((S, K)\), \(V(\lambda S, \lambda K) = \lambda V(S, K)\): a scale symmetry of the Black-Scholes equation, so the delta \(\Phi(d_1)\) depends on \(S\) and \(K\) only through the ratio \(S/K\).

1D convection-diffusion

On the real line, with constant drift \(b\) and diffusion \(D\), the terminal data is a Gaussian centered at \(x_0\):

\[ \partial_t V + b\, V_x + D\, V_{xx} = 0, \qquad V(T, x) = \exp\!\Big(-\frac{(x - x_0)^2}{2 w^2}\Big). \]

The exact solution is

\[ V(t, x) = \left(\frac{w^2}{\operatorname{var}}\right)^{1/2} \exp\!\Big(-\frac{(x + b\tau - x_0)^2}{2\,\operatorname{var}}\Big), \qquad \operatorname{var} = w^2 + 2D\tau, \qquad \tau = T - t. \]

The constant coefficients make the equation translation invariant, so the terminal Gaussian is translated by \(-b\tau\) and its variance grows to \(w^2 + 2D\tau\).

Elliptic problems

Elliptic problems have no time variable; the solution is the direct solve of

\[ -\mathcal{T} u + \rho(x)\, u = f, \]

where \(\mathcal{T}\) is the transport operator and \(\rho\) the reaction coefficient. The boundary conditions are the problem.

1D Laplace equation (Dirichlet)

On \(x \in [0,1]\):

\[ u'' = 0, \qquad u(0) = 0, \qquad u(1) = 1, \]

with exact solution \(u(x) = x\).

1D Poisson equation

On \(x \in [0,1]\):

\[ -u'' = 1, \qquad u(0) = 0, \qquad u(1) = 0, \]

with exact solution \(u(x) = \tfrac12 x(1 - x)\).

2D Poisson equation with separable forcing

On \((x, y) \in [0,1]^2\):

\[ -u_{xx} - u_{yy} = 2\pi^2\sin(\pi x)\sin(\pi y), \qquad u = 0 \text{ on the boundary}, \]

with exact solution \(u(x, y) = \sin(\pi x)\sin(\pi y)\).

1D stationary reaction-diffusion

On \(x \in [0,1]\):

\[ -u'' + a u = 0, \qquad u(0) = 0, \qquad u(1) = \sinh(\sqrt{a}), \]

with exact solution \(u(x) = \sinh(\sqrt{a}\, x)\).

1D Laplace equation (Neumann lower boundary)

On \(x \in [0,1]\):

\[ -u'' = 0, \qquad u'(0) = 0, \qquad u(1) = 1, \]

with exact solution \(u(x) = 1\).

1D Laplace equation (Robin lower boundary)

On \(x \in [0,1]\):

\[ -u'' = 0, \qquad u(0) - u'(0) = 0, \qquad u(1) = 1, \]

with exact solution \(u(x) = \tfrac12(x + 1)\).

These six elliptic problems cover constant, polynomial, separable multi-dimensional, reaction-dominated, Neumann, and Robin boundary behaviour, verifying that the solver assembles the bare elliptic operator \(-\mathcal{T} u + \rho u = f\) rather than an identity-shifted one.

Optimal Criteria

Asymptotic Approximations

Asymptotic methods [@cartea2015algorithmic, @gueant2013dealing] provide closed-form approximations to optimal quotes in limiting regimes.

Small Risk Aversion (\(\gamma \to 0\))

Certainty equivalent expansion: \(\text{CE}(X) \approx \mathbb{E}[X] - \frac{\gamma}{2}\text{Var}(X)\).

Linear quoting rule:

\[ \begin{aligned} \delta^{a*}\_\text{linear}(q) &\approx \frac{1}{k} - \frac{1}{2}\gamma\sigma^2(T - t)(2q - 1) \\ \delta^{b*}\_\text{linear}(q) &\approx \frac{1}{k} + \frac{1}{2}\gamma\sigma^2(T - t)(2q + 1) \end{aligned} \]

Infinite Horizon (\(T \to \infty\))

Stationary quotes from the principal eigenvector of the generator (see Exact Solutions in Stochastic Optimal Control for the matrix formulation):

\[\delta^{a*}_\infty(q) \approx \frac{1}{\gamma}\ln\!\left(1 + \frac{\gamma}{k}\right) - \frac{2q - 1}{2}\sqrt{\frac{\sigma^2\gamma}{2kA}\left(1 + \frac{\gamma}{k}\right)^{1 + k/\gamma}}\]

\[\delta^{b*}_\infty(q) \approx \frac{1}{\gamma}\ln\!\left(1 + \frac{\gamma}{k}\right) + \frac{2q + 1}{2}\sqrt{\frac{\sigma^2\gamma}{2kA}\left(1 + \frac{\gamma}{k}\right)^{1 + k/\gamma}}\]

Kelly Criterion and Optimal Position Sizing

The Kelly criterion maximizes the expected logarithmic growth rate of wealth [@merton1969lifetime]:

\[W^* = \arg\max_w \mathbb{E}[\log(1 + wR)]\]

where \(R\) is the per-period return and \(w\) is the fraction of wealth allocated.

Connection to CARA Utility

CARA utility \(U(x) = -e^{-\gamma x}\) with risk aversion \(\gamma\) gives the certainty equivalent \(\text{CE} = \mathbb{E}[X] - \frac{\gamma}{2}\text{Var}(X)\) for normally distributed \(X\). The optimal position maximizes CE, which for an asset with Sharpe ratio \(\text{SR} = \mu/\sigma\) is:

\[q^* = \frac{\mu}{\gamma\sigma^2}\]

For log-utility (Kelly), the certainty equivalent of \(\log(1 + wR)\) under small \(wR\) approximates \(wR - \frac{1}{2}w^2R^2\). Expanding \(\mathbb{E}[\log(1 + wR)] \approx w\mu - \frac{1}{2}w^2(\mu^2 + \sigma^2)\), the optimal fraction is the classic Kelly formula:

\[w^* \approx \frac{\mu}{\sigma^2 + \mu^2} \approx \frac{\mu}{\sigma^2} \quad (\mu \ll \sigma)\]

Comparing, the CARA optimal position \(q^* = \mu/(\gamma\sigma^2)\) maps to the Kelly fraction when \(\gamma \approx 1/q^*\) (risk aversion inversely proportional to the optimal position). In market making terms, CARA with \(\gamma = 1/\bar{q}\) approximates the growth-optimal strategy for a target inventory scale \(\bar{q}\).

Reservation Price and Target Inventory (CARA Approximation)

Note: this is a CARA-based approximation, not a rigorous Kelly HJB solution. The Avellaneda-Stoikov model with drift uses the CARA (exponential utility) ansatz \(V = -e^{-\gamma(X + qS + \theta(t,q))}\). The HJB separates because CARA factors across state variables — log-utility \(U(x) = \log x\) does not have this property.

In the CARA drift model, the optimal target inventory emerges from the quadratic drift-penalty tradeoff [@avellaneda2008high]:

\[q^* = \frac{k\mu}{2\tilde{\alpha}} = \frac{k\mu}{k\gamma\sigma^2} = \frac{\mu}{\gamma\sigma^2}\]

This matches the Kelly position to first order because CARA with normal returns and log-utility both produce mean-variance optimal positions. However, the value-function curvature (how aggressively the market maker adjusts quotes as inventory deviates from \(q^*\)) differs between the two utility functions. The CARA model uses the AS spread formula; a rigorous Kelly model would produce different spread dynamics.

Rigorous Kelly HJB

The true Kelly market-making problem maximizes expected log-terminal-wealth [@merton1969lifetime]:

\[V(t, S, q, X) = \sup_{\delta^a, \delta^b} \mathbb{E}\!\left[\log\!\left(X_T + q_T S_T\right) \;\middle\|\; S_t=S, q_t=q, X_t=X\right]\]

This objective does not admit the CARA separation ansatz. However, log-utility is homogeneous of degree 0 in \((S, X)\): \(V(t, cS, q, cX) = \log c + V(t, S, q, X)\), a scale symmetry under joint rescaling of price and cash. This admits the reduction:

\[V(t, S, q, X) = \log S + v(t, x, q), \qquad x = X/S\]

Substituting into the HJB and computing the Ito generator gives the reduced PDE for \(v(t, x, q)\):

\[\begin{aligned} 0 = \partial_t v &+ \tfrac{1}{2}\sigma^2(x^2 v_{xx} + 2x v_x - 1) \\ &+ \sup_{\delta^a} A e^{-k\delta^a}\! \big[v(t, x+1+\delta^a/S, q-1) - v(t, x, q)\big] \\ &+ \sup_{\delta^b} A e^{-k\delta^b}\! \big[v(t, x-1+\delta^b/S, q+1) - v(t, x, q)\big] \end{aligned}\]

Terminal: \(v(T, x, q) = \log(x + q)\).

The state reduces to \((x, q)\): 1 continuous dimension (wealth ratio) plus discrete inventory. The \(\delta/S\) terms in the \(x\)-jumps are O(1e-4) for typical parameter regimes and are neglected in the FD discretization (\(\Delta x = 1\) aligns grid steps with unit fill jumps).

Optimal Spreads (First-Order)

Neglecting \(\delta/S\), the FOC for each side gives:

\[\delta^{a*} = \frac{1}{k} + v(x, q) - v(x+1, q-1)\] \[\delta^{b*} = \frac{1}{k} + v(x, q) - v(x-1, q+1)\]

These are implicit in \(\delta\) when \(\delta/S\) is not neglected. The FD solver uses the first-order explicit form on the grid.

Contrast with CARA

Implemented in solver/src/models/kelly_hjb.rs and engine/src/strategies/kelly_rigorous.rs.

Engine Strategy (Approximation)

The KellyStrategy in the engine implements the CARA-based approximation: it estimates \(\mu\) and \(\sigma\) online, computes \(q^* = \mu/(\gamma\sigma^2)\), and shifts the AS reservation price toward \(q^*\) instead of zero. This is a heuristic — it uses the Kelley target with CARA spread dynamics.

Multi-Asset Kelly (Future Work)

For \(n\) correlated assets with drift vector \(\boldsymbol{\mu}\) and covariance \(\Sigma\), the Kelly-optimal fraction vector is:

\[\mathbf{w}^* = \Sigma^{-1}\boldsymbol{\mu}\]

In the market making context, this would correspond to simultaneous optimal quoting across multiple correlated instruments. The HJB framework would extend to a multi-dimensional inventory \(\mathbf{q} \in \mathbb{Z}^n\) with the risk penalty \(\frac{\gamma}{2}\mathbf{q}^\top\Sigma\mathbf{q}\). Not yet implemented.

Numerical Schemes

Finite Difference Policy Iteration

Discretize the HJB on a grid. The theory of viscosity solutions [@crandall1983viscosity, @user2013users] guarantees convergence of monotone schemes [@barles1991convergence]. Practical implementations follow [@forsyth2007numerical, @kushner2001numerical].

At each time step \(t_n\):

1. Policy evaluation: solve \(A(u^{(k)}) V^{(k)} = V^{n+1} - c(u^{(k)})\) (linear system).
2. Policy improvement [@howard1960dynamic]: \(u_i^{(k+1)} \in \arg\max_u \lbrace [A(u)V^{(k)}]_i + c(u)_i \rbrace\).
3. Iterate until \(\|u^{(k+1)} - u^{(k)}\| < \text{tol}\).

Upwind Discretization

First derivatives:

\[D_x V \approx \begin{cases} \dfrac{V(x + \Delta x) - V(x)}{\Delta x}, & \text{drift}(x, u) > 0 \\[8pt] \dfrac{V(x) - V(x - \Delta x)}{\Delta x}, & \text{drift}(x, u) < 0 \end{cases}\]

Second derivatives (central difference):

\[D_x^2 V \approx \frac{V(x + \Delta x) - 2V(x) + V(x - \Delta x)}{\Delta x^2}\]

The CFL condition [@courant1928partiellen] constrains the time step for explicit schemes. Rannacher smoothing [@rannacher1984finite] is used for Crank-Nicolson startup.

Linear solvers: SOR, Thomas (tridiagonal), LAPACK dgtsv. Schemes: Implicit, Explicit, Crank-Nicolson, Strang ADI.

Implementation: solver::numeric::finite_difference

BSDE Least-Squares Monte Carlo

Backward SDEs provide a probabilistic representation of the HJB solution [@pardoux1990adapted, @el1997backward]. The regression-based scheme follows [@gobet2005empirical].

Forward SDE (Euler-Maruyama)

The forward process is simulated under the optimal control \(u^*\):

\[dX_t = b(t, X_t, u^*_t)\,dt + \sigma(t, X_t)\,dW_t\]

Discretized as

\[X_{n+1} = X_n + b(t_n, X_n, u^*_n)\Delta t + \sigma(t_n, X_n)\Delta W_n.\]

Backward SDE

For the controlled problem with running reward \(f\) and terminal cost \(g\), the value \(Y_t = V(t, X_t)\) satisfies the nonlinear Feynman-Kac relation [@pardoux1990adapted, @el1997backward]:

\[-dY_t = f(t, X_t, u^*_t)\,dt - Z_t\,dW_t, \qquad Y_T = g(X_T),\]

with \(Z_t = \sigma(t, X_t)^\top D_x V(t, X_t)\).

Equivalently

\[V(t, x) = \mathbb{E}\!\left[g(X_T) + \int_t^T f(s, X_s, u^*_s)\,ds \;\middle\|\; X_t = x\right].\]

The backward driver is the running reward \(f\) alone, not the full HJB generator \(f + \mathcal{L}^u V\). The infinitesimal generator \(\mathcal{L}^u V = b\cdot\nabla V + \tfrac12\mathrm{tr}(\sigma\sigma^T D^2 V)\) is already accounted for by simulating the forward SDE under the control; adding it again in the backward step double-counts the drift and diffusion.

LSMC Algorithm

1. Simulate \(M\) forward paths under the optimal control.
2. Set \(Y_N^{(m)} = g(X_N^{(m)})\) (terminal condition).
3. For \(n = N-1, \dots, 0\):
   • Compute the regression target \(\text{Target}^{(m)} = Y_{n+1}^{(m)} + f(t_n, X_n^{(m)}, u^*_n)\Delta t\).
   • Regress \(\text{Target}^{(m)}\) onto basis functions \(\psi_k(X_n^{(m)})\): \[c^n = \arg\min_c \sum_{m=1}^M \left(\text{Target}^{(m)} - \sum_{k=1}^K c_k \psi_k(X_n^{(m)})\right)^2\]
   • Set \(Y_n^{(m)} = \sum_k c_k^n \psi_k(X_n^{(m)})\).

Basis types: Power, Hermite, Chebyshev, Laguerre. Feature count for degree-2: \(K = \frac{(N+1)(N+2)}{2}\) (6 for N=2, 10 for N=3).

See [@han2018solving, @han2017deep] for extensions with neural network approximation (deep BSDE).

Implementation: solver::numeric::bsde

Simulation Framework (Backtest)

Fill Probability

For a limit order at distance \(\delta\) from mid, with base intensity \(\lambda_t\):

\[p_{\text{fill}} = 1 - \exp(-\lambda_t e^{-k\delta} \Delta t)\]

Fill events are independent Bernoulli trials per side per tick.

Implementation: engine::backtest

Terminal Liquidation

At \(t = T\), remaining inventory can be liquidated by crossing the spread:

\[X_T = X_{T^-} + q_T(S_T - \delta_{\text{liq}})\]

The liquidation_half_spread parameter is passed externally (e.g. to run_backtest_with_liquidation). In solver models with TerminalCondition::LiquidationCost, the terminal value function is \(V(T,q) = -|q| \cdot \text{base_spread}\) where \(\text{base_spread} = \frac{1}{\gamma}\ln(1 + \frac{\gamma}{k})\).

Implementation: engine::backtest, solver::models

PnL Decomposition

\[ \begin{aligned} \text{PnL}\_\text{spread} &= \sum_{i=1}^{N_{\text{fills}}} (P_i - S_{t_i}) \\ \text{PnL}\_\text{dir} &= \sum_{t=1}^T q_{t-1}(S_t - S_{t-1}) \end{aligned} \]

PnL_spread uses signed edge: positive when filled passively (bought below mid / sold above mid), negative when filled aggressively (crossed the spread). Terminal liquidation cost is subtracted from PnL_spread.

Implementation: engine::backtest

Sharpe Ratio

The Sharpe ratio [@sharpe1966mutual] is the standard risk-adjusted return metric. For backtest robustness, the Sortino ratio [@sortino1994performance] (penalizing only downside deviation) is also available.

\[\text{SR} = \frac{\bar{R}}{\sigma_R}, \qquad \bar{R} = \frac{1}{K}\sum_{k=1}^K R_k, \qquad \sigma_R = \sqrt{\frac{1}{K}\sum_{k=1}^K (R_k - \bar{R})^2}\]

Hawkes Matcher Parameters

The StochasticMatcher with Hawkes uses the fluid limit for intensity evolution [@hawkes2018hawkes]:

\[\frac{d\lambda_t}{dt} = \beta(\lambda_\infty - \lambda_t) + \alpha\lambda_t\]

Bilateral variant tracks separate \(\lambda_t^b\), \(\lambda_t^a\) with independent self-excitation on fill events.

Implementation: engine::matcher::stochastic

Rigor Reference

This page catalogues every solver item and engine strategy, mapping each to the HJB equation it solves and the method used. Use it to understand what problem is being solved and with what rigour.

Rigor levels

• Rigorous: solves the stated HJB from first principles, no shortcuts. Numerical solutions (FD, BSDE) are rigorous when they discretize the exact HJB without further approximation.
• Approximation: controlled simplification with known error bounds (e.g. asymptotic expansion, small-parameter limit, stationary limit).
• Heuristic: practical rule derived from rigorous theory but without formal derivation or error bounds.

Solver: Analytical Solutions

All analytical solutions solve the CARA-utility AS HJB with \(U(x) = -e^{-\gamma x}\) and exponential fill intensities \(\lambda(\delta) = A e^{-k\delta}\). The reduced HJB and its matrix-exponential solution are derived in Exact Solutions in Stochastic Optimal Control.

Solver: Numerical Models (ControlProblem<N> trait)

All models implement ControlProblem<N> and are solved by either PolicyIterationSolver (FD on a grid) or BsdeSolver (least-squares Monte Carlo). Models solved on a diffusive grid also implement the FD-specific PdeProblem<N> transport contract.

All market-making models below use the CARA utility ansatz \(V = -e^{-\gamma(X + qS + \theta)}\). The HJB they discretize is the reduced PDE for \(\theta(t,q,\ldots)\).

Numerical method details

• FD Policy Iteration: discretizes the HJB on a tensor-product grid. At each time step: (1) policy improvement via \(\sup\) optimization on each grid node, (2) policy evaluation via solving a linear system (SOR, Thomas tridiagonal, or LAPACK dgtsv).
• BSDE: forward Euler-Maruyama simulation, backward least-squares regression onto a polynomial basis (Power, Hermite, Chebyshev, or Laguerre). Regresses continuation values, then applies the control optimiser to extract optimal \(\lambda_\pm\).
• Both solvers are generic over ControlProblem<N>. Correctness for any model reduces to correct implementation of optimize(), running_reward()/generator(), terminal(), and next_step().

Engine: Strategies

Engine strategies consume observations and emit order requests. They are separate from the solver: some consume precomputed solver output, others use closed-form formulas, and others are entirely heuristic.

HJB Derivation References

The full HJB derivation, CARA separation, and matrix-exponential solution are in Exact Solutions in Stochastic Optimal Control.

The Merton and linear-quadratic closed forms are in Exact Solutions in Stochastic Optimal Control. The Kelly criterion and its relationship to CARA are in Optimal Criteria.

Limits of Stochastic Optimal Control Models

This page catalogs the models obtained as limits of a more canonical model: a parameter tends to a boundary value and the richer model reduces to a simpler one with its own exact solution (for example the risk-neutral market maker as \(\gamma \to 0\)). It records only the limit relationships. The dynamics, HJB, and (where one exists) closed form of every model are stated once in Exact Solutions in Stochastic Optimal Control or SOC Models Without Exact Solutions and are not restated here.

The companion page Symmetries in Stochastic Optimal Control catalogs the same limits from the symmetry viewpoint; Rigor Reference records which limit is rigorous versus an approximation; and the solver index maps each model to its implementation.

Limit models

Scope

Only limits that produce a model with a stated exact solution on the exact-solutions page are listed. Asymptotic approximations (the near-maturity \(\tau \to 0\) spread) and research limits without an implemented closed form (mean-field \(m \to \infty\), the no-loss default limit) are catalogued in Rigor Reference and Symmetries in Stochastic Optimal Control respectively.

Symmetries in Stochastic Optimal Control

This appendix catalogs the symmetry structure of the stochastic optimal control (SOC) problems in this repository: which symmetries hold, which are gained or restored in a limit, and which model effect breaks each. It is a companion to Exact Solutions in Stochastic Optimal Control, which derives the solutions; this page only lists the symmetries and points to the derivations. The general framework follows [@pham2009continuous, @fleming2006controlled].

The symmetries labelled "research" below belong to the multi-asset market-making and mean-field problems, which are not yet realized by an implemented model; they are catalogued here so the symmetry structure is not lost, but they do not yet have a closed form on the exact-solutions page.

A symmetry is a transformation of the state (possibly paired with a relabeling of controls) under which the value transforms covariantly, \(V(t, \phi(x)) = \Phi(V(t, x))\) for a fixed \(\Phi\). The three covariance types that occur are multiplicative (CARA), additive (log), and power (CRRA). A symmetry is useful when it reduces the independent coordinates of \(V\), or forces a parity constraint a solver or network must satisfy. The Lie-symmetry treatment follows [@olver1993applications]; the conservation content follows [@noether1918invariante].

Symmetries of the base problems

Symmetries obtained in limits

A limit acts on a symmetry in one of three ways: it restores a symmetry removed by a finite-size or finite-horizon effect, deforms the realization of a symmetry, or exposes the deterministic skeleton by removing dissipation.

Symmetry breaking

Each model effect breaks a specific symmetry; the break is the economic correction a solver must reproduce.

What is and is not conserved

The HJB is a parabolic, dissipative equation; there is no physical energy conserved along optimal trajectories. The objects that play the role of a conserved structure are:

• Comparison principle. If \(V_1 \le V_2\) at the terminal time then \(V_1 \le V_2\) everywhere [@crandall1983viscosity]; monotone schemes preserve it [@barles1991convergence]. This is the correct law to preserve, not an energy.
• Martingale value process. Along the optimal path, the value plus the accumulated running reward is a martingale, conserved in expectation.
• Hamiltonian along characteristics. Pontryagin's Hamiltonian \(H(x, p) = \sup_u \lbrace f(x, u) + b(x, u)^\top p \rbrace\) generates a symplectic flow on \((x, p)\); a continuous symmetry of \(H\) yields a Noether-conserved quantity along characteristics [@noether1918invariante].
• Discrete parity. A reflection symmetry is a \(\mathbb{Z}_2\) charge, preserved if the terminal data and the generator preserve evenness.

Cointegration

This page states the mathematics of cointegration that the repository relies on: the definitions of integrated and cointegrated series, the exact representations, the Gaussian rank inference of the vector error-correction model, the exact theory of a mean-reverting relation with a known cointegrating vector, the filtering theory for time-varying activity, the high-dimensional spectral laws, and the optimal stopping of a mean-reverting spread. Results are stated as results, with the regime in which they hold and the reductions that validate them.

The open questions of the research programme that uses this material, none of which has a proof in the literature as stated, are in the temporal-cointegration plan, and the precise statements that must be proven to close them are on Cointegration: Proof Obligations. This page contains only established mathematics and does not restate them.

Throughout, \(p_t \in \mathbb{R}^N\) is the vector of log-prices, \(N\) is the dimension, \(T\) the number of observations, and \(\gamma = N/T\) the aspect ratio. The Ornstein-Uhlenbeck process used below is the same object defined on Stochastic Processes; the cointegrating residual is an instance of it, with the symbols of that page.

The model assumptions are Gaussian innovations, parameters that are constant within a regime, a first-order finite-state activity chain, and log-prices that are (I(1)). The technical assumptions are finite innovation variance, an irreducible and ergodic chain wherever the filter is used, and finite moments wherever a spectral limit is invoked. Assumptions specific to a single result are stated with that result. The symbol (\Delta) denotes the sampling interval, whereas (\Delta x_t = x_t - x_{t-1}) denotes the first difference; the two are distinguished by the presence of an operand.

Integrated series and the cointegrating representation

A scalar process \(x_t\) is \(I(0)\) if it has a causal linear representation \(x_t = \sum_{j \ge 0} \psi_j \varepsilon_{t-j}\) with \(\varepsilon_t\) independent, mean zero and finite variance, and \(\sum_{j \ge 0} \lvert \psi_j \rvert < \infty\); equivalently, its partial sums converge under the functional limit

\[ T^{-1/2} \sum_{t=1}^{\lfloor Tr \rfloor} x_t \;\Rightarrow\; \sigma W(r), \qquad 0 \le r \le 1, \]

where \(W\) is standard Brownian motion and \(\Rightarrow\) denotes weak convergence in \(D[0,1]\). A process is \(I(1)\) if \(\Delta x_t := x_t - x_{t-1}\) is \(I(0)\) and \(x_t\) itself is not \(I(0)\); then \(T^{-1/2} x_{\lfloor Tr \rfloor} \Rightarrow \sigma W(r)\), and \(\sigma = \psi(1) \sigma_\varepsilon\) is the long-run standard deviation.

The vector \(p_t = (p_{1t}, \dots, p_{Nt})^\top\) is cointegrated of rank \(r\) if every component is \(I(1)\) and there exist \(r\) linearly independent vectors collected in an \(N \times r\) matrix \(\beta\) such that \(\beta^\top p_t\) is \(I(0)\). The columns of \(\beta\) are the cointegrating vectors, \(\mathrm{span}(\beta)\) is the cointegrating space, and \(r\) is the rank. If \(r = 0\) the system has no cointegrating relation; if \(r = N\) every component is \(I(0)\).

The representation that all of the following uses is the exact coordinate identity. Choose any complement \(\beta_\perp \in \mathbb{R}^{N \times (N-r)}\) with \(\beta^\top \beta_\perp = 0\). Then \(\beta\) and \(\beta_\perp\) together span \(\mathbb{R}^N\), and every state decomposes into a stationary part and a trend part,

\[ p_t = \beta (\beta^\top \beta)^{-1} z_t + \beta_\perp (\beta_\perp^\top \beta_\perp)^{-1} f_t, \qquad z_t = \beta^\top p_t, \qquad f_t = \beta_\perp^\top p_t . \]

Substituting \(p_t\) into \(\beta^\top p_t\) gives \(\beta^\top p_t = z_t\) exactly, because \(\beta^\top \beta_\perp = 0\) kills the trend term. The identity is therefore exact finite-dimensional linear algebra and involves no approximation: choosing the dynamics of the stationary coordinates \(z_t\) and the trend coordinates \(f_t\) determines \(p_t\) with exact cointegration by construction.

The converse direction is the Granger representation theorem. If \(p_t\) is \(I(1)\) with cointegration rank \(r\), then it has a vector error-correction representation

\[ \Delta p_t = \alpha \beta^\top p_{t-1} + \sum_{i=1}^{k-1} \Gamma_i \Delta p_{t-i} + \varepsilon_t, \qquad \Pi := \alpha \beta^\top, \]

with \(\alpha, \beta \in \mathbb{R}^{N \times r}\) of full column rank and \(\alpha_\perp^\top \Gamma \beta_\perp\) invertible, where \(\Gamma = I_N - \sum_{i=1}^{k-1} \Gamma_i\) and \(\alpha_\perp, \beta_\perp\) are complements (Engle and Granger, 1987). The matrix \(\Pi\) has rank exactly \(r\): it is the only channel through which the levels \(p_{t-1}\) enter, and its rank is the rank of the cointegrating space.

Gaussian inference and the rank test

For inference assume the innovations are Gaussian and the dynamics are contiguous over the sample, so no parameter changes within \([0,T]\). The model is

\[ \Delta p_t = \Pi p_{t-1} + \sum_{i=1}^{k-1} \Gamma_i \Delta p_{t-i} + \varepsilon_t, \qquad \varepsilon_t \sim \mathcal{N}(0, \Omega), \]

with fixed \(N\) and \(T \to \infty\). The Gaussian likelihood is exact, and maximizing it over \(\Pi = \alpha \beta^\top\) of rank \(r\) is a reduced-rank regression (Johansen, 1988, 1991). Reduce \(\Delta p_t\) and \(p_{t-1}\) on the lagged differences \((1, \Delta p_{t-1}, \dots, \Delta p_{t-k+1})\) by least squares to obtain residuals \(R_{0t}\) and \(R_{1t}\), and set

\[ S_{ij} = \frac{1}{T} \sum_{t=1}^{T} R_{it} R_{jt}^\top, \qquad i, j \in \lbrace 0, 1 \rbrace . \]

The profile likelihood depends on \(\beta\) only through the determinant \(\lvert \beta^\top S_{11} \beta \rvert / \lvert \beta^\top \beta \rvert\), so the maximizer solves the generalized eigenvalue problem

\[ \det\!\left( \lambda S_{11} - S_{10} S_{00}^{-1} S_{01} \right) = 0 . \]

Order the eigenvalues \(\hat\lambda_1 \ge \dots \ge \hat\lambda_N\). The maximum likelihood estimate of \(\beta\) is the matrix of the \(r\) eigenvectors belonging to the \(r\) largest eigenvalues, normalized by \(\hat\beta^\top S_{11} \hat\beta = I\), and \(\hat\alpha = S_{01}\hat\beta\). This estimator is an explicit algebraic functional of the data; it is closed form in the sense of a finite eigenproblem, not an iterative search.

The likelihood-ratio statistic for the null hypothesis of rank at most \(r\) is

\[ \mathrm{LR}(r) = -T \sum_{i=r+1}^{N} \ln\!\left( 1 - \hat\lambda_i \right) . \]

Under the null its limit is not the chi-square law, because the regressor \(p_{t-1}\) is \(I(1)\) and the estimation error of \(\beta\) does not vanish fast enough. It is the trace of a functional of an \((N-r)\)-dimensional standard Brownian motion,

\[ \mathrm{LR}(r) \;\Rightarrow\; \operatorname{tr}\!\lbrace \left( \int_0^1 W \, dW^\top \right) \left( \int_0^1 W W^\top \, du \right)^{-1} \left( \int_0^1 dW \, W^\top \right) \rbrace, \]

whose percentiles are tabulated (Johansen, 1991). Two features matter for the present problem and are used repeatedly below: the distribution depends on \(N\), and it is derived under fixed \(N\) with \(T \to \infty\), so it does not apply when \(N\) and \(T\) grow together.

Exact theory of a single mean-reverting relation

If \(\beta\) is known, the relation is observed directly and the problem becomes univariate. Model the residual \(j\) as an Ornstein-Uhlenbeck process with speed \(\theta_j > 0\), long-run mean \(\mu_j\), and volatility \(\sigma_j > 0\),

\[ d z_{j,t} = \theta_j (\mu_j - z_{j,t})\, dt + \sigma_j\, dW_{j,t}, \]

the same process as on Stochastic Processes. Its stationary law and the relaxation time constant are

\[ z_j(\infty) \sim \mathcal{N}\!\left( \mu_j, \frac{\sigma_j^2}{2\theta_j} \right), \qquad h_j = \frac{\ln 2}{\theta_j}, \]

where \(h_j\) is the half-life: the time for the expected deviation from \(\mu_j\) to halve. The stationary standard deviation is \(\varsigma_j = \sigma_j / \sqrt{2\theta_j}\), and it is the natural scale of a deviation and of any trading threshold.

The exact discrete law at sampling interval \(\Delta\) follows from the integrating factor. Solving the linear SDE over one step,

\[ z_{j,t+\Delta} = \mu_j + (z_{j,t} - \mu_j) e^{-\theta_j \Delta}

• \sigma_j \int_0^{\Delta} e^{-\theta_j (\Delta - s)} \, dW_{j,s}, \]

and the Ito isometry evaluates the variance of the stochastic integral as \(\sigma_j^2 (1 - e^{-2\theta_j \Delta}) / (2\theta_j)\). Hence the sampled process is exactly an AR(1),

\[ z_{j,t+\Delta} = \mu_j (1 - \phi_j) + \phi_j z_{j,t} + \epsilon_{j,t}, \qquad \phi_j = e^{-\theta_j \Delta}, \qquad \operatorname{Var}(\epsilon_{j,t}) = \frac{\sigma_j^2}{2\theta_j} \left( 1 - \phi_j^2 \right), \]

with \(\epsilon_{j,t}\) independent and Gaussian. Everything here is exact, not an Euler approximation: the Euler scheme would replace \(\phi_j\) by \(1 - \theta_j \Delta\) and is only first order.

Two consequences of the exact law are used as design criteria. First, the resolvability condition: as \(\theta_j \Delta \to 0\) the autoregressive root satisfies \(\phi_j \to 1\), so the sampled residual converges to a random walk and cannot be distinguished from one by any test of fixed size. A non-vanishing \(\theta_j \Delta\) is therefore necessary for the relation to be identifiable at the sampling interval. Second, the observability condition: the expected deviation from the mean after a lifespan \(D\) is \(e^{-\theta_j D} = 2^{-D/h_j}\) of its initial value, so observing the relaxation to within a fraction \(\rho\) of the initial deviation requires \(D / h_j \ge \log_2(1/\rho)\). For example, \(\rho = 0.1\) requires \(D/h_j \gtrsim 3.3\), and \(\rho = 0.03\) requires \(D/h_j \gtrsim 5\).

The Gaussian estimator of \(\theta_j\) is a deterministic function of the AR(1) estimate. Conditional on the first observation, the least-squares estimator of \(\phi_j\) is the demeaned autocorrelation

\[ \hat\phi_j = \frac{\sum_{t} (z_{j,t} - \bar z_j)(z_{j,t+\Delta} - \bar z_j)}{\sum_{t} (z_{j,t} - \bar z_j)^2}, \qquad \hat\theta_j = -\frac{\ln \hat\phi_j}{\Delta}, \]

and \(\hat\phi_j\) is downward biased in finite samples, of order \(1/T\) (Kendall, 1954; Stambaugh, 1999). The bias inflates the estimated speed and understates the half-life; it is a property of the estimator, reproduced by the model, and not a numerical defect.

Time-varying activity

Cointegration that holds only on part of the sample is modelled by switching the adjustment channel, not the cointegrating space. Let \(s_{j,t} \in \lbrace 0, 1 \rbrace\) be the activity state of relation \(j\), with

\[ d z_{j,t} = s_{j,t}\, \theta_j (\mu_j - z_{j,t})\, dt + \sigma_j\, dW_{j,t}. \]

When \(s_{j,t} = 1\) the residual mean-reverts; when \(s_{j,t} = 0\) it is a driftless random walk. This is exactly the \(\alpha_j \to 0\) limit of the error-correction representation on the inactive intervals, so the local rank of the cointegrating space equals the number of active relations. The switching capability defines the class; estimating it uses the filter below.

Take the state to be a first-order Markov chain on a finite set with transition matrix \(P\), \(P_{ij} = \Pr(s_t = j \mid s_{t-1} = i)\), and let \(\eta_{t,j} = f(y_t \mid s_t = j)\) be the density of the observation under state \(j\). The Hamilton (1989) filter computes the exact state distribution by the forward recursion

\[ \xi_{t \mid t} = \frac{\xi_{t \mid t-1} \odot \eta_t}{\mathbf{1}^\top (\xi_{t \mid t-1} \odot \eta_t)}, \qquad \xi_{t+1 \mid t} = P^\top \xi_{t \mid t}, \]

where \(\odot\) is elementwise multiplication and \(\mathbf{1}\) is the vector of ones. The log-likelihood is exact and equals the sum of the normalizing constants,

\[ \ell(\psi) = \sum_{t=1}^{T} \ln \mathbf{1}^\top \!\left( \xi_{t \mid t-1} \odot \eta_t \right), \]

and the exact regime posterior is the backward-smoothed distribution \(\xi_{t \mid T}\). The objective and the posterior are exact; the maximizer over the parameters \(\psi\) is not available in closed form, so estimation is a finite-dimensional numerical optimization of an exactly computable criterion. This is the reason the switching model is classified below as exact but not closed form.

If the regime path is known, for instance from a deterministic schedule, the model is conditionally Gaussian with a known regressor partition, and the exact conditional maximum likelihood estimator is block least squares: ordinary least squares applied separately to each regime. No filter is needed, and this is the reference estimator the switching filter must reproduce.

For a single change at an unknown time, the break fraction estimator converges at rate \(T\) to an argmax of a two-sided Brownian motion with a parabolic drift (the Chernoff distribution), and the supremum of the likelihood-ratio statistic over the change date converges to a known functional of Brownian motion (Andrews, 1993; Bai, 1997). For a cointegrated system, the tests that separate a change in the adjustment \(\alpha\), in the cointegrating space \(\beta\), in the rank, and in the short-run dynamics are exactly constructed with tabulated limits (Hansen, 2003). These are the tools that distinguish the three meanings of temporary cointegration.

High-dimensional spectral theory

When \(N\) is of order \(10^3\) and \(\gamma = N/T\) is of order one, the eigenvalues of a sample covariance are governed by the Wishart ensemble. If \(X\) has \(N \times T\) independent entries and \(S = T^{-1} X X^\top\), the eigenvalues of \(S\) follow the Marchenko-Pastur law

\[ \rho_{\mathrm{MP}}(d\lambda) = \frac{\sqrt{(\lambda_+ - \lambda)(\lambda - \lambda_-)}}{2\pi \gamma \lambda} \quad \text{on } [\lambda_-, \lambda_+], \qquad \lambda_\pm = (1 \pm \sqrt{\gamma})^2 , \]

concentrated on \([\lambda_-, \lambda_+]\), with an atom at the origin when \(\gamma > 1\) (Marchenko and Pastur, 1967). The largest eigenvalue fluctuates around the upper edge \(\lambda_+\) on the scale \(T^{-2/3}\) and follows the Tracy-Widom law: \(\beta = 1\) for real entries and \(\beta = 2\) for complex entries (Tracy and Widom, 1994, 1996; Johnstone, 2001). The joint eigenvalue density is the Laguerre ensemble density, which is exact at finite \(N\), so a test built on the eigenvalues can be calibrated in finite samples and not only asymptotically.

A cointegrating relation generates a separated eigenvalue, which is the spiked covariance model. For a population covariance \(I + \ell\, v v^\top\) and aspect ratio \(\gamma\), the Baik-Ben Arous-Peche transition is sharp at

\[ \ell_c = \sqrt{\gamma} : \]

for \(\ell < \ell_c\) the spiked eigenvalue does not separate from the bulk and the sample eigenvector is asymptotically orthogonal to \(v\), whereas for \(\ell > \ell_c\) the outlier sits at \((1+\ell)(1+\gamma/\ell)\) and the squared overlap of the sample and population eigenvectors converges to \((1 - \gamma/\ell^2)/(1 + \gamma/\ell)\) (Baik, Ben Arous and Peche, 2005; Benaych-Georges and Nadakuditi, 2011). The threshold is the detection floor of the rank problem: below it no consistent eigenvector is available, however large \(T\) is.

The sample covariance is ill-conditioned at these aspect ratios, so it is replaced by a shrinkage estimator before any eigen-decomposition. The optimal linear shrinkage toward a scaled identity,

\[ \hat\Sigma = (1 - \delta) S + \delta\, \frac{\operatorname{tr}(S)}{N} I, \]

has a closed-form, Marchenko-Pastur-consistent intensity \(\delta\) (Ledoit and Wolf, 2004).

A structural caveat conditions the transfer of all of the above to cointegration. The spectral results assume independent or weakly mixing entries with a growing sample. The matrices of the rank problem are built from the levels \(p_{t-1}\), which are \(I(1)\) and therefore non-mixing: they are functionals of matrix Brownian motion, not a Wishart of independent rows. Replacing the Wishart thresholds by Marchenko-Pastur or Baik-Ben Arous-Peche thresholds in the rank problem is therefore not justified by the results as stated. What the correct threshold is in that case is open and is stated as such in the plan.

Optimal stopping of a mean-reverting spread

Once \(\beta\) is fixed, trading a detected relation is an optimal stopping problem for the spread. Centre the residual, \(x = z - \mu\), so that with a constant cost \(c\) per round trip the payoff of an entry at level \(x\) and exit at the mean is the deviation measured net of cost. The state follows

\[ d x_t = -\theta x_t\, dt + \sigma\, dW_t, \]

with infinitesimal generator \(\mathcal{L} = \tfrac{1}{2} \sigma^2 \partial_{xx} - \theta x \partial_x\). On a continuation region the value function satisfies the homogeneous equation \(\mathcal{L} u = 0\), which integrates directly. Writing \(\mathcal{L} u = 0\) as \(u'' / u' = 2\theta x / \sigma^2\) and integrating once gives \(u'(x) = c_1 \exp(\theta x^2 / \sigma^2)\), hence

\[ u(x) = c_0 + c_1 \int_0^{x} \exp\!\left( \frac{\theta s^2}{\sigma^2} \right) ds . \]

The integral is the Dawson function, equivalently expressible through the parabolic cylinder function or the confluent hypergeometric \(M(1/2, 3/2, \theta x^2/\sigma^2)\). The optimal thresholds therefore appear as the unknowns of a system of value-matching and smooth-pasting equations with right-hand sides in this family, and the thresholds are characterized by transcendental equations rather than an elementary closed form. The full transition density of the free process is the exact Gaussian kernel

\[ p(\Delta, x, y) = \frac{1}{\sqrt{2\pi v(\Delta)}} \exp\!\left( -\frac{(y - m(\Delta, x))^2}{2 v(\Delta)} \right), \qquad m = x e^{-\theta \Delta}, \quad v = \frac{\sigma^2}{2\theta} \left( 1 - e^{-2\theta \Delta} \right), \]

which is the same discrete law used in the previous section. With killing at two barriers the transition density has a spectral representation in the Ornstein-Uhlenbeck eigenfunctions, which are the Hermite functions with eigenvalues \(-n\theta\); this representation is what would make a tick-quantized observation model exactly solvable, and is stated here only as the structure of the known result.

Assumptions and rigor

The results above are classified by the project rigor levels. The model and technical assumptions are stated at the top of the page.

Validation targets

The results on this page define the checks that an implementation must pass, and the reductions that show each check is the right one. No implementation exists yet; these are targets, not passed tests.

• With \(r = 0\), the system is \(I(1)\) with no relation, and the detector must reject at its nominal size. This is the null against which the spectral thresholds, including the Ledoit-Wolf regularized ones, are calibrated.
• With \(N = 2\) and \(r = 1\), the reduced-rank regression has a hand-computable generalized eigenvalue problem, and the estimator can be checked in closed form.
• Setting \(s_{j,t} \equiv 1\) reduces the switching model to the contiguous model, and the filter likelihood reduces to the Gaussian likelihood; the two must agree.
• As \(\theta_j \Delta \to 0\) the exact AR(1) law converges to a random walk, which is the resolvability failure and must appear as a loss of power.
• As \(\sigma_j \to 0\) the residual follows the deterministic path \(\mu_j + (z_0 - \mu_j)e^{-\theta_j t}\), which fixes the half-life reduction.
• The stopping thresholds reduce to the known mean-reversion band as \(c \to 0\) and reproduce the exact Gaussian transition kernel above; the two-barrier spectral density must match the free kernel when the barriers are removed.

Cointegration: Proof Obligations

This page states, for each result the cointegration programme needs and that is not established, exactly what must be proven and why. It is the proof-target companion to Cointegration, which holds only established mathematics, and to the temporal-cointegration plan, which holds the programme and its implementation stages.

Nothing on this page is a result. Each item is an obligation: a precisely stated theorem to be proven, or refuted by an explicit counterexample, or shown to be ill-posed and replaced. An obligation is discharged in exactly one of three ways and is then migrated out of this page:

• a theorem with a complete proof, which moves to Cointegration as established mathematics;
• a counterexample with a proof, recorded as a negative result within the obligation, together with the corrected statement that replaces it;
• a well-posedness finding, showing the obligation as stated is not a well-defined mathematical question, with the reformulation.

An item is not discharged by a simulation consistent with it, by a proof sketch, or by an appeal to the stationary analogue. Where a line of the proof is routine and where the research content lies is stated for each item, so that effort is not spent on the routine part.

The setting and the symbols below are those of Cointegration and Stochastic Processes. All obligations are stated for \(N\) of order \(10^3\), aspect ratio \(\gamma = N/T\) of order one, and the observability ordering \(\Delta \ll h \ll D \le T\).

Setting and definitions

Fix a filtered probability space \((\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0}, \mathbb{P})\). Time is continuous unless a discrete index is stated; observations are recorded at integer multiples of the sampling interval \(\Delta\).

Let \(\beta \in \mathbb{R}^{N \times r}\) be fixed of full column rank, and let \(\beta_\perp\) be a complement with \(\beta^\top \beta_\perp = 0\). The residual and trend coordinates are \(z_t = \beta^\top p_t\) and \(f_t = \beta_\perp^\top p_t\), and the reconstruction identity is \(p_t = \beta(\beta^\top\beta)^{-1} z_t + \beta_\perp(\beta_\perp^\top\beta_\perp)^{-1} f_t\).

The regime process \(s_t \in \lbrace 0, 1 \rbrace^r\) is a semi-Markov process: it changes value at stopping times \(\tau_0 < \tau_1 < \dots\), the sojourn in state \(u\) has law \(F_u\) supported on \([D_{\min}, \infty)\) with finite mean \(m_u\) and \(F_u(D_{\min}) = 0\), and the embedded chain has transition matrix \(R\) with \(R_{uu} = 0\). Write \(\pi = \lim_{t\to\infty} \mathbb{E}[s_t]\) for the stationary activation vector, assumed to exist.

The residual follows the switching Ornstein-Uhlenbeck dynamics

\[ d z_t = \operatorname{diag}(s_t)\, \Theta (\mu - z_t)\, dt + \Sigma_z\, dW_t, \qquad \Theta = \theta I_r, \quad \Sigma_z \Sigma_z^\top \succ 0, \]

so that when \(s_{j,t} = 1\) the relation \(j\) mean-reverts with half-life \(h_j = \ln 2 / \theta_j\), and when \(s_{j,t} = 0\) it is a driftless random walk. The trend is a driftless random walk, \(f_t = f_0 + \Sigma_f B_t\), with \(W, B\) independent standard Brownian motions.

The observation layer maps the latent state to the data, \(y_t = \mathcal{Q}\big(p_t + \eta_t\big)\), where \(\eta_t\) is microstructure noise and \(\mathcal{Q}\) is rounding to a tick grid. During the theoretical programme the layer is switched off until an obligation concerns it.

Two definitions are needed and neither is standard. Call \(p\) interval-cointegrated with activity \(s\) and vectors \(\beta\) if, on every maximal interval on which \(s\) is constant, the restriction of \(\beta^\top p\) to that interval is stationary. The local rank at time \(t\) is \(r_t = \sum_{j=1}^{r} s_{j,t}\), the number of active relations; the ergodic rank is the rank of the ergodically averaged error-correction matrix \(\bar\Pi = \lim_{T\to\infty} T^{-1}\sum_{t} \alpha \operatorname{diag}(s_t)\beta^\top\).

P1. Well-posedness and the switching representation

Statement. Prove the following, which jointly make interval-cointegration a defensible process class.

1. For every regime path and every initial condition, the switching residual SDE has a pathwise unique strong solution, and the reconstruction \(p_t = \beta(\beta^\top\beta)^{-1} z_t + \beta_\perp(\beta_\perp^\top\beta_\perp)^{-1} f_t\) is continuous at every switching instant.
2. On each maximal active interval, the law of \(z\) converges to the Ornstein-Uhlenbeck stationary law as the interval length grows, uniformly in the initial condition, at the rate governed by \(\exp(-\theta D_{\min})\).
3. The open direction. Every \(I(1)\) process that is interval-cointegrated with activity \(s\) admits a representation of the form above, with \(\alpha_t = \alpha \operatorname{diag}(s_t)\), a short-run part, and \(\alpha_\perp^\top \Gamma \beta_\perp\) invertible on each active interval.

Why it is required. Without part 3, "interval-cointegrated" is a definition with no representation theorem behind it, and no estimator has a stated target. Part 3 is the converse of the construction, in the same position that the Granger representation theorem occupies for the contiguous case.

What would refute it. A process that is interval-cointegrated yet admits no error-correction representation on some active interval, or a process with well-defined local behaviour but no well-defined path at a switching instant.

Strategy and division of labour. Parts 1 and 2 are routine: Lipschitz drift and the explicit integrating factor give both, and part 2 is the standard Ornstein-Uhlenbeck convergence with the switching entering only through the interval length. Part 3 is the research content; it requires a Beveridge-Nelson-type decomposition adapted to a random partition of the time axis.

Depends on. Nothing.

P2. Identifiability and consistency of rank under switching

Statement. Let the regime chain be ergodic and irreducible on its state space, so \(\pi_j > 0\) for every \(j\), and let the Johansen rank estimator \(\hat r_T\) be the number of eigenvalues of \(S_{11}^{-1} S_{10} S_{00}^{-1} S_{01}\) above a threshold \(c_T \to 0\). Prove:

1. \(\operatorname{rank}(\bar\Pi) = r\); equivalently, the \(j\)-th adjustment column is attenuated by the factor \(\pi_j\), so rank is identified if and only if \(\pi_j > 0\) for all \(j\) and \(\alpha_j \neq 0\).
2. \(\hat r_T \to r\) in probability.
3. The open direction. The limiting law of \(T(\hat r_T - r)\) is a mixture of the contiguous Johansen functional, with mixing weights determined by the stationary distribution of the regime chain and the sojourn laws \(F_u\). Give the mixture explicitly.
4. The boundary. If some \(F_u\) is heavy-tailed so that active intervals have vanishing density, consistency fails. Give the precise tail condition that separates consistency from failure.

Why it is required. It decides whether the tabulated contiguous critical values may be used at all, and it is the first place where switching changes the answer rather than merely the sampling variability. Part 1 also gives the correct interpretation of an estimated loading: it is the true loading scaled by the activation probability, so a small \(\hat\alpha_j\) does not distinguish a weak relation from a rarely active one.

What would refute it. A regime law with \(\pi_j > 0\) for which \(\hat r_T\) is inconsistent, which would show that positive activation probability is not sufficient for identifiability.

Strategy and division of labour. Part 1 is routine: an ergodic theorem gives \(T^{-1}\sum_t \operatorname{diag}(s_t) \to \operatorname{diag}(\pi)\) and rank is preserved by invertible diagonal scaling. Parts 2 and 4 are the research content, requiring a functional central limit theorem for a triangular array whose regressor distribution is randomly time-varying. Part 3 requires characterizing the limit of an averaged quadratic form in matrix Brownian motion under a random time change.

Depends on. P1.

P3. Exact null law of the joint time-and-subspace detector

Statement. Define the joint detector

\[ \Lambda_T = \sup_{\tau \in [\epsilon T, (1-\epsilon) T]} \; \sup_{\beta \in G_{r,N}} \; \mathrm{LR}\big(\tau, \beta\big), \]

where \(G_{r,N}\) is the Grassmannian of \(r\)-dimensional subspaces of \(\mathbb{R}^N\) and \(\mathrm{LR}(\tau, \beta)\) is the likelihood-ratio statistic for a rank-\(r\) relation with vectors \(\beta\) active on a neighbourhood of \(\tau\). Under the null that \(p\) is \(I(1)\) with rank \(0\), derive the limiting law of \(\Lambda_T\) after the appropriate centering and scaling.

The conjectured limit object is a functional of a two-parameter process indexed by (time, subspace), the subspace parameter living on \(G_{r,N}\), obtained as the supremum of a quadratic form in matrix Brownian motion over the product of the time interval and the Grassmannian. Establish whether the two suprema commute in the limit.

Why it is required. This is the exact solution of the detection half. A detector that searches over both change times and subsets has a null distribution that is neither the Johansen law nor the single-break law; without it, the threshold is a simulated table for one parameter setting and cannot be transported.

What would refute it. A proof that \(\Lambda_T\) has no limiting law under any consistent scaling, together with the scaling that restores one, or a proof that the limit is degenerate.

Strategy and division of labour. The one-dimensional suprema each have known limits; the obligation is the interaction. The starting point is a functional central limit theorem for the score process indexed by \(\tau\) uniformly over \(G_{r,N}\), followed by a continuous-mapping argument on the product space, with care that \(G_{r,N}\) is non-compact in \(N\).

Depends on. P2, P7.

P4. Random matrix theory for the rank problem with \(I(1)\) levels

Statement. Let \(S_{ij} = T^{-1} \sum_t R_{it} R_{jt}^\top\) be the matrices of the reduced-rank regression, built from the \(I(1)\) levels. Prove, or disprove, a functional limit theorem for the spectrum of

\[ M_T = S_{11}^{-1/2} S_{10} S_{00}^{-1} S_{01} S_{11}^{-1/2} \]

as \(N, T \to \infty\) with \(N/T \to \gamma \in (0, \infty)\). Determine whether the largest eigenvalue of \(M_T\) under the null of rank zero separates from a bulk, and if so give the bulk law, the edge, and the fluctuation class. Define the cointegration threshold \(\ell_c^{\mathrm{coint}}\) as the sharp spike strength at which a rank-one alternative becomes detectable, and establish whether \(\ell_c^{\mathrm{coint}} = \sqrt{\gamma}\), or \(\ell_c^{\mathrm{coint}} > \sqrt{\gamma}\), or \(\ell_c^{\mathrm{coint}} < \sqrt{\gamma}\).

Why it is required. Marchenko-Pastur and Baik-Ben Arous-Peche require independent or weakly mixing entries with a growing sample. The matrices here are built from the levels, which are \(I(1)\) and non-mixing, so those results do not apply as stated. Until this is resolved, no threshold for the 1000-dimensional rank problem is justified, and the comparison with the stationary threshold is the quantitative measure of how much harder the non-stationary problem is.

What would refute it. An argument that no eigenvalue separates from any bulk in this scaling, which would mean that rank is not recoverable at any signal strength when \(\gamma\) is of order one. That is itself a decisive negative result.

Strategy and division of labour. The obstruction is that the residual matrices involve \(T^{-1}\sum_t p_{t-1} \varepsilon_t^\top\), a functional of matrix Brownian motion, so the standard Wishart representation fails. The candidate route is to express \(M_T\) as a functional of an \(N \times T\) array with a non-stationary, non-mixing column structure and to seek a limit theorem for such functionals; this is the least developed of the obligations.

Depends on. P7.

P5. Sparse-spike thresholds and delocalization

Statement. Let \(\beta\) have support of size \(s \ll N\), so that the alternative is a sparse spike. Prove an upper and a lower bound for the minimal spike strength at which the support of \(\beta\) is recoverable, with the aspect ratio \(\gamma\) fixed. Determine whether the threshold scales as \(\sqrt{s \log N / T}\) or as \(\sqrt{N/T}\), and characterize the region of \((s, \gamma)\) in which the estimated eigenvector is delocalized and therefore does not identify the support.

Why it is required. Sparse \(\beta\) is the only design that answers the question of which dimensions are cointegrated, and the delocalized threshold is the wrong one to use for it. This obligation also fixes what "which dimensions" can mean: if the eigenvector is delocalized, the support is not recoverable even when the relation is detected.

What would refute it. A proof that sparse structure gives no improvement over the dense threshold in this non-stationary setting.

Depends on. P4.

P6. The detection floor as a minimax lower bound

Statement. Let \(\mathcal{P}(\ell)\) be the class of laws under which the system carries a rank-one relation of spike strength \(\ell\), and \(\mathcal{P}_0\) the rank-zero null. Prove that there is a constant \(c > 0\) such that for \(\ell < \ell_c\),

\[ \inf_{\hat r} \; \sup_{P \in \mathcal{P}(\ell) \cup \mathcal{P}_0} \; P\big(\hat r \text{ incorrect}\big) \;\ge\; c, \]

the infimum over all tests, and construct an estimator whose risk vanishes for \(\ell > \ell_c\). The \(\ell_c\) is that of P4, or of P5 in the sparse case.

Why it is required. It converts "our detector failed on this instance" into "no detector can succeed on this class", which is the only statement that closes the question of how long and how strongly a relation must be active. It is the mathematical content behind the claim that there is a hard floor.

What would refute it. An estimator that beats the conjectured floor, which would invalidate P4 or P5.

Strategy and division of labour. A Le Cam or Fano argument over a packing of the parameter space by alternative cointegrating vectors, with the divergence between the resulting laws controlled by the spike strength and the effective sample size per active interval, which is where the switching enters.

Depends on. P4, P7.

P7. Sharp timescale conditions

Statement. Turn the ordering \(\Delta \ll h \ll D\) into necessary and sufficient conditions. Let \(\alpha \in (0,1)\) be the test size and \(1-\beta\) the prescribed power. Prove:

1. Resolvability. Distinguishing \(\theta = 0\) from \(\theta > 0\) at a fixed span requires \(\theta \Delta\) bounded away from zero, and the Fisher information of the discretely sampled experiment is \(\Theta\big((\theta \Delta)^2\big)\). Give the exact constant.
2. Observability. For a fixed power, \(D / h \ge c(\alpha, \beta)\); give \(c(\alpha, \beta)\) and show it is necessary, not merely sufficient.
3. Estimability. For a fixed power with detection and confirmation on disjoint sub-intervals, \(D / W \ge c'(\alpha, \beta)\), where \(W\) is the detection window.

Why it is required. These constants are what make the timescale ordering operational. Without them, every parameter choice in a generator or a study is undefended, and a failure cannot be attributed to the theory rather than to the implementation.

What would refute it. A demonstration that the conditions are sufficient but not necessary, with the correct necessary condition replacing each.

Strategy and division of labour. Local asymptotic normality for the active-interval experiment gives the Fisher information of part 1; contiguity across a switching boundary gives the loss from part 3. Part 2's constant is a power calculation under the exact Gaussian law of Cointegration and is the most tractable of the three.

Depends on. Nothing.

P8. Exact solvability class of the switching stopping problem

Statement. Consider the optimal stopping problem for the spread with regime switching, whose value on the continuation region solves the coupled system \(\mathcal{L}_k u_k = 0\) across regimes \(k\), with regime-dependent speed \(\theta_k\), volatility \(\sigma_k\), cost \(c_k\), and transition rates \(q_{kl}\). Characterize the set of \((\theta_k, \sigma_k, c_k, q_{kl})\) for which the coupled system admits a closed-form solution in the parabolic cylinder or confluent hypergeometric family, and prove that outside this set no such closure exists.

Why it is required. It defines which detected relations can be traded with an exact policy and where the boundary of the exact-solvable envelope lies. It is the same structural question asked of market making in multi-asset-market-making.md, so a shared answer is a structural result rather than a model-specific one.

What would refute it. A proof that no nontrivial switching structure closes, which fixes the envelope boundary as a negative result.

Depends on. P7.

P9. Exact likelihood under tick quantization and asynchronicity

Statement. With the observation layer active, the exact likelihood is \(\prod_t \Pr\big(y_t \in \text{bin} \mid z_{t-1}\big)\), a product of transition probabilities of the latent process between quantization bins. Prove whether these probabilities admit an exact (spectral, Hermite-function) representation in the multivariate switching case, as they do for a univariate Ornstein-Uhlenbeck with two barriers, or prove that they do not and give the correct approximation class and its error.

Why it is required. Every threshold derived under continuous synchronous observation is an upper bound on real performance, and this obligation quantifies the gap. It is also the only obligation in which the HFT observation layer, and not the latent process, is the object of study.

What would refute it. A proof of intractability, with the approximating class that replaces the exact likelihood.

Depends on. P1.

P10. Bias of the half-life estimator under switching

Statement. Let \(\hat\theta\) be the estimator obtained by applying the AR(1) estimator of Cointegration to the active sub-intervals only. Derive \(\mathbb{E}[\hat\theta] - \theta\) to leading order in \(T\), as a function of the activation vector \(\pi\), the sojourn laws \(F_u\), and the sampling ratio \(\theta \Delta\), and construct a bias-corrected estimator with vanishing bias.

Why it is required. The half-life determines tradeability. The contiguous bias is of order \(1/T\); if the switching bias is of larger order, half-life estimation is unreliable at the horizons the programme targets, and the tradeability filter would be systematically wrong.

What would refute it. A demonstration that the bias is not \(O(1/T)\) but of larger order, which would make the estimator unreliable and require a different estimator rather than a correction.

Depends on. P2.

Dependencies and discharge order

graph TD
    P7[P7 Timescale conditions] --> P3[P3 Joint null law]
    P7 --> P4[P4 RMT for I(1) levels]
    P7 --> P6[P6 Minimax floor]
    P1[P1 Representation] --> P2[P2 Rank identifiability]
    P1 --> P9[P9 Quantized likelihood]
    P2 --> P3
    P2 --> P10[P10 Half-life bias]
    P4 --> P5[P5 Sparse thresholds]
    P4 --> P6
    P7 --> P8[P8 Switching stopping]

The rule is that an obligation is started only after the obligations it depends on are discharged. The obligations that can start immediately are P1 and P7, which is why they are the first two in any schedule. P2 and P4 are the load bearing results: P2 decides whether rank is identifiable at all under switching, and P4 decides whether it is recoverable in high dimension.

What a discharged obligation looks like

For illustration of the required standard, a discharged obligation would read as a statement of the form: under assumptions A1 to Ak, for all parameter values in a specified set, a specified estimator satisfies a specified limit, with the limit object named and its distribution or value given, together with the reduction to the corresponding contiguous result as the sojourn laws degenerate to a point mass at \(D_{\min} \to \infty\). The reduction is the check that the statement is not vacuously general, and it is the analogue of the specializations listed in the validation targets of Cointegration.

Bibliography

References grouped by topic. Each citation below is linked from the relevant Mathematical Reference page. Items marked with [*] are the primary references for that topic.

Stochastic Calculus and SDEs

• Itô, K. (1944). Stochastic integral. Proc. Imperial Academy, 20(8), 519--524. [@ito1944stochastic]
• Itô, K. (1951). On a formula concerning stochastic differentials. Nagoya Math. J., 3, 55--65. [@ito1951formula] [*]
• Itô, K. (1951). On stochastic differential equations. AMS Memoirs, 4. [@ito1951stochastic]
• Stratonovich, R. L. (1966). A new representation for stochastic integrals and equations. SIAM J. Control, 4(2), 362--371. [@stratonovich1966new]
• Harrison, J. M. & Pliska, S. R. (1981). Martingales and stochastic integrals in the theory of continuous trading. Stochastic Processes and their Applications, 11(3), 215--260. [@harrison1981martingales]
• Kac, M. (1949). On distributions of certain Wiener functionals. Trans. AMS, 65(1), 1--13. [@kac1949distributions]
• Feynman, R. P. (1948). Space-time approach to non-relativistic quantum mechanics. Rev. Mod. Phys., 20(2), 367. [@feynman1948space]

Textbooks

• Shreve, S. E. (2004). Stochastic Calculus for Finance II: Continuous-Time Models. Springer. [@shreve2004stochastic] [*]
• Steele, J. M. (2001). Stochastic Calculus and Financial Applications. Springer. [@steele2001stochastic]
• Øksendal, B. (2003). Stochastic Differential Equations: An Introduction with Applications (6th ed.). Springer. [@oksendal2013stochastic] [*]
• Klebaner, F. C. (2012). Introduction to Stochastic Calculus with Applications. Springer. [@klebaner2012introduction]

Stochastic Optimal Control

Foundational Theory

• Bellman, R. (1952). On the theory of dynamic programming. PNAS, 38(8), 716--719. [@bellman1952theory] [*]
• Bellman, R. (1966). Dynamic programming. Science, 153(3731), 34--37. [@bellman1966dynamic]
• Pontryagin, L. S. (1962). Mathematical Theory of Optimal Processes. [@pontryagin2018mathematical]
• Merton, R. C. (1969). Lifetime portfolio selection under uncertainty: The continuous-time case. REStat, 247--257. [@merton1969lifetime]

Controlled Diffusion and HJB

• Krylov, N. V. (1980). Controlled Diffusion Processes. Springer. [@krylov1980controlled]
• Fleming, W. H. & Soner, H. M. (2006). Controlled Markov Processes and Viscosity Solutions. Springer. [@fleming2006controlled] [*]
• Fleming, W. H. & Rishel, R. W. (1975). Deterministic and Stochastic Optimal Control. Springer. [@fleming2012deterministic]
• Yong, J. & Zhou, X. Y. (1999). Stochastic Controls: Hamiltonian Systems and HJB Equations. Springer. [@yong1999stochastic] [*]
• Pham, H. (2009). Continuous-Time Stochastic Control and Optimization with Financial Applications. Springer. [@pham2009continuous] [*]
• Dynkin, E. B. (1965). Markov Processes. Springer. [@dynkin1965markov]
• Dynkin, E. B. & Yushkevich, A. A. (1979). Controlled Markov Processes. Springer. [@dynkin1979controlled]
• Bertsekas, D. (2012). Dynamic Programming and Optimal Control, Vol. I. [@bertsekas2012dynamic]
• Øksendal, B. & Sulem, A. (2009). Applied Stochastic Control of Jump Diffusions (3rd ed.). Springer. [@oksendal2009applied] [*]

Viscosity Solutions

• Crandall, M. G. & Lions, P.-L. (1983). Viscosity solutions of Hamilton-Jacobi equations. Trans. AMS, 277(1), 1--42. [@crandall1983viscosity] [*]
• Crandall, M. G., Ishii, H. & Lions, P.-L. (1992). User's guide to viscosity solutions of second order PDEs. Bull. AMS, 27(1), 1--67. [@user2013users] [*]

Reinforcement Learning (Context)

• Sutton, R. S. & Barto, A. G. (2018). Reinforcement Learning: An Introduction (2nd ed.). MIT Press. [@sutton2018reinforcement]
• Schulman, J. et al. (2017). Proximal policy optimization algorithms. [@schulman2017proximal]

Market Making

Foundational Models

• Avellaneda, M. & Stoikov, S. (2008). High-frequency trading in a limit order book. Quantitative Finance, 8(3), 217--224. [@avellaneda2008high] [*]
• Ho, T. & Stoll, H. R. (1981). Optimal dealer pricing under transactions and return uncertainty. J. of Financial Economics, 9(1), 47--73. [@ho1981optimal] [*]

Extensions and Analysis

• Guéant, O., Lehalle, C.-A. & Fernandez-Tapia, J. (2013). Dealing with the inventory risk: a solution to the market making problem. Math. Fin. Econ., 7(4), 477--507. [@gueant2013dealing] [*]
• Guilbaud, F. & Pham, H. (2013). Optimal high-frequency trading with limit and market orders. Quantitative Finance, 13(1), 79--94. [@guilbaud2013optimal]
• Bayraktar, E. & Ludkovski, M. (2011). Liquidation in limit order books with controlled intensity. Math. Finance, 24(4), 627--650. [@bayraktar2011liquidating]
• Cartea, Á., Jaimungal, S. & Ricci, J. (2014). Buy low, sell high: A high frequency trading perspective. SIAM J. Financial Math., 5(1), 415--444. [@cartea2014buy]
• Cartea, Á. & Sánchez-Betancourt, L. (2021). Shadow prices for optimal market making. SIAM J. Financial Math., 12(3). [@cartea2021shadow]
• Lehalle, C.-A. & Mounjid, O. (2017). Limit order strategic placement with adverse selection risk. Market Microstructure and Liquidity. [@lehalle2017limit]
• Guéant, O. (2017). Optimal market making. Applied Mathematical Finance, 24(2), 112--138. [@gueant2017optimal]
• Bergault, P., Evangelista, D., Guéant, O. & Vieira, D. (2021). Closed-form approximations in multi-asset market making. Applied Mathematical Finance, 28(2), 101--126. [@bergault2021closed] [*]
• Bergault, P. & Guéant, O. (2021). Size matters for OTC market makers: general results and dimensionality reduction techniques. Mathematical Finance, 31(3). [@bergault2021size]

Textbooks and Surveys

• Cartea, Á., Jaimungal, S. & Penalva, J. (2015). Algorithmic and High-Frequency Trading. Cambridge. [@cartea2015algorithmic] [*]
• Guéant, O. (2016). The Financial Mathematics of Market Liquidity. CRC Press. [@gueant2016financial]
• Lehalle, C.-A. & Laruelle, S. (2013). Market Microstructure in Practice. World Scientific. [@lehalle2013market]
• O'Hara, M. (1995). Market Microstructure Theory. Blackwell. [@ohara1995market]
• Menkveld, A. J. (2013). High frequency trading and the new market makers. J. Financial Markets, 16(4), 712--740. [@menkveld2013high]
• Brogaard, J., Hendershott, T. & Riordan, R. (2014). High-frequency trading and price discovery. Rev. Financial Studies, 27(8), 2267--2306. [@brogaard2014high]

Reinforcement Learning Approaches

• Spooner, T. et al. (2018). Market making via reinforcement learning. AAMAS. [@spooner2018market]
• Briola, A. et al. (2021). Deep reinforcement learning for active high frequency trading. [@briola2021deep]
• Qin, X. et al. (2023). EarnHFT: Efficient hierarchical RL for HFT. [@qin2023earnhft]

Hawkes Processes

• Hawkes, A. G. (2018). Hawkes processes and their applications to finance: a review. Quantitative Finance, 18(2), 193--198. [@hawkes2018hawkes] [*]
• Ogata, Y. (1988). Statistical models for earthquake occurrences and residual analysis for point processes. JASA, 83(401), 9--27. [@ogata1988statistical] [*]
• Rizoiu, M.-A. et al. (2018). SIR-Hawkes: Linking epidemic models and Hawkes processes. WWW. [@rizoiu2018sir]

Stochastic Volatility Models

• Heston, S. L. (1993). A closed-form solution for options with stochastic volatility. Rev. Financial Studies, 6(2), 327--343. [@heston1993closed] [*]
• Lord, R. et al. (2010). A comparison of biased simulation schemes for stochastic volatility models. Quantitative Finance, 10(2), 177--194. [@lord2010comparison]

Credit Risk and Default

Structural and Reduced-Form Models

• Merton, R. C. (1974). On the pricing of corporate debt: The risk structure of interest rates. Journal of Finance, 29(2), 449--470. [@merton1974pricing] [*]
• Jarrow, R. A., Lando, D. & Turnbull, S. M. (1997). A Markov model for the term structure of credit risk spreads. Review of Financial Studies, 10(2), 481--523. [@jarrow1997markov]
• Bielecki, T. R. & Rutkowski, M. (2004). Credit Risk: Modeling, Valuation and Hedging. Springer. [@bielecki2004credit] [*]

Optimal Investment and Contagion

• Kraft, H. & Steffensen, M. (2007). Bankruptcy, counterparty risk, and optimal investment. Finance and Stochastics, 11(1), 131--163. [@kraft2007bankruptcy]
• Capponi, A. & Figueroa-López, J. E. (2014). Dynamic portfolio optimization with a defaultable security and regime-switching. Mathematical Finance, 24(2), 207--249. [@capponi2014dynamic]
• Bo, L. & Capponi, A. (2016). Optimal investment in credit derivatives portfolio under contagion risk. Mathematical Finance, 26(4), 785--834. [@bo2016optimal]
• Sircar, R. & Zariphopoulou, T. (2010). Utility valuation of multiname credit derivatives and application to CDOs. Quantitative Finance, 10(2), 195--208. [@sircar2010utility]

Numerical Methods

Finite Difference and HJB Discretization

• Forsyth, P. A. & Labahn, G. (2007). Numerical methods for controlled Hamilton-Jacobi-Bellman PDEs in finance. J. Computational Finance, 11(2). [@forsyth2007numerical] [*]
• Achdou, Y. et al. (2013). Hamilton-Jacobi Equations: Approximations, Numerical Analysis and Applications. Springer. [@achdou2013hamilton]
• Barles, G. & Souganidis, P. E. (1991). Convergence of approximation schemes for fully nonlinear second order equations. Asymptotic Analysis, 4(3), 271--283. [@barles1991convergence]
• Kushner, H. J. & Dupuis, P. G. (2001). Numerical Methods for Stochastic Control Problems in Continuous Time. Springer. [@kushner2001numerical] [*]
• Howard, R. A. (1960). Dynamic Programming and Markov Processes. MIT Press. [@howard1960dynamic]
• Courant, R., Friedrichs, K. & Lewy, H. (1928). Über die partiellen Differenzengleichungen der mathematischen Physik. Math. Ann., 100, 32--74. [@courant1928partiellen]
• Rannacher, R. (1984). Finite element solution of diffusion problems with irregular data. Numerische Math., 43(2), 309--327. [@rannacher1984finite]

BSDE Theory and Numerics

• Pardoux, E. & Peng, S. (1990). Adapted solution of a backward stochastic differential equation. Systems & Control Letters, 14(1), 55--61. [@pardoux1990adapted] [*]
• Pardoux, E. & Peng, S. (2005). BSDEs and quasilinear parabolic PDEs. Lecture Notes in Control and Inf. Sci., 176, 200--217. [@pardoux2005backward]
• El Karoui, N., Peng, S. & Quenez, M. C. (1997). Backward stochastic differential equations in finance. Mathematical Finance, 7(1), 1--71. [@el1997backward]
• Ma, J., Protter, P. & Yong, J. (1994). Solving forward-backward SDEs explicitly---a four step scheme. Probab. Theory Rel. Fields, 98(3), 339--359. [@ma1994solving]
• Ma, J. & Yong, J. (1999). Forward-Backward Stochastic Differential Equations and Their Applications. Springer. [@ma1999forward]
• Gobet, E., Lemor, J.-P. & Warin, X. (2005). A regression-based Monte Carlo method to solve BSDEs. Annals of Applied Probability, 15(3), 2172--2202. [@gobet2005empirical] [*]

Deep Learning and High-Dimensional PDEs

• Han, J., Jentzen, A. & E, W. (2018). Solving high-dimensional PDEs using deep learning. PNAS, 115(34), 8505--8510. arXiv:1707.02568 [@han2018solving]
• Han, J. & Jentzen, A. (2017). Deep learning-based numerical methods for high-dimensional parabolic PDEs and BSDEs. arXiv:1706.04702 [@han2017deep]
• Al-Aradi, A. et al. (2022). Extensions of the deep Galerkin method. [@al2022extensions]
• Raissi, M., Perdikaris, P. & Karniadakis, G. E. (2019). Physics-informed neural networks. J. Computational Physics, 378, 686--707. arXiv:1711.10561 [@raissi2019physics]
• Karniadakis, G. E. et al. (2021). Physics-informed machine learning. Nature Reviews Physics, 3(6), 422--440. [@karniadakis2021physics]
• Song, Y. et al. (2021). Score-based generative modeling through SDEs. ICLR. arXiv:2011.13456 [@song2020score]
• Sirignano, J. & Spiliopoulos, K. (2018). DGM: A deep learning algorithm for solving partial differential equations. Journal of Computational Physics, 375, 1339--1364. arXiv:1708.07469 [@sirignano2018dgm] [*]
• Cheridito, P., Dupret, J.-L. & Hainaut, D. (2025). Deep learning for continuous-time stochastic control with jumps. NeurIPS. arXiv:2505.15602 [@cheridito2025jumps]

Neural Operators and Architectures

• Vaswani, A. et al. (2017). Attention is all you need. NeurIPS. arXiv:1706.03762 [@vaswani2017attention]
• Amos, B., Xu, L. & Kolter, J. Z. (2017). Input convex neural networks. ICML. arXiv:1609.07152 [@amos2017input]
• Zaheer, M. et al. (2017). Deep sets. NeurIPS. arXiv:1703.06114 [@zaheer2017deep]
• Lee, J. et al. (2019). Set transformer: A framework for attention-based permutation-invariant neural networks. ICML. arXiv:1810.00825 [@lee2019set]
• Li, Z. et al. (2021). Fourier neural operator for parametric partial differential equations. ICLR. arXiv:2010.08895 [@li2021fourier]
• Lu, L., Jin, P., Pang, G., Zhang, Z. & Karniadakis, G. E. (2021). Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3, 218--229. arXiv:1910.03193 [@lu2021deeponet]
• Kovachki, N. et al. (2023). Neural operator: Learning maps between function spaces with applications to PDEs. Journal of Machine Learning Research, 24(89), 1--97. arXiv:2108.08481 [@kovachki2023neural]

Group Equivariance and Hamiltonian Architectures

• Noether, E. (1918). Invariante Variationsprobleme. Nachr. Ges. Wiss. Goettingen, Math.-Phys. Kl., 235--257. [@noether1918invariante]
• Olver, P. J. (1993). Applications of Lie Groups to Differential Equations (2nd ed.). Springer. [@olver1993applications]
• Cohen, T. & Welling, M. (2016). Group equivariant convolutional networks. ICML. arXiv:1602.07576 [@cohen2016group]
• Weiler, M. & Cesa, G. (2019). General E(2)-equivariant steerable CNNs. NeurIPS. arXiv:1911.08251 [@weiler2019general]
• Thomas, N., Smidt, T., Kearnes, S., Yang, L., Li, L., Kohlhoff, K. & Riley, P. (2018). Tensor field networks: Rotation- and translation-equivariant neural networks for 3D point clouds. arXiv:1802.08219 [@thomas2018tensor]
• Fuchs, F., Worrall, D., Fischer, V. & Welling, M. (2020). SE(3)-transformers: 3D roto-translation equivariant attention networks. NeurIPS. arXiv:2006.10503 [@fuchs2020se3]
• Satorras, V. G., Hoogeboom, E. & Welling, M. (2021). E(n)-equivariant graph neural networks. ICML. arXiv:2102.09844 [@satorras2021en]
• Finzi, M., Stanton, S., Izmailov, P. & Wilson, A. G. (2020). Generalizing convolutional neural networks for equivariance to Lie groups on arbitrary continuous data. ICML. arXiv:2002.12880 [@finzi2020lie]
• Worrall, D. & Welling, M. (2019). Deep scale-spaces: Equivariance over scale. NeurIPS. arXiv:1905.11697 [@worrall2019scale]
• Greydanus, S., Dzamba, M. & Yosinski, J. (2019). Hamiltonian neural networks. NeurIPS. arXiv:1906.01563 [@greydanus2019hamiltonian]
• Cranmer, M., Greydanus, S., Hoyer, S., Battaglia, P., Sperberg-McQueen, D. & Ho, S. (2020). Lagrangian neural networks. ICLR Deep Differential Equations Workshop. arXiv:2003.04630 [@cranmer2020lagrangian]
• Bronstein, M. M., Bruna, J., Cohen, T. & Velickovic, P. (2021). Geometric deep learning: Grids, groups, graphs, geodesics, and gauges. arXiv:2104.13478 [@bronstein2021geometric] ^bronstein2021geometric
• Wang, R., Walters, R. & Yu, R. (2022). Approximately equivariant networks for imperfectly symmetric dynamics. ICML. arXiv:2201.11969 [@wang2022approximately]

Neural Operators for Optimal Control

• Hwang, R., Lee, J. Y., Shin, J. Y. & Hwang, H. J. (2021). Solving PDE-constrained control problems using operator learning. AAAI, 36(4). arXiv:2111.04941 [@hwang2021solving]
• Wang, S., Bhouri, M. A. & Perdikaris, P. (2021). Fast PDE-constrained optimization via self-supervised operator learning. arXiv:2110.13297 [@wang2021fast]
• Lanthaler, S. & Stuart, A. M. (2023). The parametric complexity of operator learning. arXiv:2306.15924 [@lanthaler2023parametric] [*]
• Lee, J. Y. & Kim, Y. (2024). Hamilton-Jacobi based policy-iteration via deep operator learning. arXiv:2406.10920 [@lee2024hamilton] [*]
• Hoischen, N., Bevanda, P., Sosnowski, S., Hirche, S. & Houska, B. (2024). Data-driven stochastic optimal control in reproducing kernel Hilbert spaces. arXiv:2407.16407 [@hoischen2024datadriven]
• Kratsios, A., Neufeld, A. & Schmocker, P. (2025). Generative neural operators of log-complexity can simultaneously solve infinitely many convex programs. arXiv:2508.14995 [@kratsios2025generative]
• Xu, W., Han, J. & Lai, R. (2025). Self-supervised amortized neural operators for optimal control: Scaling laws and applications. arXiv:2512.24897 [@xu2025amortized]
• Cohen, S. N., de Feo, F., Hebner, J. & Sirignano, J. (2026). Deep Hilbert-Galerkin methods for infinite-dimensional PDEs and optimal control. arXiv:2603.19463 [@cohen2026hilbert]
• Kratsios, A., Livieri, G. & Schmocker, P. (2026). NeuralChaos: Optimal adapted approximation of square integrable predictable processes. arXiv:2607.14361 [@kratsios2026neuralchaos]
• Gao, S., Zhou, M. & Lai, R. (2026). Self-supervised in-context operator learning for stochastic mean-field control. arXiv:2608.18282 [@gao2026selfsupervised]
• Mohanty, S. K. (2026). Explainable artificial intelligence for financial integral equations: A fixed-point neural operator approach. arXiv:2604.27127 [@mohanty2026explainable]

Monte Carlo Methods

• Glasserman, P. (2004). Monte Carlo Methods in Financial Engineering. Springer. [@glasserman2004monte]

General Numerical Analysis

• Quarteroni, A., Sacco, R. & Saleri, F. (2007). Numerical Mathematics. Springer. [@quarteroni2007numerical]

Financial Mathematics (Broader Context)

• Black, F. & Scholes, M. (1973). The pricing of options and corporate liabilities. JPE, 81(3), 637--654. [@black1973pricing]
• Schwartz, E. S. (1997). The stochastic behavior of commodity prices. J. Finance, 52(3), 923--973. [@schwartz1997stochastic]
• Vasicek, O. (1977). An equilibrium characterization of the term structure. J. Financial Economics, 5(2), 177--188. [@vasicek1977equilibrium]
• Sharpe, W. F. (1966). Mutual fund performance. J. Business, 39(1), 119--138. [@sharpe1966mutual]
• Sortino, F. A. & Price, L. N. (1994). Performance measurement in a downside risk framework. J. Investing, 3(3), 59--64. [@sortino1994performance]

Auxiliary

• Gould, J. P. (1968). Adjustment costs in the theory of investment of the firm. REStud, 35(1), 47--55. [@gould1968adjustment]
• Abel, A. B. & Eberly, J. C. (1993). A unified model of investment under uncertainty. NBER WP 4296. [@abel1993unified]
• Gode, D. K. & Sunder, S. (1993). Allocative efficiency of markets with ZI traders. JPE, 101(1), 119--137. [@gode1993allocative]
• Klabnik, S. & Nichols, C. (2018). The Rust Programming Language. No Starch Press. [@klabnik2023rust]
• Kalman, R. E. (1960). A new approach to linear filtering and prediction problems. J. Basic Eng., 82(1), 35--45. [@kalman1960new]

Plan

Navigation only. The canonical plan lives under docs/project/plan/.

• Overview
• Status
• Decisions
• Workstreams
• Overview diagram

Completed Phases

Each entry is a one-line dated note. Design rationale lives in mdbook or rustdoc pages, not here.

Phase 16: Neural Solver Workstream — 2026-08-29

Completed the neural_solver package: package scaffolding and committed reference fixtures, deep BSDE for diffusions, Rust/JAX reference agreement, jump-aware deep BSDE with deterministic and log-normal Merton jump reference models, DGM diffusion helpers, and the benchmark ladder. All six stages are Done in docs/project/plan/status.md.

Phase 15: Correlated Diffusions — 2026-08-22

Added the full symmetric Hessian (StateDerivatives::hessian_full and with_full_hessian) and the centered four-corner mixed-derivative stencil in both the BSDE and FD derivative paths. Generalized the LQ regulator generator to the full covariance C C' cross term and validated a correlated two- diffusion LQ model against its Riccati closed form (lq_correlated_fd.rs), keeping the reduced Heston rho drift correction unchanged because the spot-variance correlation is an exact reduced-form term, not a missing grid cross-derivative.

Phase 14: Arbitrary Jump Processes — 2026-08-22

Added JumpKernel<N>, JumpTransition, DimensionKind::Jump, and PdeProblem::jump_kernel with a unit-jump default, generalized the explicit-Euler FD path to sum arbitrary jump kernels, added the AvellanedaLotSize model with a lot_size = 1 reduction limit reproducing the base Avellaneda-Stoikov exact spreads, and documented the jump-kernel convention in solver_pde.md, solver_bsde.md, and the validation test matrix.

Phase 13: Exact-Reduction and Manufactured-Solution Validation — 2026-08-22

Validated every expanded market-making model's degenerate reduction to a closed-form Avellaneda-Stoikov target (Heston, Hawkes, Heston-Hawkes, bilateral Hawkes, order-flow imbalance) at high base intensity with per-component bid/ask checks, fixed the GlobalDiscrete initialization bug, and added a manufactured-solution oracle recovered by FD to machine precision and BSDE to 0.0025.

Phase 12: Stationary and Ergodic Solvers — 2026-08-15

Implemented EllipticProblem, EllipticControlProblem, StationarySolver, StationaryEigenProblem, and PerronSolver. Validated against infinite-horizon LQ (algebraic Riccati) and AvellanedaGueant.

Phase 11: Generic Control and PDE Contracts — 2026-08-15

Implemented ControlProblem<N>, PdeProblem<N>, StateDerivatives, and GridSolution. Replaced the market-making-specific Model<N> and ControlOutput contract. Validated with Merton, LQ regulator, and American put; ported the market-making models to ControlProblem.

Phase 10: Signal-Based Engine Strategies — 2026-08-02

Created SignalEngineStrategy bridging market_model::PriceStrategy to the engine's Strategy trait. Supports multi-signal voting (SignalCombinator) and configurable position sizing. 6 unit tests. Updated engine_strategies.md and implementation_matrix.md.

Phase 8: Rigorous Kelly HJB Solution — 2026-08-02

Implemented KellyHjb model solving the true log-utility HJB via FD policy iteration on a 2D [q, x] grid. Created KellyRigorousStrategy with 3D [q, x, tau] lookup tables. 11 unit tests. Updated optimal_criteria.md with full derivation.

Phase 7: Kelly Criterion Optimal Strategy (Heuristic) — 2026-08-02

Implemented KellyStrategy with online EWMA drift/vol estimation and Kelly target q* = mu/(gamma * sigma^2). 9 tests. Updated docs.

Phase 5: Edition 2024 Migration — 2026-08-01

Changed edition = "2021" to "2024" in all Cargo.toml files. Renamed gen to grid_gen. Fixed deref-pattern ambiguity and collapsed if let chains per edition 2024 clippy rules.

Phase 4: VecEnv Process Genericity — 2026-08-01

Replaced hardcoded Heston with VecEnvProcess enum (GBM, Heston, Bates). Dynamic state_dim. Updated engine_vec_env.md.

Phase 3: Unify Process Dynamics Across Crates — 2026-08-01

Added type Process and fn process() to Model<N>. Process = () for 8 models, Process = HestonProcess for Heston. Integration test confirms model-process consistency. Updated solver_ndim.md.

Phase 2: Solver Crate Structure — 2026-08-01

Moved examples, merged Matrix into solver/src/linalg/, removed a/kappa from Solver::solve_with_spreads. Restructured slow-tests.

Phase 1: Code Cleanup — 2026-08-01

Removed debug print from solver/src/models/avellaneda.rs, added all_outputs() to SimulationResult.

Phase 0: Mathematical Reference Documentation — 2026-08-01

Created docs/src/reference.md with verified equations for all processes, HJB derivations, analytical solutions, and numerical schemes.

Optimizations

Performance, memory, or allocation improvements with no behavioural change. Numbered O1, O2, ... and never part of phase numbering.

Gate requirements for every optimization:

• All existing tests pass with identical output to the pre-change baseline.
• No test assertions were relaxed.
• The relevant benchmark shows measurable improvement.

Pending

O5: Chunked Simulation Runner

SimulationRunner pre-generates all normals and keeps all paths in memory simultaneously. For large simulations (100k+ paths) this exceeds 1 GB working set.

• Add chunk_size: usize to SimulationConfig (default 0 = all at once).
• Split paths into chunks; generate normals per chunk.
• Deterministic RNG: chunk k uses seed + path_offset + 1.
• Memory benchmark: 100k paths x 1000 steps, chunk_size=0 vs 1000.

Target: Peak memory reduced by >= 5x with chunk_size=1000.

Completed

O7: Cache the LQ Riccati solution — 2026-08-29

Cached the backward Riccati trajectory (P(t), q(t), and R^-1) in an OnceLock so each time level is integrated once instead of once per grid node per time step. lq_correlated_fd.rs dropped from 15.11 s to 0.32 s in the debug profile (~47x) and 1.27 s to 0.03 s in release (~42x). Riccati values bit-identical to the post-fix baseline.

O6: Generalize Solver Output Layer — 2026-08-02

Replaced NumericalResult<N> with GridSolution<N> (generic) and SpreadResult (MM-specific). Solver::solve returns GridSolution<N>, solve_with_spreads returns SpreadResult. Updated 20+ call sites. Zero behavioural change.

O4: Process Micro-optimizations — 2026-08-02

Precomputed fields in GBM, Heston, RoughOU. GBM -12.1%, Heston -19.1%, RoughOU -12.5%. Zero behavioural change.

O3: Solver FD Throughput — 2026-08-02

Auto-route N=2 Implicit/CN to Strang ADI. 2D 101x100 500 steps: 87.9 ms (~1.5x). Zero behavioural change.

O2: Hawkes Intensity O(1) Recurrence — 2026-08-02

Recursive intensity update, periodic prune. 500x1000: 36.0 ms. <3% relative error in mean intensity. Zero behavioural change.

O1: Engine Allocation Hot Path — 2026-08-02

Stack arrays, capacity pre-allocation, rayon parallel VecEnv. 904 steps/sec for N=1024 (target >=500). Zero behavioural change.

Issues

Discovered problems out of scope for the current change. Each entry has a label, affected files, description, and fix recommendation. Resolved issues are deleted from this file, not archived.

Performance/allocation issues with no behavioural impact are tracked as optimizations (O1, O2, ...), not here.

Active

I26: Sign error in soc_exact.md AS price-diffusion intermediate step

• Files: docs/src/reference/soc_exact.md (CARA separation and reduced PDE)
• Problem: The step deriving the inventory penalty writes 1/2 sigma^2 partial_S^2 V = -1/2 gamma^2 sigma^2 q^2 V. With the ansatz V = -exp(-gamma (X + q S + theta)) the correct derivative is partial_S^2 V = +gamma^2 q^2 V, so the term is +1/2 gamma^2 sigma^2 q^2 V; the written sign is reversed. The final reduced PDE (penalty -1/2 gamma sigma^2 q^2) is correct, so the intermediate equality is internally inconsistent with it.
• Fix: Flip the sign of the intermediate equality to +1/2 gamma^2 sigma^2 q^2 V, or drop the intermediate step and state the reduced penalty directly.

I20: AvellanedaExact::value_function_theta returns a gauge-shifted theta

• Files: solver/src/numeric/ode.rs (LinearSpectralSolver::solve_tridiagonal), solver/src/analytical/avellaneda/exact.rs (value_function_theta), neural_solver/tests/data/reference.json (the avellaneda_stoikov.theta values), neural_solver/src/neural_solver/models/avellaneda_stoikov.py.
• Problem: solve_tridiagonal computes exp(M tau) y0 but returns it scaled by exp(-max_exp) for numerical stability. The scaling cancels in the spread ratios (v_q / v_{q+1}), so exact_spreads is correct, but value_function_theta applies ln to the scaled vector, so its returned theta is shifted by -max_exp / kappa. The absolute theta (and hence the value function) is therefore not the physical reduced value, even though the observable spreads are unaffected.
• Fix: Either multiply the returned vector back by exp(max_exp) in solve_tridiagonal before any ln is taken, or make value_function_theta add max_exp / kappa back, and regenerate reference.json.

I19: Heston finite-horizon BSDE diverges at the a = 140 test point

• Files: solver/tests/numerical/bsde_vs_fdm/consistency.rs (test_heston_finite_horizon_bsde_matches_fdm_spread_values), solver/src/models/heston.rs

• Problem: The decoupled BSDE regression diverged at a = 140 (bid = 10.0, ask = -5.0, fitted value ~5.6e4). Two independent root causes were traced:

  1. Reward/forward discretization inconsistency. Heston::optimize computes the fill intensity as lambda = a * exp(-kappa * delta) with no cap (unlike the Hawkes models). When a transient regression gradient drives a spread to its lower clamp (delta = -5), the intensity becomes a * exp(7.5). At a = 140 that is ~2.5e5, and running_reward = (lambda_bid + lambda_ask) / (gamma + kappa) blows up, closing a positive-feedback loop. The forward step clamps p = lambda * dt to [0, 1], but the reward used the unclamped lambda, so the backward step counted fills the forward step cannot produce. Fix: bound the reward rate at 1/dt per fill side in bsde_driver only, matching the forward probability clamp without touching the lambda returned to to_spreads. Capping the intensity itself (e.g. .min(20)) is incorrect: to_spreads inverts delta = -ln(lambda/a)/kappa, so at a = 140 even the correct symmetric lambda = 59 exceeds any absolute cap, and capping would corrupt the recovered spread and the FD path.

  2. Decoupled proxy-fill forward pass. The finite-horizon test ran BsdeSolver::new (decoupled). For a full-value market-making model the inventory dynamics depend on the optimal control, so the decoupled forward pass simulated inventory fills from the frozen base-spread proxy while the backward step recovered a skewed optimal control. This forward/backward inconsistency produced a large anti-symmetric bid/ask skew at q = 0 (bid 1.03, ask 0.13, total still matching FD). The total spread matched FD but each component was off by ~0.45. Fix: the test uses BsdeSolver::coupled (Picard iteration), which re-simulates the forward inventory under the current optimal control via next_step_controlled. With coupled mode the bid/ask components match FD to ~0.005.

  These two fixes are mathematically required, not duct tape: the reward bound restores the discretization-consistent backward driver, and coupled mode restores the forward/backward consistency of a control-dependent state process.

I9: policy_iteration.rs / operator.rs re-export shims

• Files: solver/src/numeric/policy_iteration.rs, solver/src/numeric/operator.rs
• Problem: Both are thin pub use super::finite_difference::*; shims from Phase 2. numeric/mod.rs still calls them "backward-compatible."
• Fix: Update callers to solver::numeric::finite_difference:: and delete the shims, or document them as permanent.

I11: engine::python::market::legacy duplicates simulation + engine

• File: engine/src/python/market/legacy.rs
• Problem: Defines PyEngine, PyEngineGbmAs, etc. plus run_many_* functions, all hardcoded to GBM/Heston/Bates + AS. The newer PySimEngine and run_simulation{,_batch,_summary} cover the same ground generically.
• Fix: Confirm no external consumers, then delete legacy.rs and its registrations. Otherwise mark #[deprecated].

I14: Python vec_env bindings out of sync with VecEnv

• File: engine/src/python/market/vec_env.rs
• Problem: The binding imports crate::vec_env::VecEnvProcess (no longer present), sets VecEnvConfig { process } (the field no longer exists; the config now holds Heston parameters directly), and calls VecEnv::state_dim (the method was removed). cargo check -p engine --features python fails for this reason.
• Fix: Rewrite the binding against the current VecEnv API. Either expose the fixed 6D state width or reintroduce a state_dim accessor, and drop the removed VecEnvProcess/process indirection.

I15: Orphaned convergence_bump.rs references removed solver APIs

• File: solver/tests/numerical/exact_vs_finite_diff/convergence_bump.rs
• Problem: The file is not declared in exact_vs_finite_diff/mod.rs, so it is never compiled. It still imports the removed solver::models::traits::{Gradients, Model}, calls PolicyIterationSolver::solve_grid, and calls Operator::compute_gradients.
• Fix: Either delete the file (it duplicates convergence_base.rs and convergence_market_impact.rs) or port it to the current ControlProblem/solve_grid_control API and wire it into the module.

I16: Market-making BSDE spread magnitude is ~10x the exact solution

• Files: solver/src/numeric/bsde/solver.rs, solver/src/models/control.rs, solver/src/models/avellaneda.rs (and drift/impact variants)
• Problem: The BSDE backward step used only running_reward (the reduced Hamiltonian H*) and omitted the reduced generator (the local inventory-risk source -0.5 gamma sigma^2 q^2), so the market-making value had no inventory skew and the total spread collapsed to 2 * base_spread. The configured InitializationMode was also never applied: every forward path started at the exact initial state, giving a rank-one design that could not identify the inventory gradient. The streaming normal-equation regression squared the design condition number and added a one-dimensional value bias (Merton 9.2%, LQ up to 66%). The impact model's forward proxy also omitted the permanent impact shift.
• Fix: Added a ControlProblem::bsde_driver hook (default running_reward for full value problems; reduced market-making models add the local source), applied InitializationMode to the initial forward cloud (integer lattice for discrete inventory, continuous for diffusion), switched the default regression backend to Householder QR with a minimal ridge floor for rank-deficiency, and included the impact shift in the impact forward proxy.

I18: Finite-horizon LQ BSDE value carries O(dt) discretization bias

• Files: solver/src/models/lq_regulator.rs, solver/src/numeric/bsde/solver.rs
• Problem: The LQ value estimate is biased 0.47% at dt=0.005 and 0.30% at dt=0.0025, scaling with dt (Euler-Maruyama discretization error), not with path count. The <= 0.1% high-precision target requires dt <= ~0.001 at 100_000 paths, which is a slow one-time run, not a defect. Documented here so the discretization floor is tracked separately from Monte Carlo noise.
• Fix: No code change required; record the achieved dt-limited precision in bsde_precision.md and note that the <= 0.1% target for LQ needs a finer time step.
• Stale reference: the 0.47%/0.30% figures were measured against LqRegulator::exact_value, which before the Riccati time-indexing fix returned the constant P(0) for every remaining horizon instead of the time-varying P(t). Re-measure the bias against the corrected reference before finalizing the bsde_precision.md numbers.

I21: BSDE mixed-derivative stencil is computed but never consumed

• Files: solver/src/numeric/bsde/solution.rs (state_derivatives), solver/src/models/control.rs (StateDerivatives)
• Problem: The BSDE state_derivatives centered four-corner stencil populates hessian_full[i][j] for i != j, but no model's optimize reads those off-diagonal entries. The LQ regulator's optimize returns the closed-form -R^-1 B' P(t) x, which is independent of the diffusion matrix C and ignores derivs entirely; the market-making models' optimize reads only the directional fwd/bwd entries. The Heston rho correlation is a reduced drift correction, not a hessian_full cross term. As a result the BSDE mixed derivative is dead output: it is unit-tested against a synthetic bilinear basis but never validated end to end against a closed form through a real solve. The FD mixed stencil is consumed by the LQ generator and is now guarded by cross_term_changes_fd_value.
• Fix: Add a correlated control problem whose optimizer genuinely depends on the off-diagonal Hessian (for example multiplicative/state-dependent noise, which has no simple Riccati closed form, or a cost coupling that makes the optimal control read hessian_full), then validate FD and BSDE against its exact or manufactured reference.

I23: Quadratic-generator BSDE test blows up from a wide initial range

• Files: solver/tests/numerical/bsde_solver_validation/analytical_baselines.rs (test_quadratic_generator_produces_finite_values)
• Problem: QuadraticGeneratorModel's forward step dy = y^2 dt has finite-time blow-up at t = 1/y. The solver is configured with with_initial_range(5.0), which perturbs the initial cloud over a wide interval around y_init = 0.5; paths seeded with y > 2 blow up before the horizon = 0.5 endpoint, so solve_trajectory_control returns a non-finite value and assert!(bsde_value.is_finite()) fails. The test is unseeded, so the failure is flaky.
• Fix: Narrow the initial range (for example with_initial_range(0.5)) so the whole cloud stays inside the blow-up radius for the horizon, add with_seed(...), or cap the forward state / shorten the horizon.

I24: Hawkes finite-horizon BSDE vs FDM spread test exceeds tolerance

• Files: solver/tests/numerical/bsde_vs_fdm/consistency.rs (test_hawkes_finite_horizon_bsde_matches_fdm_spread_values)
• Problem: The bid or ask spread recovered by the BSDE solver differs from the FD reference by more than tol = 0.15. The test uses 30_000 paths, 3 unseeded repeats, a Hermite(3) basis, and no fixed seed, so it is subject to Monte Carlo variance and the same forward/backward bias family documented in I16, I18, and I19. It fails intermittently rather than deterministically.
• Fix: Seed the solver with with_seed(...) and either raise num_paths / repeats or relax the tolerance after measuring the actual bias. Check whether the Hawkes forward proxy omits a local source (cf. I16) before treating the residual as pure noise.