# Novel ANN Method for Solving Ordinary and Time-Fractional Black–Scholes Equation

Authors: Saeed Bajalan; Nastaran Bajalan
Journal: Complexity
 (2021)

---

## Abstract

The main aim of this study is to introduce a 2-layered artificial neural network (ANN) for solving the Black–Scholes partial differential equation (PDE) of either fractional or ordinary orders. Firstly, a discretization method is employed to change the model into a sequence of ordinary differential equations (ODE). Subsequently, each of these ODEs is solved with the aid of an ANN. Adam optimization is employed as the learning paradigm since it can add the foreknowledge of slowing down the process of optimization when getting close to the actual optimum solution. The model also takes advantage of fine-tuning for speeding up the process and domain mapping to confront the infinite domain issue. Finally, the accuracy, speed, and convergence of the method for solving several types of the Black–Scholes model are reported.


---

## Body

## 1. Introduction

Both PDEs of the ordinary and fractional order play an important role in pricing of financial derivatives. PDEs of the ordinary order are the basis of various models proposed for pricing of different types of options. On the contrary, financial markets show fractal behavior [1–4] and fractional PDEs (FPDE) which can better reflect that the reality of them have gained a lot of interest recently. Hence, finding an accurate and efficient approach for solving both types is a critical issue in pricing.The most famous PDE in finance is the Black–Scholes (B-S) model, which is broadly adopted for option pricing. So far, studies have presented different approaches for finding the numerical solution of this model and its variation when the exact form does not exist [5–13]. By using tick-by-tick data, Cartea and del-Castillo-Negrete [14] found that the value of European-style options satisfies a FPDE containing a nonlocal operator in time-to-maturity known as the Caputo fractional derivative. Furthermore, in order to depict the fractal structure in the financial market, the classical B-S equation has been generalized by means of fractional derivatives which have the property of self-similarity and better modeling of long-range dependency [15–17]. Although the fractional B-S is a powerful tool for explaining the hereditary and memory characteristics of the financial market, it is considerably difficult to obtain an accurate solution for it, due to the memory trait of fractional derivatives [18]. As a result, numerous researchers have tried techniques for approximating such problems. Among different analytical models presented so far to solve time-fractional B-S, the most cited articles employed integral transform [15, 19–23], wavelet-based hybrid methods [24], the separation of variables [25], the homotopy analysis and homotopy perturbation methods [26–28], and Fourier Laplace transform [29]. However, due to the high computational complexity of these solutions, numerical methods are often better alternatives for solving such mathematical models.Cartea and del Castillo-Negrete used backward difference technique together with shifted Grunwald–Letnikov scheme to solve the spatial fractional FMLS process [14]. Investigation of convergence analysis and comparison of the solution of three space fractional B-S modes were provided by Marom and Momoniat [30]. Finite differences’ methods were used to provide a solution for the time-fractional B-S model in 2013 by Song and Wang [31]. Zhang et al. also proposed a second-order finite difference method in the following year [32]. Moreover, in 2017, the weighted finite difference method was utilized to find the numerical solution of the model [33]. Chen et al. investigated predictor-corrector approaches for the solution of American option in 2015 [34]. Khan and Ansari [35] were the first to use the Sumudu transform to solve the fractional model of European options. In the same year, Zhang provided an unconditionally stable implicit numerical scheme for the model by changing the Riemann–Liouville derivatives to the Caputo derivative [32]. As methods based on the radial basis function (RBF) are widely used for approximating the ordinary PDEs, the meshless method was proposed for finding the solution of the time-fractional method as well [18]. In 2020, the fractional model was numerically approximated by using Quintic and sixth order B-spline functions as the basis for a collocation method providing high accuracy for the generalized B-S mode [36, 37]. In order to provide to sixth order accuracy in solving the generalized Black–Scholes model, Roul and Goura used both of the Crank–Nicolson scheme and sextic B-spline collocation method.Regarding machine learning and ANN in particular, ANN has been traditionally used for predicting option prices [38–43]; however, to the best of our knowledge except few research studies, there is no available literature for solving the B-S differential equation by ANN [44, 45]. Even these few research studies are limited to integer order PDEs, and there is no study for the solution of the fractional order B-S. In this paper, we are going to construct a 2-layered artificial neural network (ANN) to solve the Black–Scholes model of either time-fractional or ordinary orders. The Adaptive Moment Estimation (Adam), which has been specifically created to be used by neural networks, acts as the optimizer in the ANN to add the foreknowledge of slowing down when getting close to the optimal solution. To make the process faster and increase the efficiency of the method, fine-tuning is applied to the model. Also, to overcome the problem of the infinite problem domain for approximation domain mapping is used to map the whole problem to a finite interval. The rest of paper is organized as follows: Section 2 is dedicated to problem formulation, Section 3 explains the methodology, Section 4 presents the numerical results, and finally, the conclusion is provided in Section 5.


## 2. Problem Formulation

A put option on an underlying assetS is said to follow a Geometric Brownian Motion (GBM), where σ, r, and W are the volatility, interest rate, and Brownian motion, respectively, if it obeys the stochastic differential equation as follows:(1)dS=rSdt+σSdw.Using equation (1) and risk-neutral valuation formula together with the classic Feynman–Kac formula, the Black–Scholes operator is formed as below:(2)LUS,t=−∂US,t∂t−σ22S2∂2US,t∂S2−rS∂US,t∂S+rUS,t,where US,t is the unknown function which determines the option price Robbins and Monro [46]. It has been specified that this option will have a certain payoff at a certain date in the future, depending on the value(s) taken by the stock up to that date.It is well known now that a time-fractional Black–Scholes equation with the derivative of real orderα can be obtained to describe the price of an option under several circumstances such as when the change in the underlying asset is assumed to follow a fractal transmission system. Fractional derivatives, as they are called, were introduced in option pricing in a bid to take advantage of their memory properties to capture both major jumps over small periods of time and long-range dependencies in markets. Therefore, the fractional Black–Scholes model can be formulated as follows:(3)LUS,t=−DαUS,t+γ1∂2US,t∂S2+γ2∂US,t∂S+γ3US,t+fS,t,where D denotes the fractional derivative, α is a real number, and γi,i=1,2,and3, and fS,t are functions dependent on the values of σ, r, and α written using these notations for simplicity. It can be shown that the time derivative DαUS,t appearing in 3 equals the α-order Caputo fractional derivative by Chen et al. [47].Based on the type of the option, the corresponding condition set is as follows:(4)US,T=gS,Ua,t=Mat,Ub,t=Mbt.In some cases, e.g., European,Mbt moves toward infinity, thus the problem domain is semi-infinite. Some possible strategies are defined in Section 3.4 to overcome this obstacle when pricing. The corresponding parameters and their definitions are addressed in Table 1.Table 1
Option parameters and their definitions.
ParameterDefinitionTMaturity (years)KStrike priceSUnderlying asset (stock price)rInterest rateDFractional derivativeαReal-valued derivative orderσVolatilitygSPayoff function


## 3. Methodology

In this section, four main concepts employed in the present approach are explained. Then, the necessary steps to be taken for finding the solution are combined to form the proposed method.


### 3.1. Time Discretization

Solving multivariable equations increases the time complexity and the risk of producing inconsistent answers by computational software. In summary, time discretization methods are useful tools for converting such models into a series of ordinary differential equations (ODE). This approach is the result of applying finite difference methods on one dimension of an equation to approximately calculate the value of derivatives with respect to that dimension.Since here fractional equations are also investigated, ordinary and fractional time discretization approaches are discussed. Suppose that the PDE to be solved is defined onS,t∈a,b∗0,T. Based on this method, US,t in each time step is defined as UiS=US,iΔt, where US,t is the answer, Δ=T−0/N, and N is the number of time steps.


#### 3.1.1. Ordinary Time Discretization

Consider the PDE to be solved as follows:(5)∂US,t∂t=Ω∂US,t∂S,∂2US,t∂S2,US,t,S,t,where US,t is the answer and Ω is a linear or nonlinear function. Instead of solving this problem on the two dimensions, it can be converted to a series of dependent ODEs.By using the definedUiS, the following equation is constructed:(6)∂UiS∂t=θΩ∂UiS∂S,∂2UiS∂S2,UiS,S,iΔt+1−θΩ∂Ui−1S∂S,∂2Ui−1S∂S2,Ui−1S,S,i−1Δt.The method is implicit, ifθ=1. In this case, only posterior time step is used. The method is called explicit if θ=0 where only computing time step is utilized to approximate the solution. If the value of θ is equal to 1/2, the method is the common Crank–Nicolson method which is unconditionally stable and of the second order in time, and it uses both posterior and computing time steps for approximating the solution of the model. Due to the nonsmoothness of the payoff function and the activation functions in our ANN, the Crank–Nicolson cannot reach its second-order convergence. It can also cause extra inconsistencies because of the same problem. Hence, from this point on, θ=0 is considered.


#### 3.1.2. Fractional Time Discretization

Suppose that the time-fractional PDE (FPDE) to be solved is as follows:(7)DtαUS,t=Ω∂US,t∂S,∂2US,t∂S2,US,t,S,t,where α denotes the Caputo derivatives of the function and 0<α<1. As the first step for such FPDEs, the Caputo derivative should be discretized (readers are advised to see Hadian Rasanan et al. [48] for the preliminaries and thorough information about this type of derivative). Consider the following theorem.Theorem 1.
Suppose that0,T is divided to N parts with step size of Δt=T/N, 0<α<1, and qt∈C20,tk where tk=kΔt, and the following holds for this interval:(8)1Γ1−α∫0tkq′tt−tkαΔt−Δt−αΓ2−αb0qx,tk−∑m=1k−1bk−m−1−bk−mqtm−bk−1qt0≤1Γ2−α1−α12+22−α2−α−1+2−αmax0≤t≤tkq″tΔt2−α,where bm=m+11−α−m1−α.Proof.
See the proof in the work of Suna and Wub [49].According to the above theorem, equation (3) can be discretized in the following form:(9)DtαUn+1S≈1Γ2−α∑m=0nU^n+1−jS−U^n−jSΔtαbm=HU^n+1S.Now, the unknown function,Ui,i=0,…,n, should be approximated in each time step such that it satisfies equation (3) and also its initial and boundary conditions in equation (4). In this regard, the boundary conditions can be satisfied by considering UNS=gS in computations. On the contrary, to satisfy the boundary condition, the sum of least square error methods is used in Section 4.The remaining step is satisfying equation (3), so for this purpose, the cost function is chosen as below:(10)CostS,W=12Nr∑i=1r∑n=0N−11Γ2−α∑m=0nU^Si,tn+1−m−U^Sm,tn−mdtαbm−HU^Si,tn+12,where S=S1,S2,…,Sr and Si is the ith training data. To find the optimum weights for the network, this cost function should be minimized subject to W, so the following nonlinear least square problem is obtained:(11)minW12Nr∑i=1r∑n=0N−11Γ2−α∑m=0nU^n+1−mSi−U^n−mSmdtαbm−HU^n+1Si2.It is noteworthy that whenα=1, the ordinary time discretization method will be used.It should be noted that the method is based on the discretization of Caputo-type derivatives, and the result based on other types is only possible if the discretization is possible for other types.


### 3.2. Function Approximation

According to the universal approximation theorem, every continuous function can be approximated by a feedforward neural network [50]. This theorem states that any linear function can be approximated by an ANN without any hidden layers. But for functions of higher orders, the approximation can be well-established, if the ANN has at least one hidden layer. To calculate the price of an option based on equations (2) and (3), the value of option at each time step using a 2-layered network should be approximated as follows:(12)NW,S=ΨVφWS+B0+B1,where n is the number of neurons in the hidden layer, B0=b1,b2,…,bn, B1=β1, and φi,i∈1,2,…,n, are the activation functions in the hidden layer, and Ψ is the activation function for the output layer. The above formula can be seen in Figure 1.Figure 1
The topology of the network used for solving the Black–Scholes model.The nodes commensurate with the edgesbi and β1 in this network are called the biases whose inputs are unchangeable 1s, and their weights are additional parameters as means of adjusting the output of the next layer. In other words, they help the network fit best for the given data.The most famous activation function in deep learning and neural networks is sigmoid.The two main reasons for the broad application of the sigmoid function are the straight-forward calculations of its derivatives and its bounded value. The sigmoid function is defined as follows, and its values are in0,1:(13)sigmoidx=11+exp−x.These features make it a perfect candidate for problems that produce probabilities and for the Black–Scholes model. The values of options are nonnegative, and the effect of other parameters will not increase the calculations as they are multiplied by smaller values produced by sigmoid. On the contrary, the derivative of this well-known activator, its slope, is easily calculable between any two arbitrary points:(14)sigmoid′x=exp−x1+exp−x2.Although linear functions such as identity are not favorable for hidden layers as they take away the chance to generalize and adapt from the network, it is possible to use them as the activator of the output layer; hence, the hidden layers are present and are directly interacting with inputs.


### 3.3. Fine-Tuning

One of the key factors in the present study is the possibility of applying the fine-tuning methods during the training process. Fine-tuning is employing a previously trained neural network to find the solution of a new similar task. This process is normally applied to datasets related to images and voices. However, following the same approach, it is possible to increase the accuracy and speed of the network in this work.Building and validating an ANN from scratch can be a huge task in its own right, depending on what data being trained on it, and many parameters such as the number of layers and the number of nodes in hidden layers, the proper activation functions, and learning rate should be found through trial and error. If a trained model that already does one task well exist and that task is similar to ours in at least some remote way, then everything the model has already learned can be taken advantage of and applied to the new specific task. If the task is completely similar, like what we are facing when solving the problem at different time steps, the exact weights can be used as the initial values. If the models are somewhat similar, still some knowledge exists on the previous network which is notable for speeding up the process of building/modifying and training the network for the new task. Then, the only job remains for the network is learning the new features and properties that were not available in the former task. Here, once the network is trained for the first time steps, the obtained parameters, weights, and biases can be effectively reemployed for training the data fed to the network in other time steps. The approximated solution and its convergence rate are compared inSection 5.


### 3.4. Domain Mapping

Considering the vulnerability of neural networks due to the bounded domain of their activation functions in calculations on infinite domains, the domain mapping approach is utilized to shift the problem from its semi-infinite domain to a finite interval. This helps to prevent the error caused by common solutions such as truncation of the domain.On finite domains,S=x will be considered, but in semi-infinite domains, transformation formulas should be used for shifting the problem to a desired finite one. Here, xS=2/πarctanS/L is used to shift the problem’s domain which is 0,∞ to 0,1, in which L is the characteristic length of the mapping Boyd [51].Here,L is chosen in a way that 60% of all training points stand before the mapped strike price because the price significantly differs from zero in 0,K (Rad et al. [52]). It means that, by defining l as an indicator for 0.6, in the view of the fact that these points are equidistant, L is computed as follows:(15)L=Ktanπ/2l.First, let us introduce the following notations:(16)US,t=U˜x,t,S=Ltanπ2x,ϒ≜∂S∂x=Lπ2cos2π/2x,Θ≜∂Υ−1∂x=−2cosπ/2xsinπ/2xL.Hence, the derivatives needed for the calculations according toSection 2 are(17)∂US,t∂S=∂U˜x,t∂x∂x∂S=1Υ∂U˜x,t∂x,∂2US,t∂S2=∂∂S∂US,t∂S=1Υ2∂2U^x,t∂x2+ΘΥ∂U^x,t∂x.By substituting equation (17) in equations (3) and (2), the domain of the obtained Black–Scholes model is 0,1.Since the transformation is applied to the whole problem, the payoff function and boundary conditions of equation (4) should be changed as well.


#### 3.4.1. Discussion

Using domain mapping helps to approximate the answer on the whole interval of the problems. On the contrary, the number of training data points needed for solving the model on smaller intervals is significantly less as can be seen in Section5. However, it can decrease the accuracy since the whole domain is being compacted into one small domain and loss of information may occur, meaning that truncating the domain causes a perfect approximation on a subdomain of problems but with acceptance of a bit of loss in accuracy, the whole domain can be covered. So, there lays a trade-off between these two methods.


### 3.5. Adaptive Moment Estimation Learning

Regression modeling is used to determine coefficients of mathematical functions based on empirical data. The method of least squares determines the coefficients such that the sum of the squares of the deviations between the data and the curve fit is minimized. Finding a satisfactory solution to nonlinear least square problems is one of the famous topics among scientists who work on nonlinear systems of equations. For minimizing a vector function,Λx, that is, matrix A is defined as Λ:ℝn⟶ℝm, and m≥n with respect to a predefined x=x1,x2,…,xn, that is to say, x∗∈ℝn is found in a way that(18)ϒx=12∑i=1mΛix2=12Λx2=12fTx⋅fx,in which(19)x∗=minxΥx.Several methods have been introduced for solving this nonlinear least square model so far. As we can see the same vector can be constructed when solving differential equations, similarly, for the Black–Scholes models, in each time step, the final goal is finding the proper network weights for solving the following system of equations:(20)eq1=US,ti−U^S,ti=0,i=0,…,M,eq2=U^0,ti−fS=0,i=0,…,M,eq3=limx⟶boundaryU^S,ti−gS=0,i=0,…,M.The best possible solution to this system of equations is calculated when the sum of the squares of these equations gets smaller. Therefore, the problem is converted into an optimization problem. By minimizing the following equation, the appropriate weights and biases for our network are found:(21)objectiveU^S,t=minW,V,B1,B2eq12+eq22+eq32.Gradient descentGD is the most famous iterative algorithm employed as a learning paradigm to solve regression problems. In GD, after initializing the weights, the gradients, G, of the cost function is calculated. The cost function is the sum square error of the output based on the desired output for each member of the training dataset. Then, based on G, the weights become updated W=W−ηG. This process is repeated by considering the new values as the initial ones until the cost function is desirably minimized. η is known as the learning parameter and is used to balance the rate of increase or decrease in each iteration. It should be noted that the only constraint on this value is 0<η<1. Various methods and theorems are introduced to find the boundaries of this value according to the problem. Some suggest using a big value and decreasing its value as the result approaches the correct value. In contrast, some others suggest choosing a very small value and then increasing it exponentially in time when the correct direction is found. But generally, they all prefer making this parameter a function of time, while none of them propose a single formula to calculate the exact amount of it for the best approximation [46, 53, 54].GD family has different optimizers such as Stochastic Gradient Descent (SGD), Adaptive Moment Estimation (Adam), Root Mean Square Propagation (RMSprop), and Nesterov Accelerated Gradient (NAG) which are mostly used in deep learning because of their speed and strength according to various control parameters such as the size of the training datasets and the pattern in which the training data is scattered.In GD for updating only one parameter, all available samples in the dataset should be visited; however, in SGD [55], minibatches, which are small subsets of the whole dataset, are used to update a single parameter. For relatively large datasets, this causes the algorithm to converge faster. GD is an actual optimizer trying to find the exact gradients while in SGD; as explained, the algorithm only approximates the gradients and not the precise value. Since SGD fluctuates a lot, due to frequent updates with high variance, it shows a paradoxical behavior. It can explore new and potential directions to find the minimum but, at the same time, this behavior puts the network in danger of completely missing the local or global minimum.A solution was proposed by the father of propagation, Geoffrey Hinton Tieleman and Hinton [56]. This thorough study which has not been academically submitted or published gained a lot of attention. The proposed algorithm fights the possibilities of vanishing or exploding the gradients. In other words, RMSprop normalizes the gradient using a moving average of squared gradients. This normalization balances the step size, reducing the step size for large gradients to avoid exploding and increasing it for small ones to avoid vanishing. Since this approach uses the exponentially decaying average, it is related to the most recent gradients, so the past gradient would not play a great role in updating the parameters. This leads to slow changes in the learning rate; however, it is relatively faster than GD.So far, it can be seen that RMSProp and SGD are the best options. Adam is an adaptive algorithm which is generally considered as the combination of these two paradigms with momentum. This methodology has been specifically created to be used by neural networks.Like RMSprop, Adam employs squared gradients to modify the learning rate. Also, the first and second moments are utilized using the moving average of the gradient such as SGD. However, a specific learning rate is calculated for each network parameter (weights) using two hyperparameters. Here, a summary of how Adam optimization works for the present model is provided. For the complete explanation, readers are encouraged to study Kingma and Ba [57]. The convergence of the method has been described in several great papers. But finally, all of the studies confirm the convergence proof provided in the first papers [58].All computations are done using autograd package of Python. In this package, Adam optimization is implemented as follows:def adam (grad,w, callback = None, num_iters = 100,step_size = 0.001,b1 = 0.9, b2 = 0.999, eps = 10∗∗−8):m=np⋅zeroslenwv=np⋅zeroslenwfori in range (num_iters):g=gradw,iif callback: callbackw,i,gm=1−b1∗g+b1∗m# First moment estimate.v=1−b2∗g∗∗2+b2∗v# Second moment estimate.mhat=m/1−b1∗∗i+1# Bias correction.vhat=v/1−b2∗∗i+1x=x−step_size∗mhat/np⋅sqrtvhat+eps.returnxAs the name of the method describes, it is derived from adaptive moment estimation.nth moment of a random variable is defined as the expected value of that variable to the power of n:(22)mn=Ewn,where m shows the nth moment and w is a random variable. The gradient of the cost function of the neural network can be considered a random variable since it usually evaluated on some small random batch of data. The first moment is mean, and the second moment is uncentered variance. To estimate the moments, Adam utilizes exponentially moving averages and computed on the gradient evaluated on a current minibatch:(23)mi=β1mi−1+1−β1gi,vi=β2vi−1+1−β2gi2,where m and v denote the moving averages and g is the gradient of the current data presented to the network. According to Kingma and Ba [57], which is also mentioned in the above snippet from autograd package, the values of hyperparameters β1 and β2 have two default values of 0.9 and 0.999, respectively. While the authors did not discuss the choosing process of these two variables, all studies reported very promising and, in most cases, perfect estimations using these two default values (see also [59]). The vectors of moving averages are initialized with zeros at the first iteration.The remaining problem with these moments was being biased towards zero sincemi and vi are initialized as vectors of 0’s. In other words, especially during the initial epochs and when the decay rates are small (i.e., β1 and β2 are close to 1), the values of mi and vi will not change significantly or even at all. So, the authors proposed the following bias corrections in order to surmount this obstacle:(24)m^i=m1−β1i,v^i=v1−β2i.Now, for each of the parameters (weights), a specific updating rule can be created:(25)w=w−ηm^v^+ϵ,where ϵ is an control parameter preventing the fractional part from producing a division by zero error.Different scientific studies have shown that Adam outperforms other methods. According to empirical practices, this method has better performance and accuracy. This is also discussed inSection 5. One problem that is stated by many studies is the convergence of the method. However, Kingma and Ba [57] provided the analysis for the convex problems; other papers argued the convergence of the method on a few nonconvex problems. And, with some modification, they finally agreed on its usability.In [60], the full analysis of the convexity of the Black–Scholes model is proposed. Due to differences such as the failure of put-call parity in real markets instead of theory, this paper proves that, for all American options, they preserve their convexity in bubbled markets as well as nonbubbled ones. They showed that European options are convexity preserving only for bounded payoffs. Thus, in this respect, the prices of American options are more robust than their European counterparts. In the same study, it is shown that models for bubbles are convexity preserving for bounded contracts. More precisely, consider x,t∈0,∞×0,T, and let u1x,t and u2x,t be the option prices such that their corresponding volatilities are nonnegative α1 and α2 which satisfy α1x,t≤α2x,t.Theorem 2.
Assume thatg is concave. Then, ux,t is concave in x for any t∈0,T. Moreover, the option price is decreasing in the volatility, that is, u1x,t≥u2x,t for all x,t∈0,∞×0,T. Similarly, if g is convex and bounded, then ux,t is convex in x for any t∈0,T. Moreover, the option price is increasing in the volatility, that is, u1x,t≤u2x,t for all x,t∈0,∞×0,T.The full proof is available in [60]. As they mentioned, the proof is valid under the assumption of the uniqueness of the result for such an option, which is proved in their thorough study on properties of Black–Scholes models in more realistic markets as well.


## 3.1. Time Discretization

Solving multivariable equations increases the time complexity and the risk of producing inconsistent answers by computational software. In summary, time discretization methods are useful tools for converting such models into a series of ordinary differential equations (ODE). This approach is the result of applying finite difference methods on one dimension of an equation to approximately calculate the value of derivatives with respect to that dimension.Since here fractional equations are also investigated, ordinary and fractional time discretization approaches are discussed. Suppose that the PDE to be solved is defined onS,t∈a,b∗0,T. Based on this method, US,t in each time step is defined as UiS=US,iΔt, where US,t is the answer, Δ=T−0/N, and N is the number of time steps.


### 3.1.1. Ordinary Time Discretization

Consider the PDE to be solved as follows:(5)∂US,t∂t=Ω∂US,t∂S,∂2US,t∂S2,US,t,S,t,where US,t is the answer and Ω is a linear or nonlinear function. Instead of solving this problem on the two dimensions, it can be converted to a series of dependent ODEs.By using the definedUiS, the following equation is constructed:(6)∂UiS∂t=θΩ∂UiS∂S,∂2UiS∂S2,UiS,S,iΔt+1−θΩ∂Ui−1S∂S,∂2Ui−1S∂S2,Ui−1S,S,i−1Δt.The method is implicit, ifθ=1. In this case, only posterior time step is used. The method is called explicit if θ=0 where only computing time step is utilized to approximate the solution. If the value of θ is equal to 1/2, the method is the common Crank–Nicolson method which is unconditionally stable and of the second order in time, and it uses both posterior and computing time steps for approximating the solution of the model. Due to the nonsmoothness of the payoff function and the activation functions in our ANN, the Crank–Nicolson cannot reach its second-order convergence. It can also cause extra inconsistencies because of the same problem. Hence, from this point on, θ=0 is considered.


### 3.1.2. Fractional Time Discretization

Suppose that the time-fractional PDE (FPDE) to be solved is as follows:(7)DtαUS,t=Ω∂US,t∂S,∂2US,t∂S2,US,t,S,t,where α denotes the Caputo derivatives of the function and 0<α<1. As the first step for such FPDEs, the Caputo derivative should be discretized (readers are advised to see Hadian Rasanan et al. [48] for the preliminaries and thorough information about this type of derivative). Consider the following theorem.Theorem 1.
Suppose that0,T is divided to N parts with step size of Δt=T/N, 0<α<1, and qt∈C20,tk where tk=kΔt, and the following holds for this interval:(8)1Γ1−α∫0tkq′tt−tkαΔt−Δt−αΓ2−αb0qx,tk−∑m=1k−1bk−m−1−bk−mqtm−bk−1qt0≤1Γ2−α1−α12+22−α2−α−1+2−αmax0≤t≤tkq″tΔt2−α,where bm=m+11−α−m1−α.Proof.
See the proof in the work of Suna and Wub [49].According to the above theorem, equation (3) can be discretized in the following form:(9)DtαUn+1S≈1Γ2−α∑m=0nU^n+1−jS−U^n−jSΔtαbm=HU^n+1S.Now, the unknown function,Ui,i=0,…,n, should be approximated in each time step such that it satisfies equation (3) and also its initial and boundary conditions in equation (4). In this regard, the boundary conditions can be satisfied by considering UNS=gS in computations. On the contrary, to satisfy the boundary condition, the sum of least square error methods is used in Section 4.The remaining step is satisfying equation (3), so for this purpose, the cost function is chosen as below:(10)CostS,W=12Nr∑i=1r∑n=0N−11Γ2−α∑m=0nU^Si,tn+1−m−U^Sm,tn−mdtαbm−HU^Si,tn+12,where S=S1,S2,…,Sr and Si is the ith training data. To find the optimum weights for the network, this cost function should be minimized subject to W, so the following nonlinear least square problem is obtained:(11)minW12Nr∑i=1r∑n=0N−11Γ2−α∑m=0nU^n+1−mSi−U^n−mSmdtαbm−HU^n+1Si2.It is noteworthy that whenα=1, the ordinary time discretization method will be used.It should be noted that the method is based on the discretization of Caputo-type derivatives, and the result based on other types is only possible if the discretization is possible for other types.


## 3.1.1. Ordinary Time Discretization

Consider the PDE to be solved as follows:(5)∂US,t∂t=Ω∂US,t∂S,∂2US,t∂S2,US,t,S,t,where US,t is the answer and Ω is a linear or nonlinear function. Instead of solving this problem on the two dimensions, it can be converted to a series of dependent ODEs.By using the definedUiS, the following equation is constructed:(6)∂UiS∂t=θΩ∂UiS∂S,∂2UiS∂S2,UiS,S,iΔt+1−θΩ∂Ui−1S∂S,∂2Ui−1S∂S2,Ui−1S,S,i−1Δt.The method is implicit, ifθ=1. In this case, only posterior time step is used. The method is called explicit if θ=0 where only computing time step is utilized to approximate the solution. If the value of θ is equal to 1/2, the method is the common Crank–Nicolson method which is unconditionally stable and of the second order in time, and it uses both posterior and computing time steps for approximating the solution of the model. Due to the nonsmoothness of the payoff function and the activation functions in our ANN, the Crank–Nicolson cannot reach its second-order convergence. It can also cause extra inconsistencies because of the same problem. Hence, from this point on, θ=0 is considered.


## 3.1.2. Fractional Time Discretization

Suppose that the time-fractional PDE (FPDE) to be solved is as follows:(7)DtαUS,t=Ω∂US,t∂S,∂2US,t∂S2,US,t,S,t,where α denotes the Caputo derivatives of the function and 0<α<1. As the first step for such FPDEs, the Caputo derivative should be discretized (readers are advised to see Hadian Rasanan et al. [48] for the preliminaries and thorough information about this type of derivative). Consider the following theorem.Theorem 1.
Suppose that0,T is divided to N parts with step size of Δt=T/N, 0<α<1, and qt∈C20,tk where tk=kΔt, and the following holds for this interval:(8)1Γ1−α∫0tkq′tt−tkαΔt−Δt−αΓ2−αb0qx,tk−∑m=1k−1bk−m−1−bk−mqtm−bk−1qt0≤1Γ2−α1−α12+22−α2−α−1+2−αmax0≤t≤tkq″tΔt2−α,where bm=m+11−α−m1−α.Proof.
See the proof in the work of Suna and Wub [49].According to the above theorem, equation (3) can be discretized in the following form:(9)DtαUn+1S≈1Γ2−α∑m=0nU^n+1−jS−U^n−jSΔtαbm=HU^n+1S.Now, the unknown function,Ui,i=0,…,n, should be approximated in each time step such that it satisfies equation (3) and also its initial and boundary conditions in equation (4). In this regard, the boundary conditions can be satisfied by considering UNS=gS in computations. On the contrary, to satisfy the boundary condition, the sum of least square error methods is used in Section 4.The remaining step is satisfying equation (3), so for this purpose, the cost function is chosen as below:(10)CostS,W=12Nr∑i=1r∑n=0N−11Γ2−α∑m=0nU^Si,tn+1−m−U^Sm,tn−mdtαbm−HU^Si,tn+12,where S=S1,S2,…,Sr and Si is the ith training data. To find the optimum weights for the network, this cost function should be minimized subject to W, so the following nonlinear least square problem is obtained:(11)minW12Nr∑i=1r∑n=0N−11Γ2−α∑m=0nU^n+1−mSi−U^n−mSmdtαbm−HU^n+1Si2.It is noteworthy that whenα=1, the ordinary time discretization method will be used.It should be noted that the method is based on the discretization of Caputo-type derivatives, and the result based on other types is only possible if the discretization is possible for other types.


## 3.2. Function Approximation

According to the universal approximation theorem, every continuous function can be approximated by a feedforward neural network [50]. This theorem states that any linear function can be approximated by an ANN without any hidden layers. But for functions of higher orders, the approximation can be well-established, if the ANN has at least one hidden layer. To calculate the price of an option based on equations (2) and (3), the value of option at each time step using a 2-layered network should be approximated as follows:(12)NW,S=ΨVφWS+B0+B1,where n is the number of neurons in the hidden layer, B0=b1,b2,…,bn, B1=β1, and φi,i∈1,2,…,n, are the activation functions in the hidden layer, and Ψ is the activation function for the output layer. The above formula can be seen in Figure 1.Figure 1
The topology of the network used for solving the Black–Scholes model.The nodes commensurate with the edgesbi and β1 in this network are called the biases whose inputs are unchangeable 1s, and their weights are additional parameters as means of adjusting the output of the next layer. In other words, they help the network fit best for the given data.The most famous activation function in deep learning and neural networks is sigmoid.The two main reasons for the broad application of the sigmoid function are the straight-forward calculations of its derivatives and its bounded value. The sigmoid function is defined as follows, and its values are in0,1:(13)sigmoidx=11+exp−x.These features make it a perfect candidate for problems that produce probabilities and for the Black–Scholes model. The values of options are nonnegative, and the effect of other parameters will not increase the calculations as they are multiplied by smaller values produced by sigmoid. On the contrary, the derivative of this well-known activator, its slope, is easily calculable between any two arbitrary points:(14)sigmoid′x=exp−x1+exp−x2.Although linear functions such as identity are not favorable for hidden layers as they take away the chance to generalize and adapt from the network, it is possible to use them as the activator of the output layer; hence, the hidden layers are present and are directly interacting with inputs.


## 3.3. Fine-Tuning

One of the key factors in the present study is the possibility of applying the fine-tuning methods during the training process. Fine-tuning is employing a previously trained neural network to find the solution of a new similar task. This process is normally applied to datasets related to images and voices. However, following the same approach, it is possible to increase the accuracy and speed of the network in this work.Building and validating an ANN from scratch can be a huge task in its own right, depending on what data being trained on it, and many parameters such as the number of layers and the number of nodes in hidden layers, the proper activation functions, and learning rate should be found through trial and error. If a trained model that already does one task well exist and that task is similar to ours in at least some remote way, then everything the model has already learned can be taken advantage of and applied to the new specific task. If the task is completely similar, like what we are facing when solving the problem at different time steps, the exact weights can be used as the initial values. If the models are somewhat similar, still some knowledge exists on the previous network which is notable for speeding up the process of building/modifying and training the network for the new task. Then, the only job remains for the network is learning the new features and properties that were not available in the former task. Here, once the network is trained for the first time steps, the obtained parameters, weights, and biases can be effectively reemployed for training the data fed to the network in other time steps. The approximated solution and its convergence rate are compared inSection 5.


## 3.4. Domain Mapping

Considering the vulnerability of neural networks due to the bounded domain of their activation functions in calculations on infinite domains, the domain mapping approach is utilized to shift the problem from its semi-infinite domain to a finite interval. This helps to prevent the error caused by common solutions such as truncation of the domain.On finite domains,S=x will be considered, but in semi-infinite domains, transformation formulas should be used for shifting the problem to a desired finite one. Here, xS=2/πarctanS/L is used to shift the problem’s domain which is 0,∞ to 0,1, in which L is the characteristic length of the mapping Boyd [51].Here,L is chosen in a way that 60% of all training points stand before the mapped strike price because the price significantly differs from zero in 0,K (Rad et al. [52]). It means that, by defining l as an indicator for 0.6, in the view of the fact that these points are equidistant, L is computed as follows:(15)L=Ktanπ/2l.First, let us introduce the following notations:(16)US,t=U˜x,t,S=Ltanπ2x,ϒ≜∂S∂x=Lπ2cos2π/2x,Θ≜∂Υ−1∂x=−2cosπ/2xsinπ/2xL.Hence, the derivatives needed for the calculations according toSection 2 are(17)∂US,t∂S=∂U˜x,t∂x∂x∂S=1Υ∂U˜x,t∂x,∂2US,t∂S2=∂∂S∂US,t∂S=1Υ2∂2U^x,t∂x2+ΘΥ∂U^x,t∂x.By substituting equation (17) in equations (3) and (2), the domain of the obtained Black–Scholes model is 0,1.Since the transformation is applied to the whole problem, the payoff function and boundary conditions of equation (4) should be changed as well.


### 3.4.1. Discussion

Using domain mapping helps to approximate the answer on the whole interval of the problems. On the contrary, the number of training data points needed for solving the model on smaller intervals is significantly less as can be seen in Section5. However, it can decrease the accuracy since the whole domain is being compacted into one small domain and loss of information may occur, meaning that truncating the domain causes a perfect approximation on a subdomain of problems but with acceptance of a bit of loss in accuracy, the whole domain can be covered. So, there lays a trade-off between these two methods.


## 3.4.1. Discussion

Using domain mapping helps to approximate the answer on the whole interval of the problems. On the contrary, the number of training data points needed for solving the model on smaller intervals is significantly less as can be seen in Section5. However, it can decrease the accuracy since the whole domain is being compacted into one small domain and loss of information may occur, meaning that truncating the domain causes a perfect approximation on a subdomain of problems but with acceptance of a bit of loss in accuracy, the whole domain can be covered. So, there lays a trade-off between these two methods.


## 3.5. Adaptive Moment Estimation Learning

Regression modeling is used to determine coefficients of mathematical functions based on empirical data. The method of least squares determines the coefficients such that the sum of the squares of the deviations between the data and the curve fit is minimized. Finding a satisfactory solution to nonlinear least square problems is one of the famous topics among scientists who work on nonlinear systems of equations. For minimizing a vector function,Λx, that is, matrix A is defined as Λ:ℝn⟶ℝm, and m≥n with respect to a predefined x=x1,x2,…,xn, that is to say, x∗∈ℝn is found in a way that(18)ϒx=12∑i=1mΛix2=12Λx2=12fTx⋅fx,in which(19)x∗=minxΥx.Several methods have been introduced for solving this nonlinear least square model so far. As we can see the same vector can be constructed when solving differential equations, similarly, for the Black–Scholes models, in each time step, the final goal is finding the proper network weights for solving the following system of equations:(20)eq1=US,ti−U^S,ti=0,i=0,…,M,eq2=U^0,ti−fS=0,i=0,…,M,eq3=limx⟶boundaryU^S,ti−gS=0,i=0,…,M.The best possible solution to this system of equations is calculated when the sum of the squares of these equations gets smaller. Therefore, the problem is converted into an optimization problem. By minimizing the following equation, the appropriate weights and biases for our network are found:(21)objectiveU^S,t=minW,V,B1,B2eq12+eq22+eq32.Gradient descentGD is the most famous iterative algorithm employed as a learning paradigm to solve regression problems. In GD, after initializing the weights, the gradients, G, of the cost function is calculated. The cost function is the sum square error of the output based on the desired output for each member of the training dataset. Then, based on G, the weights become updated W=W−ηG. This process is repeated by considering the new values as the initial ones until the cost function is desirably minimized. η is known as the learning parameter and is used to balance the rate of increase or decrease in each iteration. It should be noted that the only constraint on this value is 0<η<1. Various methods and theorems are introduced to find the boundaries of this value according to the problem. Some suggest using a big value and decreasing its value as the result approaches the correct value. In contrast, some others suggest choosing a very small value and then increasing it exponentially in time when the correct direction is found. But generally, they all prefer making this parameter a function of time, while none of them propose a single formula to calculate the exact amount of it for the best approximation [46, 53, 54].GD family has different optimizers such as Stochastic Gradient Descent (SGD), Adaptive Moment Estimation (Adam), Root Mean Square Propagation (RMSprop), and Nesterov Accelerated Gradient (NAG) which are mostly used in deep learning because of their speed and strength according to various control parameters such as the size of the training datasets and the pattern in which the training data is scattered.In GD for updating only one parameter, all available samples in the dataset should be visited; however, in SGD [55], minibatches, which are small subsets of the whole dataset, are used to update a single parameter. For relatively large datasets, this causes the algorithm to converge faster. GD is an actual optimizer trying to find the exact gradients while in SGD; as explained, the algorithm only approximates the gradients and not the precise value. Since SGD fluctuates a lot, due to frequent updates with high variance, it shows a paradoxical behavior. It can explore new and potential directions to find the minimum but, at the same time, this behavior puts the network in danger of completely missing the local or global minimum.A solution was proposed by the father of propagation, Geoffrey Hinton Tieleman and Hinton [56]. This thorough study which has not been academically submitted or published gained a lot of attention. The proposed algorithm fights the possibilities of vanishing or exploding the gradients. In other words, RMSprop normalizes the gradient using a moving average of squared gradients. This normalization balances the step size, reducing the step size for large gradients to avoid exploding and increasing it for small ones to avoid vanishing. Since this approach uses the exponentially decaying average, it is related to the most recent gradients, so the past gradient would not play a great role in updating the parameters. This leads to slow changes in the learning rate; however, it is relatively faster than GD.So far, it can be seen that RMSProp and SGD are the best options. Adam is an adaptive algorithm which is generally considered as the combination of these two paradigms with momentum. This methodology has been specifically created to be used by neural networks.Like RMSprop, Adam employs squared gradients to modify the learning rate. Also, the first and second moments are utilized using the moving average of the gradient such as SGD. However, a specific learning rate is calculated for each network parameter (weights) using two hyperparameters. Here, a summary of how Adam optimization works for the present model is provided. For the complete explanation, readers are encouraged to study Kingma and Ba [57]. The convergence of the method has been described in several great papers. But finally, all of the studies confirm the convergence proof provided in the first papers [58].All computations are done using autograd package of Python. In this package, Adam optimization is implemented as follows:def adam (grad,w, callback = None, num_iters = 100,step_size = 0.001,b1 = 0.9, b2 = 0.999, eps = 10∗∗−8):m=np⋅zeroslenwv=np⋅zeroslenwfori in range (num_iters):g=gradw,iif callback: callbackw,i,gm=1−b1∗g+b1∗m# First moment estimate.v=1−b2∗g∗∗2+b2∗v# Second moment estimate.mhat=m/1−b1∗∗i+1# Bias correction.vhat=v/1−b2∗∗i+1x=x−step_size∗mhat/np⋅sqrtvhat+eps.returnxAs the name of the method describes, it is derived from adaptive moment estimation.nth moment of a random variable is defined as the expected value of that variable to the power of n:(22)mn=Ewn,where m shows the nth moment and w is a random variable. The gradient of the cost function of the neural network can be considered a random variable since it usually evaluated on some small random batch of data. The first moment is mean, and the second moment is uncentered variance. To estimate the moments, Adam utilizes exponentially moving averages and computed on the gradient evaluated on a current minibatch:(23)mi=β1mi−1+1−β1gi,vi=β2vi−1+1−β2gi2,where m and v denote the moving averages and g is the gradient of the current data presented to the network. According to Kingma and Ba [57], which is also mentioned in the above snippet from autograd package, the values of hyperparameters β1 and β2 have two default values of 0.9 and 0.999, respectively. While the authors did not discuss the choosing process of these two variables, all studies reported very promising and, in most cases, perfect estimations using these two default values (see also [59]). The vectors of moving averages are initialized with zeros at the first iteration.The remaining problem with these moments was being biased towards zero sincemi and vi are initialized as vectors of 0’s. In other words, especially during the initial epochs and when the decay rates are small (i.e., β1 and β2 are close to 1), the values of mi and vi will not change significantly or even at all. So, the authors proposed the following bias corrections in order to surmount this obstacle:(24)m^i=m1−β1i,v^i=v1−β2i.Now, for each of the parameters (weights), a specific updating rule can be created:(25)w=w−ηm^v^+ϵ,where ϵ is an control parameter preventing the fractional part from producing a division by zero error.Different scientific studies have shown that Adam outperforms other methods. According to empirical practices, this method has better performance and accuracy. This is also discussed inSection 5. One problem that is stated by many studies is the convergence of the method. However, Kingma and Ba [57] provided the analysis for the convex problems; other papers argued the convergence of the method on a few nonconvex problems. And, with some modification, they finally agreed on its usability.In [60], the full analysis of the convexity of the Black–Scholes model is proposed. Due to differences such as the failure of put-call parity in real markets instead of theory, this paper proves that, for all American options, they preserve their convexity in bubbled markets as well as nonbubbled ones. They showed that European options are convexity preserving only for bounded payoffs. Thus, in this respect, the prices of American options are more robust than their European counterparts. In the same study, it is shown that models for bubbles are convexity preserving for bounded contracts. More precisely, consider x,t∈0,∞×0,T, and let u1x,t and u2x,t be the option prices such that their corresponding volatilities are nonnegative α1 and α2 which satisfy α1x,t≤α2x,t.Theorem 2.
Assume thatg is concave. Then, ux,t is concave in x for any t∈0,T. Moreover, the option price is decreasing in the volatility, that is, u1x,t≥u2x,t for all x,t∈0,∞×0,T. Similarly, if g is convex and bounded, then ux,t is convex in x for any t∈0,T. Moreover, the option price is increasing in the volatility, that is, u1x,t≤u2x,t for all x,t∈0,∞×0,T.The full proof is available in [60]. As they mentioned, the proof is valid under the assumption of the uniqueness of the result for such an option, which is proved in their thorough study on properties of Black–Scholes models in more realistic markets as well.


## 4. Numerical Results and Discussion

In this section, three test examples with exact solutions are chosen according to the previous works for examining the accuracy and efficiency of the proposed ANN. According to what is mentioned inSection 3.5, this approach is applicable to other kinds of options such as barriers and American options. All computations are performed using Python.3.7 software on a 2.7GHz Intel Core i7 CPU machine with 16GB of memory. Only one hidden layer with 20 hidden neurons is used in all of the samples. The initial weights are scattered in −0.01,0.01. The number of epochs for the first time stamp is 5000, and for the rest of the steps, thanks to fine-tuning, decreasing this number to 1200 provides very promising results. All values in 2000,7000 for the first time-step iterations and all values in 600,7000 do not lead to overfitting/underfitting, but the best results in our experiment achieved using the stated values. If the learning rate is low, then training is more reliable, but optimization will take a lot of time because steps towards the minimum of the loss function are small, and on the contrary, if the learning rate is high, then training may not converge or may even diverge. Weight changes can be so big that the optimizer overshoots the minimum and makes the loss worse. In this research, best values for the learning rate are found using the genetic algorithm.Example 1.
Let us consider a European call option, in which its interest rate, volatility, and strike price are 0.05, 0.2, and 10, respectively. The governing equation is similar to equation (31), and the boundary conditions set is as follows:(26)US,T=maxS−K,0,U0,t=0,limS⟶∞US,t=S−Kexp−rt.
With the maturity of 1, in years, the approximate price att=0 is shown in Figures 2(a) and 2(b). The exact solution of this call option can be obtained using the following analytical solution denoted by UexactS,t:(27)d1=lnS/K+r+1/2σ2T−tσT−t,d2=d1−σT−t,NS=12π∫−∞Sexp−12y2dy,UexactS,t=SNd1−Kexp−rT−tNd2.
When truncating the domain at 15, the cardinal of the training dataset is 150 containing equidistant points scattered between0,15 and η=0.03. The number of time steps, N, is 20, and the average calculation time per each epoch on the abovementioned configuration is 0.098s. Figure 2(a) demonstrates that when the problem domain is truncated, the approximate price is behaving fine until it reaches the truncation point which is 15 in these examples.Figure 2
Plots of the approximated solutions of Example1 when the number of hidden neurons is 20 and dt=1/20. (a) US,0. (b) US,t, t∈0,T.
(a)(b)To solve this issue, the problems are mapped to0,1 using Section 3.4; then, Ux,t is computed using the proposed ANN and sigmoid functions as the activation functions; then, Ux,t is reverted to the original model’s domain using the inverse mapping function so that US,t is calculated. Only 10 equidistant points are used as training points in this case, and the logarithmic absolute errors obtained from two approaches are compared in Figure 3. It should be noted that, after the truncation point, the error increases rapidly for the first approach, but when truncating the domain, the overall error is higher at the beginning of the interval but it remains steady and even falls at the end of the domain. Since the mapping function converges to infinity on x=1, the numerical calculation on software such as Python will not be able to perform the calculations. So, these comparisons are done using a very big value for x=0.9999999. Figure 4 confirms the fact that the Adam optimizer performs better than the other two optimizers, as it starts to converge and moves towards the answer in earlier epochs for the first time step. The average calculation time for SGD and RMSprop are, respectively, 0.45s and 0.63s per epoch. It is noteworthy that RMSprop crashed due to overflow encounters, and the depicted figure is just for comparing, with the learning rate of 0.01 instead of 0.03 which somehow might make the comparison unreliable. But the point is observed, and this method fails in comparison to the other methods for solving this type of the Black–Scholes model.Figure 3
Plots of the logarithmic absolute error for the solution of Example1 when the number of hidden neurons is 20 and dt=1/20, the logarithmic error (a) using the truncating approach and (b) using the mapping function.
(a)(b)Figure 4
Comparison of networks’ convergence toward the exact solution of Example1 using (a) Adam optimizer, (b) SGD optimizer, and (c) RMSprop η=0.01
(a)(b)(c)Figures5(a) and 5(b) illustrate the incredible influence of fine-tuning on the objective convergence.Figure 5
The effect of fine-tuning on the convergence rate of the method. Solving Example1 with (a) fine-tuning and (b) without fine-tuning.
(a)(b)Example 2.
Consider the following fractional model of a European option with homogeneous boundary conditions as follows:(28)DαUs,t=a∂2Us,t∂t2+b∂Us,t∂t−cUs,t+fS,t=0,U0,t=t,U1,t=0,US,0=S2S−1.
The interest rate and volatility are 0.05 and 0.25, respectively. Other variables are calculated based on these two variable,a=1/2sigma2, b=r−a, and c=r. The fractional order of the equation is α=0.7. Also,(29)fS,t=2t2−αΓ3−α+2t1−αΓ2−αS21−S−t+12a2−6S+b2S−3S2−cS21−S.
The exact solution for this equation is formulated as below:(30)UexactS,t=t+12S21−S.
The number of time steps,N, is 10, the number of hidden neurons in this example is 6, and the average calculation time per each epoch on the abovementioned configuration is 0.0012s. Also, η=0.03. The approximate solution calculated using the Adam neural network is plotted in Figure 6. However, smaller number of training data points, 20 and 40, still produced errors in 0,10−2, and 60 equidistant points are used as training points in this case to reach the logarithmic absolute errors illustrated in Figure 7(a), giving us the freedom to increase the accuracy even more. Because the number of epochs in this example is very small, the absolute errors obtained from SGD, RMSprop, and Adam are compared in Figure 7. Here, it can be seen that, with a small number of neurons and training points, the accuracy of the model is more promising than the other optimizers. The average calculation time for SGD and RMSprop are, respectively, 0.0059s and 0.17s per epoch.Figure 6
Plots of the approximated solutions of Example2 when the number of hidden neurons is 60 and dt=1/20. (a) US,0 (the red dotted plot shows the approximate solution). (b) US,t, t∈0,T.
(a)(b)Figure 7
Comparison of networks’ convergence toward the exact solution of Example2 using (a) Adam optimizer, (b) SGD optimizer, and (c) RMSprop.
(a)(b)(c)Figure8 illustrates the incredible influence of fine-tuning on the objective convergence. When fine-tuning is used, the objective function starts with a very small value, and hence, it converges rapidly even for very small values of the cost 10−4.Figure 8
The effect of fine-tuning on the convergence rate of the method. Solving Example2 with (a) fine-tuning and (b) without fine-tuning.
(a)(b)Example 3.
Let us consider a European put option, in which its interest rate, volatility, and strike price are 0.05, 0.2, and 10, respectively:(31)LUS,t=∂US,t∂t−σ22S2∂2US,t∂S2−rS∂US,t∂S+rUS,t,(32)VS,T=maxK−S,0,V0,t=Kexp−rt,limS⟶∞VS,t=0.
With the maturity of 1, in years, the achieved result is shown in Figure9. The exact solution of this put option can be obtained using the following analytical solution denoted by UexactS,t:(33)d1=lnS/K+r+1/2σ2T−tσT−t,d2=lnS/K+r−1/2σ2T−tσT−t,NS=12π∫−∞Sexp−12y2dy,UexactS,t=−SN−d1+Kexp−rT−tN−d2.
Since the problem domain is unbounded according to the boundary conditions, when truncating the domain at 15, the cardinal of the training dataset is 110, containing the equidistant point scattered between0,15 and η=0.2. The number of time steps, N, is 10, and the average calculation time per each epoch on the abovementioned configuration is 0.032s. Increasing the number of data points will increase the accuracy for this configuration slightly (other parameters might need to be adjusted as well); however, we preferred this size to reduce complexity and memory usage.
In Figure9(a), it is shown that when the problem domain is truncated, unlike Example 3, the approximate price is behaving fine throughout the whole unbounded interval. But this does not state that this behavior is the expected behavior of the option considering that the boundary condition makes the option price move towards zero. To make it clearer, the errors for truncated and mapped approximate solutions are compared in Figure 10. In Figure 10(a), the absolute error before the truncation point is relatively better than Figure 10(b). However, after the truncation point, it starts to increase and then again flattens out which is predictable according to the boundary condition of this specific function. In other words, this behavior cannot be generalized to other options as well because farther points are outside the training dataset and the network cannot learn their values. So, the preferred way is employing a mapped domain with lower accuracy but more stable behavior. Besides, only 10 equidistant points are used as training points in this case.
Figure11 shows the superiority of the Adam optimizer as it starts to converge and moves towards the answer in earlier epochs for the first time step. The average calculation time for SGD and RMSprop are, respectively, 0.038s and 0.045s per epoch. Figures 12(a) and 12(b) illustrate the influence of fine-tuning on the objective convergence.Figure 9
Plots of the approximated solutions of Example3 when the number of hidden neurons is 60 and dt=1/20. (a) US,0 (the red dotted plot is the approximate solution). (b) US,t, t∈0,T.
(a)(b)Figure 10
Plots of the convergence and logarithmic absolute error for the solution of Example3 when the number of hidden neurons is 60 and dt=1/30. (a) Cost function in each ANN iteration. (b) UexactS,T−UappS,0.
(a)(b)Figure 11
Comparison of networks’ convergence toward the exact solution of Example3 using (a) Adam optimizer, (b) SGD optimizer, and (c) RMSprop.
(a)(b)(c)Figure 12
The effect of fine-tuning on the convergence rate of the method. Solving Example3 with (a) fine-tuning and (b) without fine-tuning.
(a)(b)Remark for future work: this paper did not cover the application of the method on space fractional or time-space fractional equations [61]. As the sigmoid functions are not fractional, the end result will lack the proper accuracy; therefore, readers are encouraged to use other basis functions such as fractional Chebyshev functions to approximate space fractional models.


## 5. Conclusion

This study investigates neural networks with the famous Adam optimizer for solving financial Black–Scholes equations. Converting the PDE into a series of time-dependent ODEs using the backward Euler finite difference method and then solving each of these equations using the proposed model confirm the satisfactory result and fast calculation of the method. The speed of the method is caused by the parallel computations in the neural network for each independent neuron, the straightforward calculations of sigmoid activation functions that do not add to the complexity of the model, and also the small number of training points and hidden neurons for achieving very promising accuracy. The neural network outperforms other methodologies regarding the consistency and accuracy of the model in confrontation with machines or calculation mistakes because of its fault tolerance. Fine-tuning plays a significant role in this method by reducing the building, validation, and calculation time. It also helps the method converge faster by finding the appropriate direction for gradients as depicted in three examples inSection 4. Domain mapping, which has not been used in ANNs before to the best of found knowledge, is employed to make calculations possible on bigger or infinite problem domains. As a result of combining these approaches into one single ANN, reliable, fast, and accurate results were calculated. The methodology is applicable to other types of options priced by either ordinary or fractional models as well as other partial differential equations in any other field of study that can be solved using this network.



---

*Source: 5511396-2021-07-30.xml*