Study of Pricing of High-Dimensional Financial Derivatives Based on Deep Learning

Abstract

Many problems in the fields of finance and actuarial science can be transformed into the problem of solving backward stochastic differential equations (BSDE) and partial differential equations (PDEs) with jumps, which are often difficult to solve in high-dimensional cases. To solve this problem, this paper applies the deep learning algorithm to solve a class of high-dimensional nonlinear partial differential equations with jump terms and their corresponding backward stochastic differential equations (BSDEs) with jump terms. Using the nonlinear Feynman-Kac formula, the problem of solving this kind of PDE is transformed into the problem of solving the corresponding backward stochastic differential equations with jump terms, and the numerical solution problem is turned into a stochastic control problem. At the same time, the gradient and jump process of the unknown solution are separately regarded as the strategy function, and they are approximated, respectively, by using two multilayer neural networks as function approximators. Thus, the deep learning-based method is used to overcome the “curse of dimensionality” caused by high-dimensional PDE with jump, and the numerical solution is obtained. In addition, this paper proposes a new optimization algorithm based on the existing neural network random optimization algorithm, and compares the results with the traditional optimization algorithm, and achieves good results. Finally, the proposed method is applied to three practical high-dimensional problems: Hamilton-Jacobi-Bellman equation, bond pricing under the jump Vasicek model and option pricing under the jump diffusion model. The proposed numerical method has obtained satisfactory accuracy and efficiency. The method has important application value and practical significance in investment decision-making, option pricing, insurance and other fields.

1. Introduction

High-dimensional (dimension ≥ 3) nonlinear partial differential equations (PDEs) are among the subjects that draw a lot of interest and are widely used in a variety of fields. Many practical problems require the use of a high-dimensional PDE, for example: the Schrodinger equation of the quantum many-body problem, whose PDE dimension is about three times that of the electron or quantum (particle) in the system; the Black Scholes equation used to price financial derivatives, where the dimension of the PDE is the number of relevant financial assets under consideration; and the Hamilton-Jacobi-Bellman equation in dynamic programming. High-dimensional nonlinear partial differential equations are very practical, but because there are few explicit solutions or the analytical expressions are too complicated, it is frequently necessary to use some numerical techniques to solve them. However, in practical applications, high-dimensional nonlinear partial differential equations are usually very difficult to solve and remain one of the challenging topics in the academic community. The difficulty is that, due to the “curse of dimensionality” [1], the time complexity of traditional numerical solutions will increase exponentially with the increase of dimensions, thus requiring extensive computing resources. Nevertheless, we urgently need to approximate the numerical solution of these high-dimensional nonlinear partial differential equations because these equations can resolve a variety of real-world issues.

In the field of finance, the method of solving high-dimensional PDEs is frequently employed. However, the “heavy tail” of financial markets has been demonstrated by numerous empirical research studies. “Abnormal” events such as new inventions, economic policies, wars or other news, as well as various arbitrage and hedging activities, can lead to sudden and intermittent changes in asset prices. These sudden changes in price and price fluctuations intuitively reflect the jumping effect that exists in financial markets. Although jumping behavior has a low likelihood of happening, when it does happen, it causes significant harm and often results in investors suffering enormous losses or even going bankrupt. How to build adequate models to depict the jumping behavior of asset prices has become a hot topic for academics as financial market research has become more sophisticated and miniaturized. Therefore, it is necessary to incorporate common jumping behaviors in financial problems, such as buying (selling) property defaults, corporate bankruptcies, operational failures, insurance events, etc., into high-dimensional PDE numerical solutions.

Backward stochastic differential equations (BSDEs) are widely used in stochastic control and financial mathematics. One important reason is that a BSDE can be expressed as the stochastic counterpart of a PDE through the Feynman-Kac formula. El Karoui N et al. [2] provided a detailed introduction to the different properties of BSDEs and their applications to finance. Zhang J [3] proposed a numerical scheme for a class of BSDEs with possible path-dependent terminal values and proved its convergence. Recently, Barigou K and Delong Ł [4] also examined the pricing of equity-linked life insurance contracts by framing the price as the solution of a system of non-linear PDEs, reformulating the problem as a BSDE with jumps and solving it numerically using efficient neural networks.

In recent years, deep learning algorithms have rapidly gained popularity and success in various application sectors. As a result, several scientists have had success using deep learning methods to handle high-dimensional PDE issues. In 2017, E W. Han et al. [5,6] first proposed the deep BSDE method and systematically applied deep learning to general high-dimensional PDEs. Following that, Han et al. (2020) [7] continued to apply the deep BSDE method to the study of random games and also produced promising outcomes. The method was expanded to include the case of 2BSDE and the corresponding fully nonlinear PDE by Beck et al. (2019) [8]. Raissi et al. (2017) [9] use the latest development of probabilistic machine learning to infer the control equation represented by a parametric linear operator and modify the prior value of the Gaussian process according to the special form of the operator, which is then used to infer the parameters of the linear equation from the scarce and possibly noisy observations. To approximate the solution of a high-dimensional PDE, Sirignano et al. (2018) [10] suggested using a deep neural network technique. The network was trained to satisfy boundary conditions, initial conditions and differential operators in a batch of randomly sampled time and space. The achievements listed above show that it is feasible to solve high-dimensional PDEs with a deep learning-based method.

In this paper, we mainly apply the deep learning algorithm to a special class of high-dimensional nonlinear partial differential equations with jumps and obtain a numerical solution of the equation. Specifically, through mathematical derivation and equivalent formula expression of a high-dimensional PDE and backward stochastic differential equation, the problem of solving a partial differential equation is equivalent to the problem of solving a BSDE, and then it is transformed into a stochastic control problem. Then, a deep neural network framework is created to address the issue. Therefore, this method can be used to solve high-dimensional PDEs and corresponding backward stochastic differential equations simultaneously.

The method introduced in this paper to solve a kind of high-dimensional PDE with jump by using deep learning is mainly applied in the financial field, and it has important applications in financial derivatives pricing, insurance investment decision-making, small and micro enterprise financing, risk measurement mechanism and other issues.

2. Background Knowledge

2.1. A Class of PDE

We consider a class of semilinear parabolic PDEs with jump term in the following form:

Let $T\in \left(0,\infty \right)$, $d\in \mathbb{N}$, and $f:ℝ\times {ℝ}^d\to ℝ$ and $g:{ℝ}^d\to ℝ$ be continuous functions. For all $t\in \left[0,T\right]$, $x\in {ℝ}^d$, the function to be solved $u=u\left(t,x\right)\in C^{1,2}\left(\left[0,T\right]\times {ℝ}^d,ℝ\right)$ has a terminal value condition $u\left(T,x\right)=g\left(x\right)$ and satisfies:

$$\frac{\partial u}{\partial t}\left(t,x\right)-Lu\left(t,x\right)-f\left(t,x,u\left(t,x\right),{\sigma }^T\left(t,x\right){\nabla }_{x}u\left(t,x\right),Bu\left(t,x\right)\right)=0$$

(1)

where we introduce two operators:

$$Lu\left(t,x\right)={\nabla }_{x}u\left(t,x\right)\left[\mu \left(x\right)-\lambda \beta \left(x\right)\right]+\frac{1}{2}Tr\left(\sigma {\sigma }^T\left(t,x\right)\left(Hes{s}_{x}u\right)\left(t,x\right)\right)+\lambda \left\{u\left(t,x+\beta \left(x\right)\right)-u\left(t,x\right)\right\}$$

(2)

$$Bu\left(t,x\right)=u\left(t,x+\beta \left(x\right)\right)-u\left(t,x\right)$$

(3)

Here $t$ is the time variable, $x$ is the $d$-dimensional space variable, $f$ is the nonlinear part of the equation, ${\nabla }_{x}u$ represents the gradient of $u$ with respect to $x$, $Hes{s}_{x}u$ represents the Hessian matrix of $u$ with respect to $x$. In particular, this equation can be regarded as a special case of partial integro-differential equation. What we are interested in is the solution at $t=0$, $x=\xi \in {ℝ}^d$, which is $u\left(0,\xi \right)$ [11,12].

2.2. Backward Stochastic Differential Equations with Jumps

Let $\left(\mathsf{\Omega },\mathcal{F},P\right)$ be a complete probability space, $W:\left[0,T\right]\times \mathsf{\Omega }\to {ℝ}^d$ be the d-dimensional standard Brownian motion in this probability space, ${\left\{{\mathcal{F}}_t\right\}}_{t\in \left[0,T\right]}$ be the normal filtration generated by $W$ in space $\left(\mathsf{\Omega },\mathcal{F},P\right)$. ${\left\{X_t\right\}}_{0\le t\le T}$, ${\left\{Y_t\right\}}_{0\le t\le T}$, ${\left\{Z_t\right\}}_{0\le t\le T}$, ${\left\{U_t\right\}}_{0\le t\le T}$ are integrable $\mathbb{F}$-adapted stochastic processes. For a class of forward backward stochastic differential equations with jump terms, it has the following form:

$$\left\{\begin{array}{c}{X}_t=X_0+{{\displaystyle \int }}_0^t\mu \left(s,X_s\right)ds+{{\displaystyle \int }}_0^t\sigma \left(s,X_s\right)d{W}_s+{{\displaystyle \int }}_0^t\beta \left(s ,X_s\right)d{\tilde{N}}_s\\ Y_t=g\left(X_T\right)+{{\displaystyle \int }}_t^{T}f\left(s,X_s, Y_s,Z_s,U_s\right)ds-{{\displaystyle \int }}_t^T{Z}_{s}d{W}_s-{{\displaystyle \int }}_t^T{U}_{s}d{\tilde{N}}_s\end{array}\right.$$

(4)

where $d{\tilde{N}}_t$ is the compensated Poisson measure, ${\tilde{N}}_t$ is a centralized Poisson process with the compensator intensity $\lambda t$ such that ${\tilde{N}}_t=N_t-\lambda t$, $N_t$ is the Poisson process in the probability space $\left(\mathsf{\Omega },\mathcal{F},P\right)$ such that $E\left(N_t\right)=\lambda t$.

Under the standard Lipschitz assumptions on the coefficients $\mu$, $\sigma$, $\beta$, $f$, $g$, the existence and uniqueness of the solution have been proved [13].

2.3. The Generalized Nonlinear Feynman-Kac Formula

Under the appropriate regularization assumption, and if $u\left(t,x\right)\in C\left(\left[0,T\right]\times R^d\right)$ satisfies Equation (1) and linear growth condition $\left|u\left(t,x\right)\right|\le K\left(1+\left|x\right|\right)$, $\left|{\nabla }_{x}u\left(t,x\right)\right|\le K\left(1+\left|x\right|\right)$, $\left(t,x\right)\in C\left(\left[0,T\right]\times R^d\right)$, the following relationships are established almost everywhere:

$$Y_t=u\left(t,X_t\right)\in R$$

(5)

$$Z_t={\nabla }_{x}u\left(t,X_t\right)\sigma \left(t,X_t\right)\in R^d$$

(6)

$$U_t=u\left(t,X_t+\beta \left(t,X_t\right)\right)-u\left(t,X_t\right)\in R$$

(7)

$\left(Y_t,Z_t,U_t\right)$ is the only solution of the BSDE [11,12].

2.4. Improvement of Neural Network Parameter Optimization Algorithm

The Adam (Adaptive Moment Estimation) algorithm is an algorithm that combines the RMSProp algorithm with classical momentum in physics. It uses the first-order moment estimate and the second-order moment estimation of gradient to dynamically modify the learning rate of each parameter. The Adam optimizer is one of the most popular classical optimizers in deep learning, and it also shows excellent performance in practice [14].

Although Adam combines RMSprop with momentum, the adaptive moment estimation with Nesterov acceleration is often superior to momentum. Therefore, we consider including the Nesterov acceleration effect [15] into the Adam algorithm, that is, using the Nadam (Nesterov-accelerated Adaptive Moment Estimation) optimization algorithm. The calculation formula is as follows:

$${\widehat{g}}_t=\frac{g_t}{1-{{\displaystyle \prod }}_{i=1}^t{\mu }_i}$$

(8)

$$m_t=\mu \ast m_{t-1}+\left(1-\mu \right)\ast g_t$$

(9)

$${\widehat{m}}_t=\frac{m_t}{1-{{\displaystyle \prod }}_{i=1}^{t+1}{\mu }_i}$$

(10)

$$n_t=\nu \ast n_{t-1}+\left(1-\nu \right)\ast g_t^2$$

(11)

$${\widehat{n}}_t=\frac{n_t}{1-{\nu }^t}$$

(12)

$${\overline{m}}_t=\left(1-{\mu }_t\right){\widehat{g}}_t+{\mu }_{t+1}{\widehat{m}}_t$$

(13)

$${\theta }_t={\theta }_{t-1}-\eta \frac{{\overline{m}}_t}{\sqrt{{\widehat{n}}_t}+\epsilon }$$

(14)

For parameters of the neural network ${\theta }_t$, $g_t$ is the gradient of ${\theta }_t$. Besides, $m_t$ and $n_t$ are the first order moment estimate and the second order moment estimate of the gradient, respectively, which can be regarded as the estimation of the expectation $E\left|g_t\right|$ and $E\left|g_t^2\right|$, and here $\mu$ and $\nu$ are their attenuation rates, respectively. Moreover, ${\widehat{m}}_t$ and ${\widehat{n}}_t$ are the correction for $m_t$ and $n_t$, which can be approximated as an unbiased estimate of the expectation. It can be seen that NAdam has a stronger constraint on the learning rate and more direct influence on the gradient update. As a result, we attempt to apply the new optimization algorithm to our method and compare the outcomes with the traditional Adam algorithm.

3. Main Theorem

3.1. Basic Ideas

In this paper we propose a deep learning-based PDE numerical solution for nonlinear PDEs with jump terms in the form of Equation (1). The basic ideas of the algorithm are as follows:

(1)

By using the generalized nonlinear Feynman-Kac formula, the PDEs to be solved can be equivalently constructed using BSDEs with jumps.

(2)

Taking the gradient operator and the jump term of the solution of the function to be solved as policy functions, the problem of solving the numerical solution of BSDEs can be considered as a stochastic control problem, which can be further considered as a reinforcement learning problem.

(3)

Two different deep neural networks are used to approximate this pair of high-dimensional strategy functions, and the neural networks are trained through a deep learning method to obtain the numerical solution of the original equation.

3.2. Transforming the Nonlinear PDE Numerical Solution Problem into a Stochastic Control Problem

It is well known that the generalized nonlinear Feynman-Kac formula can connect a PDE and a BSDE under appropriate assumptions. Thus, the solution $Y_t\left(x\right)$ of the above equation is equivalent to the viscous solution $u\left(t,x\right)$ of the semilinear parabolic partial differential equation with jump term in Equation (1) [16,17]. We can estimate the solution $u\left(t,x\right)$ of Equation (1) by estimating the solution of Equation (4) because the numerical solution of the traditional partial differential equation does not perform well in high dimensions. Since BSDEs developed from the study of stochastic control problems and have common characteristics and internal relations, the nonlinear PDE problem can also be transformed into a stochastic control problem through the nonlinear Feynman-Kac formula.

As for the stochastic control problem related to Equation (4), the following formula holds almost everywhere for all $\left(t,x\right)\in C\left(\left[0,T\right]\times R^d\right)$:

$$\left\{\begin{array}{c}{X}_s^{t,x}=x+{{\displaystyle \int }}_{\mathrm{t}}^s\mu \left(r,X_r^{t,x}\right)dr+{{\displaystyle \int }}_t^s\sigma \left(r,X_r^{t,x}\right)d{W}_r+{{\displaystyle \int }}_t^s\beta \left(r,X_r^{t,x}\right)d{\tilde{N}}_r\\ Y_s^{t,x}=g\left(X_T^{t,x}\right)+{{\displaystyle \int }}_s^{T}f\left(r,X_r^{t,x},Y_r^{t,x},Z_r^{t,x},U_r^{t,x}\right)ds-{{\displaystyle \int }}_s^T{Z}_r^{t,x}d{W}_r-{{\displaystyle \int }}_s^T{U}_r^{t,x}d{\tilde{N}}_r\end{array}\right.$$

(15)

Among them, $s\in \left[t,T\right]$ [18].

Under the appropriate regularity assumption, for the nonlinear function $f$, it holds that a group of solution functions composed of $u\left(0,\xi \right)\in ℝ$, ${\nabla }_{x}u\left(t,X_t\right)\sigma \left(t,X_t\right)\in \mathcal{A}$ and $u\left(t,X_t+\beta \left(t,X_t\right)\right)-u\left(t,X_t\right)\in ℝ$ is the unique global minimum of the following functions $\left(y,Z,U\right)\in ℝ\times \mathcal{A}\times ℝ$:

$$\left(y,Z,U\right)⟼\mathbb{E}\left[{\left|Y_T^{y,Z,U}-g\left(X_T^{t,x}\right)\right|}^2\right]\in \left[0,\infty \right]$$

(16)

We regard $Z$ and $U$ as the policy functions of a deep learning problem and use DNN to approximate $Z$ and $U$. In this way, the stochastic process $u\left(t,X_t\right)$ corresponds to the solution function of the stochastic control problem and can be solved by using the policy functions $Z$ and $U$, and thus we can transform the numerical solution problem of the nonlinear partial differential equation into a stochastic control problem.

3.3. Forward Discretization of the Backward Stochastic Differential Equations with Jumps

For the following BSDE:

$$\left\{\begin{array}{c}{X}_t=x+{{\displaystyle \int }}_0^t\mu \left(s,X_s\right)ds+{{\displaystyle \int }}_0^t\sigma \left(s,X_s\right)d{W}_s+{{\displaystyle \int }}_0^t\beta \left(s ,X_s\right)d{\tilde{N}}_s\\ Y_t=g\left(X_T\right)+{{\displaystyle \int }}_t^{T}f\left(s,X_s, Y_s,Z_s\right)ds-{{\displaystyle \int }}_t^T{Z}_{s}d{W}_s-{{\displaystyle \int }}_t^T{U}_{s}d{\tilde{N}}_s\end{array}\right.$$

(17)

For all $t_1,t_2\in \left[0,T\right]$ and $t_1\le t_2$, the following equations hold almost everywhere:

$$X_{t_2}=X_{t_1}+{{\displaystyle \int }}_0^t\mu \left(s,X_s\right)ds+{{\displaystyle \int }}_0^t\sigma \left(s,X_s\right)d{W}_s+{{\displaystyle \int }}_0^t\beta \left(s ,X_s\right)d{\tilde{N}}_s$$

(18)

$$Y_{t_2}=Y_{t_1}-{{\displaystyle \int }}_{t_1}^{t_2}f\left(s,X_s, Y_s,Z_s\right)ds+{{\displaystyle \int }}_{t_1}^{t_2}{Z}_{s}d{W}_s+{{\displaystyle \int }}_{t_1}^{t_2}{U}_{s}d{\tilde{N}}_s$$

(19)

By substituting (6) and (7) into (19), we can obtain:

$$\begin{gathered} Y_{t_2}=Y_{t_1}-{{\displaystyle \int }}_{t_1}^{t_2}f\left(s,X_s, Y_s,\left({\nabla }_{x}u\sigma \right)\left(s,X_{\mathrm{s}}\right)\right)ds+{{\displaystyle \int }}_{t_1}^{t_2}\langle \left({\nabla }_{x}u\sigma \right)\left(s,X_{\mathrm{s}}\right),d{W}_s\rangle {}_{R^d} \\ +{{\displaystyle \int }}_{t_1}^{t_2}\langle u\left(t,X_s+\beta \left(s,X_s\right)\right)-u\left(s,X_s\right),d{\tilde{N}}_s\rangle {}_R \end{gathered}$$

(20)

Now, we will discretize them in the time dimension and divide the time interval $\left[0,T\right]$ into N partitions, so that $t_0,t_1,\dots ,t_N\in \left[0,T\right]$ and $0=t_0<t_1<\dots <t_N=T$. For $N\in \mathbb{N}$ sufficiently large, according to (18) and (19) we have:

$$X_{t_{n+1}}-X_{t_n}=\mu \left(t_n,X_{t_n}\right)\Delta t_n+\sigma \left(t_n,X_{t_n}\right)\Delta W_{t_n}+\beta \left(t_n,X_{t_n}\right)\Delta {\tilde{N}}_{t_n}$$

(21)

$$Y_{t_{n+1}}=Y_{t_n}-f\left(t_n,X_{t_n},Y_{t_n},\left({\nabla }_{x}u\sigma \right)\left(t_n,X_{t_n}\right)\right)\Delta t_n+\langle \left({\nabla }_{x}u\sigma \right)\left(t_n,X_{t_n}\right),W_{t_{n+1}}-W_{t_n}\rangle +\langle U\left(t_n,X_{t_n}\right),{\tilde{N}}_{t_{n+1}}-{\tilde{N}}_{t_n}\rangle$$

(22)

where $\Delta t_n=t_{n+1}-t_n$, $\Delta W_{t_n}=W_{t_{n+1}}-W_{t_n}$, $\Delta {\tilde{N}}_{t_n}={\tilde{N}}_{t_{n+1}}-{\tilde{N}}_{t_n}$. [5,18].

The above is the Euler scheme for discretization. In this way, we can start from the initial value $X_0$ of (4) and finally obtain the approximation of the Euler scheme of (4) at N partitions.

3.4. General Framework for Neural Networks

Two feedforward neural networks are established at each time step $t=t_n$ in (21) and (22). One is to approximate the gradient of the unknown solution, which means to approximate the function $X_{t_n}\to Z_{t_n}: {\nabla }_{x}u\left(t_n,X_{t_n}\right)\sigma \left(t_n,X_{t_n}\right)$, and this neural network is recorded as ${\mathrm{NN}}_{{\mathsf{\theta }}_{{\mathrm{z}}_{\mathrm{n}}}}\left(\mathrm{x}\right)$, such that ${\mathsf{\theta }}_{{\mathrm{z}}_{\mathrm{n}}}$ represents all parameters of this neural network and $Z_{\mathrm{n}}$ indicates that this neural network is used to approximate $Z_{t_n}$ at time $t_n$. Another is to approximate the jump process of the unknown solution, which means to approximate the function $X_{t_n}\to U_{t_n}:u\left(t_n,X_{t_n}+\beta \left(t_n,X_{t_n}\right)\right)-u\left(t_n,X_{t_n}\right)$, and this neural network is recorded as ${\mathrm{NN}}_{{\mathsf{\theta }}_{{\mathrm{U}}_{\mathrm{n}}}}\left(\mathrm{x}\right)$, such that ${\mathsf{\theta }}_{{\mathrm{U}}_{\mathrm{n}}}$ represents all parameters of this neural network, and ${\mathrm{U}}_{\mathrm{n}}$ indicates that this network is used to approximate $U_{t_n}$ at time $t_n$. For the convenience of expression, we can suppose ${\mathsf{\theta }}_Z=\left\{{\mathsf{\theta }}_{{\mathrm{Z}}_1},\dots ,{\mathsf{\theta }}_{{\mathrm{Z}}_{\mathrm{N}-1}}\right\}$,${\mathsf{\theta }}_U=\left\{{\mathsf{\theta }}_{{\mathrm{U}}_1},\dots ,{\mathsf{\theta }}_{{\mathrm{U}}_{\mathrm{N}-1}}\right\}$, $\theta =\left\{{\mathsf{\theta }}_Z,{\mathsf{\theta }}_{\mathrm{U}}\right\}$. As shown in Figure 1, all sub-neural networks are stacked together to form a complete neural network.

$X_{t_0}=\xi \in {ℝ}^d$ in the initial layer is a random variable. $u\left(t_0,X_{t_0}\right)$, ${\nabla }_{x}u\left(t_0,X_{t_0}\right)$ and $u\left(t_0,X_{t_0}+\beta \left(t_0,X_{t_0}\right)\right)-u\left(t_0,X_{t_0}\right)$ in the initial layer are unknown, and we treat them as parameters of the neural network. $u\left(t_0,X_{t_0}\right)$, ${\nabla }_{x}u\left(t_0,X_{t_0}\right)$ and $u\left(t_0,X_{t_0}+\beta \left(t_0,X_{t_0}\right)\right)-u\left(t_0,X_{t_0}\right)$ correspond to the value of the BSDE with jumps as follows:

$$Y_0=u\left(t_0,\xi \right)$$

(23)

$$Z_0={\nabla }_{x}u\left(t_0,\xi \right)$$

(24)

and

$$U_0=u\left(t_0,\xi +\beta \left(t_0,\xi \right)\right)-u\left(t_0,\xi \right)$$

(25)

In this neural network, the $X_{t_n}$ of the current layer is related to the $X_{t_{n-1}}$ of the previous layer, and at the same time the $u\left(t_n,X_{t_n}\right)$ of the current layer is related to the $X_{t_{n-1}}$, $u\left(t_{n-1},X_{t_{n-1}}\right)$, ${\nabla }_{x}u\left(t_{n-1},X_{t_{n-1}}\right)$ and $u\left(t_{n-1},X_{t_{n-1}}+\beta \left(t_{n-1},X_{t_{n-1}}\right)\right)-u\left(t_{n-1},X_{t_{n-1}}\right)$ of the previous layer. However, there are no ${\nabla }_{x}u\left(t_{n-1},X_{t_{n-1}}\right)$ and $u\left(t_{n-1},X_{t_{n-1}}+\beta \left(t_n,X_{t_n}\right)\right)-u\left(t_{n-1},X_{t_{n-1}}\right)$ in the current layer. Therefore, as shown in Figure 1, our solution is to start from the $X_{t_n}$ of the current layer to build two sub-neural networks to represent these two values. In addition, the construction of the final loss function can be obtained from a given terminal value condition $u\left(T,x\right)$, that is, $u\left(t_N,X_{t_N}\right)$ in the neural network, which also corresponds to $g\left(X_T\right)$ in the nonlinear Feynman-Kac formula.

Specifically, for $n=\left\{0,1,\dots ,N-1\right\}$, we can set ${\widehat{Y}}_0$, ${\widehat{Z}}_0$ as parameters to obtain the appropriate Euler approximation $\widehat{Y}$ of the process forward:

$${\widehat{Y}}_{t_{n+1}}={\widehat{Y}}_{t_n}-f\left(t_n,X_{t_n},{\widehat{Y}}_{t_n},{\widehat{Z}}_{t_n}\right)\Delta t_n+\langle {\widehat{Z}}_{t_n},\Delta W_{t_n}\rangle +\langle {\widehat{U}}_{t_n},\Delta {\tilde{N}}_{t_n}\rangle$$

(26)

For all appropriate ${\mathsf{\theta }}_{{\mathrm{z}}_{\mathrm{n}}}$, we have ${\widehat{Z}}_{t_n}={\mathrm{NN}}_{{\mathsf{\theta }}_{{\mathrm{z}}_{\mathrm{n}}}}\left(X_{t_n}\right)\approx \left({\nabla }_{x}u\sigma \right)\left(t_n,X_{t_n}\right)$. For all appropriate ${\mathsf{\theta }}_{{\mathrm{U}}_{\mathrm{n}}}$, we have ${\widehat{U}}_{t_n}={\mathrm{NN}}_{{\mathsf{\theta }}_{{\mathrm{U}}_{\mathrm{n}}}}\left(X_{t_n}\right)\approx u\left(t_n,X_{t_n}+\beta \left(t_n,X_{t_n}\right)\right)-u\left(t_n,X_{t_n}\right)$. Then, we can obtain a suitable approximation of $u\left(0,\xi \right)$:

$${\widehat{Y}}_0\approx u\left(0,\xi \right)$$

(27)

The mean square error between the final output ${\widehat{Y}}_{t_N}$ of neural network and the true value $g\left(X_{t_N}\right)$ at the terminal time is selected as the loss function:

$$\theta ⟼loss\left(\theta \right)=\mathbb{E}\left[{\left|{\widehat{Y}}_{t_N}-g\left(X_{t_N}\right)\right|}^2\right]\in \left[0,\infty \right]$$

(28)

$\theta$ is the set of all training parameters of the above system, such that $\theta =\left\{{\widehat{Y}}_0,{\widehat{Z}}_0,{\widehat{U}}_0,{\mathsf{\theta }}_{{\mathrm{z}}_1},\dots {\mathsf{\theta }}_{{\mathrm{z}}_{\mathrm{N}-1}},{\mathsf{\theta }}_{{\mathrm{U}}_1},\dots ,{\mathsf{\theta }}_{{\mathrm{U}}_{\mathrm{N}-1}}\right\}$. This loss function is used since $Y_T=g\left(X_T\right)$. The expectation in the loss function in (28) is the expectation for all sample paths, but due to the infinite number of sample paths, it is not possible to traverse the entire training set. Therefore, we adopt an optimization method based on stochastic gradient descent. In each iteration of gradient descent, only a portion of samples is selected to estimate the loss function, thereby obtaining the gradient of the loss function over all parameters and obtaining a one-step neural network parameter update.

Like this, the back propagation of this neural network uses an optimizer to update the parameters layer by layer. When the DNN is trained, ${\widehat{Y}}_0$ is fixed as the parameter. Take this parameter out and it is the required value.

3.5. Details of the Algorithms

The detailed steps of our proposed algorithm based on deep learning to numerically solve the BSDE with jumps is presented as follows.

According to our algorithm, which uses a deep learning-based neural network solver, the BSDE with jumps in the form of Equation (4) can be numerically solved in the following ways:

• Simulate sample paths using standard Monte Carlo methods
• Use a deep neural network (DNN) to approximate $Z_t={\nabla }_{x}u\left(t,X_t\right)\sigma \left(t,X_t\right)$ and $U_t=u\left(t,X_t+\beta \left(t,X_t\right)\right)-u\left(t,X_t\right)$, and then plug them into the BSDE with jumps to perform a forward iterative operation in time

For simplicity, here we use the one-dimensional case as an example, and the high-dimensional case is similar. We divide the time interval $\left[0,T\right]$ into N partitions, so that $t_0,t_1,\dots ,t_N\in \left[0,T\right]$ and $0=t_0<t_1<\dots <t_N=T$, assuming $h_i=t_{i+1}-t_i$, $d{W}_i=W_{i+1}-W_i$ and $d{\tilde{N}}_i={\tilde{N}}_{i+1}-{\tilde{N}}_i$. The detailed steps are as follows:

(1)

N Monte Carlo paths of the diffusion process $\left\{X_i^{\left(j\right)}:i=0,1,\dots ,N; j=1,2,\dots ,M; X_0^{\left(j\right)}\equiv X_0\right\}$ are sampled by the Euler scheme:

$$X_{i+1}^{\left(j\right)}=X_i^{\left(j\right)}+\mu \left(t_{i },X_i^{\left(j\right)} \right)h_i+\sigma \left(t_{i },X_i^{\left(j\right)} \right)d{W}_i^{\left(j\right)}+\beta \left(t_i,X_i^{\left(j\right)}\right)d{\tilde{N}}_i{}^{\left(j\right)}$$

(29)

This step is the same as the standard Monte Carlo method. Other discretization schemes can also be used, such as logarithmic Euler discretization or Milstein discretization.

(2)

At $t_0=0$, the initial value $Y_0$, $Z_0=\frac{\partial Y}{\partial X}\left(0, X_0\right)\sigma \left(0, X_0\right)$ and $U_0=u\left(0,X_0+\beta \left(0,X_0\right)\right)-u\left(0,X_0\right)$ are randomly selected as parts of the neural network parameters, and $Y_0$, $Z_0$ and $U_0$ are all constant random numbers selected from an empirical range.

(3)

At each time step $t_i$, given $\left\{X_i^{\left(j\right)}: j=1,2,\dots ,M\right\}$, $\left\{Z_i^{\left(j\right)}: j=1,2,\dots ,M\right\}$ and $\left\{U_i^{\left(j\right)}: j=1,2,\dots ,M\right\}$ are approximated using deep neural networks, respectively. Note that every time $t=t_i$, two sub-neural networks are used for all Monte Carlo paths. In our one-dimensional case, as described in Section 3.4, we introduce two sub-deep neural networks: ${\mathsf{\theta }}_{{\mathrm{z}}_{\mathrm{i}}}: ℝ\to ℝ$, such that $\left\{Z_i^{\left(j\right)}={\mathsf{\theta }}_{{\mathrm{Z}}_{\mathrm{i}}}\left(X_i^{\left(j\right)}\right): j=1,2,\dots ,M\right\}$ and ${\mathsf{\theta }}_{{\mathrm{U}}_{\mathrm{i}}}: ℝ\to ℝ$, such that $\left\{U_i^{\left(j\right)}={\mathsf{\theta }}_{{\mathrm{U}}_{\mathrm{i}}}\left(X_i^{\left(j\right)}\right): j=1,2,\dots ,M\right\}$. ${\mathsf{\theta }}_Z=\left\{{\mathsf{\theta }}_{{\mathrm{Z}}_1},\dots ,{\mathsf{\theta }}_{{\mathrm{Z}}_{\mathrm{N}-1}}\right\}$, ${\mathsf{\theta }}_U=\left\{{\mathsf{\theta }}_{{\mathrm{U}}_1},\dots ,{\mathsf{\theta }}_{{\mathrm{U}}_{\mathrm{N}-1}}\right\}$, $\theta =\left\{{\mathsf{\theta }}_Z,{\mathsf{\theta }}_{\mathrm{U}}\right\}$.

(4)

For $t_i\in \left[0,T\right]$, we have

$$Y_{i+1}^{\left(j\right)}=Y_i^{\left(j\right)}-f\left(t_i,X_i^{\left(j\right)},Y_i^{\left(j\right)},Z_i^{\left(j\right)},U_i^{\left(j\right)}\right)h_i+\langle Z_i^{\left(j\right)},d{W}_i^{\left(j\right)}\rangle +\langle U_i^{\left(j\right)},d{\tilde{N}}_i{}^{\left(j\right)}\rangle$$

(30)

According to this formula, we can directly calculate the $Y_{i+1}^{\left(j\right)}$ at the next time step. This process does not require any parameters to be optimized. In this way, the BSDE with jumps propagates forward in time direction from $t_i$ to $t_{i+1}$. Along each Monte Carlo path, as the BSDE with jumps propagates forward in time from 0 to T, $Y_N^{\left(j\right)}$ is estimated as $Y_N^{\left(j\right)}\left(Y_0, Z_0,U_0,{\theta }^{\prime }\right)$, where ${\theta }^{\prime }=\left({\theta }_Z,{\theta }_U\right)$ is all the hyper-parameters for the N − 1 layers neural networks.

(5)

Calculate the following loss function:

$$L_{Forward}=\frac{1}{M}{\displaystyle \sum }_{j=1}^M{\left(Y_N^{\left(j\right)}\left(Y_0, Z_0,U_0,\theta \right)-g\left(X_N^{\left(j\right)}\right)\right)}^2$$

(31)

where g is the terminal function.

(6)

Use stochastic optimization algorithms to minimize the loss function:

$$\left(\widetilde{Y_0}, \widetilde{z_0},\tilde{\theta }\right)=\underset{\left(Y_0, z_0,\theta \right)}{\mathrm{argmin}}\frac{1}{M}{\displaystyle \sum }_{j=1}^M{\left(Y_N^{\left(j\right)}\left(Y_0, z_0,\theta \right)-g\left(X_N^{\left(j\right)}\right)\right)}^2$$

(32)

The estimated value $\widetilde{Y_0}$ is the desired initial value at time t = 0.

4. Numerical Results

In this section, we will use the deep neural network to code the theoretical framework of the proposed algorithm. Many problems in the financial field can be solved through numerical solutions using high-dimensional PDEs and corresponding BSDEs. Among them, the jump process can depict the sudden impact of the outside world in its process, so that it can accurately depict a class of uncertain things, and this situation is most common in financial markets. So, we extend the numerical solution problem of a high-dimensional PDE to the case of jump diffusion. In this paper, three classical high-dimensional PDEs with important applications in finance-related fields are selected for numerical simulation. They are financial derivative pricing under the jump-diffusion model, Bond pricing under the Vasicek model with jump, and the Hamilton-Jacobi-Bellman equation.

The numerical experiments in this paper use a 64-bit Windows 10 operating system, based on the PyTorch deep learning framework in Python. All examples in this section are calculated based on 1000 sample paths, each chosen for a time partition of N = 100. We ran each example 10 times independently to calculate the average result. All neural network parameters are initialized by uniform distribution. Each sub-neural network of each time node contains 4 layers: one d-dimensional input layer, two (d + 10)-dimensional hidden layers and one d-dimensional output layer. We use the rectifier function (ReLU) as the activation function. Batch Normalization (BN) techniques are used after each linear transformation and before activation.

4.1. Pricing of Financial Derivatives under Jump Diffusion Model

In this section, we will apply the proposed method to the pricing of derivatives related to a 100-dimensional jump diffusion model. In this model, the stock price $X_t$ satisfies the following jump diffusion model [19]:

$$d{X}_t=\left(r-\lambda k\right)X_{t}dt+\sigma X_{t}d{W}_t+\left(V-1\right)X_{t}d{\tilde{N}}_t$$

(33)

where $r$ is the constant discount rate, $\lambda$ is the average number of jumps per unit time of the stock price and $\sigma$ is the constant volatility. $V$ represents the jump magnitude. Assuming that the jump amplitude is fixed and V is constant, we can let $k=V-1$ and $k>-1$. This equation can be written as $d{X}_t=r{X}_{t}dt+\sigma X_{t}d{W}_t+k{X}_{t}d{\tilde{N}}_t$.

It is known that the European call option with the stock as the underlying asset satisfies the following partial differential equation:

$$\frac{\partial u}{\partial t}\left(t,x\right)+\frac{1}{2}{\sigma }^2{x}^2{\Delta }_{x}u\left(t,x\right)+\left(r-\lambda k\right)x{\nabla }_{x}u\left(t,x\right)+\lambda \left[u\left(t,xV\right)-u\left(t,x\right)\right]-ru\left(t,x\right)=0$$

$$u\left(T,x\right)={\left(x-K\right)}^{+}$$

(34)

For all $t\in \left[0,T\right]$, $x,\omega \in {ℝ}^d$, $y\in ℝ$, $z\in {ℝ}^{1\times d}$. Suppose $d=100$, $T=1$, $N=100$, $\lambda =1$, $X_0=\left(1,\dots ,1\right)$, $\mu \left(t,x\right)=r{X}_t$, $\sigma \left(t,x\right)=\sigma X_t$, $\beta \left(t,x\right)=\left(V-1\right)X_t=k{X}_t$, $f\left(t,x,y,z\right)=-rY$, $g\left(x\right)={\left(\underset{1\le i\le 100}{\max }{x}_i-5\right)}^{+}$, here $r^i$, ${\sigma }^i$, $k^i$ are randomly selected in $\left[0,1\right]$, $i=1,\dots 100$.

After calculation,

Table 1. Numerical results with Adam optimizer.

Number of Iteration Steps n	$\mathbf{Mean} \mathbf{of} \mathit{u}\left(0,{\mathit{X}}_0\right)$	$\mathbf{Standard} \mathbf{Deviation} \mathbf{of} \mathit{u}\left(0,{\mathit{X}}_0\right)$	Mean of the Loss Function	Standard Deviation of the Loss Function	Runtime in Second
1000	0.5989	0.1719	1.2743	7.7509	474
2000	0.4664	0.1235	1.2532	4.3095	891
3000	0.4183	0.1623	1.0894	3.7375	1271
4000	0.3940	0.1465	0.6269	2.4771	1676

and

Table 2. Numerical results with NAdam optimizer.

Number of Iteration Steps n	$\mathbf{Mean} \mathbf{of} \mathit{u}\left(0,{\mathit{X}}_0\right)$	$\mathbf{Standard} \mathbf{Deviation} \mathbf{of} \mathit{u}\left(0,{\mathit{X}}_0\right)$	Mean of the Loss Function	Standard Deviation of the Loss Function	Runtime in Second
1000	0.2821	0.0625	1.1528	2.7450	549
2000	0.3019	0.0484	0.6747	2.0031	828
3000	0.3084	0.0406	0.4911	1.6566	1206
4000	0.3117	0.0356	0.4134	1.4535	1651

show the important numerical results of solving the related derivative pricing problem of the 100-dimensional jump diffusion model with Adam and NAdam optimizers, respectively, including that with the change of iteration steps, mean and standard deviation of loss function $Y_0$, mean and standard deviation of loss function and running time. Only the numerical results of iteration steps $n\in \left\{1000, 2000, 3000, 4000\right\}$ are selected as typical examples for display. We set the initial value $u\left(0,X_0\right)=0.3124$ obtained from the classical Monte Carlo method as the exact solution and compare it with the numerical results.

It can be seen from Table 1 and Table 2 that when the iteration steps are the same, the calculation time required for using the two optimizers is not very different. As the iteration progresses, the loss function value of the NAdam optimizer is smaller than that of the Adam optimizer at the same iteration steps and the final loss function value is also smaller.

Figure 2 and Figure 3 show the changes of the initial value estimate and loss function with the number of iteration steps when Adam and NAdam optimizers are used for solving. It can be seen that the initial value estimates and loss functions of the numerical solutions of the two methods converge rapidly. However, the convergence speed of the latter is obviously faster than that of the former, whether it is the initial value estimate or the loss function. Among them, the former tends to converge after about 1500 iterations, while the latter tends to converge after about 700 iterations. When the iteration is 4000 times, the numerical solution obtained with the Adam algorithm is 0.32143775, while the numerical solution obtained with the NAdam algorithm is 0.32273737, which is similar and not significantly different from the exact value simulated by the Monte Carlo method.

4.2. Bond Pricing under the Jumping Vasicek Model

In this section, we use the proposed method to solve the pricing problem of a class of bonds with interest rates subject to the Vasicek model with jumps [20,21,22]. In this model, the short-term interest rate $X_t$ obeys the following stochastic differential equations:

$$d{X}_t=a\left(b-X_t\right)dt+\sigma d{W}_t+kd{\tilde{N}}_t$$

(35)

i.e., each part of $X_t$ follows:

$$d{X}_t^i=a^i\left(b^i-X_t^i\right)dt+{\sigma }^{i}d{W}_t+k^{i}d{\tilde{N}}_t$$

(36)

For all $t\in \left[0,T\right]$, $x,\omega \in {ℝ}^d$, $y\in ℝ$, $z\in {ℝ}^{1\times d}$. Suppose $d=100$, $T=1$, $N=100$, $\lambda =1$, $X_0=\left(1,\dots ,1\right)$, $\mu \left(t,x\right)=\left({\mu }^1,\dots ,{\mu }^i,\dots ,{\mu }^d\right)$, ${\mu }^i=a^i\left(b^i-x^i\right)$, $\sigma \left(t,x\right)=\left({\sigma }^1,\dots ,{\sigma }^i,\dots ,{\sigma }^d\right)$, $\beta \left(t,x\right)=\left({\beta }^1,\dots ,{\beta }^i,\dots ,{\beta }^d\right)$, $a^i,b^i,{\sigma }^i,$ ${\beta }^i$ are randomly selected in $\left[0,1\right]$, $f\left(t,x,y,z\right)=-\left(\underset{1\le i\le n}{\max }{X}^i\right)Y$, $g\left(x\right)=1$. Then, for the zero-coupon bond price $u$ paying 1 at maturity T under the above jump Vasicek model, it satisfies that $u\left(T,x\right)=1$ for all $t\in \left[0,T\right]$, $x\in {ℝ}^d$, and the following partial differential equation holds:

$$\begin{gathered} \frac{\partial u}{\partial t}\left(t,x\right)+{\displaystyle \sum }_{i=1}^n\left[a^i\left(b^i-x^i\right)-{\beta }^i\right]\frac{\partial u}{\partial x^i}+\frac{1}{2}{\displaystyle \sum }_{1\le i,j\le n}{\sigma }^i{\sigma }^j\sqrt{x^i{x}^j}\frac{{\partial }^{2}u}{\partial x^i\partial x^j}+u\left(t,x+\beta \left(x\right)\right)-u\left(t,x\right)-\left(\underset{1\le i\le n}{\max }{X}^i\right)u \\ =0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \end{gathered}$$

(37)

After calculation,

Table 3. Numerical results with Adam optimizer.

Number of Iteration Steps n	$\mathbf{Mean} \mathbf{of} \mathit{u}\left(0,{\mathit{X}}_0\right)$	$\mathbf{Standard} \mathbf{Deviation} \mathbf{of} \mathit{u}\left(0,{\mathit{X}}_0\right)$	Mean of the Loss Function	Standard Deviation of the Loss Function	Runtime in Second
1000	0.2096	0.0435	0.1405	0.4825	538
2000	0.2300	0.0370	0.0820	0.3464	924
3000	0.2396	0.0331	0.0585	0.2849	1296
4000	0.2450	0.0302	0.0463	0.2476	1779

and

Table 4. Numerical results with NAdam optimizer.

Number of Iteration Steps n	$\mathbf{Mean} \mathbf{of} \mathit{u}\left(0,{\mathit{X}}_0\right)$	$\mathbf{Standard} \mathbf{Deviation} \mathbf{of} \mathit{u}\left(0,{\mathit{X}}_0\right)$	Mean of the Loss Function	Standard Deviation of the Loss Function	Runtime in Second
1000	0.2828	0.0293	0.0323	0.1310	652
2000	0.2772	0.0215	0.0224	0.1033	1326
3000	0.2757	0.0177	0.0172	0.0855	1957
4000	0.2750	0.0154	0.0150	0.0753	2581

show the important numerical results of solving a class of bond pricing problems with interest rates subject to the Vasicek model with jumps using the Adam and NAdam optimizers, respectively, including that with the change of iteration steps, mean and standard deviation of loss function $Y_0$, mean and standard deviation of loss function and running time. Only the numerical results of the iteration steps $n\in \left\{1000, 2000, 3000, 4000\right\}$ are selected as typical examples for display. We set the initial value $u\left(0,X_0\right)=0.2702$ obtained from the classical Monte Carlo method as the exact solution and compare it with the numerical results.

It can be seen from Table 3 and Table 4 that with the iteration, the loss function value of the NAdam optimizer is smaller than that of the Adam optimizer at the same iteration steps, and the final loss function value is smaller, but the cost is that the operation time becomes longer.

Figure 4 and Figure 5 show the changes of the initial value estimate and loss function with the number of iteration steps when Adam and NAdam optimizers are used for solving. It can be seen that the initial value estimates and loss functions of the numerical solutions of the two methods converge rapidly. However, the convergence speed of the latter is obviously faster and more stable than the former, whether it is the initial value estimate or the loss function. When the iteration is 4000 times, the numerical solution obtained with the Adam algorithm is 0.25530547, while the numerical solution obtained with the NAdam algorithm is 0.27403888, which is similar and not significantly different from the exact value simulated by the Monte Carlo method.

4.3. Hamilton-Jacobi-Bellman (HJB) Equation

In fields such as finance, investment, and risk management, optimization and control problems are often involved, and these problems are often represented by stochastic optimal control models. One way to solve this kind of control problem is to obtain and solve the Hamilton-Jacobi-Bellman equation (HJB equation for short) of the corresponding control problem according to the principle of dynamic programming. In this section, we use the proposed method to solve a class of 100-dimensional HJB equations [12]:

For all $t\in \left[0,T\right]$, $x,\omega \in {ℝ}^d$, $y\in ℝ$, $z\in {ℝ}^{1\times d}$. Suppose $d=100$, $T=1$, $N=100$, $\lambda =1$, $X_0=\left(1,\dots ,1\right)$, $\mu \left(t,x\right)=0$, $\sigma \left(t,x\right)=\sqrt{2}$, $\beta \left(t,x\right)={\beta }^{T}x$, where ${\beta }^T=\left({\beta }^1,\dots ,{\beta }^i,\dots ,{\beta }^d\right)$, ${\beta }^i$ are randomly selected in $\left[0,1\right]$, $f\left(t,x,y,z\right)=-{\Vert z\Vert }_{{ℝ}^{1\times d}}^2$, $g\left(x\right)=ln\left(\left(1+{\Vert x\Vert }^2\right)/2\right)$, it satisfies that $u\left(T,x\right)=ln\left(\left(1+{\Vert x\Vert }^2\right)/2\right)$ for all $t\in \left[0,T\right]$, $x\in {ℝ}^d$, and the following partial differential equation holds:

$$\frac{\partial u}{\partial t}\left(t,x\right)+{\Delta }_{x}u\left(t,x\right)+u\left(t,x+\beta \left(x\right)\right)-u\left(t,x\right)-{\nabla }_{x}u\left(t,x\right)\beta \left(x\right)=\Vert {\nabla }_{x}u{\left(t,x\right)\Vert }^2$$

(38)

After calculation, the important numerical results of solving a class of HJB equations with Adam and NAdam optimizers are shown in

Table 5. Numerical results with Adam optimizer.

Number of Iteration Steps n	$\mathbf{Mean} \mathbf{of} \mathit{u}\left(0,{\mathit{X}}_0\right)$	$\mathbf{Standard} \mathbf{Deviation} \mathbf{of} \mathit{u}\left(0,{\mathit{X}}_0\right)$	Mean of the Loss Function	Standard Deviation of the Loss Function	Runtime in Second
2000	1.2504	0.1521	1.3142	1.4709	857
4000	1.4850	0.2820	0.9090	1.1166	1770
6000	1.7704	0.4768	0.7447	0.9412	2906
8000	2.0139	0.5904	0.6314	0.8386	3864

and

Table 6. Numerical results with NAdam optimizer.

Number of Iteration Steps n	$\mathbf{Mean} \mathbf{of} \mathit{u}\left(0,{\mathit{X}}_0\right)$	$\mathbf{Standard} \mathbf{Deviation} \mathbf{of} \mathit{u}\left(0,{\mathit{X}}_0\right)$	Mean of the Loss Function	Standard Deviation of the Loss Function	Runtime in Second
2000	0.8762	0.1941	1.2294	1.5767	1649
4000	1.2456	0.4419	0.8682	1.1733	3104
6000	1.6858	0.7241	0.6715	0.9980	4527
8000	1.9266	0.7531	0.5606	0.8855	5948

, respectively, including that with the change of iteration steps, mean and standard deviation of loss function $Y_0$, mean and standard deviation of loss function, and running time. Only the numerical results of iteration steps $n\in \left\{2000, 4000, 6000, 8000\right\}$ are selected as typical examples for display. We set the initial value $u\left(0,X_0\right)=2.6119$ obtained from the classical Monte Carlo method as the exact solution and compare it with the numerical results.

It can be seen from Table 5 and Table 6 that with the iteration, the loss function value of the NAdam optimizer is smaller than that of the Adam optimizer at the same iteration steps, and the final loss function value is smaller, but the operation time is longer.

Figure 6 and Figure 7 show the changes of the initial value estimate and loss function with the number of iteration steps when Adam and NAdam optimizers are used for solving. It can be seen that the convergence speed of the latter is obviously faster than that of the former, whether it is the initial value estimate or the loss function. Among them, the former tends to converge after about 7000 iterations, while the latter tends to converge after about 5500 iterations. After 8000 iterations, the numerical solution obtained with the Adam algorithm is 2.73136686, while the numerical solution obtained with the NAdam algorithm is 2.6456454, which is similar and not significantly different from the exact value simulated by the Monte Carlo method.

5. Conclusions

In this paper, we propose an algorithm that can be used to solve a class of partial differential equations with jump terms and their corresponding backward stochastic differential equations. Through the nonlinear Feynman-Kac formula, the above-mentioned high-dimensional nonlinear PDE with jumps and its corresponding BSDE can be expressed equivalently, and the numerical solution problem can be regarded as a stochastic control problem. Next, we treat the gradient and jump process of the unknown solution separately as policy functions and use two neural networks at each time division to approximate the gradient and jump process of the unknown solution respectively. In this way, we can use deep learning to overcome the “curse of dimensionality” caused by high-dimensional PDE with jumps and obtain numerical solutions. In addition, we attempt to replace the traditional Adam algorithm with the new stochastic optimization approximation algorithm and apply it to our algorithm.

Next, focusing on the financial field, we applied our algorithm to solve three common high-dimensional problems in finance-related fields and then compared the results. We concluded that our algorithm performs well in numerical simulation. It is also concluded that after applying the new optimizer to the deep learning algorithm, the convergence speed of the model is mostly faster, and the generalization ability is significantly improved, at the cost of the operation time possibly being longer.

The algorithm proposed in this paper can be used for a wide variety of problems in finance, insurance, actuarial modeling and other fields, such as the pricing problem of some special path-dependent financial derivatives, the bond pricing problem under the jump diffusion model of interest rates, the pricing of equity-linked life insurance contracts and so on.

6. Discussion and Limitations

The algorithm proposed in this article only targets a special class of nonlinear parabolic PDEs with jump terms, and its application objects have certain limitations. We can consider improving the algorithm to extend the range of application objects, such as other special forms of PDE. For example, the improved deep learning-based algorithm can be extended to solve reflected PDEs, using the penalty function method to approximate reflected PDEs to a general form of PDEs. We can also consider delayed PDEs, restricted PDEs, and so on. In addition, this paper only considers the most common Poisson process, not the general jump process. We hope to apply this algorithm to other jump processes in the future. In addition, this paper does not focus on the changes in the intensity of the jump process. We hope to extend this algorithm to the general integro-partial differential equations in the future.

This paper mainly focuses on numerical simulations of examples in finance. In the future, more cases can be considered to test the practical effect of the proposed algorithm in other jump-diffusion instances. In addition, this paper only introduces the specific steps of the algorithm but does not analyze the convergence, error estimation and robustness of stochastic simulation, and thus, we hope to improve it in future research.