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.