Deep BSDE for high-dimensional Black-Scholes

A clean-room PyTorch replication of Han, Jentzen, E (PNAS 2018) — "Solving high-dimensional partial differential equations using deep learning" — applied to the multi-asset Black-Scholes PDE.

The headline: Monte Carlo's pricing variance scales as $1/\sqrt{N}$ per dimension and finite-difference schemes die at $d \approx 4$, but a deep BSDE network prices a 100-asset basket European call directly from the PDE with convergence behaviour that is essentially flat in dimension and agrees with a $10^{7}$-path MC oracle to 0.22 % at $d = 100$.

Headline numbers

d	Deep BSDE $Y_0$	MC oracle (10⁷ paths)	rel. err	train time	device
1	10.4513	10.4524 ± 0.0064	0.011 %	1555 s	CPU
10	5.6387	5.6457 ± 0.0015	0.123 %	790 s	CPU
50	4.9134	4.9221 ± 0.0004	0.178 %	461 s	GPU
100	4.8705	4.8814 ± 0.0003	0.224 %	552 s	GPU

Against the brief's success criteria:

• 1-D Deep BSDE matches closed-form Black-Scholes to 0.011 % (target <0.5 %). ✔
• 100-D basket call matches 10⁷-path MC to 0.22 % (target <1 %). ✔
• python eval/convergence_vs_dim.py renders the hero figure above. ✔
• Variance across three independent training seeds at $d = 10$: σ/μ = 0.03 %. ✔
• Delta hedge using learned $Z$ reduces terminal P&L variance by 42–204 × across dimensions (see below). ✔

Approach in one paragraph

Under the risk-neutral measure the underlyings $S_t = (S^1_t, \dots, S^d_t)$ follow correlated geometric Brownian motions. The option price $u(t, S)$ satisfies the $d$-dimensional Black-Scholes PDE with terminal payoff $g(S_T)$. The nonlinear Feynman-Kac theorem reformulates this as a backward stochastic differential equation (BSDE):

$$ Y_t ;=; g(X_T) ;+; \int_t^T f(s, X_s, Y_s, Z_s),\mathrm ds ;-; \int_t^T Z_s^{\top},\mathrm dW_s, \qquad Y_t = u(t, X_t),; Z_t = \sigma(t, X_t)^{\top}\nabla u(t, X_t). $$

The Deep BSDE algorithm parameterises the initial value $Y_0$ as a free scalar, the initial gradient proxy $Z_0$ as a free vector, and each $Z_{t_n}$ for $n \ge 1$ as a neural network $\mathcal Z_n(X_{t_n})$. It then time-marches the discretised BSDE forward along simulated paths and minimises $\mathbb E[,|Y_N - g(X_N)|^2,]$ over $(Y_0, Z_0, {\theta_n})$. The learned $Y_0$ is the price.

Full derivation with the Itô step and the martingale-representation argument: docs/derivation.md.

Results

Convergence vs dimension (hero)

Deep BSDE tracks the MC oracle to a fraction of 1 % across $d \in {1, 10, 50, 100}$ — no visible dimension dependence.

1-D price-function reconstruction

Each marker is an independent Deep BSDE training run at the given initial spot $S_0$; the solid curve is the closed-form Black-Scholes price. The learned points sit on the curve across $S_0 \in [70, 130]$ — deep OTM ($S_0 = 70$) is the hardest regime for the algorithm and reaches ~1.8 % relative error after tightening the $Y_0$ init range; near-ATM and ITM runs agree with the closed form to well under 0.5 %.

Training dynamics

Loss drops ~3 orders of magnitude over training. The learned $Y_0$ curves converge onto the MC-oracle lines (dotted) within roughly 1000 iterations at every dimension.

Animated version of the 100-D run: docs/animations/training_100d.gif.

Variance across training seeds

Three independent seeds at $d = 10$ give a relative standard deviation of 0.03 % — the learned price is stable to three basis points under re-initialisation and mini-batch resampling.

Delta-hedge validation of the learned $Z$ process

The convergence table only checks $Y_0$. The Deep-BSDE algorithm also learns $Z_t = \sigma(t, X_t)^{\top} \nabla u(t, X_t)$, which — inverted — is the per-asset delta $\partial u / \partial S_i$. To validate that $Z$ is meaningful (not just a nuisance that makes the loss finite), a seller runs a discrete-time self-financing delta hedge with $\Delta_t$ read off the learned $Z_t$ and the residual placed in a risk-free bank account; the terminal P&L $V_T - g(S_T)$ is the hedge residual.

d	$\sum_i \Delta^0_i$ (basket delta)	un-hedged $\sigma$	hedged $\sigma$	variance reduction
1	0.637 (= $N(d_1)$, closed-form 0.637)	15.55	1.09	× 204
10	0.791	5.60	0.45	× 156
100	0.996	2.12	0.33	× 42

In 1-D the learned time-0 delta agrees with the closed-form Black-Scholes delta $N(d_1)$ to four decimal places — i.e. the Deep-BSDE algorithm does not just learn the scalar price $Y_0$, it learns a genuine gradient process usable as a hedge ratio. The basket-delta pattern across $d$ is what you'd expect from the effective-volatility argument $\sigma_{\text{eff}} = \sigma / \sqrt{d}$ for an equally-weighted, uncorrelated basket: $d = 100$ asymptotes to near-1 because the basket is nearly deterministic at its expectation, pushing the call deep ITM in effective-vol units. Variance reduction is 156× at $d = 10$ and 42× at $d = 100$; the fall-off at large $d$ is consistent with the $Z$-subnetworks being harder to train than the scalar $Y_0$, and is a candidate for further tuning.

Hyper-parameter sensitivity (d = 10)

Ablation around the default schedule — learning rate, number of time steps, mini-batch size. Deep BSDE is not balanced on a knife-edge: across all three axes the relative error sits in a 0.1–0.3 % band, except the one genuinely bad setting (lr = 1e-3, which is simply too small for 4 000 iterations at $d = 10$ to reach convergence). Larger batch size and finer time discretisation both help monotonically, consistent with standard Euler-Maruyama intuition.

Second PDE family — Allen-Cahn

The whole point of the docs/derivation.md translation is that only three callables change between PDEs: drift $\mu$, diffusion $\sigma$, and the generator $f$. To prove the point, pde/allen_cahn.py + the generic solver in bsde/generic_solver.py tackle the Allen-Cahn benchmark from Han-Jentzen-E 2018 sec. 4.3:

$$ \partial_t u + \Delta u + u - u^3 = 0, \qquad u(T, x) = \frac{1}{2 + \tfrac{2}{5}|x|^2}, \qquad T = 0.3,\ d = 100. $$

The nonlinear term flips sign on transit to the BSDE generator ($f(y) = y^3 - y$, not $y - y^3$ — see the docstring in pde/allen_cahn.py). Learned $Y_0 = 0.0529$ against the paper's branching-diffusion reference $\approx 0.0528$ — relative error 0.14 % (paper's own Deep BSDE entry: 0.30 %). Wall clock: ~8 min on a single GPU.

Two practical lessons, worth naming for the next person:

1. Generator sign. The BSDE form of the PDE puts the nonlinearity on the right-hand side of ∂_t u + ½ Tr(σσᵀ Hess u) + ∇u·μ = f. For Allen-Cahn this sends u - u³ over to u³ - u. Getting this backwards converges the solver to a physically wrong but loss-minimal solution (see commit history for the debug trail).
2. BN in the Z-subnet. The paper's BN-heavy Z-subnet ZSubnet is fine for Black-Scholes (linear $f$) but is too unstable for $y^3 - y$: BN at the output forces $Z \sim \mathcal{N}(0, I)$ per coordinate, which pushes $Y_n$ outside the stable region $|y| &lt; 1$ on the first rollout step. SimpleZSubnet (no BN, Tanh activations) + tight gradient clipping stabilises training. Included as an opt-in subnet="plain" on GenericSolverConfig.

PINN comparison on 1-D Black-Scholes

Side-by-side with a canonical Physics-Informed Neural Network (Raissi, Perdikaris, Karniadakis 2019) on the same 1-D BS problem — a $2 \to 50 \to 50 \to 50 \to 50 \to 1$ MLP with $\tanh$ activations, trained by jointly minimising the PDE residual and the terminal condition, inputs normalised to $[-1, 1]$.

With normalised inputs, PINN converges to 0.02 % relative error and is actually faster in wall-clock than Deep BSDE on this problem (~2 min vs ~26 min). This is the honest finding: on a smooth 1-D PDE with an easy terminal, PINN is hard to beat on raw wall time. The documented fragility of PINN shows up elsewhere — scale mismatch in inputs (without the $[-1, 1]$ normalisation the PINN above stalls at ~25 % error), higher dimensions, stiff or discontinuous payoffs, $\lambda_\text{PDE}/\lambda_\text{BC}$ balancing. Deep BSDE side-steps all of these by design; its headline advantage is not "faster at $d = 1$", it is "works at $d = 100$ where PINN gradient signals through a 100-D second derivative essentially vanish".

Wall-clock / compute study

Deep BSDE training is dominated by per-iteration cost, not by a cost that blows up with dimension. The right-hand panel shows the number of MC paths needed to match each Deep BSDE relative error — the basket average reduces per-path variance with $d$, so plain MC is actually cheap for this payoff; the point of the paper is not that Deep BSDE beats MC for vanilla basket calls, but that the same algorithm extends to nonlinear, path-dependent, or high-dimensional-without-variance-reduction PDEs where MC's constant factor explodes.

Limitations

What this replication does not do, by design, and where it's weakest:

• $Y_0$ is a point estimate, not a pricing surface. Each training run learns $u(0, S_0)$ at a single initial spot. The 1-D reconstruction plot above is five independent trainings — there is no amortised evaluation of $u(t, x)$ for arbitrary $(t, x)$. A mesh-free extension (e.g. training with batched initial conditions) would fix this but is out of scope.
• Deep OTM is the hardest regime. At $S_0 = 70$ (ATM strike 100) the relative error sits at 1.76 % even after tightening the $Y_0$ init range — the terminal-residual loss has little signal on $Y_0$ when $g(X_T)$ is zero on most paths. Deep BSDE prices OTM by path mass, and rare-payoff paths are under-represented in a 64–256-sample batch.
• Plain MC is very hard to beat on this payoff. An arithmetic basket averages out per-path noise, so MC's variance shrinks with $d$ rather than growing — see the right panel of the wall-clock plot. The method's true advantage lies on payoffs where MC's variance blows up (American options, path-dependent claims, nonlinear generators) — those are scope extensions, not covered here.
• Z-quality at large $d$ is not as tight as $Y_0$-quality. The delta-hedge study below validates the learned $Z$ against closed-form $N(d_1)$ in 1-D (agreement to 4 decimal places) and shows 156× variance reduction at $d = 10$, but the $d = 100$ hedge only reaches 42× variance reduction — the $Z$-subnetworks there are near-converged but still imperfect. Plausible improvements: larger batch size (currently 64 at $d = 100$), longer training after the $Y_0$ plateau, or finer time discretisation.
• Batch-norm during inference. The $Z$-subnetworks are trained in train() mode with BN statistics; evaluation of $Y_0$ itself is a read of a learnable parameter and is unaffected, but using the learned $Z_n$ on new inputs would require careful BN-eval handling.
• Architecture and schedules are not tuned. Hidden width $d+10$ and the LR milestones follow the paper; no search has been performed. The headline errors are likely sub-optimal by a factor of 2–3 with mild tuning.

Layout

deep-bsde-black-scholes/
├── pde/           Black-Scholes PDE: coefficients, payoffs, 1D closed form
├── bsde/          Euler-Maruyama forward SDE + Deep BSDE network
├── baselines/     Vanilla Monte-Carlo pricer (exact log-Euler, antithetic)
├── eval/          Convergence / variance / wall-clock / reconstruction
├── docs/          PDE → BSDE derivation, training-curve GIF
├── figures/       Rendered plots (hero figure, results table, etc.)
├── runs/          Cached training artifacts (JSON + state_dict)
├── tests/         Pytest suite (28 tests, all passing)
└── train.py       CLI training entry point

Reproduction

python -m venv .venv
# Windows bash:  source .venv/Scripts/activate
# macOS / Linux: source .venv/bin/activate
pip install -e ".[dev]"

# Unit tests — closed-form, SDE moments, MC, Deep BSDE loss / integration
pytest

# Single training run (1-D sanity)
python train.py --d 1 --iterations 2000 --output-dir runs/d1_seed0

# Regenerate every figure in the README from cached artifacts
# (auto-trains anything missing; ~1 h full fresh clone with a CUDA GPU)
python eval/make_all.py

Each training run is fully reproducible from --seed: the torch global RNG, the per-batch path sampler, and the $Y_0$/$Z_0$ initial draws are all seeded off a single integer.

Citations

• Han, J.; Jentzen, A.; E, W. (2018). Solving high-dimensional partial differential equations using deep learning. PNAS 115(34), 8505–8510.
• E, W.; Han, J.; Jentzen, A. (2017). Deep learning-based numerical methods for high-dimensional parabolic PDEs and BSDEs. arXiv:1706.04702.
• Black, F.; Scholes, M. (1973). The pricing of options and corporate liabilities. J. Polit. Econ. 81(3), 637–654.
• Pardoux, É.; Peng, S. (1990). Adapted solution of a backward stochastic differential equation. Systems Control Lett. 14, 55–61.
• El Karoui, N.; Peng, S.; Quenez, M. C. (1997). Backward stochastic differential equations in finance. Math. Finance 7(1), 1–71.

License

MIT. See LICENSE.