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