Deep Learning of Transition Probability Densities for Stochastic Asset Models with Applications in Option Pricing

1. Introduction

In computational finance, such common tasks as option pricing and estimating sensitivities (Andricopoulos et al. 2003, 2007; Glasserman 2003; Chen et al. 2014; Su and Newton 2020; and references therein) and simulating likelihood estimators (Aït-Sahalia 2002, Aït-Sahalia and Kimmel 2007, Yu 2007, Giesecke and Schwenkler 2019) greatly benefit from access to fast and accurate evaluation of transition probability density functions (TPDFs) for stochastic differential equations (SDEs) modeling asset prices.

It is only in a few cases that we are able to find a TPDF in closed form; thus, we rarely have an exact formula which does not require numerical methods to calculate the TPDF. The well-known case of geometric Brownian motion, which is behind the Black-Scholes equation for European options, is one such with a closed form TPDF. In the absence of a closed form solution, sometimes we are able to obtain a semiclosed form for the TPDF. These cases require additional calculations such as use of Bessel functions or calculations of integrals in the complex plane. For example, the CEV process suggested by Cox (1996) has a semiclosed form TPDF (Chen et al. 2014). Lewis (2016) presented the semiclosed form of the joint transition density for the stochastic volatility process of Heston (1993) and a 3/2-stochastic volatility process, a variation of the Heston process (Heston 1997, Lewis 2016). In certain cases, TPDFs are not known in closed form but their characteristic functions are. Densities can subsequently be calculated using an inverse Fourier transform (O’Sullivan 2005, Lord et al. 2008, Su et al. 2017). For the majority of processes used for modeling asset prices, where we have no other ways to calculate density, we can turn to density approximations, as introduced by Chen et al. (2014).

There are several directions for obtaining density approximations described in the literature. First is the class of expansion methods: Aït-Sahalia (2002, 2008) and Aït-Sahalia and Kimmel (2007) rely on small-time asymptotic expansions to solve multivariate diffusion processes. This direction is subsequently extended in Yu (2007) to solve multivariate jump-diffusions; Henry-Labordère (2008) and Henry-Labordere et al. (2017) use different expansion approaches for multivariate diffusion processes that do not rely on small-time expansion; and the polynomial expansion approach suggested by Filipović et al. (2013) is able to evaluate the density for affine jump-diffusion. Another class of density approximation explores exact sampling based on Monte Carlo simulation, first developed by Beskos and Roberts (2005) and Chen and Huang (2013) and then extended in a number of papers including Giesecke and Smelov (2013) and Giesecke and Schwenkler (2019). The latest development in this direction, from Guay and Schwenkler (2021), gives an unbiased density estimator that can be used to tackle nonaffine multivariate jump diffusion processes. In addition to these two major classes of approximation techniques, there are other techniques available; for example, the saddlepoint approximation approach suggested by Aït-Sahalia and Yu (2006) that is based on the characteristic function.

Nevertheless, density approximations have limitations, an important one being that they cannot universally provide approximations for all underlying processes. For example, the approximation techniques used in Chen et al. (2014) cannot solve for jump-diffusion processes and the density approximation is only accurate when the time step is very small. On top of this, density approximations can be cumbersome and may require complex computations. Thus, universality across any of the models of asset prices is available for those willing to put in the necessary effort but its practicality is more limited than we would wish. It is true that conventional numerical methods, particularly finite differences, can also be used to solve for the TPDF (Floc’h and Kennedy 2014, Hagan et al. 2014, Su and Newton 2020), but they are computationally expensive. More importantly, these only allow calculation for a single set of model parameters at a time, so they are far from ideal. A faster and more generic way of calculating density is needed for practical universal application. This leads us to deep learning.

The use of deep learning in option pricing has a long history, dating back at least as far as the early 1990s (Malliaris and Salchenberger 1993, Hutchinson et al. 1994) and has built a large technical literature of applications and improvements aimed at creating more reliable setups linking underlying data directly with option prices in the markets and in calibrating implied volatilities. In contrast, we apply deep learning to very accurately approximate the TPDF, which can then be used (among other applications) in option pricing. The difficulties in extracting densities from market data are considerable, particularly in the tails (Figlewski 2009, 2018) but that is not our aim. Rather, we use deep learning in place of other mathematical techniques in the calculation of TPDFs, training artificial neural networks (NNs) on simulated data according to various models of assets underlying the options. This approach offers the potential to obtain TPDFs to a high accuracy more quickly and easily than previously, via pretrained networks. A novel advance in our approach is in the parametric solution of the related partial differential equations (PDEs). It is usual in option pricing to solve a PDE for fixed parameters but then to have to recalculate when changes are made (e.g., within a calibration procedure). Here we create an engine that does not require us to solve the PDE again.

Deep learning algorithms offer an effective way of solving ordinary differential equations and PDEs. Prototypes of algorithms as differential equation solvers may be found in the work of Lee and Kang (1990) and Dissanayake and Phan-Thien (1994). Subsequently, similar feedforward approaches have been introduced by van Milligen et al. (1995), Lagaris et al. (1998), and Lagaris et al. (2000). Schmidhuber (2015) and Yadav et al. (2015) review the area. Calculations are very fast once the learning networks are trained but they suffer the drawback of needing to sample the data using a mesh. Mesh-free approaches to solving PDE problems based on deep learning are presented by Sirignano and Spiliopoulos (2018) (see also Rackauckas et al. (2020) and references therein). Here, we adapt the “deep Galerkin method” (DGM) deep learning algorithm developed by Sirignano and Spiliopoulos for computing densities of underlying processes. Through this algorithm we arrive at an approximate PDE solution represented by a deep neural network, which is trained to satisfy the differential operator, the initial condition, and the boundary conditions. The key for this work is that deep learning algorithms for PDEs can be extended to parametric PDE problems (Khoo et al. 2017, Kutyniok et al. 2019, Geist et al. 2020). By parametric problems, we understand that we are interested not just in solving a PDE in a domain of its independent variables (time and space), as standard PDE solvers do, but also in having a representation of a PDE’s solution for a range of its parameters’ values. There are several types of methods under the umbrella of model order reduction capable of dealing with parametric PDE problems (Antoulas et al. (2015) and references therein). The type used in this paper, based on deep learning networks, has been shown over recent years to be very robust and universal (Khoo et al. 2017, Kutyniok et al. 2019, Geist et al. 2020; also its use in this paper). Parametric PDE problems arise naturally in finance. Indeed, in finance applications, we need to be able to find option prices, densities for the underlying price process, etc. not only as functions of spot price and maturity but also as functions of volatility and other parameters of the underlying process to satisfy the practical need for frequent recalibration of models. Frequent recalibration imposes the requirement for valuations to be very fast for a range of parameter values. Deep learning applied to parametric PDE problems can achieve this, as we shall demonstrate. We emphasize that standard PDE solvers (e.g. finite differences) give approximations of PDE solutions only for given, fixed values of parameters and need to be run online each time a parameter changes, making them uncompetitive for the parametric problem described above. An important feature of using deep learning to solve parametric problems, as applied to the evaluation of TPDFs, is that it is “single solve” such that, once trained, the algorithm delivers results for any set of parameters without the need to rerun calculations.

These advantages are the prize. Nevertheless, formulating the use of deep learning techniques to calculate transition probability densities is a challenge. We know that the solution of the Fokker-Planck (Kolmogorov forward) equation is the transition probability density. Its initial condition is a Dirac delta function, which has zero value everywhere except at one point where it is infinite. The presence of the delta function in the initial condition makes it impractical to use deep learning methods for the Fokker-Plank equation. Al-Aradi et al. (2018) make an attempt to solve the Fokker-Planck equation with DGM for the Ornstein-Uhlenbeck process by assuming a Gaussian initial condition. Pleasing as this is, it cannot be generalized to other underlying processes simply because other processes do not possess Gaussian shapes. Here, we circumvent this difficulty by considering instead the backward Kolmogorov equation with the terminal condition being a step function, whose solution is the cumulative distribution function. Consequently, we propose a different direction by employing deep learning to solve the governing backward Kolmogorov equation for the cumulative probability function of an underlying stochastic process, thereby retrieving the transition density through differentiating the approximate cumulative probability. The governing backward Kolmogorov equation exists for any Markovian stochastic process meaning that, in theory, the deep learning solver is able to tackle any asset price model described by SDEs.

To illustrate the computational efficiency and accuracy of the proposed neural TPDF generators, we use them as a new “engine” for QUAD methods for option pricing. QUAD calculates option prices through integrating the product of the payoff function and the TPDF and does not require massive calculation between the key points in time that specify the option type (Andricopoulos et al. 2003, 2007). Lack of available closed form or accurate approximations of TDPFs limited the range of underlying processes that could be tackled by QUAD methods for derivatives pricing. Chen et al. (2014) finessed this through insertion of interchangeable approximation techniques. This works well but a limitation is the accuracy of the TPDF approximation engines, leading to a computational bottleneck for the more practically interesting models of the underlying. The neural TDPF generators developed in our paper not only remove the bottleneck but, more significantly, innovate the deep learning route by being “single solve,” across all option pricing techniques using the TPDF: There is no need to recalculate densities and the network rapidly returns results for any fresh inputs. We demonstrate the new approach starting with geometric Brownian motion (the classical Black-Scholes model), then the Heston model, the SABR model (a nonaffine case), Kou’s double exponential jump-diffusion model, and stochastic volatility jump diffusion. In the case of jump-diffusion models, we show how to use the deep learning approach to parametrically solve partial integral differential equations (PIDEs) to evaluate the corresponding TPDF. By doing so, we simultaneously contribute to the literature on deep PDE solvers by extending it to solving PIDEs. We also advance the use of deep learning for solving parametric PDE/PIDE problems by extensive computational investigation of the choices of loss functions, NN hyperparameters and NN architecture, as well as testing the performance of various GPU cards.

Our paper introduces a novel numerical approach for solving TPDFs, offering a fast and universal approximation method. We pioneer the utilization of deep learning in successfully addressing this issue. Furthermore, we comprehensively assess the accuracy of neural TPDFs through option pricing results with QUAD, thereby enhancing our understanding of the implementation of deep learning techniques for solving PDEs.

2. Methods

We set the scene by expressing option prices via TPDFs of the underliers because we will test the performances of NN approximated densities in the application setting of option pricing. For clarity of exposition, we begin with a simple one-dimensional model. In Section 2.2, we show how to formulate PDE problems to find TPDFs. In Section 2.2.2, we apply deep learning techniques to approximating parametric PDEs and using this to solve backward Kolmogorov equations to find parametric TPDFs in a computationally effective manner.

2.1. Option Pricing via the Transition Density of an Underlier

Before moving on to financially more interesting cases, we consider the simplest model where a single underlying, $S(s),$ follows a stochastic differential equation (SDE) written under the forward measure:

Here, $W(s)$ is a one-dimensional standard Wiener process and the volatility $\sigma (s,x)$ is a deterministic function. For example, the famous Black-Scholes option pricing formula for European style options rests on the model in which the asset’s price follows geometric Brownian motion:

Neglecting the discounting factor, the European option price $V(0,S_0)$ at time $t=0$ on an underlying $S_{0,S_0}(s)$ with the spot price $S_{0,S_0}(0)=S_0$ and with a payoff function $f(x)$ and maturity T is expressed by

Because the solution of SDE (1) is a Markov process, we can write

where $p(t,S_0;T,y)$ is the TPDF for the random variable $S_{t,S_0}(T)$ that is a solution of (1) at time T starting from the point $S_0$ at time t.

Computationally, although working in the asset price space is viable, it is preferable to use the log-asset price space since this change of variables allows us to stretch the phase space, that is, the left bound of the integration becomes $-\infty$. Consequently, this change of variables is beneficial for training a NN for the backward Kolmogorov equation, as we will discuss later (see Section 2.2.2). By Ito’s formula, we obtain the SDE in the log-space, that is, for $X(s)=\ln S(s)$:

After this change of variable, the pricing Equation (4) becomes

where $p(t,x;T,y)$ is the TPDF for the solution $X_{t,x}(T)$ of the SDE (5) at time T starting from x at time t. When the volatility $\sigma$ is time independent, without loss of generality, we can start the solution $X(s)$ at time $t=0$ and use the simpler notation $p(x;s,y)$ or in short $p(s,y)$ instead of $p(0,x;s,y)$ that can also be used when we aim to find an option price at time $t=0$.

Thus, in Equation (6), if we know the TPDF $p(t,x;T,y)$, we can readily calculate the price of the option, $V(0,S_0)$. To this end, the integral in Equation (6) can be approximated via quadrature after appropriate truncation of the integration range $[-\infty ,+\infty ]$. Details of implementation for a given density can be found in the early papers in the QUAD series (Andricopoulos et al. 2003, 2007).

2.2. Preliminaries on Transition Densities

It is straightforward to write the Fokker-Planck (forward Kolmogorov) PDE for a TPDF. However, its initial condition contains the Dirac delta function that hinders the effective application of numerical techniques including the deep learning approach we present in this paper. We circumvent this problem by proposing a novel approach through solving the backward Kolmogorov PDE for the cumulative distribution function (CDF) and then differentiating it to obtain the TPDF.

2.2.1. Forward Equation for Transition Density.

The Fokker-Planck (forward Kolmogorov) equation governs the time evolution of the transition probability for SDEs’ solutions (Gichman and Skorochod 1972, Freidlin 1985, Gardiner 2004). In the case of the SDE (5), it takes the form

with initial condition for a fixed xwhere $\delta (·)$ is the Dirac delta function. Recall that the full notation for the TPDF is $p(s,y)=p(0,x;s,y)$ and note that the PDE problem (Equations (7)–(8)) is with respect to the variables $(s,y)$, whereas x is a parameter. Also, the Dirac delta function $\delta (y-x)$ is a function with zero value everywhere except at x:and ${\int }_{-\infty }^{+\infty }\delta (y-x)dy=1$. Because the function is infinite at the point $y=x$, it is difficult to solve the problem (7)–(8) numerically. This limitation extends to approximation via deep learning. Hence, we need to find a way to deal with it. To this end, instead of solving for TPDF directly, we first solve for the CDF and then use the CDF to calculate the TPDF. This method, as we will show, is efficient in dealing with the evaluation of TPDFs.

2.2.2. Backward Kolmogorov Equation for Cumulative Distribution Function.

We define the CDF of the log-stock price process $X_{t,x}(T)$ at terminal time T and the starting point $(t,x)$ as

where $\text{Prob}(X_{t,x}(T)\le y)$ is the probability of the random variable $X_{t,x}(T)$ not being greater than y.

In the case of the SDE (5), the function $C(t,x)=C(t,x;T,y)$ satisfies the backward Kolmogorov equation (Gichman and Skorochod 1972, Freidlin 1985, Gardiner 2004):

with the terminal condition for a fixed ywhere $C(T,x)=C(T,x;T,y)$. The problem (11)–(12) is written with negative direction in time, whereas the Fokker-Plank Equation (7) has positive direction in time. The PDE problem (11)–(12) is with respect to the variables $(t,x)$, and T and y are parameters. It is easier to solve numerically Equations (11)–(12) than Equations (7)–(8) because the terminal condition (12) contains the bounded function.

Using a deep NN, we approximate $C(t,x;T,y)$ and then, using an automatic differentiation built-in function from a deep learning library, such as TensorFlow, or finite differences, we evaluate the transition probability density $p(t,x;T,y)$ by differentiating $C(t,x;T,y)$ with respect to y; that is

For asset price models of practical interest, the CDF $C(t,x;T,y)$ is a smooth function for $t<T$; hence, the differentiation in Equation (13) is well defined. Consequently, numerical differentiation of an approximate $C(t,x;T,y)$ can be done accurately. We discuss this in more detail in Section 3.

The solution to the problem (11)–(12) depends not only on the independent variables $(t,x)$ but also on the set of parameters T, y and the volatility $\sigma (t,x)$, which can be parametrically encoded. A universal option pricing engine (as well as evaluation of sensitivities) requires an engine which can quickly and sufficiently accurately produce the TPDF for a range of values of the parameters so that we do not need to solve a PDE again, for example, after recalibration of the volatility. Thus, we need to solve a parametric PDE problem, not just for fixed values of parameters (and the variables $(t,x)$) but for a set of them. Traditional numerical methods (finite difference, finite element, Monte Carlo, etc.) are designed to solve PDE problems with fixed parameters. In contrast, model order reduction methods for PDEs (Antoulas et al. 2015), and the most modern and universal of them, based on deep learning (Khoo et al. 2017, Kutyniok et al. 2019, Geist et al. 2020), are aimed at parametric PDE problems.

2.3. Deep Learning as a PDE Solver

Recent developments in deep learning allow us to solve parametric PDE problems effectively. In Section 2.3.1, we present the key ingredients of the DGM proposed in Sirignano and Spiliopoulos (2018) for solving PDEs with fixed parameters. In Section 2.3.2, we recall the fundamentals of deep learning and explain how to extend DGM to solving parametric problems. Implementation of DGM methods requires us to have a NN of suitable architecture. We find that even a feedforward network, the simplest NN architecture, can be used successfully for neural TPDF generators, although here, for consistency of the presentation, we prefer to use the NN architecture introduced together with DGM in Sirignano and Spiliopoulos (2018) (henceforth we refer to it as the DGM NN), which we find effective for our purposes (see Section 2.3.3). In Online Appendix B, we compare the performance of the feedback forward NN (its simplest version: the multilayer perceptron) with the DGM NN.

2.3.1. Constructing the Loss Function.

We start with a formal description of DGM used to solve a PDE problem under fixed parameters, which we later extend into solving a parametric PDE problem. Let $G\subset {\mathbb{R}}^d$ be a domain with the boundary $\partial G$ and $\mathcal{L}$ be an elliptic operator. Consider the parabolic PDE

with terminal conditionand boundary condition

We assume that this problem has a classical solution $u(t,\mathit{x})$ (Ladyzhenskaya et al. 1968).

Equations (11)–(12) are of the form of Equations (14)–(15) and the backward Kolmogorov equation for any SDE driven by Wiener process has this form.

Next, we would like to use deep learning algorithms to “learn” the function $u(t,\mathit{x})$ and, hence, approximate the solution to the PDE problem (14)–(15). In other words, we wish to train a deep NN $f(t,\mathit{x};\mathit{\theta })$ to approximate $u(t,\mathit{x})$, where $\mathit{\theta }$ are trainable parameters of the NN.

To train the NN, we need to construct a loss (objective) function. In machine learning and deep learning, the loss function quantifies the performance of the approximation. In the PDE context, the loss function measures the fit (in other words, the residual) of the differential operator (i.e., the left-hand side of the Equation (14)), the terminal condition (15) and the boundary condition (16). The loss function $L(f)$ used in this paper has the form

where $L_1(f)$ is a loss in the differential operator term (the residual) and $L_2(f)$ is a loss in the terminal condition term. Here ${\Vert f(y)\Vert }_{\mathcal{Y},\nu }^2$ is the weighted $L^2$-norm defined as ${\Vert f(y)\Vert }_{\mathcal{Y},\nu }^2={\int }_{\mathcal{Y}}{|f(y)|}^2\nu (y)dy,$ where $\nu (y)$ is a positive probability density on a domain $\mathcal{Y}$. In the context of this paper’s application, we choose ${\nu }_1(y)$ as the probability density for uniform distribution on the time-price domain $\mathcal{Y}=[0,T]\times G$ and ${\nu }_2(y)$ as the probability density for uniform distribution on the price domain $\mathcal{Y}=G$.

Our loss function (19) has two differences compared with the loss function of the original DGM of Sirignano and Spiliopoulos (2018). First, we have a new hyperparameter ${\lambda }_L$ in the loss function, which controls the relative significance of the differential operator term in the NN training. We take ${\lambda }_L\ge 1$ because the term $L_1(f)$ is more related to the accuracy of the approximation (see further discussion in Section 3). Second, we do not include a term in the loss function corresponding to the boundary condition (16). The reason is that, for most of the underlying processes that we encounter in practice, we solve the Cauchy problem (11)–(12) which does not have a boundary condition. As is typical for all solvers of the Cauchy problem, we solve the PDE problem in a sufficiently large computational space domain G, so that the domain of interest (e.g. for pricing) is smaller than G to avoid imposing boundary conditions. However, having said that, in this paper we also consider a special case of the underlying process, a constant elasticity of variance (CEV)-like process, which has a boundary condition: we illustrate this using the example of the SABR model in Section 4.2.

2.3.2. Overview of Neural Network Training.

Next, we briefly recall the deep learning terminology. As discussed in Section 2.3.1, our goal is to train the NN $f(t,\mathit{x};\mathit{\theta })$ (with $\mathit{\theta }$ being the trainable NN parameters) to accurately approximate $u(t,\mathit{x})$ so that the loss function $L(f)$ is minimized. Hence, we have an optimization problem: we need to find the NN parameters $\mathit{\theta }$ that make the loss the smallest. To capture the complexity of the PDE solution, a complex deep NN has to be used to represent the solution, and finding optimal NN parameters could become a laborious task. This task is dealt with much more efficiently using a GPU rather than a CPU, predominately thanks to GPUs having a large number of so-called tensor cores working in parallel and aimed to accelerate matrix operations.

The NN training starts with initialization of the parameters $\mathit{\theta }$ and with generating randomized points $(t,\mathit{x})$ in $[0,T]\times G$. These points are fed into the input layer of the NN and propagated to evaluate $f(t,\mathit{x};\mathit{\theta })$ (the forward propagation step) and to calculate the loss $L(f)$. We then use an optimization algorithm (or optimizer) to update the trainable parameter $\mathit{\theta }$ to minimize the loss function. One of the simplest examples of optimization algorithms is gradient descent. The optimizer requires efficient calculation of gradients of the loss function with respect to the NN parameters that is a back-propagation algorithm (Goodfellow et al. 2016).

Let us now summarize. The idea of solving PDE problems via a NN is that, by applying the optimization algorithm of choice and back propagation from the deep learning (Goodfellow et al. 2016), one can efficiently find the set of NN parameters $\mathit{\theta }$ which minimizes the loss function $L(f)$. The NN solves PDE problems without building a mesh. We also remark that the deep learning approach allows us to break away from the “curse of dimensionality” and efficiently solve high-dimensional PDEs, including high-dimensional parametric PDEs as shown in this paper.

The algorithm adapted for our setup begins by initializing its parameters $\mathit{\theta }$ with some ${\mathit{\theta }}_0$, followed by the loss optimization loops. The nth loop is as follows:

1. For the differential operator term, generate random points $(t_n,{\mathit{x}}_n)$ from $[0,T]\times G$ according to the uniform distribution density ${\nu }_1$ on $[0,T]\times G$ ; for the terminal condition term, generate the random points ${\mathit{w}}_n$ from G according to the uniform distribution density ${\nu }_2$ on G.

2. Calculate the mean-square error $G({\mathit{\theta }}_n,{\mathit{s}}_n)$ at the M randomly sampled points ${\mathit{s}}_n=\{{((t_{ni},{\mathit{x}}_{ni}),{\mathit{w}}_{ni})}_{i=1}^M\}$:

   $$G({\mathit{\theta }}_n,{\mathit{s}}_n)=\frac{{\lambda }_L}{M}{\displaystyle \sum _{i=1}^M{\left(\frac{\partial }{\partial t}f(t_{ni},{\mathit{x}}_{ni};{\mathit{\theta }}_n)+\mathcal{L}f(t_{ni},{\mathit{x}}_{ni};{\mathit{\theta }}_n)\right)}^2}+\frac{1}{M}{\displaystyle \sum _{i=1}^M{(f(T,{\mathit{w}}_{ni};{\mathit{\theta }}_n)-u_T({\mathit{w}}_{ni}))}^2}.$$(20)

3. Perform the optimization algorithm to update the parameters $\mathit{\theta }$ at the random point ${\mathit{s}}_n$. In this paper, we use the adaptive movement estimation algorithm (ADAM) introduced by Kingma and Ba (2014), which is one of the most popular and efficient stochastic optimization algorithms, to perform optimization. In ADAM we choose the learning rate ${\alpha }_n$ ourselves (see further details in Section 2.3.3) while we use the default values for the rest of its parameters as in Kingma and Ba (2014).

4. Stop iterations if a convergence criterion is satisfied. The convergence criterion could be $G({\mathit{\theta }}_n,{\mathit{s}}_n)$ smaller than some predetermined threshold or ${\Vert {\mathit{\theta }}_{n+1}-{\mathit{\theta }}_n\Vert }_2$ smaller than some predetermined threshold. Our choice of convergence criteria in this paper is to train for a fixed number of epochs (as described in later sections) and choose the outcome (i.e., the parameters $\mathit{\theta }$) of the one with the least loss.

We recall that an epoch is a single pass of the entire data set through the NN during training whereby updating its parameters.

As it was pointed out in Sirignano and Spiliopoulos (2018), ${\mathit{\theta }}_n$ will converge to a critical point of the loss function $L(f(·;\mathit{\theta }))$ as $n\to \infty$; that is

However, it is likely that ${\mathit{\theta }}_n$ only converges to a local minimum rather than a global minimum since the NN $f(t,\mathit{x};\mathit{\theta })$ is nonconvex. It is for this reason that the hyperparameter ${\lambda }_L\ge 1$ in Equation (20) is needed to improve the accuracy of the approximation.

The algorithm can solve the PDE problem (11)–(12) with varying time, t, and spot price, x, with other input parameters being fixed. Its extension to solving parametric PDE problems is explained and illustrated in later sections.

We can use different NN architectures to represent PDE solutions, even the simplest feedforward neural network architecture known as multilayer perceptrons (MLP). In this paper, we mostly use the DGM NN architecture proposed by Sirignano and Spiliopoulos (2018) that we found effective for the PDE and PIDE problems under consideration, and we present its architecture details in the following section. We briefly discuss the comparison results of MLP versus DGM in Online Appendix B.

2.3.3. DGM Neural Network.

It is known (Sirignano and Spiliopoulos 2018) that the speed of NN training is problem-dependent and that it is important to choose an appropriate, problem-specific NN structure to achieve reasonable training speed for a problem under consideration. Thanks to the activities of the machine learning community in recent years, we can gain access to numerous NN structures for PDE applications. The DGM network, used in this paper, can be considered as a variant of the widely used long short-term memory networks (LSTMs) introduced by Hochreiter and Schmidhuber (1997). LSTMs are a special kind of recurrent NN that are useful for modeling sequential data such as time series prediction or speech recognition problems. It is shown by Sirignano and Spiliopoulos (2018) that an LSTM-type architecture works well for PDE problems with “sharp turns” in the initial or terminal conditions, which is exactly the situation we have here, because we have a step function in the terminal condition of the backward Kolmogorov equation that we approximate.

Next, we introduce the DGM NN architecture:

where $\mathbf{x}$ denotes the input layer, $\vartheta (·)$ is the activation function, $\odot$ denotes element-wise multiplication, the number of hidden layers is $L+1$ (there are 1 initial hidden layer and L LSTM like layers), and $\mathbf{U},\mathbf{W},\mathbf{b}$ are the parameters of the network which together form $\mathit{\theta }$. Particularly, $\mathbf{x}$ is a $D\times 1$ vector, where D is the number of trainable features; $\mathbf{U}$ and $\mathbf{W}$ are $M\times D$ matrices, where M is the number of nodes in the neural network; $\mathbf{b}$ is a $M\times 1$ vector. The input $\mathbf{x}$ can be just $(t,\mathit{x})$ from $[0,T]\times G$ for solving a PDE problem with fixed parameters as in the previous subsection or can also include parameters of the PDE when we solve parametric PDE problems—as in the following sections. We use the same letter f in slightly different contexts, but this should not lead to any confusion.

For completeness of the presentation, we visualize the NN used for this paper in Figure 1. The input layer $\mathbf{x}$ represents the input data and in the case of parametric PDEs that includes randomly sampled time points t and all the space and parameters points $\mathit{x}$. The input layer $\mathbf{x}$ only contains the simulated data, and we use different simulated data for every five epochs. In the NN considered, $\mathbf{x}$ is not just the input of the initial layer but also the input of all the hidden layers. If we zoom in on the LSTM-like hidden layer (Figure 1(b)), then we can see that apart from the initial data $\mathbf{x}$, the input of the current hidden layer also includes the output of the previous hidden layer ${\mathbf{S}}^{\mathcal{\ell }}$. A hidden layer with this setting could strengthen the long-term memory (${\mathbf{S}}^{\mathcal{\ell }}$, $\mathcal{\ell }=1,\dots ,L+1$) of the network (e.g., the architecture for LSTM networks (Hochreiter and Schmidhuber 1997) and highway networks (Srivastava et al. 2015)), which, as suggested in Sirignano and Spiliopoulos (2018), could then subsequently translate to the improved performance of capturing the nonlinearity and “sharp turn” in solving PDEs as presented in this paper. Once we have calculated the output f, we then use it to construct the loss function before applying an efficient optimization algorithm (like ADAM we are using here) to update the NN parameters with the aim of minimizing the loss function.

A NN with this LSTM-like architecture is quite complex, consisting of a chain of hidden layers with each hidden layer containing sublayers. Goodfellow et al. (2016) suggest that empirically a “deep” NN consisting of more layers and more nodes per layer performs better than a “shallow” NN. However, a deeper NN requires more training time with potentially only a marginal increase in accuracy. Therefore, one needs to find an optimal set of hyperparameters (the set of parameters that controls the NN architecture and the learning process) for training. Sirignano and Spiliopoulos (2018) suggest using the following hyperparameters: four hidden layers in the network ($L=3$), 50 nodes per layer ($M=50)$, using the hyperbolic tangent function as the activation function $\vartheta (·)$, which guarantees that $f(\mathbf{x};\mathit{\theta })$ is smooth and, hence, is suitable for approximating classical solutions of PDEs. Xavier initialization is used to initialize parameters and, as noted earlier, the ADAM optimization algorithm is used to update the parameters. We use the default settings for ADAM in TensorFlow but adjust the learning rate. We find that this set of hyperparameters also works well for NN training in this paper. Here, we use a fixed learning rate ${\alpha }_n=0.0001$, and we find that this ${\alpha }_n$ value gives an optimal convergence result for our problems. After training, we choose the NN with parameters such that it has the smallest total loss. We implement this algorithm in TensorFlow, which is an open-source software library developed by Google and designed specifically for deep learning. TensorFlow has a built-in function for parameter initialization and optimization algorithms. This allows great simplification of the NN development process.¹

For the numerical illustrations in this paper, we use 5,000 sample points per epoch. An epoch means passing the entire 5,000 randomly generated sample points backward and forward through the NN. We note that the numbers of sample points used here are greater than the 1,000 sampling points used in Sirignano and Spiliopoulos (2018) because we focus on solving parametric PDEs and greater numbers of sample points are needed to get to the required accuracy. To increase the training efficiency, each newly generated data set is trained on five epochs. Depending on their complexity, different models require different numbers of epochs to train to achieve the required accuracy. In theory, a more accurate NN approximation of the PDE/PIDE solution can be obtained with more training. However, at the same time, the cost of improvement becomes higher, and users face a situation of diminishing returns when more training is used. In Sirignano and Spiliopoulos (2018), six GPU nodes are used along with the asynchronous stochastic gradient descent method for parallel training. The training setup presented here used a single GPU. For illustration purposes, we performed deep learning training via the Google Colab Pro server, where we used a single NVIDIA Tesla P100 GPU or, alternatively, a single NVIDIA RTX 2080 GPU in a gaming laptop. In other words, even without access to high-powered computers, we are able to use the techniques developed in this paper, in a reasonable amount of time, to train the networks. We compare training time using various GPUs in Online Appendix A.1. Moreover, once trained, these networks are used in option pricing calculations that require only minor computational power (see Online Appendix A.2).

3. Tuning and Performance Analysis Under Geometric Brownian Motion

We first illustrate our implementation using geometric Brownian motion (GBM) as in the derivation of the Black-Scholes-Merton formula (Black and Scholes 1973, Merton 1973). Being one of the simplest models for the underlying asset, the GBM model is also one of the easiest to train. It is a lot less costly than others when it comes to tuning the hyperparameter ${\lambda }_L$ and hence useful in understanding the performance of the proposed approach.

Under the forward measure of the log asset price space, the backward Kolmogorov equation for the CDF $C(t,x)=C(t,x;T,y)$ in the GBM case is

with terminal condition

To find the TPDF from the cumulative distribution, we must first solve the backward Kolmogorov equation in a parametric manner for a range of $y\in [y_{\text{min}},y_{\text{max}}]$, given a minimum value of the range $y_{\text{min}}$ and a maximum value $y_{\text{max}}$, so that after taking the derivative of $C(t,x;T,y)$ with respect to y, we obtain the TPDF $p(t,x;T,y)$. Also, from the asset price process, we have that the CDF is a function of the volatility, $\sigma$. Thus, the problem needs to be solved in the $(t,x,y,\sigma )$-domain $Q\subset {\mathbb{R}}^4$, which contains not only the independent variables $(t,x)\in [0,T]\times G\subset {\mathbb{R}}^2$ with respect to which the PDE is formulated, but also the parameters $(y,\sigma )\in {\mathbb{R}}^2$. Using the NN $f(\mathbf{x};\mathit{\theta })=f(t,x,y,\sigma ;\mathit{\theta })$ described in the previous section, we can find a solution to this parametric problem for the backward Kolmogorov equation.

The optimization process starts with initializing the network parameters $\mathit{\theta }$. Then, for each epoch, we uniformly sample points from Q before feeding them into the DGM network and using an optimization algorithm to minimize the loss function. The process is repeated until the required accuracy is reached. The loss function in this example is defined as

The first and second partial derivatives terms $\partial f/\partial t$, $\partial f/\partial x$, and ${\partial }^{2}f/\partial x^2$ can either be computed directly using the gradient function in TensorFlow (or similar derivatives calculation functions in other deep learning libraries) or approximated via numerical differentiation with finite differences.

The domain Q in which we train the network is $x,y\in [-2.3,2.3]$, $\sigma \in [0,0.6]$, $t\in [0,1.2].$ We note that the choice of $\pm 2.3$ for the range of x and y corresponds approximately to $\pm 3.5{\sigma }_{\text{max}}{t}_{\text{max}}$, where ${\sigma }_{\text{max}}$ and $t_{\text{max}}$ are the upper bounds of $\sigma$ and t parameters, respectively. The time horizon considered for the network is $T=1.2$. We also note that $T-t$ is the time to maturity, which in future we will often denote as T again, but this should not cause any confusion. Although we fix the log-spot price $x_0=0$, we still need to train the network in x as the PDE is written in x. We train 2 million epochs. For every five epochs, we use 5,000 random points. Thus, in total, we use 2 billion random space time points to train this network. A single NVIDIA P100 GPU takes 0.022 seconds to complete 1 epoch, which means it takes about 12.5 hours to complete NN training using 2 million epochs. Once the offline training of the NN has been completed, it can be used online for the pricing of any option with this underlying model (within the range of volatility $\sigma$ and maturity T for which the network was trained) and for other computational finance tasks. This is unlike conventional numerical methods, which have to be run over again (i.e., “online”) for each set of the parameters’ values separately: considerably slower than taking the required values from a trained NN. In Online Appendices A.1 and A.2, we report the speed of both offline training (across different GPUs) and online access time.

Once the cumulative distribution function $C(t,x;T,y)=C(t,x;T,y,\sigma )$ has been approximated (via the NN training) by $f(\mathbf{x};\mathit{\theta })=f(t,x,y,\sigma ;\mathit{\theta })$, we use numerical differentiation to approximate the TPDF $p(t,x;T,y)=p(0,x;T-t,y,\sigma )$ by $\tilde{p}(0,x;T-t,y)=\tilde{p}(0,x;T-t,y;\mathit{\theta })$:

for some $\mathrm{\Delta }>0$. In our experiments we chose $\mathrm{\Delta }=0.005$. The choice of $\mathrm{\Delta }$ should be neither too small nor too large: if it is too small, that could introduce oscillations; if it is too large, the approximation accuracy around the peak area of the density is insufficient. We find that a good choice for $\mathrm{\Delta }$ is in the range of [0.001, 0.01].

3.1. ${\lambda }_L$ Tuning

We use most of the NN hyperparameters as set in Sirignano and Spiliopoulos (2018), and we find their modeling setup effective (apart from a few alterations we have adapted for this application; see Section 2.3.3). We introduced a new hyperparameter ${\lambda }_L$ in the loss function, and in this section, we use numerical experiments to understand the impact of this parameter. This hyperparameter gives relative weighting to the two loss terms. Without the hyperparameter ${\lambda }_L$, the loss terms $L_1$ and $L_2$ bear the same weight. When ${\lambda }_L>1$, this means that the loss term with the differential operator weights more than the terminal condition term; likewise, when ${\lambda }_L<1$, it indicates that the terminal condition needs to be satisfied more compared with the differential operator. Intuitively, the loss term $L_1$ needs to have a higher weight than $L_2$ because the error $L_1$ is expected to have a higher influence on the accuracy of the NN. We tune the hyperparameter ${\lambda }_L$ and validate our assertion. To this end, we train the NN for different values of ${\lambda }_L=0.1,1,10,100,1,000$. We compare the obtained approximation $\tilde{p}$ with the exact density p, which is Gaussian in the example considered:

For each ${\lambda }_L$, we train up to 2,000,000 epochs and record the best training results (i.e., with the smallest loss) within checkpoints set at 5,000, 25,000, 50,000, 100,000, 250,000, 500,000, 1,000,000, 1,500,000, and 2,000,000. Figure 2 demonstrates the tuning. Across all the considered ${\lambda }_L$, the training accuracy generally increases as more epochs (or sampling points) are used. However, after a certain amount of training, for example, after 500,000 epochs, the root mean square error (RMSE) does not go down as fast as at the start of the training: This is especially apparent when ${\lambda }_L=1$ and 10, for which the algorithm cannot find better NN parameters beyond 1 million epochs (as shown by the dashed lines). This is because of the vanishing gradient problem in deep learning, meaning that it is becoming increasingly difficult to find optimal NN parameters to reduce the loss. When ${\lambda }_L=10$, the vanishing gradient problem is alleviated; ${\lambda }_L=100$ helps to improve the performance further; but the performance gets worse when ${\lambda }_L=1,000$. When the number of epochs used is greater than 500,000, we see that ${\lambda }_L=100$ produces the best training results for the given time domain ($t=0.25,0.5,0.75,1.0)$. Overall, we find that ${\lambda }_L=100$ produces accurate results for our NN approximation framework.

3.2. Performance Across Maturities

Figures 2 and 3 demonstrate that the NN approximation can closely approximate the exact TPDF in terms of RMSE. However, when considering the accuracy of the relative approximation, the log scale plots in Figure 3 reveal that the approximation is accurate around the center but not for the tails. Specifically, the approximation cannot yield highly accurate results for density values smaller than ${10}^{-3}$. This limitation may be attributed to the error measure used during training. As has been discussed in Section 2.3.1, the training is designed to minimize the $L^2$ loss, that is, the mean square error. Hence, it is unsurprising that the results are more accurate in terms of RMSE than relative error and that the approximations of the tails yield large percentage errors, which is a limitation of our approach.

We also see from Figure 2 that when maturity t is decreased, the approximation accuracy decreases. This observation is not surprising because when t is made smaller, the TPDF concentrates around the initial state (i.e., $x=0$ in this example) and the derivatives of the CDF C with respect to y tend to infinity when t goes to zero (see Equations (13) and (23)), affecting the accuracy of the numerical differentiation. If the training range in x is too wide ($[-2.3,2.3]$ in this example), then the narrow range around the zero initial state may not acquire sufficient data points, resulting in less accurate approximations compared with those obtained for larger t. If we did not train specifically for t being close to zero, the approximation accuracy for this region could be low. Moreover, because we are dealing with the backward Kolmogorov equation with a terminal condition, when translating the time to maturity, t, to the time variable used in the loss function, it becomes $T-t$, where T is the terminal time. That means, as $t\to 0$, the time variable in the training $T-t\to T$.

As discussed earlier, when ${\lambda }_L=100$, the emphasis is on the differential operator $L_1$ term in the loss function L rather than on the terminal condition $L_2$ term, which increases the overall accuracy of the approximated results while simultaneously it exacerbates the inaccuracy for $t\to 0$. In Figure 4, we repeat the same tuning as in Figure 2 but this time we show the plot for various small times to maturity $t=0.025$, 0.05, 0.075, and 0.1. We see that for these small times to maturity, ${\lambda }_L=1$ is the clear winner. As t increases, the accuracy gap narrows, and at $t=0.1$ the differences between the results obtained with ${\lambda }_L=1$, 10, and 100 are very small. The reason why ${\lambda }_L=1$ performs best for small times to maturity is that equal emphasis is put on the terminal condition term and the differential operator term, resulting in higher accuracy for the region close to T (or time to maturity close to zero). It should be noted that when ${\lambda }_L=0.1$, although more training emphasis is placed on the terminal condition term, performance is much worse due to the lack of emphasis on the main body of the approximated solution (i.e., on the differential operator term).

Transfer learning allows us to address the additional complexity arising from NN training across different maturities as discussed above. The use of transfer learning is a common strategy to train a new NN for a problem using the knowledge from an already trained NN for a related problem (Goodfellow et al. 2016). Here, we transfer the NN we already have from one time range to another. Specifically, if the TPDF with smaller T is what we require, we can use an already trained NN on a larger time range to train for smaller T in order to fine tune the NN for a small time-to-maturity domain. Transfer learning is much more efficient than training a NN from scratch. For instance, here we use the NN trained on the domain as defined above (training done with ${\lambda }_L=100)$ and then change the domain Q to the new domain $x,y\in [-0.75,0.75]$, $\sigma \in [0,0.6]$, $t\in [1.08,1.2].$ We only train a NN using another 2,500 epochs and we select the NN with the smallest total loss. Consequently, we achieve an accurate NN approximation that is dedicated to smaller maturities T. The performance is shown in Figure 4. We see that the performance of the NN obtained via transfer learning is better than that of the NN trained with ${\lambda }_L=1$. We fine tune the NN for the new domain using only a small number of epochs because we have altered the weight in the loss function from ${\lambda }_L=100$ to 1; training with more epochs would make the NN lose the traits of overall accuracy brought by the NN trained with ${\lambda }_L=100$ and, instead, converge to the NN trained with ${\lambda }_L=1$ as more epochs are used. Figure 5 illustrates the TPDF before and after tuning and the improvement is visible. Similarly, when we need a result that is out of the training domain range, we can use transfer learning to obtain the required NN efficiently. This transfer learning idea is applied to all the stochastic models considered in this paper.

3.3. Pricing Performance

In this section, we test the accuracy of the neural TPDF by using it in option pricing. For illustrative purposes, we begin with a simple example. Table 1 shows the pricing performance for put options for a given set of parameters under the GBM process. Overall, the pricing is quite accurate in terms of absolute error, with all absolute errors less than $5\times {10}^{-4}.$ However, the relative errors for deep out-of-the-money options are large. This is because the prices of deep out-of-the-money closed form put options are very close to zero. If the approximate TPDF does not reach the required level of accuracy, it will inevitably fail to accurately price the out-of-the-money put options. In addition, for deep out-of-the-money options, when the prices are very small, the approximate TPDF could give negative prices, although the negative prices are very close to zero. This can be seen as a limitation of the NN approximated density.

Table 1. Put Option Prices Calculated from GBM’s Closed Form TPDF (CF) and Neural Network Approximate TPDF (NN)

An important feature of the NN approximation approach is that, once the model is trained in the parametric space, it can quickly produce the TPDF given a set of parameters. To assess its performance, we compare the pricing errors of the neural TPDF against the closed form solutions over the randomly generated parameter sets across different moneyness positions and maturities. Zhou (2022) defines moneyness ranges from a comprehensive data set of all US-listed equity call and put options from 1996 to 2019, and we base our definition of moneyness ranges for put options on this: Deep in-the-money (DITM) (1.20–1.40); In-the-money (ITM) (1.05–1.20); At-the-money (ATM) (0.95–1.05); Out-of-the-money (OTM) (0.80–0.95); Deep out-of-the-money (DOTM) (0.60–0.80). We acknowledge that this classification can only be applied to normal market conditions and not to extreme market conditions. We do not take into account the extreme strikes at either end, where pricing accuracy is lower, typically for any generic approximation method for a density. Here we use the price percentage error (PrPCTE) and the price root mean square error (PrRMSE) to evaluate the pricing performance with the NN approximated TPDF. The results shown in Table 2 confirm our observation from the example in Table 1. Pricing is accurate in terms of RMSE but less accurate in terms of absolute error. We recall that, in the training of the Black-Scholes-Merton CDF, the training loss reaches ${10}^{-3}$ and the RMSE of the density reaches $2\times {10}^{-3}$. It is not surprising that the RMSE of the prices also reaches a similar level of accuracy. However, this is more of a problem with using an approximated density to price options rather than a problem specific to the NN approach.

Table 2. Pricing Errors of the NN Estimated Density of the Black-Scholes-Merton Model, Benchmarked Against BSM’s Closed Form Put Pricing Solutions

We have shown in Section 3.2 that the neural TPDF accuracy decreases as the time to maturity decreases, and the option price evaluations in Table 2 confirm this observation in general but the impact in terms of RMSE is still small. We emphasize that, although the pricing error is of the order ${10}^{-4}$, the relative error can still be large because put prices are close to zero when strikes and volatility are small. This can be seen from the extreme PrPCTE values in Table 2. The extreme PrPCTE is produced because the closed form solutions can give results with high accuracy, whereas the results using the NN approximated TPDF can only give a certain level of accuracy, depending on how small the loss achieved in training is. For example, if the closed form price is ${10}^{-16}$ but its approximation is ${10}^{-4}$, and the percentage error is about ${10}^{12}$, whereas the absolute error is still ${10}^{-4}$. This is a general limitation of using an approximate TPDF for pricing options regardless of the method used to approximate the TPDF. In terms of pricing performance across different volatilities, there is no significant difference in accuracy even when the volatility, $\sigma$, is as small as 0.05, and thus, we do not present the results here.

Once the NN is trained, retrieving the neural TPDF is very quick, and this brings us to the important feature of replacing previous calculation engines for TPDF with one developed through deep learning: accessing a trained NN to retrieve the density takes very little time, meaning that calculating a single option price using neural TPDF with QUAD is as fast as using the closed form density with QUAD. This preserves the immense speed advantage of not just the original QUAD shown in Andricopoulos et al. (2003) but for more complex models of the underlying which had previously slowed the QUAD calculations (Chen et al. 2014).

In the next section, we provide further illustrations of the proposed approach to approximating TPDFs using deep learning by applying it to several stochastic volatility and jump diffusion models.

4. Neural TPDFs for Stochastic Volatility and Jump Diffusion Models

In the previous section, we introduced our approach using the simplest model: the geometric Brownian motion, a.k.a., the original Black-Scholes-Merton setup. In this section, we apply the proposed deep learning approach to a number of more complex, representative models of underliers. We first consider the standard Heston model (we also successfully applied our methodology to a Heston model with time-dependent parameters, but the corresponding results are not included in the paper). Next, considering the SABR model, we show that nonaffine models can be tackled by our approach. To further confirm the universality of our approach, we add jumps to models for underliers and demonstrate how partial integral differential equations (PIDEs) can be approximated efficiently by NNs.

4.1. Heston Model

Under the forward measure, the Heston model (Heston 1993) written for the log stock price, X, and variance, V, of the stock price has the form

where $W_1$ and $W_2$ are Wiener processes with correlation $\rho$, $\omega$ represents the long term variance, $\kappa$ is the rate at which V reverts to $\omega$, and $\xi$ is the volatility of the volatility (vol of vol). We require that $2\kappa \omega \ge {\xi }^2$ to ensure that the variance process is strictly positive at any finite time.

The backward Kolmogorov equation for the joint CDF $C(t,x,v)=C(t,x,v;T,y,z)$ is

with terminal condition

This PDE problem is of the form (14)–(15).

We build a NN, f, to approximate the CDF, C. For the Heston model, the neural network f is trained in the $(t,x,y,v,z,\kappa ,\omega ,\xi ,\rho )$-domain $Q\subset {\mathbb{R}}^9$, which is high-dimensional. Traditional methods aimed at solving PDEs under fixed parameters are uncompetitive in this challenging parametric PDE setting in comparison with the deep learning algorithm. The high dimensionality here is of the parametric space. Recall that we need to solve high-dimensional parametric PDEs to find TPDFs for a range of values in the time-state and parametric spaces to effectively use the TPDFs for financial engineering tasks such as pricing options. The problem (14)–(15) is higher dimensional than the geometric Brownian motion case in Section 3. Therefore, it can be expected that the number of epochs and sample points for NN training are larger. We follow the same training process as in the previous section with natural adjustments.

Similar to the geometric Brownian motion case, we take the spot log-price $x_0=0$. Furthermore, we choose Q so that $x,y\in [-3.5,3.5]$, $v,z\in [0,1.0]$, $\kappa \in [0.8,1.2]$, $\omega \in [0.1,0.5]$, $\xi \in [0,0.5]$, $\rho \in [-0.5,0.5]$, and $t\in [0,1.2]$. We compare the approximate joint TPDF with the semiclosed form of the TPDF from Lewis (2016).

We first tune the hyperparameter ${\lambda }_L$. Although we ultimately calculate the joint TPDF $p(t,x,v;T,y,z)$ for the Heston model, it is convenient to perform tuning based on the accuracy of the NN approximation of the marginal density $p_X(t,x;T,y)$. The marginal density is calculated by integrating the joint density $p(t,x,v;T,y,z)$ from 0 to $\infty$ with respect to z, that is,

The integration is done using the trapezium rule. The NN training is performed by running 500,000 epochs. Figure 6 displays the joint TPDF and the corresponding marginal density. We find that the approximate TPDF is close to the true TPDF. It is difficult to make a direct visual comparison using the two-dimensional joint TPDF plots; the marginal density, on the other hand, gives a better visualization and shows that the neural TPDF provides a very good approximation of the true TPDF obtained from the semiclosed formula as shown in Figure 6: The two marginal density plots are indistinguishable. In addition, the marginal density itself is more closely related to European option pricing and extracting the volatility smile. Thus, here we tune ${\lambda }_L$ according to the marginal density. The tuning results are summarized in Table 3. We find that ${\lambda }_L=1$ to 100 produce similar results, although ${\lambda }_L=100$ performs slightly better, whereas ${\lambda }_L=1,000$ performs worse than the rest. Moreover, ${\lambda }_L=100$ has an advantage compared with ${\lambda }_L=1$ or 10. Indeed, from Figure 7, one can see that when ${\lambda }_L=100$ the $L_2$ loss still shows a stronger downward sloping trend and is less volatile whereas for ${\lambda }_L=1$ and 10 the $L_2$ losses are already flattening and more volatile. This means that further training for ${\lambda }_L=100$ can improve the $L_2$ loss, whereas its $L_1$ loss is already smaller than for the other values of $\lambda$. Thus, we again observe that the introduction of ${\lambda }_L>1$ helps to mitigate the vanishing gradient problem commonly encountered in machine learning training. This tuning result agrees with the one obtained for the GBM model in Section 3.1. We tested 500,000 epochs here, and the performance of ${\lambda }_L=100$ is already better than the others. It can be assumed that with more training, the difference in accuracy will become larger, as has been shown in the GBM model case. Thus, we recommend the setting ${\lambda }_L=100$ for training all NNs.

Table 3. Tuning ${\lambda }_L$

Following the previous discussion, we choose ${\lambda }_L=100$ and, to gain higher accuracy, we train 2 million epochs and 5,000 points every 5 epochs, which is equivalent to using 2 billion sample points to train this NN. Using a single NVIDIA Tesla P100 GPU, each epoch takes about 0.038 seconds, which is slightly more than for the one-dimensional Black-Scholes-Merton example in the previous section. Completing the whole 2 million epochs of training takes about 21 hours. This time can be greatly reduced by using more powerful GPUs. Far more important, we emphasize that this is a training time that is not part of the computation times for later use of the neural TPDF, for example, for option pricing calculations, just as the time to write computer code for options calculations is not part of the later option pricing calculation times. NN training is done separately (offline). Once the NN is trained, it is used as a highly convenient engine (online) to obtain the approximate TPDF quickly. Indeed, next, we illustrate by showing calculations of option prices using the neural TPDF that are faster than when the semiclosed form of the TPDF is used.

We compare the pricing performance incorporated into QUAD. As can be seen in Table 4, option pricing using the neural joint TPDF produces sufficiently accurate results when compared with the semiclosed formulas. At this point, we highlight a considerable advantage of the NN over the semiclosed form of TPDF. Based on QUAD calculation with 51 mesh points per direction, QUAD with the neural TPDF is 7 times faster than with the semiclosed form TPDF. If we increase the number of mesh points per direction to 101, that is 10,201 points in the QUAD mesh, the speed advantage is about 11 times. The speed difference is larger with more points because the computing time in accessing the NN is independent of the number of points used (see Online Appendix A.2). This is a good property brought by the NN method. In practice, QUAD with 51 mesh points per direction can produce very accurate results. The neural TPDF is very fast to evaluate and it is a perfect match with QUAD.

Table 4. Put Option Prices and Option Implied Volatility (IV) Calculated from Semiclosed Form TPDF (CF) and Neural TPDF (NN)

In this illustration, we calculate put prices and then obtain the volatility smile by inverting the Black-Scholes-Merton formula. We use an 11-point calculation and interpolate the smile. Figure 8 shows the implied volatility smiles using the semiclosed form TPDF and the neural TPDF. The results are also presented in Table 4. The neural TPDF can reproduce a smile close to the semiclosed form pricing formula for the Heston model. The smile is accurate for both deep in-the-money options and out-of-the-money options.

Similar to Section 3.3, to further investigate the accuracy in the parametric space, we compare the pricing errors of the NN approximated density against the semiclosed form solutions over the random parameter sets across different moneyness positions and maturities.

Table 5 shows the results of pricing errors for the Heston model. The pricing performance with the NN approximated density is good overall, showing that the NN approximated density is sufficiently accurate to produce reasonable prices and implied volatilities. Similar to the results shown for the Black-Scholes-Merton model, a number of features can be observed from these results. First, the pricing performance is better at larger t, because the CDF has less rapid changes, making it easier for the NN to approximate. Second, the approximation is only accurate to a certain extent. For example, if the training is accurate up to ${10}^{-3}$, then this will translate to a similar level of accuracy in the density approximation. Therefore, when the approximated density is subsequently used to price options, the accuracy may be lower due to the loss of accuracy incurred in either the density calculation step of the CDF or the QUAD method. Thus, it is seen that when pricing a deep out-of-the-money option with a small maturity, although the RMSE of the price is small, the percentage errors of the prices can still be quite large because the prices of deep out-of-the-money put options are usually very small, making it difficult to obtain prices with low relative (percentage) errors. This is a problem with using approximate density to price options rather than a problem specific to the NN approach.

Table 5. Pricing Errors of the NN Estimated Density of the Heston Model, Benchmarked Against Heston’s Semiclosed Form Put Pricing Solutions

4.2. SABR Model

The Black-Scholes-Merton model and the standard Heston model are affine. For these, closed forms (or semiclosed forms) of European option prices and the corresponding TPDFs are known. Nonaffine models, on the other hand, have neither closed form European prices nor TPDFs. The popular SABR model (Hagan et al. 2002) is nonaffine when the parameter $\beta <1$, and it does not have a closed form or semiclosed form TPDF. Hagan et al. (2014) subsequently derived an approximate effective forward equation for SABR’s marginal density over price using an asymptotic expansion and the marginal density can then be calculated numerically using finite difference techniques. The joint TPDF for the SABR model, however, is still difficult to evaluate, even with numerical techniques. In this section, we show that our method can tackle the SABR model and can approximate the joint TPDF successfully.

The SABR model for forward price F and volatility $\sigma$ is written as (Hagan et al. 2002)

Here we consider the more challenging case of this model when $0<\beta <1$. It is known (Rebonato et al. 2009, Lewis 2016) that, for $0<\beta <1$, $F=0$ is an attainable barrier for the process $F(t)$. Because for $1/2\le \beta <1$ the solution $(F(s)$, $\sigma (s))$ of (33)–(34) is unique, $F=0$ is an absorbing boundary. For $0<\beta <1/2$, the absorbing boundary condition is prescribed at $F=0$ to satisfy the no arbitrage condition.

For training purposes, it is easier to train if we transform the volatility from $\sigma$. Let $\omega =\ln \sigma$, then Equation (34) becomes

Consequently, Equation (33) takes the form

To be specific, consider the pricing of a European option with payoff $f(x)$ and maturity T. Its price at time t, with spot price x and spot log-volatility v can be expressed as

where ${\tau }_{t,x,v}$ is the Markov time when $F_{t,x,v}({\tau }_{t,x,v})$ from (35)–(36) hits the value zero.

To compute the first expectation in (37), we need to find the transition probability function $C(t,x,v)=C(t,x,v;T,y,z)$:

where $p(t,x,v;T,y^{\prime },z^{\prime })$ is the corresponding density. The function $C(t,x,v)$ satisfies the backward Kolmogorov equationwith the terminal conditionand the Dirichlet boundary condition

To compute the second expectation in (37), we need to find the probability

which satisfies the same backward Kolmogorov equationwith the terminal conditionand the Dirichlet boundary condition

The solutions to each of the Dirichlet problems need to be approximated by their own appropriately trained NN. The problem (41)–(43) is easier to solve as it has fewer parameters. In order not to make the presentation overly intricate, we will price call options for which $f(0)=0$, so that only the problem (38)–(40) needs to be solved.

In contrast with the models considered earlier, here we have to prescribe the boundary condition for the function $C(t,x,v)$. Consequently, we need to add an extra term corresponding to the boundary condition in the loss function $L(f)$ (here f is a NN as usual):

where $L_1$ loss is for the differential operator term Equation (38), and $L_2$ loss is for the terminal condition term Equation (39), and $L_3$ loss is for the Dirichlet boundary condition Equation (40) with ${\nu }_3(y)$ being the probability density for the sampling distribution on the time-price domain $\mathcal{Y}=[0,T]\times \partial G.$ As before, ${\lambda }_L$ is the hyperparameter for the loss function. We assume that both the terminal condition and the Dirichlet boundary condition are of equal importance and thus have the same weight in the loss function.

For the SABR model, the NN f is trained in the $(t,x,v,y,z,\beta ,\rho )$-domain $Q\subset {\mathbb{R}}^7$. We simplify the training by removing the parameter $\nu$ since we can choose $\nu =1$ in (38) and recover the solution for arbitrary $\nu >0$ through the change of variables:

We choose the domain Q so that $t\in [0,1.2]$, $x,y\in$ [0.0, 5.0], $v\in [-5.0,0.0]$, $z\in [-7.0,3.0]$, $\beta \in [0.5,1.0]$, $\rho =[-1.0,0.0]$. The range for the terminal volatility value, z, is higher in the SABR case because the marginal density of the volatility process has a fat tail, observed in the lognormal distribution. As with the Heston model, to attain the required accuracy, we choose ${\lambda }_1=100$ and ${\lambda }_2={\lambda }_3=1$, and we train 2.5 million epochs and generate 5,000 independent points every five epochs, which is equivalent to using 2.5 billion sample points to train this NN. We select the NN f according to the epoch with the lowest loss. Using a single NVIDIA Tesla P100 GPU, each epoch takes about 0.043 seconds, slightly longer than for the Heston model due to the increased complexity. It takes about 30 hours to complete 2.5 million epochs of training.

Because there are no closed-form or semiclosed form solutions for the SABR model, as the benchmark, we use option prices calculated numerically by the simplest random walk algorithm for solving Dirichlet problems from (Milstein and Tretyakov 2003) (see also Milstein and Tretyakov (2021)) combined with the Monte Carlo technique; we take the very small time step 0.0005 and ${10}^8$ Monte Carlo runs. We also compare NN results for the TPDF and option prices with the ones obtained by the finite difference method (implemented in QuantLib² with 200 points across the time dimension, 800 points across the price dimension, and 200 points across the volatility dimension), by the approximate marginal density for price calculated using the finite difference method (Hagan et al. 2014), and by the approximate implied volatility formula from Hagan et al. (2002, 2014).

We see in Figure 9 that the marginal density for the stock price produced by the NN is very close to the results obtained using the finite difference methods. Hagan et al. (2014) report that the implied volatility formula can produce negative densities as price $F\to 0$, which violates the no arbitrage condition, and the same result is also observed here. The NN approximation has no such deficiency.

Figure 10 shows that the approximate SABR implied volatility formula, although very convenient to use, does not produce accurate results, compared with results using the other methods. The implied volatility smile obtained through the neural TPDF is accurate, and it is very close to the ones generated by the Monte Carlo method and the finite difference method (Table 6; Figure 10).

Table 6. Call Option Prices and Option Implied Volatility Calculated from the Three Numerical Approaches: Monte Carlo Method (MC), Finite Difference Method (FD), and QUAD with the Neural TPDF (NN)

Table 7 shows the pricing errors for the SABR model. Similar to the Heston model, the SABR performance shows that the NN approximated density can be used to obtain option prices up to a certain level of accuracy. Although the performance in terms of absolute pricing accuracy (PrRMSEs) is consistent across different option moneyness positions and maturities, the relative pricing accuracy (PrPCTE) is lower for OTM options, especially at shorter maturities, because the values of OTM options at shorter maturities are minuscule and prices obtained using the NN approximated density have difficulty in accurately capturing those minuscule prices.

Table 7. Pricing Errors of the NN Estimated Density of the SABR Model, Benchmarked Against Monte Carlo’s Call Pricing Solutions (with an Average Approximation Error $2\times {10}^{-4}$)

4.3. Adding Jumps

When it comes to computing TPDFs, jump-diffusion processes pose a unique challenge: The Kolmogorov backward equations for jump-diffusion processes involve the integral term due to the presence of jumps in the models and so the PDEs occurring in the case of the diffusion processes considered in the previous sections become PIDEs. Although solving PIDEs may be awkward by other numerical techniques such as finite difference methods, the deep learning approach presented in this paper can solve them efficiently. The use of deep learning to solve (especially parametric) PIDE problems has rarely been explored in the literature. Here we illustrate the deep learning approach to finding approximate TPDFs arising from jump-diffusion models by considering the two well-known jump-diffusion models: Kou’s double exponential jump-diffusion model (Kou 2002) and the stochastic volatility jump-diffusion model (Bates 1996).

4.3.1. Kou’s Double Exponential Jump Diffusion.

Kou’s model for the asset price S under a forward measure can be written in the following form (Kou 2002):

where $W(t)$ is a standard Wiener process, $N(t)$ is a Poisson process with intensity $\lambda >0$, and the jumps $J_1,J_2,\dots$ are such thatwith ${\xi }_i$ being independent and identically distributed double-exponential random variables with density

Here the parameters $p,q\ge 0,$ $p+q=1$, ${\eta }_1>1$, ${\eta }_2>0$, and ${\mathrm{𝟙}}_{\{A\}}$ is the indicator function of a set A. The jumps with random sizes $J_1,J_2,\dots$ occur at jump times ${\tau }_1,{\tau }_2,\dots$ at which the Poisson process $N(t)$ increases by one. The constant

It is assumed that the Wiener process $W(t),$ the Poisson process $N(t),$ and the jumps $J_i$ are independent.

Let $X(t)=\ln S(t)$ and rewrite the SDE (48) accordingly (we do not present this SDE here). The CDF $C(t,x)=C(t,x;T,y)$ for the $X(T)$ satisfies the backward Kolmogorov equation (Applebaum 2009):

with terminal condition

Note that (50) is a PIDE.

We rewrite the integral term of (50) as

where the integrals on the right-hand side can be approximated using a Gauss-Laguerre quadrature in the implementation of the NN training. However, the range of abscissas in a Gauss-Laguerre quadrature rule is too wide for NN training purposes. For instance, the seven-point Gauss-Laguerre rule needs to have abscissa up to 19.40 and, as a result, $z/{\eta }_1$ should have a much wider range than is required for the price range of interest if e.g. ${\eta }_1=2$. In comparison, the seven-point Gauss-Hermite rule has an abscissa range of $\pm 2.65$, which is well within the price range in the NN training. Therefore, we further transform the integral term (52):

The previous integrals cannot be approximated by standard Gaussian-Hermite quadratures because the lower bound is zero rather than infinity. Instead, we use Gaussian quadrature weights and abscissas suggested in Steen et al. (1969) to approximate ${\int }_0^{\infty }f(z)e^{-z^2}dz$. We find that using seven quadrature points suffices in terms of accuracy and efficiency.³

The domain Q used for training this model is $x,y\in [-5.0,5.0]$, $\sigma \in [0,0.5]$, $\lambda \in [0.0,2.0]$, $p\in [0.0,1.0]$, ${\eta }_1\in [1.1,20.0]$, ${\eta }_2\in [0.1,20.0]$, and $t\in [0,1.2]$. We choose ${\lambda }_L=100$ and, to get higher accuracy, we train the NN using 2 million epochs with 5,000 points every 5 epochs. We choose the NN with the smallest total loss. The training time per epoch using an NVIDIA P100 GPU is about 0.047 seconds, and it takes about 26 hours to complete the training of 2 million epochs.

Because Kou’s model does not have a closed form TPDF, we use the inverse Fourier transform of its characteristic function to derive the true TPDF in the log price space, before using it to compare with the neural TPDF. Figure 11 shows the comparisons. As can be seen from this figure, the neural TPDF provides a close approximation to the true density. Similar to the discussion in Section 3.3, the neural TPDF for Kou’s double exponential jump diffusion is accurate in terms of RMSE. However, when it comes to relative error, it performs well around the center, but its accuracy decreases toward the tails where the density values are low.

We also measure the performance of the trained NN using option prices and the implied volatility smile against the semiclosed form option prices. Table 8 gives the results. We observe that option prices obtained through the neural TPDF are accurate according to the benchmarking against the semiclosed form solutions; the implied volatility error is not larger than 1.7% (see also the implied volatility smile in Figure 12). The highlight here is the immense speed advantage of the combined QUAD plus the neural TPDF over the semiclosed form solution: the semiclosed form solutions require a number of expensive numerical procedures and take over 100 seconds to evaluate 11 options, whereas our deep learning parametric approach delivers very high-speed option valuations, costing only about 0.15 second to calculate 11 options.

Table 8. Put Option Prices and Option Implied Volatility Calculated from Semiclosed Form Option Prices (Kou 2002) and the Neural TPDF with QUAD in the Case of Kou’s Double Exponential Jump Diffusion Model

We further investigate pricing accuracy in the parametric space. To this end, we train Kou’s model in the parametric space and compare the pricing errors of the NN approximated density against the semiclosed form solutions over the random parameter sets across different moneyness positions and maturities. We use the same moneyness definition as in Section 4.1.

From Table 9, we see that the pricing performance is consistent in terms of RMSE (reaching ${10}^{-3})$. However, if we switch to the relative error, we observe that the performance drops notably for deep out-of-the-money options and small maturities. This is caused by the prices for deep out-of-the-money options under Kou’s model being very small (about ${10}^{-3}$ or less), and hence option prices calculated using QUAD with approximated density will inevitably produce high percentage errors, depending on the approximation accuracy of the density.

Table 9. Pricing Errors of the NN Approximated Density of Kou’s Double Exponential Jump Diffusion Model Compared with Kou’s Semiclosed Form Put Pricing Solutions

Table 10. Put Option Prices and Option Implied Volatility Calculated Using the One-Dimensional QUAD Method with Semiclosed Form Marginal Density of the SVJ Model and the Two-Dimensional QUAD with Neural Joint TPDF

4.3.2. Stochastic Volatility Jump Diffusion.

In this section, we consider a stochastic volatility model with jumps. Jumps can be added to the underlying alone or to both the underlying and the volatility. Gatheral (2011) points out that stochastic volatility with jumps or SVJ (i.e., jumps in the underlying only), suggested by Bates (1996), performs much better empirically than stochastic volatility with simultaneous jumps in both the price and volatility (SVJJ). The main reason, as the author notes, is that the SVJJ model has more parameters, making it harder to fit to observed option prices. Here we focus on Bates’ SVJ model, for which the TPDF is available in a semiclosed form Bates 1996 and hence allows us to perform insightful benchmarking of our approach to illustrate its accuracy and speed.

Bates’ SVJ model written for the log stock price, X, and the variance, V, of the stock price under a forward measure takes the form

where $W_1(t)$ and $W_2(t)$ are correlated Wiener processes, $N(t)$ is a Poisson process with intensity $\lambda >0$, and the independent jumps $J_1,J_2,\dots$ are so that

Here, ${\mu }_J$ is the mean logarithmic jump amplitude and ${\sigma }_J$ is the jump volatility, and the density of $\eta$ is

The constant $k=E(J)=\exp ({\mu }_J+{\sigma }^2/2)-1.$ The rest of the notation in (54) is the same as in the Heston model (29).

The backward Kolmogorov equation for the joint CDF $C(t,x,v)=C(t,x,v;T,y,z)$ is

with terminal condition

We approximate the integral term in Equation (57) by Gauss-Hermite quadrature with seven points. The domain Q used here for NN training is $x,y\in [-3.5,3.5]$, $v,z\in [0,0.6]$, $\kappa \in [0.8,1.2]$, $\omega \in [0.1,0.3]$, $\xi \in [0,0.4]$, $\rho \in [-0.5,0.5]$, $\lambda \in [0.0,2.0]$, ${\mu }_J\in [-0.5,0.5],$ ${\sigma }_J\in [0.0,0.5]$, and $t\in [0,1.2]$. We choose ${\lambda }_L=100$, and to get higher accuracy, we train 2 million epochs and sample 5,000 points every 5 epochs. As usual, we choose the NN with the smallest total loss. The training time per epoch using an NVIDIA P100 GPU is about 0.051 seconds, and it takes about 28.5 hours to complete the NN training of 2 million epochs.

For Bates’ SVJ model, the marginal density for X is known in semiclosed form (Bates 1996). In Figure 13, we plot the marginal density obtained through the neural joint TPDF against the semiclosed form marginal density. We see that the marginal density calculated from neural joint TPDF gives a very accurate result benchmarked against the semiclosed form marginal density (the RMSE is $1.72\times {10}^{-3}$).

We further analyze the NN performance using option pricing and implied volatility smile (see Table 10 and Figure 14). Once again, in terms of speed, we see that the two-dimensional QUAD method with neural TPDF is remarkably faster than the one-dimensional QUAD method with semiclosed form marginal density even when using fewer quadrature points. In terms of accuracy, the neural TPDF produces accurate option prices and implied volatility, with percentage implied volatility error smaller than 0.5%.

We train the stochastic volatility with jump model in the parametric space and compare the pricing errors of the NN approximated density against the semiclosed form solutions over the random parameter sets across different moneyness positions and maturities. We use the same moneyness definition as in Section 4.1. From Table 11, we see that the pricing errors for the SVJ model are similar to the case of the Heston model. In fact, the performance of our approach may even be slightly better in the case of the SVJ model. This could be because with the introduction of the jump component, the density for a given maturity is more spread out, that is, the turns are smoother for CDF, making it even easier for NN to approximate the density.

Table 11. Pricing Errors of the NN Approximated Density of Stochastic Volatility Jump Diffusion Model Compared with SVJ’s Semiclosed Form Put Pricing Solutions

Remark 1.

It is known (Baron 1994) that, for problems satisfying certain assumptions, increasing the size of the NN and the number of training points causes the NN approximation to tend asymptotically toward the exact solution of the problem. However, in contrast with conventional numerical methods (e.g., finite difference or Monte Carlo), there is currently no theory to guarantee the accuracy of a NN approximation. In this paper, we aim to illustrate the feasibility of the proposed computational approach and, consequently, we mostly consider models for underlyings where closed form or semiclosed form TPDFs or option prices are known. When such knowledge is not available, the following benchmarking is suggested. We train a neural TPDF for a given model of underlying and compute option prices using, for example, QUAD for a small set of parameters. Then for the same set of parameters we compute option prices accurately using conventional numerical methods where there is a well-established theory on controlling the error of approximations. If the option prices produced by the neural TPDF are close to accurate results obtained by the conventional methods, then we declare that the neural TPDF is sufficiently accurate to use as a neural TPDF generator. The generator is available for all points in the parametric space on which the NN is trained, whereas the expensive sanity check of accuracy is done using traditional numerical methods for a small set of parameters just once and before the NN is considered finally trained (“offline”). Such a sanity check was illustrated on the SABR model in Section 4.2. Within our approach, a neural TPDF itself and the results produced based on it do not have a statistical error because we fix the NN parameters that the training suggests are best and use these afterward as a neural TPDF generator posttraining (“online”); in other words, the error of a neural TPDF is deterministic. This is in contrast to how deep learning can be used to speed up Monte Carlo simulations in option pricing (Hinds and Tretyakov 2023).

Remark 2.

Here, we illustrate the use of neural TPDFs in conjunction with QUAD to evaluate option prices. Because QUAD is quadrature based, it is only effective for relatively low-dimensional situations (i.e., a small number of underlyings, demonstrated under geometric Brownian motion for barrier options with 5 underlyings by Andricopoulos et al. (2007) who note that it is overtaken by Monte Carlo at around 15 underlyings). To use neural TPDFs for option prices in higher-dimensional situations, one can exploit sparse grid methods (Smolyak 1963, Garcke and Griebel 2013).

5. Conclusion

This paper is the first in the literature to use a deep learning approach in approximating TPDF for any model of the underlying in a parametric fashion. Our approach not only provides the required TPDF approximations but also its “single solve” capability across all parameters means that effort is expended only once in a single precomputation, delivering a parametric neural TPDF generator that does not require recalibration.

Instead of aiming at directly solving the forward Kolmogorov (Fokker-Planck) equation for the TPDF, which would require handling the initial condition (a “nasty” Dirac delta function), we use an ingenious way of first solving the backward Kolmogorov equation for the cumulative probability function which involves the terminal condition, a more “tamable” step function, before calculating the TPDF by differentiating the cumulative probability. We train each artificial NN by optimizing a purposefully constructed cost function, thereby solving the parametric partial differential equation problem for the cumulative distribution function of any given underlying price process. For this, we can use as many data points as needed to train the NN such that it sufficiently accurately approximates the TPDF. Once the NN is trained, we have an ultra-fast neural TPDF generator, tailored to the underlying model and portable to other calculation setups and computers.

We illustrate the approach first on the simplest Black-Scholes-Merton setup under geometric Brownian motion, noting that TPDF approximations for other one-dimensional models of underlyings can likewise be obtained. Next, we use the model of Heston (1993) as an example of stochastic volatility models which can be handled by the deep learning approach. We also show that this approach can deal with other scenarios including nonaffine processes such as the SABR model of the underlying (benchmarked against Monte Carlo simulation for option prices). Finally, we demonstrate how neural TPDFs can be obtained for jump diffusion models (Kou’s double exponential jump-diffusion (Kou 2002) and stochastic volatility jump diffusion (Bates 1996)) by solving in a parametric fashion the backward Kolmogorov equations taking the form of partial integro-differential equations.

The applications of learning the TPDF using deep NNs go beyond those illustrated, such as simulating likelihood estimators and default probability calculations in credit risk management (Filipović et al. 2013). Options are ubiquitous in finance, from the large range of directly available options on exchanges and over the counter, to credit risk and Real Options. During the five decades since the founding framework of Black, Scholes and Merton, practitioners and academics have developed ways of dealing with option features (such as early exercise) in combination with models of the underlying but often either the models are too simple or the solutions too hard. Although parallel computing, using GPUs, speeds up conventional approaches (indeed, investment banks may use hundreds or even thousands of these to accelerate parallel calculations) nevertheless, massive calculations still need to be repeated in frequent recalibrations. The deep learning approach we have introduced, resting on finding parametric neural TPDFs, avoids this.

Our paper also makes contributions to the currently fast developing area of deep learning algorithms used for solving differential equations, which includes demonstrating how to computationally efficiently find solutions to the Fokker-Planck equations, both for diffusion and jump-diffusion models in a parametric fashion; extensive testing to demonstrate how to choose cost functions for training NNs with higher accuracy; comparing the performance of different NN architectures and the suitability of different GPUs for the task.

There is room for improvement in our method, as demonstrated by the analysis of pricing errors. A limitation of the approach presented here is its accuracy in terms of absolute errors, which means it may struggle with precise approximation of values very close to 0.0, resulting in high relative errors (percentage errors). Future research should address this issue. Additionally, the use of the DGM NN (Sirignano and Spiliopoulos 2018) in our approach does not necessarily imply that it is the optimal network architecture for solving the underlying problems (see Online Appendix B for a brief comparison between DGM and MLP networks for our approach). There is potential for further investigation into the search for or development of a better and more effective network architecture. Moreover, an interesting avenue for future research involves exploring model calibration modified through using market data.

The introduction of targeted deep learning for TPDF, as initiated in this paper, holds promise for the development of comprehensive and efficient programs capable of addressing a wide range of models and applications.