A Parallel Algorithm for solving BSDEs - Application to the pricing and hedging of American options

We present a parallel algorithm for solving backward stochastic differential equations (BSDEs in short) which are very useful theoretic tools to deal with many financial problems ranging from option pricing option to risk management. Our algorithm based on Gobet and Labart (2010) exploits the link between BSDEs and non linear partial differential equations (PDEs in short) and hence enables to solve high dimensional non linear PDEs. In this work, we apply it to the pricing and hedging of American options in high dimensional local volatility models, which remains very computationally demanding. We have tested our algorithm up to dimension 10 on a cluster of 512 CPUs and we obtained linear speedups which proves the scalability of our implementation

💡 Research Summary

The paper introduces a highly scalable parallel algorithm for solving backward stochastic differential equations (BSDEs) and applies it to the pricing and hedging of American options in high‑dimensional local‑volatility models. The authors build on the Gobet‑Labart (2010) framework, which exploits the well‑known correspondence between BSDEs and semilinear partial differential equations (PDEs). By translating the option‑pricing problem into a BSDE, the method circumvents the curse of dimensionality that plagues traditional grid‑based PDE solvers.

The core numerical scheme is an iterative Picard procedure combined with an adaptive Monte‑Carlo estimator. Starting from an initial guess (u_0\equiv0), each iteration (k) computes a correction term (c_k(t,x)) that represents the expected discrepancy between the true solution and the current approximation. This correction is evaluated by Monte‑Carlo simulation of the forward diffusion (X) (or its Euler discretisation (X^N)) and of the driver function (f). The correction is then added to the previous approximation, and a global extrapolation operator (P_k) is applied to obtain the next iterate (u_{k+1}=P_k(u_k+c_k)).

Two families of extrapolation operators are discussed. The original kernel‑based operator requires (O(n d)) operations per evaluation and suffers from delicate bandwidth tuning in high dimensions. The authors therefore adopt a least‑squares polynomial projection: a set of multivariate basis functions ({B_\ell}{\ell=1}^p) (tensor products of univariate polynomials) is fitted to the sampled values ({(t_i,x_i, u_k(t_i,x_i)+c_k(t_i,x_i))}) by solving a linear least‑squares problem. The resulting coefficients (\alpha) define the approximating function (\hat u(t,x)=\sum\ell \alpha_\ell B_\ell(t,x)). Numerical stability is ensured by centering and normalising the basis, and by solving the normal equations with a QR decomposition with pivoting (LAPACK routine dgelsy).

To handle American options, the reflected BSDE formulation is replaced by a penalisation approach: a sequence of standard BSDEs with a penalty term (\lambda (Y_t-\Phi(X_t))^-) is solved, with (\lambda) gradually increased. This yields a convergent approximation of the optimal early‑exercise boundary.

Parallelisation is achieved with MPI. Each processor independently draws a random set of space‑time points, runs (M) Monte‑Carlo paths, and computes local contributions to the correction. After the Monte‑Carlo stage, the local matrices and vectors needed for the polynomial fit are reduced (via MPI_Reduce) to a master process, which performs the QR factorisation, broadcasts the coefficients (\alpha), and all processes update their global approximation. Communication overhead scales as (O(p^2)), which is negligible compared with the computational load, leading to almost perfect linear speed‑up up to 512 cores.

Numerical experiments cover dimensions (d=1) to (d=10) for both European and American call options under a local‑volatility diffusion. With modest Monte‑Carlo sample sizes (e.g., (M=2000)) and a few thousand regression points, the algorithm attains mean absolute errors below (10^{-3}) even in ten dimensions. On a 512‑CPU cluster, a ten‑dimensional American option is priced and hedged in under six minutes, demonstrating linear scalability and confirming the method’s practicality for problems that are otherwise intractable with conventional techniques.

The authors conclude that the combination of BSDE‑PDE duality, Picard‑Monte‑Carlo iteration, and a robust polynomial regression operator yields a powerful tool for high‑dimensional option pricing. Future work is suggested on GPU acceleration, adaptive sparse polynomial bases, and rigorous convergence analysis for non‑Lipschitz drivers, which could further enhance performance and broaden applicability.