A backward differential deep learning-based algorithm for solving high-dimensional nonlinear backward stochastic differential equations

Abstract

1. Introduction

In this paper we are concerned with the numerical solution of the decoupled forward-backward stochastic differential equation (FBSDE) of the form

where |$W_{t} = \left ( W_{t}^{1}, \ldots , W_{t}^{d} \right )^\top $| is a |$d$|-dimensional Brownian motion, |$a: [0, T] \times \mathbb{R}^{d} \to \mathbb{R}^{d}$|⁠, |$b: [0, T] \times \mathbb{R}^{d} \to \mathbb{R}^{d \times d}$|⁠, |$f:\left [0,T\right ]\times \mathbb{R}^{d}\times \mathbb{R}\times \mathbb{R}^{1\times d} \to \mathbb{R}$| is the driver function and |$g: \mathbb{R}^{d} \to \mathbb{R}$| is the terminal condition that depends on the final value |$X_{T}$| of the forward stochastic differential equation (SDE). Hence, the randomness in the backward stochastic differential equation (BSDE) is driven by the forward SDE. Usually, the coupled FBSDE is referred to as a FBSDE. Hence, to avoid confusion, we refer to the decoupled FBSDE (1.1) as a BSDE. We shall work under the standard well-posedness assumptions of Pardoux & Peng (1990) to ensure the existence of a unique solution pair of (1.1).

The main motivation for studying BSDEs lies in their significance as essential tools for modeling problems across various scientific domains, including finance, economics, physics, etc., due to their connection to partial differential equations (PDEs) through the well-known (nonlinear) Feynman–Kac formula. As an illustrative example of their applications in finance it was demonstrated in Karoui et al. (1997) that the price and delta hedging of an option can be represented by a BSDE. Such an approach via a BSDE has a couple of advantages when compared with the usual one of considering the associated PDE. First, the delta hedging strategy is inclusive in the BSDE solution. Secondly, many market models can be presented in terms of BSDEs, ranging from the Black–Scholes model to more advanced ones such as local volatility models (Labart & Lelong, 2011), stochastic volatility models (Fahim et al., 2011), jump-diffusion models (Eyraud-Loisel, 2005), defaultable options (Ankirchner et al., 2010) and many others. Thirdly, BSDEs can also be used in incomplete markets (Karoui et al., 1997). Furthermore, using BSDEs eliminates the need to switch to the so-called risk-neutral measure. Therefore, BSDEs represent a more intuitive and understandable approach for option pricing and hedging.

Under the Black–Scholes framework such a BSDE is linear and the solution is given in a closed form. However, in most practical scenarios, BSDEs cannot be explicitly solved. For instance, the Black–Scholes model under different interest rates for lending and borrowing (Bergman, 1995) leads to a nonlinear BSDE for which finding an analytical solution becomes challenging. Hence, advanced numerical techniques to approximate their solutions become desired. In recent years, various numerical methods have been proposed for solving BSDEs, e.g. (Bouchard & Touzi, 2004; Zhang, 2004; Gobet et al., 2005; Lemor et al., 2006; Zhao et al., 2006; Bender & Zhang, 2008; Ma et al., 2008; Gobet & Labart, 2010; Zhao et al., 2010; Crisan & Manolarakis, 2012; Zhao et al., 2014; Ruijter & Oosterlee, 2015, 2016; Teng et al., 2020; Teng & Zhao, 2021) and many others. However, most of them are not suitable for tackling high-dimensional BSDEs due to the well-recognized challenge known as the ‘curse of dimensionality’. The computational cost associated with solving high-dimensional BSDEs grows exponentially with the increase in dimensionality. Some of the most important equations are naturally formulated in high dimensions. For instance, the Black–Scholes equation for option pricing exhibits the dimensionality of the BSDE with the number of underlying financial assets under consideration. Some techniques such as parallel computing using GPU computing (Gobet et al., 2016; Kapllani & Teng, 2022) or sparse grid methods (Zhang, 2013; Fu et al., 2017; Chassagneux et al., 2023) have proven effective in solving only moderately dimensional BSDEs within reasonable computation time.

In recent years machine learning models have demonstrated remarkable success in the field of artificial intelligence, inspiring applications in other domains where the curse of dimensionality has been a persistent challenge. Consequently, different approaches using machine learning have been proposed to solve high-dimensional BSDEs: the deep learning-based methods using deep neural networks (DNNs) and the regression tree-based methods (Teng, 2021, 2022). The first deep learning-based scheme called the deep BSDE (we refer to it as the DBSDE scheme) was introduced in E et al. (2017); Han et al. (2018). The authors conducted numerical experiments with various examples, demonstrating the effectiveness of their proposed algorithm in high-dimensional settings. It proved proficient in delivering both accurate approximations of the solution and computational efficiency. Therefore, the method opened the door to solving BSDEs in hundreds of dimensions in a reasonable amount of time. Several articles have been published after the original publication of the DBSDE method, some adjusting, reformulating or extending the algorithm (Fujii et al., 2019; Huré et al., 2020; Ji et al., 2020; Kremsner et al., 2020; Beck et al., 2021; Chen & Wan, 2021; Ji et al., 2021; Liang et al., 2021; Pham et al., 2021; Abbas-Turki et al., 2022; Germain et al., 2022; Gnoatto et al., 2022; Ji et al., 2022; Takahashi et al., 2022; Andersson et al., 2023; Gnoatto et al., 2023; Kapllani & Teng, 2024; Negyesi et al., 2024; Raissi, 2024), while others focused on error analysis (Han & Long, 2020; Jiang & Li, 2021) and uncertainty quantification (Kapllani et al., 2025). It has been pointed out in the literature that the DBSDE method suffers from different issues such as convergence to an approximation far from the solution or even divergence when the problem has a complex structure and a long terminal time. To tackle these drawbacks many alternative methods have been proposed; we refer to, e.g., (Huré et al., 2020; Teng, 2022; Andersson et al., 2023; Chassagneux et al., 2023; Kapllani & Teng, 2024). High-accurate gradient approximations are of great significance, especially in financial applications, where the process |$Z$| represents the hedging strategy for an option contract. Except the works in Kapllani & Teng (2024); Negyesi et al. (2024) other deep learning-based schemes do not discuss in detail the approximations for |$Z$| in high-dimensional spaces, as they are generally more challenging than approximating |$Y$| for BSDEs. In this work we develop a novel algorithm that ensures high accuracy, not only for the process |$Y$|⁠, but also for the process |$Z$|⁠.

The authors in Huré et al. (2020) approximate the unknown solution pair of (1.1) using DNNs. The network parameters are optimized at each time step through the minimization of loss functions defined recursively via backward induction. More precisely, the loss is formulated from the Euler–Maruyama discretization of the BSDE at each time interval. The method is referred to as the deep backward dynamic programming (DBDP) scheme. Such formulation gives an implicit approximation of the process |$Z$|⁠. Hence, the stochastic gradient descent (SGD) algorithm lacks explicit information about |$Z$|⁠, which impacts its approximation accuracy. To address this we enhance the SGD algorithm by providing it with additional information to achieve accurate approximations of |$Z$|⁠. We make use of differential deep learning (Huge & Savine, 2020), a general extension of supervised deep learning. In this framework the DNN model is trained, not only on inputs and labels, but also on differentials of labels with respect to (w.r.t.) inputs. Differential deep learning offers an efficient approximation, not only of the labels, but also of their derivatives when compared with traditional supervised deep learning. We use Malliavin calculus to formulate the BSDE problem as a differential deep learning problem. The Malliavin derivatives of the BSDE solution pair |$(Y, Z)$| themselves satisfy another BSDE, resulting in a system of BSDEs. This formulation also requires estimating the Hessian matrix of the solution. In the context of option pricing this matrix corresponds to |$\varGamma $| sensitivity, which can be used to indicate a potential acceleration in changes in the option’s value.

Our method works as follows. First, we discretize the system of BSDEs using the Euler–Maruyama method. Subsequently, we utilize DNNs to approximate the unknown solution of these BSDEs, requiring the estimation of the triple of the processes |$\left (Y, Z, \varGamma \right )$|⁠. The network parameters are optimized backwardly at each time step by minimizing a loss function defined as a weighted sum of the dynamics of the discretized BSDE system. Through this way SGD is equipped with explicit information about the dynamics of the process |$Z$|⁠. As a result our method can yield more accurate approximations than the scheme proposed in Huré et al. (2020), not only for the process |$Z$|⁠, but also for the process |$\varGamma $|⁠. The computation time of our scheme is higher compared with that of Huré et al. (2020). Note that the authors in Negyesi et al. (2024) also used the Malliavin derivative to improve the accuracy of |$Z$|⁠. However, their method significantly differs from ours, as they only employ supervised deep learning. Their approach requires training the BSDE system separately, which gives a higher computational cost compared with our method. This is demonstrated in our numerical experiments. Furthermore, our approach using differential deep learning can be straightforwardly extended, not only to Huré et al. (2020), which operates backward in time through local optimization at each discrete time step, but also to other deep learning-based schemes (E et al., 2017; Kapllani & Teng, 2024; Raissi, 2024) formulated forward in time as a global optimization problem (this is part of our ongoing research). In contrast, the scheme presented in Negyesi et al. (2024) cannot be integrated into such methodologies, as it cannot be formulated as a global optimization problem. To the best of our knowledge only Lefebvre et al. (2023) apply differential deep learning to solve high-dimensional PDEs, where the authors consider the associated dual stochastic control problem instead of working with BSDEs.

The outline of the paper is organized as follows. In the next section we recall some of the well-known results concerning BSDEs. In Section 3 DNNs and differential deep learning techniques are described. Our backward differential deep learning-based algorithm is presented in Section 4. Section 5 is devoted to the convergence analysis of our algorithm. The numerical experiments presented in Section 6 confirm the theoretical results and show high accuracy of the solution, its gradient and the Hessian matrix of the solution over different option pricing problems. Finally, Section 7 concludes this work.

2. Preliminaries

2.1 Spaces and notation

Let |$\left (\varOmega ,\mathscr{F},\mathbb{P},\{\mathscr{F}_{t}\}_{0\le t \le T}\right )$| be a complete, filtered probability space. In this space a standard |$d$|-dimensional Brownian motion |$\{W_{t}\}_{\leq t \leq T}$| is defined, such that the filtration |$\{\mathscr{F}_{t}\}_{0\le t\le T}$| is the natural filtration of |$W_{t}.$| As usual, we identify random variables that are equal |$\mathbb{P}$|-a.s. and, accordingly, understand equalities and inequalities between them in the |$\mathbb{P}$|-a.s. sense. For the expectation we omit the superscript |$\mathbb{P}$| if it is meant under probability measure |$\mathbb{P}$| (unless stated otherwise). We denote further

• |$x \in \mathbb{R}^{d}$| as a column vector. |$x \in \mathbb{R}^{1\times d}$| as a row vector.

• |$| x |$| for the Frobenius norm of any |$x \in \mathbb{R}^{d \times{\mathfrak{q}}}$|⁠. In the case of scalar and vector inputs these coincide with the standard Euclidian norm.

• |$\mathbb{S}^{2}\left ([0, T] \times \varOmega ; \mathbb{R}^{d\times{\mathfrak{q}}} \right )$| for the space of continuous and progressively measurable stochastic processes |$X: [0, T] \times \varOmega \to \mathbb{R}^{d\times{\mathfrak{q}}} $| such that |$\mathbb{E}\bigl [\sup _{ 0 \leq t\leq T}\left |X_{t}\right |^{2}\bigr ] < \infty $|⁠.

• |$\mathbb{H}^{2}\left ([0, T] \times \varOmega ; \mathbb{R}^{d\times{\mathfrak{q}}} \right )$| for the space of progressively measurable stochastic processes |$Z: [0, T] \times \varOmega \to \mathbb{R}^{d\times{\mathfrak{q}}} $| such that |$\mathbb{E}\left [ \int _{0}^{T} \left |Z_{t}\right |^{2} \, {\text{d}}t \right ] < \infty $|⁠.

• |$\mathbb{L}^{2}_{\mathscr{F}_{t}}\left (\varOmega ; \mathbb{R}^{d \times{\mathfrak{q}}} \right )$| for the space of |$\mathscr{F}_{t}$|-measurable random variable |$\xi : \varOmega \to \mathbb{R}^{d \times{\mathfrak{q}}} $| such that |$\mathbb{E}\bigl [\left |\xi \right |^{2} \bigr ] < \infty $|⁠.

• |$L^{2}\left ( [0, T] ; \mathbb{R}^{{\mathfrak{q}}}\right )$| for the Hilbert space of deterministic functions |$h: [0, T] \to \mathbb{R}^{{\mathfrak{q}}}$| such that |$\int _{0}^{T} \left | h\left (t\right ) \right |^{2} {\text{d}}t < \infty $|⁠.

• |$\nabla _{x} f:= \left ( \frac{\partial f}{\partial x_{1}}, \ldots , \frac{\partial f}{\partial x_{d}} \right ) \in \mathbb{R}^{1 \times d}$| for the gradient of scalar-valued multivariate function |$f\left (t, x, y, z\right )$| w.r.t. |$x \in \mathbb{R}^{d}$|⁠, and analogously for |$\nabla _{y} f \in \mathbb{R}$| and |$\nabla _{z} f \in \mathbb{R}^{1 \times d}$| w.r.t. |$y \in \mathbb{R}$| and |$z \in \mathbb{R}^{1 \times d}$|⁠, respectively. Similarly, we denote the Jacobian matrix of a vector-valued function |$u: \mathbb{R}^{d} \to \mathbb{R}^{{\mathfrak{q}}}$| by |$\nabla _{x} u \in \mathbb{R}^{{\mathfrak{q}} \times d}$|⁠.

• |$\operatorname{Hess}_{x} u \in \mathbb{R}^{d \times d}$| the Hessian matrix of a function |$u: \mathbb{R}^{d} \to \mathbb{R}$|⁠.

• |$C^{{\mathfrak{l}}}_{{\mathfrak{b}}}\left ( \mathbb{R}^{d}; \mathbb{R}^{{\mathfrak{q}}} \right )$| and |$C^{{\mathfrak{l}}}_{{\mathfrak{p}}}\left ( \mathbb{R}^{d}; \mathbb{R}^{{\mathfrak{q}}} \right )$| for the set of |${\mathfrak{l}}$|-times continuously differentiable functions |$\varphi : \mathbb{R}^{d} \to \mathbb{R}^{{\mathfrak{q}}}$| such that all partial derivatives up to order |${\mathfrak{l}}$| are bounded or have polynomial growth, respectively.

• |$\varDelta = \{t_{0}, t_{1}, \ldots , t_{N}\}$| is the time discretization of |$[0, T]$| with |$t_{0} = 0 < t_{1} < \ldots < t_{N} = T$|⁠, |$\varDelta t_{n} = t_{n+1} - t_{n}$| and |$| \varDelta |:= \max _{ 0\leq n \leq N-1 } t_{n+1} - t_{n}$|⁠.

• |$\mathbb{E}_{n}\left [ Y \right ]:=\mathbb{E}\left [ Y | \mathscr{F}_{t_{n}} \right ]$| for the conditional expectation w.r.t. the natural filtration, given the time partition |$\varDelta $|⁠.

• |$x^{\top }\in \mathbb{R}^{{\mathfrak{q}} \times d}$| for the transpose of any |$x \in \mathbb{R}^{d \times{\mathfrak{q}}}$|⁠.

• |$\operatorname{Tr}\left [x\right ]$| for the trace of any |$x \in \mathbb{R}^{d \times d}$|⁠.

• |$\mathbf{0}_{d,d}$|⁠, |$\mathbf{1}_{d,d}$| for |$\mathbb{R}^{d \times d}$| matrices of all zeros and ones, respectively.

2.2 Malliavin calculus

We shall use techniques of the stochastic calculus of variations. To this end we use the following notation. For more details we refer the reader to Nualart (2006). Let |$\mathscr{S}$| be the space of smooth random variables of the form

where |$\varphi \in C^{\infty }_{{\mathfrak{p}}}\left (\mathbb{R}^{d};\mathbb{R}\right )$|⁠, |$h_{1}, \ldots , h_{d} \in L^{2}\left ( [0, T] ; \mathbb{R}^{{\mathfrak{q}}}\right )$|⁠. The Malliavin derivative of smooth random variable |$\xi \in \mathscr{S}$| is the |$\mathbb{R}^{1 \times{\mathfrak{q}}}$|-valued stochastic process given by

We define the domain of |$D$| in |$\mathbb{L}^{2}_{\mathscr{F}_{T}}$| as |$\mathbb{D}^{1,2}\left ( \varOmega ;\mathbb{R} \right )$|⁠, meaning that |$\mathbb{D}^{1,2}$| is the closure of the class of smooth random variables |$\mathscr{S}$| w.r.t. the norm

Note that in case of vector-valued Malliavin differentiable random variables |$\xi = \left ( \xi _{1}, \ldots , \xi _{{\mathfrak{q}}} \right )$|⁠, |$\xi \in \mathbb{D}^{1,2}\left (\varOmega ; \mathbb{R}^{{\mathfrak{q}}} \right )$|⁠, its Malliavin derivative |$D_{s} \xi \in \mathbb{R}^{{\mathfrak{q}} \times{\mathfrak{q}}}$| is the matrix-valued stochastic process.

The following lemma represents the Malliavin chain rule, which can be extended to Lipschitz continuous functions.

2.3 Some results on BSDEs

We recall some results on BSDE known from the literature that are relevant for this work. For the functions in BSDE (1.1) we hierarchically structure the properties that they are assumed to fulfil.

AX1. The initial condition |$x_{0} \in \mathbb{L}^{2}_{\mathscr{F}_{0}}\left (\varOmega ; \mathbb{R}^{d}\right )$| and |$a, b$| satisfy a linear growth condition in |$x$|⁠, i.e.,

|$\forall \, t \in [0, T], x \in \mathbb{R}^{d}$| and some constant |$C>0$|⁠. Furthermore, |$a, b$| are uniformly Lipschitz continuous in the spatial variable, i.e.,

|$\forall \, t \in [0, T], x_{1}, x_{2} \in \mathbb{R}^{d}$|⁠, for some constant |$L_{a,b}>0$|⁠.

AX2. Assumption AX1 holds. Moreover, |$a(t, 0)$|⁠, |$b(t, 0)$| are uniformly bounded |$\forall $|  |$0 \leq t \leq T$| and |$a \in C_{{\mathfrak{b}}}^{0, 1}\left ( [0, T] \times \mathbb{R}^{d}; \mathbb{R}^{d} \right )$|⁠, |$b \in C_{{\mathfrak{b}}}^{0, 1}\left ( [0, T] \times \mathbb{R}^{d}; \mathbb{R}^{d \times d} \right )$|⁠.

AX3. Assumption AX2 holds. Moreover, |$a \in C_{{\mathfrak{b}}}^{0, 2}\left ( [0, T] \times \mathbb{R}^{d}; \mathbb{R}^{d} \right )$|⁠, |$b \in C_{{\mathfrak{b}}}^{0, 2}\left ( [0, T] \times \mathbb{R}^{d}; \mathbb{R}^{d \times d} \right )$| and there exist a positive constant |$C>0$| such that

AY1. The function |$f(t,x,y,z)$| is uniformly Lipschitz continuous w.r.t. |$y$| and |$z$|⁠, i.e.,

|$\forall \, (t, x, y_{1}, z_{1})$| and |$(t, x, y_{2}, z_{2}) \in [0, T] \times \mathbb{R}^{d} \times \mathbb{R} \times \mathbb{R}^{1\times d}$|⁠, for some constant |$L_{f}>0$|⁠. Moreover, |$f, g$| satisfy a quadratic growth condition in |$x$|⁠, i.e.,

|$\forall \, (t, x, y, z) \in [0, T] \times \mathbb{R}^{d} \times \mathbb{R} \times \mathbb{R}^{1\times d}$| for some constant |$C>0$|⁠.

AY2. Assumption AY1 holds. Moreover, |$f \in C_{{\mathfrak{b}}}^{0,1,1,1}\left ( [0, T] \times \mathbb{R}^{d} \times \mathbb{R} \times \mathbb{R}^{1 \times d}; \mathbb{R}\right )$| and |$g \in C_{{\mathfrak{b}}}^{1}\left ( \mathbb{R}^{d}; \mathbb{R} \right )$|⁠.

AY3. Assumption AY2 holds. Moreover, |$f \in C_{{\mathfrak{b}}}^{0,2,2,2}\left ( [0, T] \times \mathbb{R}^{d} \times \mathbb{R} \times \mathbb{R}^{1 \times d}; \mathbb{R}\right )$| and |$g \in C_{{\mathfrak{b}}}^{2}\left ( \mathbb{R}^{d}; \mathbb{R} \right )$|⁠.

In the following theorem we state the well-known result on SDEs.

The well-posedness of the BSDE (1.1) is guaranteed by Assumption AY1. The following theorem guarantees the existence of a unique solution triple of (1.1).

An important property of BSDEs is that they provide a probabilistic representation for the solution of a specific class of PDEs given by the nonlinear Feynman–Kac formula. Consider the semilinear parabolic PDE

for all |$(t, x) \in ([0, T]\times \mathbb{R}^{d})$| and the terminal condition |$u(T,x)=g(x)$|⁠. Assume that (2.1) has a classical solution |$u \in C^{1,2}_{{\mathfrak{b}}}\left ([0, T] \times \mathbb{R}^{d}; \mathbb{R}\right )$| and the aforementioned standard Lipschitz assumptions of (1.1) are satisfied. Then the solution of (1.1) can be represented |$\mathbb{P}$|-a.s. by

Next, we collect some Malliavin differentiability results on BSDEs, as we are interested in BSDEs such that their solution triple |$\{X_{t}, Y_{t}, Z_{t} \}_{0\leq t \leq T}$| is differentiable in the Malliavin sense. The results are stated in the following theorems.

The final important result that is relevant for this work is the path regularity result of the processes |$Y$| and |$Z$|⁠, which we state in the following theorem.

3. Differential deep learning

In this section we discuss differential machine learning in the context of DNNs, specifically differential deep learning, which plays a crucial role in formulating our algorithm. We start by describing DNNs, which are designed to approximate unknown or a large class of functions.

3.1 Deep neural networks

Let |$d_{0}, d_{1}\in \mathbb{N}$| be the input and output dimensions, respectively. We fix the global number of layers as |$L+2$|⁠, |$L \in \mathbb{N}$| the number of hidden layers each with |$\eta \in \mathbb{N}$| neurons. The first layer is the input layer with |$d_{0}$| neurons and the last layer is the output layer with |$d_{1}$| neurons. A DNN is a function |$\phi (\cdot ; \theta ): \mathbb{R}^{d_{0}} \to \mathbb{R}^{d_{1}}$| composed of a sequence of simple functions, which can be expressed in the following form:

where |$\theta :=\left ( \theta (1), \ldots , \theta (L+1) \right ) \in \mathbb{R}^{P}$| and |$P$| is the total number of network parameters; |$x \in \mathbb{R}^{d_{0}}$| is called an input vector. Moreover, |$A_{l}(\cdot ; \theta (l)), l = 1, 2, \ldots , L+1$| are affine transformations: |$A_{1}(\cdot ;\theta (1)): \mathbb{R}^{d_{0}} \to \mathbb{R}^{\eta }$|⁠, |$A_{l}(\cdot ;\theta (l)), l = 2, \ldots , L: \mathbb{R}^{\eta } \to \mathbb{R}^{\eta }$| and |$A_{L+1}(\cdot ;\theta (L+1)): \mathbb{R}^{\eta } \to \mathbb{R}^{d_{1}}$|⁠, represented by

where |$\theta (l):=\left (\mathscr{W}_{l}, \mathscr{B}_{l}\right )$|⁠, |$\mathscr{W}_{l} \in \mathbb{R}^{\eta _{l} \times \eta _{l-1}}$| is the weight matrix and |$\mathscr{B}_{l} \in \mathbb{R}^{\eta _{l}}$| is the bias vector with |$\eta _{0} = d_{0}, \eta _{L+1} = d_{1}, \eta _{l} = \eta $| for |$l = 1, \ldots , L$| and |$\varrho : \mathbb{R} \to \mathbb{R}$| is a nonlinear function (called the activation function), and applied component-wise on the outputs of |$A_{l}(\cdot ;\theta (l))$|⁠. Common choices are |$\tanh (\cdot ), \sin (\cdot ), \max (0,\cdot )$|⁠, etc. All these matrices |$\mathscr{W}_{l}$| and vectors |$\mathscr{B}_{l}$| form the parameters |$\theta $| of the DNN and they have the dimension

for fixed |$d_{0}, d_{1}, L$| and |$\eta $|⁠. We denote by |$\varTheta $| the set of possible parameters for the DNN |$\phi (\cdot ; \theta )$| with |$\theta \in \varTheta $|⁠. The Universal Approximation Theorem (UAT) (Cybenko, 1989; Hornik et al., 1989) justifies the use of DNNs as function approximators.

3.2 Training of DNNs using supervised deep learning

Once the DNN architecture is defined what determines the mapping of a certain input to an output are the parameters |$\theta $| incorporated in the DNN model. These parameters need to be optimized such that the DNN approximates the unknown function, which is called the training of the DNN. The loss function acts as the objective function to be minimized during the training procedure, in which the DNNs optimal set of parameters is searched.

Consider the training data sampled from some (unknown) multivariate joint distribution |$(\mathscr{X}, \mathscr{Y}) \sim \mathscr{P}$|⁠, where the random variable |$\mathscr{X} \in \mathbb{R}^{d}$| is referred to as the input and the random variable |$\mathscr{Y} \in \mathbb{R}$| as the label. The goal (in a regression setting) is then to approximate the deterministic function |$F(x):= \mathbb{E}^{\mathscr{P}} \left [ \mathscr{Y} |\mathscr{X} = x \right ]$| by DNN |$\phi (\mathscr{X};\theta )$| using |$(\mathscr{X}, \mathscr{Y}) \sim \mathscr{P}$|⁠. The loss function measures how well the current approximation of the DNN is compared with the label. A common choice is the expected squared error, which is given as

Then, the optimal parameters |$\theta ^{*}$| in (3.1) are given as

which can be estimated by using SGD-type algorithms.

3.3 Training of DNNs using differential deep learning

One of the biggest challenges w.r.t. finding the optimal parameter set of the DNN is to avoid learning training data-specific patterns, namely overfitting, and rather enforcing better generalization of the fitted models. Hence, regularization approaches have been developed for DNNs to avoid overfitting and thus improve the performance of the model. Such approaches penalize certain norms of the parameters |$\theta $|⁠, expressing a preference for |$\theta $|⁠. Differential deep learning (Huge & Savine, 2020) has the same motivation as regularization, namely to improve the accuracy of the model. This is achieved by not expressing a preference, but correctness, in particular enforcing differential correctness. It assumes that the derivative of the label w.r.t. input is known. Let us consider the function |$F_{x}(x)=\nabla _{x} F(x)$| and the random variable |$\mathscr{Z}:= F_{x}(\mathscr{X}) \in \mathbb{R}^{1\times d}$|⁠. The goal in differential deep learning is to approximate the label function |$F(x)$| by DNN |$\phi (\mathscr{X};\theta )$| using data |$( \mathscr{X}, \mathscr{Y}, \mathscr{Z}) \sim \mathscr{P}$| and minimizing the extended loss function (3.1) given as

where |$\nabla _{x}\phi $| is calculated using AD and |$\lambda \in \mathbb{R}_{+}$|⁠. Our numerical experiments indicated that approximating the derivatives using AD resulted in worse performance compared with utilizing a separate DNN. This is consistent with the results in Huré et al. (2020). Therefore, we chose to employ a separate DNN for the derivatives, namely we consider a slightly different formulation of differential deep learning compared with Huge & Savine (2020). We use one DNN |$\phi ^{y}\left ( \mathscr{X} ; \theta ^{y}\right )$| to approximate the function |$F(x)$| and another |$\phi ^{z}\left ( \mathscr{X} ; \theta ^{z}\right )$| for |$F_{x}(x)$|⁠, and rewrite the loss function (3.2) as

where |$\theta = \left ( \theta ^{y}, \theta ^{z} \right )$|⁠, |$\omega _{1}, \omega _{2} \in [0, 1]$| and |$\omega _{1} +\omega _{2} = 1$|⁠. Then, the optimal parameters |$\theta ^{*}$| in (3.3) are given as

estimated using an SGD method. Since the derivatives are integrated in the loss function (3.3) as an additional term we consider this modification to remain within the framework of differential deep learning.

4. A backward differential deep learning-based scheme for BSDEs

In this section we introduce the proposed backward differential deep learning-based method. In order to formulate BSDE as a differential learning problem we firstly discretize the integrals in the resulting BSDE system given as

where we introduced the notations |$\mathbf{X}_{t}:= \left ( X_{t}, Y_{t}, Z_{t}\right )$|⁠, |$\mathbf{D}_{s}\mathbf{X}_{t}:= \left ( D_{s} X_{t}, D_{s} Y_{t}, D_{s} Z_{t}\right )$| and |$f_{D}\left (t, \mathbf{X}_{t}, \mathbf{D}_{s}\mathbf{X}_{t} \right ):= \nabla _{x} f\left ( t, \mathbf{X}_{t} \right ) D_{s} X_{t} + \nabla _{y} f\left ( t, \mathbf{X}_{t}\right ) D_{s} Y_{t} + f\left ( t, \mathbf{X}_{t}\right ) D_{s} Z_{t}$|  |$\forall \, 0 \leq s, t \leq T$|⁠. Note that the solution of the above BSDE system is a pair of triples of stochastic processes |$\left\{\left (X_{t}, Y_{t}, Z_{t}\right )\right\}_{0\leq t \leq T}$| and |$\left\{\left (D_{s} X_{t}, D_{s} Y_{t}, D_{s} Z_{t}\right )\right\}_{0\leq s, t \leq T}$| such that (4.1)–(4.4) hold |$\mathbb{P}$|-a.s.

Let us consider the time discretization |$\varDelta $|⁠. For notational convenience we write |$\varDelta W_{n} = W_{t_{n+1}} - W_{t_{n}}$|⁠, |$(X_{n}, Y_{n}, Z_{n}) = (X_{t_{n}}, Y_{t_{n}}, Z_{t_{n}})$|⁠, |$(D_{n} X_{m}, D_{n} Y_{m}, D_{n} Z_{m}) = (D_{t_{n}} X_{t_{m}}, D_{t_{n}} Y_{t_{m}}, D_{t_{n}} Z_{t_{m}})$|⁠, and |$\left (X^{\varDelta }_{n}, Y^{\varDelta }_{n}, Z^{\varDelta }_{n}\right )$|⁠, |$\left (D_{n}X^{\varDelta }_{m}, D_{n}Y^{\varDelta }_{m}, D_{n} Z^{\varDelta }_{m}\right )$| for the approximations, where |$0 \leq n, m \leq N$|⁠. The forward SDE (4.1) is approximated by the Euler–Maruyama scheme, i.e.,

for |$n = 0, 1, \ldots , N-1$|⁠, where |$X^{\varDelta }_{0} = x_{0}$|⁠.

Next, we apply the Euler–Maruyama scheme to (4.2). For the time interval |$\left [t_{n}, t_{n+1}\right ]$| we have

Applying the Euler–Maruyama scheme in (4.6) one obtains

for |$n = N-1, N-2, \ldots , 0,$| where |$\mathbf{X}_{n}^{\varDelta }:= \left ( X_{n}^{\varDelta }, Y_{n}^{\varDelta }, Z_{n}^{\varDelta }\right )$| and |$Y^{\varDelta }_{N} = g\left (X^{\varDelta }_{N}\right )$|⁠.

Next, we discretize the BSDE for the Malliavin derivatives, i.e., (4.3)–(4.4) in a similar manner. The Malliavin derivative (4.3) approximated by the Euler–Maruyama method gives the estimates

From |$\left [t_{n}, t_{n+1}\right ]$| (4.4) is given as

Using Euler–Maruyama scheme in (4.9) we get

with |$\mathbf{D}_{n}\mathbf{X}_{n}^{\varDelta }:= \left ( D_{n} X_{n}^{\varDelta }, D_{n} Y_{n}^{\varDelta }, D_{n} Z_{n}^{\varDelta }\right )$|⁠. Given the Markov property of the underlying processes the Malliavin chain rule (Lemma 2.1) implies that

for some deterministic functions |$y\!: [0, T] \times \mathbb{R}^{d} \to \mathbb{R}$| and |$z\!: [0, T] \times \mathbb{R}^{d} \to \mathbb{R}^{1 \times d}$|⁠, where |$\gamma \!: [0, T] \times \mathbb{R}^{d} \to \mathbb{R}^{d \times d} $| is the Jacobian matrix of |$z\left (t_{m}, X_{m}\right )$|⁠. Note that from the Feynman–Kac relation (2.2) we have that |$z\left (t_{m}, X_{m}\right ) = \nabla _{x} y\left (t_{m}, X_{m}\right ) b\left (t_{m}, X_{m}\right )$|⁠. Hence, one can write that |$D_{n} Y_{m} = z\left (t_{m}, X_{m}\right ) b^{-1} \left (t_{m}, X_{m}\right ) D_{n} X_{m}$|⁠. Using Theorem 2.4 we have that (4.11) becomes

where due to the aforementioned relations |$f_{D}\left (t, \mathbf{X}_{n}^{\varDelta }, \mathbf{D}_{n} \mathbf{X}_{n}^{\varDelta } \right )= \nabla _{x} f\left ( t_{n}, \mathbf{X}_{n}^{\varDelta } \right ) D_{n} X_{n}^{\varDelta } + \nabla _{y} f\left ( t_{n}, \mathbf{X}_{n}^{\varDelta }\right ) Z_{n}^{\varDelta } + \nabla _{z} f\left ( t_{n}, \mathbf{X}_{n}^{\varDelta }\right ) \varGamma _{n}^{\varDelta } D_{n} X_{n}^{\varDelta }.$|

After discretizing the integrals our scheme is made fully implementable at each discrete time point |$t_{n}$| by an appropriate function approximator to estimate the discrete unknown processes |$\left (Y^{\varDelta }_{n}, Z^{\varDelta }_{n}, \varGamma ^{\varDelta }_{n}\right )$| in (4.7) and (4.12). We estimate these unknown processes using DNNs and propose the following scheme:

• Generate approximations |$X^{\varDelta }_{n+1}$| for |$n = 0, 1, \ldots , N-1$| of SDE (4.1) via (4.5) and its discrete Malliavin derivative |$D_{n} X_{n}^{\varDelta }$|⁠, |$D_{n} X_{n+1}^{\varDelta }$| using (4.8).

• Set

  $$ \begin{align*} &Y_N^{\varDelta, \hat{\theta}}:= g(X_N^{\varDelta}), \quad Z_N^{\varDelta, \hat{\theta}}:= \nabla_x g(X_N^{\varDelta}) b(t_N, X_N^{\varDelta}), \quad \varGamma_N^{\varDelta, \hat{\theta}}:= \left[\nabla_x (\nabla_x g\, b)\right](t_N, X_N^{\varDelta}).\end{align*} $$
• For each discrete time point |$t_{n}$|⁠, |$n = N-1, N-2, \ldots , 0$| we use three independent DNNs |$\phi ^{y}_{n}(\cdot ; \theta ^{y}_{n}): \mathbb{R}^{d} \to \mathbb{R}$|⁠, |$\phi ^{z}_{n}(\cdot ; \theta ^{z}_{n}): \mathbb{R}^{d} \to \mathbb{R}^{1 \times d}$| and |$\phi ^{\gamma }_{n}(\cdot ; \theta ^{\gamma }_{n}): \mathbb{R}^{d} \to \mathbb{R}^{d \times d}$| to approximate the discrete processes |$\left (Y_{n}^{\varDelta }, Z_{n}^{\varDelta }, \varGamma _{n}^{\varDelta }\right )$|⁠, respectively. Train the parameter set |$\theta _{n} = \left ( \theta ^{y}_{n}, \theta ^{z}_{n}, \theta ^{\gamma }_{n}\right )$| using the differential learning approach by constructing a loss function—as in (3.3)—such that the dynamics of the discretized process |$Y$| and |$Z$| given by (4.7) and (4.12) are fulfilled, namely

  $$ \begin{align} \mathbf{L}_{n}^{\varDelta}\left( \theta_{n} \right) &:= \omega_{1} \mathbf{L}^{y,\varDelta}_{n}\left( \theta_{n} \right) + \omega_{2} \mathbf{L}^{z,\varDelta}_{n}\left( \theta_{n} \right)\!, \nonumber \\ \mathbf{L}^{y,\varDelta}_{n}\left( \theta_{n} \right)&:=\mathbb{E}\left[ \left\vert Y^{\varDelta, \hat{\theta}}_{n+1} - \phi^{y}_{n}\left( X^{\varDelta}_{n}; \theta^{y}_{n} \right) + f\left(t_{n}, \mathbf{X}^{\varDelta, \theta}_{n}\right) \varDelta t_{n} - \phi^{z}_{n}\left( X^{\varDelta}_{n}; \theta^{z}_{n} \right) \varDelta W_{n} \right\vert^{2} \right]\!, \nonumber \\ \mathbf{L}^{z,\varDelta}_{n}\left( \theta_{n} \right) &:= \mathbb{E}\left[ \left\vert\vphantom{\left(\left(\phi^{\gamma}_{n}\left( X_{n}^{\varDelta}; \theta^{\gamma}_{n}\right) D_{n} X_{n}^{\varDelta}\right)^\top \varDelta W_{n}\right)^\top} Z^{\varDelta, \hat{\theta}}_{n+1} b^{-1}\left( t_{n+1}, X_{n+1}^{\varDelta} \right) D_{n} X^{\varDelta}_{n+1} - \phi^{z}_{n}\left( X^{\varDelta}_{n}; \theta^{z}_{n} \right) \right. \right. \nonumber \\ & \quad \left. \left. +\, f_{D}\left(t_{n}, \mathbf{X}^{\varDelta, \theta}_{n}, \mathbf{D}_{n}\mathbf{X}_{n}^{\varDelta, \theta}\right)\varDelta t_{n} - \left(\left(\phi^{\gamma}_{n}\left( X_{n}^{\varDelta}; \theta^{\gamma}_{n}\right) D_{n} X_{n}^{\varDelta}\right)^\top \varDelta W_{n}\right)^\top \right\vert^{2}\right]\!, \end{align} $$

  (4.13)

  where for notational convenience |$\mathbf{X}^{\varDelta , \theta }_{n}:=\left ( X_{n}^{\varDelta }, \phi ^{y}_{n}\left ( X^{\varDelta }_{n}; \theta ^{y}_{n} \right ), \phi ^{z}_{n}\left ( X^{\varDelta }_{n}; \theta ^{z}_{n} \right ) \right )$| and |$\mathbf{D}_{n}\mathbf{X}^{\varDelta , \theta }_{n}:=\left ( D_{n} X_{n}^{\varDelta }, \phi ^{z}_{n}\left ( X^{\varDelta }_{n}; \theta ^{y}_{n} \right ), \phi ^{\gamma }_{n}\left ( X^{\varDelta }_{n}; \theta ^{z}_{n} \right ) D_{n} X_{n}^{\varDelta } \right )$|⁠. We approximate the optimal parameters |$\theta ^{*}_{n} \in \operatorname{arg\,min}_{\theta _{n} \in \varTheta _{n}} \mathbf{L}^{\varDelta }_{n}\left ( \theta _{n} \right )$| using an SGD method and receive the estimated parameters |$\hat{\theta }_{n} = \left ( \hat{\theta }^{y}_{n}, \hat{\theta }^{z}_{n}, \hat{\theta }^{\gamma }_{n} \right )$|⁠. Then, we define

  $$ \begin{align*} &Y_n^{\varDelta, \hat{\theta}}:= \phi^y_n\left( X^{\varDelta}_n; \hat{\theta}^y_n \right), \quad Z_n^{\varDelta, \hat{\theta}}:= \phi^z_n\left( X^{\varDelta}_n; \hat{\theta}^z_n \right), \quad \varGamma_n^{\varDelta, \hat{\theta}}:= \phi^{\gamma}_n \left( X^{\varDelta}_n; \hat{\theta}^{\gamma}_n \right).\end{align*} $$

We refer to our scheme as differential learning backward dynamic programming (DLBDP) scheme, where |$\omega _{1} = \frac{1}{d+1}$| and |$\omega _{2} = \frac{d}{d+1}$| is considered due to dimensionality of the processes |$Y$| and |$Z$|⁠.

Note that the DBDP scheme from Huré et al. (2020) (specifically DBDP1) can be considered as a special case of our scheme by choosing |$\omega _{1} = 1$|⁠, |$\omega _{2} = 0$|⁠, and using AD for approximating the process |$\varGamma $|⁠. It can be formulated as follows:

• Generate approximations |$X^{\varDelta }_{n+1}$| for |$n = 0, 1, \ldots , N-1$| using (4.5).

• Set

  $$ \begin{align*} &Y_N^{\varDelta, \hat{\theta}} = g(X_N^{\varDelta}), \quad Z_N^{\varDelta, \hat{\theta}} = \nabla_x g(X_N^{\varDelta}) b(t_N, X_N^{\varDelta}), \quad \varGamma_N^{\varDelta, \hat{\theta}}= \left[\nabla_x (\nabla_x g\, b)\right](t_N, X_N^{\varDelta}).\end{align*} $$
• For each discrete time point |$t_{n}$|⁠, |$n = N-1, N-2, \ldots , 0$| we use two independent DNNs |$\phi ^{y}_{n}(\cdot ; \theta ^{y}_{n}): \mathbb{R}^{d} \to \mathbb{R}$| and |$\phi ^{z}_{n}(\cdot ; \theta ^{z}_{n}): \mathbb{R}^{d} \to \mathbb{R}^{1 \times d}$| to approximate the discrete processes |$\left (Y_{n}^{\varDelta }, Z_{n}^{\varDelta }\right )$|⁠, respectively. We train the parameter set |$\theta _{n} = \left ( \theta ^{y}_{n}, \theta ^{z}_{n}\right )$| by constructing a loss function such that the dynamics of the discretized process |$Y$| given by (4.7) are fulfilled, namely

  $$ \begin{align*}& \mathbf{L}^{y,\varDelta}_{n}\left( \theta_{n} \right) =\mathbb{E}\left[ \left\vert Y^{\varDelta, \hat{\theta}}_{n+1} - \phi^{y}_{n}\left( X^{\varDelta}_{n}; \theta^{y}_{n} \right) + f\left(t_{n}, \mathbf{X}^{\varDelta, \theta}_{n}\right) \varDelta t_{n} - \phi^{z}_{n}\left( X^{\varDelta}_{n}; \theta^{z}_{n} \right) \varDelta W_{n} \right\vert^{2} \right]. \end{align*} $$

  Approximate the optimal parameters |$\theta ^{*}_{n} \in \operatorname{arg\,min}_{\theta _{n} \in \varTheta _{n}} \mathbf{L}^{y,\varDelta }_{n}\left ( \theta _{n} \right )$| using an SGD method and receive the estimated parameters |$\hat{\theta }_{n} = \left ( \hat{\theta }^{y}_{n}, \hat{\theta }^{z}_{n}\right )$|⁠. Estimate the discrete process |$\varGamma _{n}^{\varDelta }$| using AD. Then,

  $$ \begin{align*} &Y_n^{\varDelta, \hat{\theta}} = \phi^y_n\left( X^{\varDelta}_n; \hat{\theta}^y_n \right), \quad Z_n^{\varDelta, \hat{\theta}} = \phi^z_n\left( X^{\varDelta}_n; \hat{\theta}^z_n \right), \quad \varGamma_n^{\varDelta, \hat{\theta}} = \nabla_x \phi^z_n(x;\hat{\theta}^z_n)\Bigr|_{x = X_n^{\varDelta}}.\end{align*} $$

Our scheme offers several advantages over the DBDP scheme and other well-known deep learning-based approaches (E et al., 2017; Germain et al., 2022; Kapllani & Teng, 2024; Raissi, 2024):

• (i)

  By explicitly incorporating the dynamics of the process |$Z$| via the BSDE (4.4) in the loss function (4.13) we enhance the accuracy of |$Z$| approximations through the SGD method.

• (ii)

  Additionally, the inclusion of the process |$\varGamma $| in the loss function through BSDE (4.4) allows for better estimation of |$\varGamma $| within the DLBDP scheme compared with the deep learning-based schemes, where AD is required for approximation of |$\varGamma $|⁠.

The scheme in Negyesi et al. (2024)—called the one-step Malliavin (OSM) scheme—also uses the Malliavin derivative to improve the accuracy of |$Z$| in the DBDP method. Hence, in the numerical experiments we compare our approach with both the schemes DBDP and OSM. The latter can be formulated as follows:

• Generate approximations |$X^{\varDelta }_{n+1}$| for |$n = 0, 1, \ldots , N-1$| of SDE (4.1) via (4.5) and its discrete Malliavin derivative |$D_{n} X_{n}^{\varDelta }$|⁠, |$D_{n} X_{n+1}^{\varDelta }$| using (4.8).

• Set

  $$ \begin{align*} &Y_N^{\varDelta, \hat{\theta}} = g(X_N^{\varDelta}), \quad Z_N^{\varDelta, \hat{\theta}} = \nabla_x g(X_N^{\varDelta}) b(t_N, X_N^{\varDelta}), \quad \varGamma_N^{\varDelta, \hat{\theta}} = \left[\nabla_x (\nabla_x g\, b)\right](t_N, X_N^{\varDelta}).\end{align*} $$
• For each discrete time point |$t_{n}$|⁠, |$n = N-1, N-2, \ldots , 0$| we consider two optimization problems. In the first one we use two independent DNNs |$\phi ^{z}_{n}(\cdot ; \theta ^{z}_{n}): \mathbb{R}^{d} \to \mathbb{R}^{1 \times d}$| and |$\phi ^{\gamma }_{n}(\cdot ; \theta ^{\gamma }_{n}): \mathbb{R}^{d} \to \mathbb{R}^{d \times d}$| to approximate the discrete processes |$\left (Z_{n}^{\varDelta }, \varGamma _{n}^{\varDelta }\right )$|⁠, respectively. We train the parameter set |$\theta _{n} = \left (\theta ^{z}_{n}, \theta ^{\gamma }_{n}\right )$| using a loss function such that the dynamics of the discretized process |$Z$| given by (4.12) (with the Malliavin derivative of the driver function evaluated at time points |$t_{n}$| and |$t_{n+1}$| (Negyesi et al., 2024)) are fulfilled, namely

  $$ \begin{align*} \mathbf{L}^{z,\varDelta}_{n}\left( \theta_{n} \right) & = \mathbb{E}\left[ \vphantom{\left(\left(\phi^{\sum\gamma}_{n}\left( X_{n}^{\varDelta}; \theta^{\gamma}_{n}\right) D_{n} X_{n}^{\varDelta}\right)^\top \varDelta W_{n}\right)^\top}\Bigg\vert Z^{\varDelta, \hat{\theta}}_{n+1} b^{-1}\Big( t_{n+1}, X_{n+1}^{\varDelta} \Big) D_{n} X^{\varDelta}_{n+1} - \phi^{z}_{n}\left( X^{\varDelta}_{n}; \theta^{z}_{n} \right) \right. \\ & \quad \left. \left. +\, f_{D}\left(t_{n+1}, \mathbf{X}^{\varDelta, \hat{\theta}}_{n+1}, \mathbf{D}_{n}\mathbf{X}_{n+1,n}^{\varDelta, \hat{\theta}}\right)\varDelta t_{n} - \left(\left(\phi^{\gamma}_{n}\left( X_{n}^{\varDelta}; \theta^{\gamma}_{n}\right) D_{n} X_{n}^{\varDelta}\right)^\top \varDelta W_{n}\right)^\top \right\vert^{2}\right], \end{align*} $$

  where |$\mathbf{X}^{\varDelta , \hat{\theta }}_{n+1}=\left ( X_{n+1}^{\varDelta }, Y^{\varDelta , \hat{\theta }}_{n+1}, Z^{\varDelta , \hat{\theta }}_{n+1} \right )$| and |$\mathbf{D}_{n}\mathbf{X}^{\varDelta , \hat{\theta }}_{n+1,n}:=\left ( D_{n} X_{n+1}^{\varDelta }, D_{n} Y^{\varDelta , \hat{\theta }}_{n+1}, \varGamma ^{\varDelta , \hat{\theta }}_{n} D_{n} X_{n}^{\varDelta } \right )$|⁠. Approximate the optimal parameters |$\theta ^{*}_{n} \in \operatorname{arg\,min}_{\theta _{n} \in \varTheta _{n}} \mathbf{L}^{z, \varDelta }_{n}\left ( \theta _{n} \right )$| using an SGD method and receive the estimated parameters |$\hat{\theta }_{n} = \left ( \hat{\theta }^{z}_{n}, \hat{\theta }^{\gamma }_{n} \right )$|⁠. Then, we define

  $$ \begin{align*} &Z_n^{\varDelta, \hat{\theta}} = \phi^z_n\left( X^{\varDelta}_n; \hat{\theta}^z_n \right), \quad \varGamma_n^{\varDelta, \hat{\theta}} = \phi^{\gamma}_n \left( X^{\varDelta}_n; \hat{\theta}^{\gamma}_n \right).\end{align*} $$

  For the second optimization problem use another DNN |$\phi ^{y}_{n}(\cdot ; \theta ^{y}_{n}): \mathbb{R}^{d} \to \mathbb{R}^{1 \times d}$| to approximate the discrete processes |$Y_{n}^{\varDelta }$|⁠. Train the parameters |$\theta ^{y}_{n}$| using a loss function such that the dynamics of the discretized process |$Y$| given by (4.7) are fulfilled, namely

  $$ \begin{align*}& \mathbf{L}^{y,\varDelta}_{n}\left( \theta_{n}^{y} \right) =\mathbb{E}\left[ \left\vert Y^{\varDelta, \hat{\theta}}_{n+1} - \phi^{y}_{n}\left( X^{\varDelta}_{n}; \theta^{y}_{n} \right) + f\left(t_{n}, X^{\varDelta}_{n}, \phi^{y}_{n}\left( X^{\varDelta}_{n}; \theta^{y}_{n} \right), Z_{n}^{\varDelta, \hat{\theta}}\right) \varDelta t_{n} - Z_{n}^{\varDelta, \hat{\theta}} \varDelta W_{n} \right\vert^{2} \right]. \end{align*} $$

  Approximate the optimal parameters |$\theta ^{*,y}_{n} \in \operatorname{arg\,min}_{\theta _{n}^{y} \in \varTheta _{n}^{y}} \mathbf{L}^{y, \varDelta }_{n}\left ( \theta _{n}^{y} \right )$| using an SGD method and receive the estimated parameters |$\hat{\theta }^{y}_{n}$|⁠. Then, we define

  $$ \begin{align*} &Y_n^{\varDelta, \hat{\theta}} = \phi^y_n\left( X^{\varDelta}_n; \hat{\theta}^y_n \right).\end{align*} $$

In comparison with the OSM scheme our approach demonstrates the following advantages:

• (i)

  Since the OSM scheme employs supervised deep learning it requires solving two optimization problems per time step—one for BSDE (4.2) and another for BSDE (4.4)—to approximate the unknown processes |$\left (Y, Z, \varGamma \right )$|⁠. Consequently, the computational cost of the OSM scheme is up to twice as high as that of our scheme. This is demonstrated in our numerical experiments.

• (ii)

  Our scheme can be seamlessly extended, not only to DBDP scheme, which is formulated backward in time through local optimizations at each discrete time step, but also to other supervised deep learning-based approaches, such as E et al. (2017); Kapllani & Teng (2024); Raissi (2024), which are formulated forward in time as a global optimization problem. This is part of our ongoing research. The OSM approach cannot be integrated to such schemes, as it cannot be formulated as a global optimization problem.

5. Convergence analysis

The main goal of this section is to prove the convergence of the DLBDP scheme towards the solution |$\left (Y, Z, \varGamma \right )$| of the BSDE system (4.1)–(4.4), and provide a rate of convergence that depends on the discretization error from the Euler–Maruyama scheme and the approximation or model error by the DNNs.

For the functions figuring in the BSDE system (4.1)–(4.4) the following assumptions are in place.

AX4. Assumption AX3 holds, with the Malliavin derivative |$\left \vert D_{s} b(t, X_{t}) \right \vert \leq C$|  |$\mathbb{P}$|-a.s. for |$0\leq s \leq t \leq T$|⁠. The functions |$a(t,x)$| and |$b(t,x)$| are |$1/2$|-Hölder continuous in time.

AY4. Assumption AY3 holds. Moreover, |$g \in C^{2+{\mathfrak{l}}}_{{\mathfrak{b}}}\left ( \mathbb{R}^{d}; \mathbb{R} \right )$|⁠, |${\mathfrak{l}}>0$|⁠. The function |$f(t,x,y,z)$| and its partial derivatives |$\nabla _{x} f$|⁠, |$\nabla _{y} f$| and |$\nabla _{z} f$| are all |$1/2$|-Hölder continuous in time.

We emphasize that our assumption for the SDE is less restrictive than that in Negyesi et al. (2024), where arithmetic Brownian motion is assumed. When pricing and hedging options usually the stock dynamics are modeled by the geometric Brownian motion (GBM). To ensure the applicability of our convergence analysis in such cases we consider in the numerical section the ln-transformation of stock prices. Consequently, we obtain a drift and diffusion function that satisfy Assumption AX4, thereby ensuring that our theoretical analysis holds in the numerical experiments. Moreover, in case of more advanced models than the GBM, if the Malliavin derivative of |$b(t, X_{t})$| is bounded, our analysis still holds.

The following lemma is a consequence of the considered assumptions.

From the mean-value theorem, for |$f \in C^{0,2,2,2}_{{\mathfrak{b}}}\left ( [0, T]\times \mathbb{R}^{d} \times \mathbb{R} \times \mathbb{R}^{1 \times d}; \mathbb{R} \right )$|⁠, we have that |$f$| and all its first-order derivatives in |$(x, y, z)$| are Lipschitz continuous. Therefore, the following holds (using also Assumption AY4 and Lemma 5.1):

with |$\mathbf{x}_{{\mathfrak{i}}} = \left ( x_{{\mathfrak{i}}}, y_{{\mathfrak{i}}}, z_{{\mathfrak{i}}} \right )$|⁠, |$\mathbf{Dx}_{{\mathfrak{i}}} = \left ( Dx_{{\mathfrak{i}}}, Dy_{{\mathfrak{i}}}, Dz_{{\mathfrak{i}}} \right )$| and |$t_{{\mathfrak{i}}} \in [0,T]$|⁠, |$x_{{\mathfrak{i}}} \in \mathbb{R}^{d}$|⁠, |$y_{{\mathfrak{i}}} \in \mathbb{R}$|⁠, |$z_{{\mathfrak{i}}}, Dy_{{\mathfrak{i}}} \in \mathbb{R}^{1 \times d}$|⁠, |$Dx_{{\mathfrak{i}}}, Dz_{{\mathfrak{i}}} \in \mathbb{R}^{d \times d}$|⁠, where |$L_{f}, L_{f_{D}}>0$| and |${{\mathfrak{i}}} = 1, 2$|⁠.

Under Assumptions AX4 and AY4, using Theorems 2.1, 2.3 and 2.5, we have that the processes |$(X, Y, Z, DX, DY )$| are all mean-squared continuous in time, more specific, there exists some constant |$C>0$| such that |$\forall \, s, r, t \in [0, T]$|

From Assumptions AX4 and AY4 and Lemma 5.1 we also see for |$0\leq s \leq t \leq T$| that

Moreover, we have the well-known error estimate that the Euler–Maruyama approximations in (4.5) admit to

under Assumption AX1 and the Hölder continuity assumption in AX4 (see Zhang, 2017, theorem 5.3.1), where the notation |$\mathscr{O}\left ( \left | \varDelta \right | \right )$| means that |$\limsup _{\left | \varDelta \right | \to 0 } \left | \varDelta \right |^{-1} \mathscr{O}\left ( \left | \varDelta \right | \right ) < \infty $|⁠. Note that under Assumption AX2 and the Hölder continuity assumption in AX4 it can be showed that the Euler–Maruyama Malliavin derivative approximations |$D_{n} X_{n+1}^{\varDelta }$| in (4.8) admit to similar error estimates as in (5.4)

Let us introduce the |$\mathbb{L}^{2}$|-regularity of |$DZ$|⁠:

with

the |$\mathbb{L}^{2}$|-projection of the corresponding Malliavin derivative w.r.t. the |$\mathscr{F}_{t_{n}}$|  |$\sigma $|-algebra. Based on relations in (4.11) we have that

Subsequently, using Lemma 5.1, the mean-squared continuity in time of |$DX$| given by Theorem 2.3 and that the terminal condition of the Malliavin BSDE (4.4) is Lipschitz continuous (due to Assumption AY4) we have that

after applying (Zhang, 2004, theorem 3.1).

We now define

for |$n=0, \ldots , N-1$|⁠, where |$\hat{\mathbf{X}}_{n}:= \left ( X_{n}^{\varDelta }, \hat{Y}_{n}^{\varDelta }, \hat{Z}_{n}^{\varDelta }\right )$| and |$\mathbf{D}_{n} \hat{\mathbf{X}}_{n}:= \left ( D_{n} X_{n}^{\varDelta }, \hat{Z}_{n}^{\varDelta }, \hat{\varGamma }_{n}^{\varDelta } b(t_{n}, X_{n}^{\varDelta })\right )$|⁠. Note that |$\hat{Y}_{n}^{\varDelta }$| and |$\hat{Z}_{n}^{\varDelta }$| in (5.7) are calculated by taking |$\mathbb{E}_{n}[\cdot ]$| in (4.7) and (4.12), where |$\mathbb{E}_{n}\left [\hat{Z}_{n}^{\varDelta } \varDelta W_{n}\right ] = 0$| and |$\mathbb{E}_{n}\left [\hat{\varGamma }_{n}^{\varDelta } b(t_{n}, X_{n}^{\varDelta })\varDelta W_{n}\right ] = 0$|⁠. Moreover, |$\hat{\varGamma }_{n}^{\varDelta }$| in (5.7) is calculated by multiplying both sides of (4.12) with |$\varDelta W_{n}$|⁠, where |$\mathbb{E}_{n}\left [\varDelta W_{n} f_{D}\left (t_{n},\hat{\mathbf{X}}_{n}^{\varDelta }, \mathbf{D}_{n} \hat{\mathbf{X}}_{n}^{\varDelta }\right )\right ] = 0$|⁠. Finally, applying the Itô isometry gives |$\hat{\varGamma }_{n}^{\varDelta }$| in (5.7).

By the Markov property of the underlying processes there exist some deterministic functions |$\hat{y}_{n}$|⁠, |$\hat{z}_{n}$| and |$\hat{\gamma }_{n}$| such that

Moreover, by the martingale representation theorem, there exists an |$\mathbb{R}^{d \times d}$|-valued square integrable process |$D_{n} \hat{Z}_{t}$| such that

and by Itô isometry we have

Hence, |$D \hat{Z}^{\varDelta }$| is an |$\mathbb{L}^{2}$|-projection of |$D \hat{Z}$|⁠. Moreover, |$\hat{Z}^{\varDelta }$| is an |$\mathbb{L}^{2}$|-projection of |$\hat{Z}$| such that

Finally, we define the approximation errors of |$\hat{y}_{n}$|⁠, |$\hat{z}_{n}$| and |$\hat{\gamma }_{n}$| by the DNNs |$\phi ^{y}_{n}$|⁠, |$\phi ^{z}_{n}$| and |$\phi ^{\gamma }_{n}$| defined as

for |$n = 0, \ldots , N-1.$| The goal is now to find an upper bound of the total approximation error of the DLBDP scheme defined as

in terms of the discretization error (from the Euler–Maruyama scheme) and the approximation errors (5.11) by the DNNs, where |$D_{n} Z_{s} - D_{n} Z^{\varDelta , \hat{\theta }}_{n} = \varGamma _{s} b(s) - \varGamma _{n}^{\varDelta , \hat{\theta }} b(t_{n})$| due to relations (4.11) and Assumption AX4.

According to Theorem 5.1 the total approximation error of the DLBDP scheme consists of four terms. The first term corresponds to

• (i)

  the strong approximation of the terminal condition and its gradient, depending on the Euler–Maruyama scheme and the functions |$\left (g(x), \nabla _{x} g(x)\right )$|⁠,

• (ii)

  the strong approximation of the Euler–Maruyama scheme and the path regularity of the processes |$\left (Y, Z\right )$|⁠, see Theorem 2.5.

The second term represents the |$\mathbb{L}^{2}$|-regularity of |$DZ$|⁠. All the aforementioned terms converge to zero as |$|\varDelta |$| goes to zero, with a rate of |$|\varDelta |$| when Assumptions AX4 and AY4 are satisfied. For the last two terms the better the DNNs are able to estimate the functions (5.8), the smaller is their contribution in the total approximation error. Note that from the UAT (Cybenko, 1989; Hornik et al., 1989) the approximation error from the DNNs can be made arbitrarily small for a sufficiently large number of hidden neurons. It is crucial noting that, in contrast to both the DBDP scheme and the method outlined in Negyesi et al. (2024), the DLBDP scheme provides a means to manage the impact of the DNN’s approximation error. This is accomplished by selecting the values of |$\omega _{1}$| and |$\omega _{2}$|⁠, resulting in improved accuracy for the processes |$\left (Y, Z, \varGamma \right )$|⁠, as we demonstrate in the next section.

6. Numerical results

In this section we illustrate the improved performance of the DLBDP scheme compared with the DBDP scheme, not only when approximating the solution, but also its gradient and the Hessian matrix. Moreover, we show that our scheme achieves similar accuracy compared with the OSM scheme for less computation time. As high-accurate gradient approximations are of great importance in finance we consider linear and nonlinear option pricing examples. All the experiments below were run in PYTHON using TensorFlow on the PLEIADES cluster (no parallelization), which consists of 268 workernodes and additionally five GPU nodes with eight NVidia HGX A100 GPUs (128 cores each, 2 TB memory and 16 GB per thread). We run the algorithms on the GPU nodes. For more information, see PLEIADES documentation.¹

6.1 Experimental set-up

In all the following examples we consider the same hyperparameters for our scheme and both the DBDP and OSM schemes for a fair comparison. For the DNNs we choose |$L = 2$| hidden layers and |$\eta = 100 + d$| neurons per hidden layer. The input is normalized based on the true moments. The input is not normalized at discrete time point |$t_{0},$| as the standard deviation is zero. A hyperbolic tangent activation |$\tanh (\cdot )$| is applied on each hidden layer. It is crucial to mention that one cannot apply batch normalization for the hidden layers as AD is required to approximate the process |$\varGamma $| in the DBDP scheme. This is because using batch normalization creates dependence for the gradients in the batch, since it normalizes across the batch dimension. When using tf.GradientTape.batch_jacobian to approximate |$\varGamma $|⁠, and if the DNN approximating |$Z$| involves tf.keras.layers.BatchNormalization layers, the resulting output has the expected shape, but its contents have an unclear meaning (see TensorFlow documentation,² batch Jacobian section). Therefore, batch normalization is emitted, not only in the DBDP scheme, but also in our scheme to ensure a fair comparison. For the SGD iterations we use the Adam optimizer with a stepwise learning rate decay approach. We choose a batch size of |$B=1024$| for each of |$\kappa $| optimization steps. At the discrete time point |$t_{N-1}$| we consider |${\mathfrak{K}} = 24000$| optimization steps, where the learning rate |$\alpha $| is adjusted as follows:

For the next discrete time points (i.e., |$t_{N-2}, \ldots , t_{0}$|⁠) we make use of the transfer learning approach, and reduce the number of optimization steps to |${\mathfrak{K}} = 10000$|⁠, and use the following learning rates:

The gradient of the driver function |$f$| w.r.t. each variable |$(x, y, z)$| and the function |$g$| w.r.t. to variable |$x$| are calculated by using AD, namely tf.GradientTape in TensorFlow. For the gradient of the function representing |$Z_{t}$| (when available) in (2.2) w.r.t. to variable |$x$|⁠, tf.GradientTape.batch_jacobian is used. Note that we consider a uniform time discretization |$\varDelta $| of |$[0, T]$|⁠. The DLBDP algorithm (without ln-transformation) calculating the final estimates |$\left (Y^{\varDelta ,\hat{\theta }}_{n}, Z^{\varDelta ,\hat{\theta }}_{n}\right )$| for |$n=N-1,\ldots ,1,0$| is given in Algorithm 1 when using the aforementioned learning rate decay and transfer learning approaches. The parameters |$\hat{\theta }$| are an estimation of |$\theta ^{*}$| due to the optimization error resulting from the Adam optimization algorithm and the estimation error from the empirical version of loss (4.13) given as

for a batch size |$B$|⁠.

We define the following mean squared errors as performance metrics for a sample with the size |$B$|⁠:

for each process. To account for the stochasticity of the underlying Brownian motion and the Adam optimizer we conduct |$Q = 10$| independent runs (training’s) of the algorithms and define, e.g.,

as the mean MSE for the process |$Y$|⁠, and similarly for the other processes. Note that as a relative measure of the MSE we consider, e.g.,

for the process |$Y$|⁠, and similarly for the other processes. We choose a testing sample of size |$B = 1024.$| The computation time (runtime) for one run of the algorithms is defined as |$\tau $|⁠, and the average computation time over |$Q=10$| runs as |$\overline{{\tau }}: = \frac{1}{Q} \sum _{q=1}^{Q} \tau _{q}.$|

6.2 The Black–Scholes BSDE

We start with a linear BSDE—the Black–Scholes BSDE—that is used for pricing of European options.

where |$c_{k}>0$| and |$\sum _{k=1}^{d} c_{k} = 1$|⁠. Note that |$a_{k}$| represents the expected return of the stock |$X_{t}^{k}$|⁠, |$b_{k}$| the volatility of the stock returns, |$\delta _{k}$| is its dividend rate and |$x_{0}^{k}$| is the price of the stock at |$t =0$|⁠. Moreover, |$X_{T}$| is the price of the stocks at time |$T$|⁠, which denotes the maturity of the option contract. The value |$K$| represents the contract’s strike price. Finally, |$R$| corresponds to the risk-free interest rate. The analytic solution (the option price |$Y_{t}$| and its delta hedging strategy |$Z_{t}$|⁠) is given by the Black–Scholes formula:

where |$\varPhi \left (\cdot \right )$| is the standard normal cumulative distribution function. The analytical solution |$\varGamma _{t} = \nabla _{x} \left ( \nabla _{x} u\left (t, X_{t}\right ) b\left (t, X_{t}\right ) \right )$| is calculated by using AD. As we mentioned in Section 5, when dealing with a forward SDE represented by the GBM, we apply the ln-transformation to ensure that the theoretical analysis is applicable to our numerical experiments. We define |$\check{X}_{t}:=\ln \left (X_{t}\right )$| and |$\check{u}(t, \check{X}_{t}):= u(t, X_{t})$|⁠. Using the Feynman–Kac formula we write the Black–Scholes BSDE in the ln-domain

The ln-transformation simplifies the Malliavin derivatives as |$D_{n} X_{n}^{k} = b_{k} X_{n}^{k}$|⁠, |$D_{n} X_{n+1}^{k} = b_{k} X_{n+1}^{k}$| and |$D_{n} \check{X}_{n}^{k} = D_{n} \check{X}_{n+1}^{k} = b_{k}$| for |$k=1,\ldots ,d$|⁠. Note that |$\left ( \check{Y}_{t}, \check{Z}_{t} \right ) = \left ( Y_{t}, Z_{t} \right )$| since |$\check{Y}_{t} = \check{u}(t, \check{X}_{t}) = u(t, X_{t}) = Y_{t}$| and |$\check{Z}_{t}^{k} = \frac{\partial \check{u}}{\partial \check{x}_{k}} b_{k} = \frac{\partial u}{\partial x_{k}} b_{k} X_{t}^{k} = Z_{t}^{k}$| for |$k = 1, \ldots , d$|⁠. Hence, we can compare the approximated solution of (6.2) in the ln-domain with the exact solution of Example 1 given in (6.1). In case of the process |$\varGamma $| we have that |$\check{\varGamma }_{t}^{k_{1}, k_{2}} \frac{1}{X_{t}^{k_{2}}} = \varGamma _{t}^{k_{1}, k_{2}}$| for |$k_{1}, k_{2}= 1, \ldots , d$|⁠. In the following tests, for |$k=1,\ldots , d$|⁠, we set |$x_{0}^{k} = 100$|⁠, |$a_{k} = 0.05$|⁠, |$b_{k} = 0.2$|⁠, |$R = 0.03$|⁠, |$c_{k} = \frac{1}{d}$| and |$\delta _{k} = 0$|⁠. Moreover, we set |$K = 100$|⁠, |$T = 1$| and |$d \in \{1, 10, 50\}$|⁠.

To provide a comparison of the approximation of each process across the discrete domain |$\varDelta $| (using the testing sample) we visualize in Fig. 1 the mean MSE values for |$d \in \{1, 10, 50\}$| from all schemes. The STD of the MSE values is given in the shaded area.

Fig. 1.

Mean MSE values of the processes |$\left (Y, Z, \varGamma \right )$| from DBDP, OSM and DLBDP schemes over the discrete time points |$\{t_{n}\}_{n=0}^{N-1}$| using the testing sample in Example 1, for |$d \in \{1, 10, 50\}$| and |$N = 64$|⁠. The STD of MSE values is given in the shaded area.

Open in new tabDownload slide

First, we compare our scheme with the DBDP scheme. For the case of |$d=1$|⁠, Fig. 1(c) clearly shows a substantial improvement in approximating the process |$\varGamma $| across the discrete time points |$\{t_{n}\}_{n=0}^{N-1}$| achieved by our scheme. Furthermore, Fig. 1(b) demonstrates that the DLBDP scheme also outperforms in approximating the process |$Z$|⁠. However, there is no improvement achieved with our scheme for the process |$Y$|⁠, as shown in Fig. 1(a). As the dimension increases to |$d=10$| and |$d=50$| our scheme further exhibits a higher accuracy for approximating the processes |$\left (Z, \varGamma \right )$|⁠. Moreover, an improvement in approximating the process |$Y$| is evident for |$d=50$| from the DLBDP scheme compared with DBDP scheme, as displayed in Fig. 1(g). The approximations from our scheme and the OSM scheme are comparable. Specifically, both schemes yield similar approximations for the process |$\varGamma $|⁠, while the OSM scheme performs better for the process |$Z$|⁠, and our scheme gives higher accuracy for the process |$Y$|⁠.

Next, we report in Table 1 the mean relative MSE of each process at |$t_{0}$| while varying |$N$| for |$d \in \{1, 10, 50 \}$| along with the average computation time from the DBDP, OSM and DLBDP schemes. The STD of the relative MSE values at |$t_{0}$| is given in the brackets.

The mean relative MSE of |$\left (Y_{0}, Z_{0}\right )$| decreases as |$N$| increases for each dimension in all schemes. This trend is also observed for |$\varGamma _{0}$| in the OSM and DLBDP schemes, but not in the DBDP scheme, which actually diverges. Note that the mean relative MSE values start to flatten out for |$N=64$|⁠, indicating that the overall contribution of the approximation error from the DNNs increases for higher |$N$| and becomes larger than the discretization error. This is consistent with the error analysis in Section 5 (see theorem 4.1 for the DBDP scheme (Huré et al., 2020) and theorem 5.2 for the OSM scheme (Negyesi et al., 2024)). Compared with the DBDP scheme our approach consistently yields the smallest mean relative MSE for each process, especially as the dimension increases. Both the OSM and DLBDP schemes provide overall comparable approximations. The average computation time of the DLBDP algorithm is higher compared with that of the DBDP algorithm. Note that we compare the computational time of all schemes including the computation of |$\varGamma $| at each optimization step. In Negyesi et al. (2024) it is mentioned that the runtime of their algorithm is roughly double of the DBDP one, as it requires solving two optimization problems per discrete time step. Since in the second optimization problem only the parameters of the DNN for the process |$Y$| are optimized one can reasonably infer that our algorithm may be up to twice as fast as the one proposed in Negyesi et al. (2024). This is observed in Table 1 when comparing the computation of the OSM and DLBDP schemes, especially as |$d$| and |$N$| increase (the algorithm’s complexity grows due to the higher number of network parameters with increasing dimensionality and the increased number of optimization problems with larger |$N$|⁠).

To train the algorithms we set a high number of optimization steps (and a high number of hidden neurons) as described in Section 6.1 such that the same hyperparameters are used for each example. However, the computation time of the algorithms can be reduced, e.g., by reducing the number of optimization steps. This can be seen in Fig. 2, where we display the mean loss and MSE values of each process for all the algorithms using a validation sample |$B=1024$|⁠, at discrete time points |$\left (t_{32}, t_{63} \right )$| in case of |$d=50$| and using |$N =64$|⁠. The mean loss is defined as |$\overline{\tilde{\mathbf{L}}}_{n}^{\varDelta }\left ( \hat{\theta }_{n} \right ):= \frac{1}{Q} \sum _{q=1}^{Q} \tilde{\mathbf{L}}_{n, q}^{\varDelta }\left ( \hat{\theta }_{n} \right )$|⁠. The STD of the loss and MSE values is given in the shaded area.

Fig. 2.

Mean loss and MSE values of the process |$\left (Y, Z, \varGamma \right )$| from DBDP, OSM and DLBDP schemes at discrete time points |$\left (t_{32}, t_{63} \right ) = \left (0.5000, 0.9844\right )$| using the validation sample in Example 1, for |$d=50$| and |$N=64$|⁠. The STD of the loss and MSE values is given in the shaded area.

Open in new tabDownload slide

By choosing for instance |${\mathfrak{K}} = 16000$| at |$t_{63}$| and |${\mathfrak{K}} = 5000$| at other discrete time points the runtime of the algorithms is substantially reduced with almost an insignificant loss of accuracy.

6.3 Option pricing with different interest rates

We now consider a pricing problem involving a European option in a financial market where the different interest rates for borrowing and lending are different. This model, originally introduced in Bergman (1995), and has been addressed in e.g., E et al. (2017, 2019); Teng (2021, 2022) is represented by a nonlinear BSDE.

where |$R_{1}$| and |$R_{2}$| are the interest rates for lending and borrowing, respectively. Note that instead of solving the above BSDE directly we solve the transformed BSDE in the ln-domain.

We test all schemes in the case of |$d=50$|⁠, using |$K_{1}=120$|⁠, |$K_{2} = 150$|⁠, |$T=0.5$|⁠, |$x_{0}=100 \mathbf{1}_{50}$|⁠, |$a = 0.06$|⁠, |$b = 0.2$|⁠, |$R_{1} = 0.04$| and |$R_{2} = 0.06$|⁠. The benchmark value is |$Y_{0} \doteq 17.9743$|⁠, which is computed using the multilevel Monte Carlo approach (E et al., 2019) with seven Picard iterations and |$Q=10$| independent runs. For |$N \in \{2, 8, 32, 64\}$|⁠, we show in Table 2 the approximation for |$Y_{0}$| (the reference results for |$Z_{0}$| are not available) from all algorithms and their average runtime. More precisely, we report the mean approximation of |$Y_{0}$| defined as |$\overline{Y}_{0}^{\varDelta , \hat{\theta }}:= \frac{1}{Q} \sum _{q=1}^{Q} Y_{0,q}^{\varDelta , \hat{\theta }}$|⁠, with the mean relative MSE and its STD given in the brackets.

We observe that our scheme consistently provides higher accurate approximations of |$Y_{0}$| for the |$50$|-dimensional nonlinear BSDE in Example 2 compared with the other schemes, resulting in smaller relative MSE value. The DBDP scheme achieves the shortest computation time, while our scheme is faster than the OSM scheme. Note that the mean relative MSE can be further reduced by increasing the number of hidden neurons or layers provided that the optimization error is sufficiently small.

6.4 The Black–Scholes extended with local volatility

The next example is taken from Ruijter & Oosterlee (2016) in order to demonstrate the effectiveness of our scheme in case of a time-dependent diffusion function. Consider an European call option as in Example 1, where each underlying asset follows a GBM with time-dependent drift and diffusion.

where for |$a(t)$| and |$b(t)$| we choose the following periodic functions:

The exact solution of this local volatility model is given by the Black–Scholes formula with volatility parameter |$ \bar{b} = \sqrt{\frac{1}{T-t} \int _{t}^{T} b(s)^{2} \,{\text{d}}s}$|⁠. More precisely, the exact solution is given by (6.1) with

We apply the ln-transformation in this example, which is similar as in the case of Example 1. Moreover, we set |$T = 0.25$|⁠, |$d=50$| and the other following parameter values

Using |$N=32$|⁠, the mean MSE values for each process over discrete domain |$\varDelta $| are visualized in Fig. 3 for the testing sample. The STD of the MSE values is displayed in the shaded area.

Fig. 3.

Mean MSE values of the processes |$\left (Y, Z, \varGamma \right )$| from DBDP, OSM and DLBDP schemes over the discrete time points |$\{t_{n}\}_{n=0}^{N-1}$| using the testing sample in Example 3, for |$d=50$| and |$N = 32$|⁠. The STD of MSE values is given in the shaded area.

Open in new tabDownload slide

Compared with previous examples we notice significant improvements from our scheme, not only in approximating the process |$Z$|⁠, but also the process |$Y$| compared with the DBDP scheme. In the case of the process |$\varGamma $| such improvements are evident only near |$t_{0}$|⁠. Interestingly, the DLBDP scheme outperforms the OSM scheme in this example for the processes |$Y$| and |$\varGamma $| while providing comparable approximations of the process |$Z$|⁠.

In Table 3 we report the mean relative MSE values at |$t_{0}$| for each process from all schemes, using |$N \in \{2, 8, 16, 32\}$|⁠. The corresponding STD is given in the brackets. The average runtime of the algorithms is also included.

Our scheme gives overall the smallest relative MSE values. In this example the improvement in approximating |$Y_{0}$| is more evident than in previous examples.

6.5 BSDE with nonadditive diffusion

We now consider the nonsymmetric example in Negyesi et al. (2024) to demonstrate the performance of our scheme when the noise in the forward SDE is nonadditive.

where |$c_{1}, c_{2} \in \mathbb{R}_{+}.$| The analytical solution is given by

We choose |$d = 50$|⁠, |$T = 10$|⁠, |$c_{1} = 10 d$|⁠, |$c_{2} = 1$| and |$x_{0} = \mathbf{1}_{d}$|⁠. In Fig. 4 we display the mean MSE values for each process over discrete domain |$\varDelta $| using the testing sample and |$N=64$|⁠, where the STD of the MSE values is visualized in the shaded area. Note that for |$N=64$| the approximations from the OSM scheme are not available, because the scheduled scripts in the GPU nodes of PLEIADES cluster have a time limit of 3 days. Therefore, only the approximations from the DBDP and DLBDP schemes are displayed. Our scheme clearly outperforms the DBDP scheme in approximating each process during the entire discrete time domain.

Fig. 4.

Mean MSE values of the processes |$\left (Y, Z, \varGamma \right )$| from DBDP and DLBDP schemes over the discrete time points |$\{t_{n}\}_{n=0}^{N-1}$| using the testing sample in Example 4, for |$d=50$| and |$N = 64$|⁠. The STD of MSE values is given in the shaded area.

Open in new tabDownload slide

We report in Table 4 the mean MSE values (due to small magnitude of the exact solution) at |$t_{0}$| for each process and the algorithm average runtime using |$N \in \{ 2, 8, 32, 64\}$|⁠. The STD of the relative MSE values at |$t_{0}$| is given in the brackets. The same conclusions can be drawn with our scheme when compared with the DBDP and OSM schemes even in the case of a more general diffusion term.

6.6 The Black–Scholes BSDE with correlated noise

Finally, we test all the schemes using an example with correlated noise. Specifically, we consider a European max call option within the Black–Scholes framework for a basket of stocks with distinct parameters (expected return, volatility and correlation). The dynamics of the stocks are therefore given as

By applying the Cholesky decomposition to the correlation matrix |$\left (\rho _{k,j}\right )_{k, j = 1,\ldots , d}$| and transforming the stock dynamics into the ln-domain the corresponding BSDE is given as follows.

where |$\check{W_{t}}$| are |$d$| independent Brownian motions, |$\check{\rho }_{k,j}$| represents the elements of the lower triangular matrix from Cholesky decomposition of |$\left (\rho _{k,j}\right )_{k, j = 1,\ldots , d}$| and

We set |$d = 20$|⁠, |$T = 0.5$|⁠, |$R = 0.05$|⁠, |$K=100$| and |$\delta _{k} = 0$| for |$k=1, \ldots , d$|⁠. The expected returns, volatilities and correlation matrix are generated randomly. Specifically, |$x_{0}^{k} \sim \mathscr{U}\,\, [K-0.05 K, K + 0.05K]$|⁠, |$a^{k} \sim \mathscr{U}\,\, [0.01, 0.1]$| and |$b^{k} \sim \mathscr{U}\,\, [0.05, 0.3]$| for |$k=1, \ldots , d$|⁠. The correlation matrix is sampled from |$\mathscr{U}\,\, [-1, 1]$|⁠, ensuring that it is symmetric, positive definite and has diagonal elements equal to one. To compute a benchmark value of |$Y_{0}$| we use the Monte Carlo method (under the exact solution of the stock dynamics in the ln-domain) with |$10^{7}$| Brownian motion samples and |$50$| independent runs. Table 5 reports the mean approximation of |$Y_{0}$|⁠, the mean relative MSE values and the average runtime for all schemes using |$N \in \{2, 8, 32, 64\}$|⁠. Standard deviations are provided in parentheses.