Numerical Solution of Time Fractional Black–Scholes Model Based on Legendre Wavelet Neural Network with Extreme Learning Machine

: In this paper, the Legendre wavelet neural network with extreme learning machine is proposed for the numerical solution of the time fractional Black–Scholes model. In this way, the operational matrix of the fractional derivative based on the two-dimensional Legendre wavelet is derived and employed to solve the European options pricing problem. This scheme converts this problem into the calculation of a set of algebraic equations. The Legendre wavelet neural network is constructed; meanwhile, the extreme learning machine algorithm is adopted to speed up the learning rate and avoid the over-ﬁtting problem. In order to evaluate the performance of this scheme, a comparative study with the implicit differential method is constructed to validate its feasibility and effectiveness. Experimental results illustrate that this scheme offers a satisfactory numerical solution compared to the benchmark method. data


Introduction
The options pricing model plays an important role in modern financial theory, which has been employed in cost budgets [1][2][3], asset valuations [4], resource allocation [5,6], and so on. It is of great significance to study options pricing and its modeling problem. The mathematical study of the financial market dates back to 1900 [7]. In 1973, the classic Black-Scholes model was proposed by Black and Scholes [8,9], which serves as a milestone of fair market options pricing. Over several decades, the B-S model has been extensively accepted because it can effectively model the options price and provides a mechanism for extracting implied volatility. Although the options price obtained from the B-S model can be well approximated to the observed one, it still has some deficiencies such as the failure to capture jump-diffusion over a small time frame in the financial market [10].
Being a generalization and extension of classic calculus, fractional calculus dates back to a letter from Leibniz to l'Hôpital in 1665 [11]. However, the theory and foundation of fractional calculus were developed by Liouville in 1832. After several decades, Grunwald and Letnikov defined the fractional derivative with differential quotient approach based on a fractional binomial expansion proposed by Newton. In 1967, Caputo defined Caputo's fractional derivative, which requires the function to be differentiable [12]. In 2006, Jumarie made some improvements to the Riemann-Liouville derivative by subtracting a constant from the integrand to avoid the derivative of the constant being nonzero [13]. Guariglia and Silvestrov introduced the Ortigueira-Caputo fractional derivative operator to describe a complex function in the distribution sense [14]. Indeed, since some remarkable work appearing in [15], fractional calculus has become the key invariant, providing rich information about the level of complexity that a certain system presents. With plenty of research focusing on the fractional calculus field, it becomes imperative and necessary to solve the can be applied in image analysis and computer vision, such as video segmentation [56], speech recognition, and data compression [57]. Besides, some improved neural-networkbased methods have recently emerged; however, these methods lack a strong theoretical foundation and require massive amounts of training data. Shi et al. proposed the deep scattering network (DSN), a variant of the deep convolutional neural network (DCNN), where the data-driven linear filters are replaced by predefined fixed multi-scale wavelet filters; for the non-stationary image textures, the fractional wavelet transform (FRWT) is employed because it can be regarded as a bank of linear translation-variant multi-scale filters [58]. Aiming at poor resolution in the high-fractional-frequency domain, Shi and his partners revealed a novel fractional wavelet package transform (FRWPT) and its recursive algorithm [59].
Moreover, previous literature demonstrated the reliability and stability of the WNN in the learning procedure [60,61]. Due to the decomposition capacity of the WNN in the time and frequency domain, this paper presents a novel Legendre wavelet neural network (LWNN) with extreme learning machine (ELM) to solve the time fractional B-S model governed by the European option. Taking the final and boundary conditions into consideration, this time fractional B-S model can be converted into a group of algebraic equations. In the end, numerical experiments illustrate the superiority of the proposed method compared to the implicit differential method.
The objective of this paper is to solve a time fractional B-S model of the European option, where ELM is used for the output layer weight learning of the LWNN. The advantages of the proposed scheme are given as follows: • It is a single hidden layer feed-forward network; we only should train the weights of the output layer; the input layer weights can be randomly selected. • ELM can be classified as unsupervised learning; the optimization scheme is unnecessary; besides, it has a fast learning rate and can avoid over-fitting. • It converts the European options pricing problem to a group of algebraic equations, which can substantially facilitate analysis.
The structure of this paper is as follows. Section 2 gives a description of the time fractional B-S model governed by the European option with double barriers. Section 3 introduces the properties of the Legendre wavelet. Section 4 deduces the operational matrix of the fractional derivative for the two-dimensional Legendre wavelet. Section 5 presents the topology of the LWNN and the algorithm of ELM. Section 6 validates the feasibility and effectiveness of the proposed method by numerical experiments. In the end, some conclusions and directions for improvement are presented.

Time Fractional B-S Model
Liang [62] proposed a bi-fractional B-S model in the time and spatial domain, where the underlying asset price follows a fractional Ito process and options price variation is a fractal transmission system. As a result, the time fractional B-S model can be regarded as a special case of the bi-fractional B-S model. Let S, t, r, σ, and q be the underlying asset price, time, risk-free interest rate, volatility, and dividend payment, respectively, then options price V(S, t) should satisfy ∂V ∂t According to the arguments in [62][63][64], the variation of the options price regarding time can be assumed as a fractal transmission system, which implies that the average options price C(S, t) from time t to expiration time T should satisfy where F(t) and d f represent the transmission function and Hausdorff dimension of the fractional transmission system, respectively. As pointed out in [62], Equation (2) is a conservation equation containing an explicit reference to the diffusion process of the options price based on a fractal structure. It can be assumed that the diffusion sets are the underlying fractal, where F(t) = A α t −α /Γ(1 − α) is selected as the transmission function with A α and α a constant and the transmission exponent, respectively. Making partial differentiation regarding t on both sides of Equation (2) yields On the other hand, because of the B-S model, we have substituting Equation (3) into (4) yields when α approaches 1, this yields It is clear that the time fractional B-S model is consistent with the classic B-S model when A α = d f = 1. In this paper, A α = d f = 1 are selected. As a matter of fact, the solution procedure is applicable to other values of A α and d f .
It should be mentioned that the time fractional order model is a simplification compared with that proposed by Liang [62]. In this paper, we assume that the underlying asset follows standard Brownian motion, only considering the options price variation with time as a fractal transmission system; this is why the fractional derivative only appears in the time domain. Actually, Jumaire [65,66] defined the stock exchange fractional dynamics as fractional exponential growth driven by Gaussian noise; it can be found that the derived equation is similar to Equation (5). However, the derived equation has a time variable coefficient in front of the second derivative, which is different from the time fractional B-S model mentioned above.

Double Barrier Option
The barrier option is probably the oldest exotic option, which was sporadically traded on the U.S. market before the establishment of the Chicago Board of Options Exchange. The barrier option is essentially a conditional option, depending on whether the barrier can be triggered or breached within its expiration; meanwhile, it is considerably less expensive than the corresponding vanilla option. Obviously, a double barrier option features two barriers, and the payoff to the holder depends on the breaching behaviors of the underlying asset process with respect to these two barriers. Let V(S, t) be the price of a European option with double barriers; obviously, if the price variation regarding time is considered to be a fractal transmission system, then V(S, t) should satisfy where P(t) and Q(t) denote the discount paid when the corresponding barrier is breached and K and T are the strike price and expiration time, respectively. In general, it is complicated to obtain the solution of this partial derivative equation; the following variables and functions are introduced to non-dimensionalize this model: , and Q(t) = KQ(τ); as a consequence, the European Equation (7) can be rewritten as is consistent with the modified Riemann-Liouville derivative in [13]. Substituting Equation (8) into Equation (7) yields

Legendre Wavelet and Its Properties
Legendre polynomials are constructed by the polynomials', {1, x, · · · , x n , · · · }, orthogonalization with respect to weighting function w(x) = 1. In 1814, Rodrigul presented the simple expression as In [67], the recurrence formula of Legendre polynomials was given as For the practical implementation of Legendre polynomials on domain [l d , l u ], it is feasible to shift the domain by means of linear transformation x = [2x − (l u + l d )]/(l u − l d ). Hence, the shifted Legendre polynomials L * n (x) can be obtained as L * n (x) = L n ( x). Specifically, the orthogonality characteristic of shifted Legendre polynomials on domain [0, 1] is Moreover, the analytical form of shifted Legendre polynomials L * n (x) can be expressed as [67] Furthermore, the wavelet function is generally constructed by the dilation and translation of the mother wavelet. When the dilation and translation parameters vary continuously, we can obtain the continuous wavelet: where a and b are the dilation and translation parameters, respectively. When they are constricted to only discrete values a = a k where ψ k,n (x) represents the basis of L 2 (R). In particular, ψ k,n (x) constructs an orthogonal basis when a 0 = 2 and b 0 = 1. Legendre wavelet ψ n,m (t) = ψ(k, n, m, t) is composed of four arguments, n = 2n + 1, n = 0, 1, 2, · · · , 2 k−1 , m = 0, 1, · · · , m 0 − 1 is the degree of the Legendre polynomials, and t is normalized time [68]. The Legendre wavelet on domain [0, 1] is defined as

Theorem 1. [[69]]
Supposing that H and Y ⊂ H are the inner product and complete space, respectively, we define {e 0 , e 1 , · · · , e n } as an orthogonal basis for H. For f ∈ H, the best approximation in Y: where (, ) denotes the inner product, that is f 2 = ( f , f ).  (17) where c n,m = ( f (x), ψ n,m (x)). If Equation (17) is substituted with finite series expansion, then f (x) can be written as where T represents the matrix transposition operator and C and Ψ(x) are the 2 k m 0 -column vector shown as below. (18) can be rewritten as where Similarly, an arbitrary function with two variables f (x, y) ∈ L 2 (R × R) defined on domain [0, 1] × [0, 1] can be approximated with the Legendre wavelet: where
Here, for a better understanding of the calculation procedure, B (1.5) and F (1.5) are provided as below, where m 0 = 4, k = 2, and α = 1.5 are selected.
In order to solve the time fractional B-S model of the European option, U(x, τ) in Equation (9) can be approximated with where U = [u i,j ] M×M is an unknown matrix, which needs to be determined; from Equation (30), we can obtain From Definition 1, Equation (32) can be rewritten as and I M is an M-dimensional identity matrix. Then, Equation (9) can be converted to Define and let ∆x = B u − B d N s , x i = B d + i∆x, i = 0, 1, 2, · · · , N s ∆τ = T/N t , τ j = j∆τ, j = 0, 1, 2, · · · , N t be the step size of stock price and time, respectively. Denoting which can be expressed as Let

ELM Algorithm for LWNN Training
ELM has a fast convergence rate owing to its single hidden layer feed-forward structure. It has many advantages over traditional BP methods, because traditional BP methods can easily become stuck in local optima and have a slow convergence rate owing to the continuous weight updating. The parameters of ELM only depend on the analytical resolution of the output layer weights and the random selection of input and hidden layer parameters (weights and biases); as a consequence, many issues, such as the learning rate, local optima and learning epochs, can be avoided. It can solve the slow convergence rate and over-fitting issues of gradient descent methods based on neural networks. Here, the ELM algorithm is illustrated in Algorithm 1. The Neural network structure of LWNN-ELM is showed in Figure 1.

Input:
a i and b i denote the weights and bias of the input layer, respectively. N i , N h , N o represent the nodes of the input, hidden, and output layer, respectively.
Step 1: Attribute random parameters of the hidden layer nodes, weights, and biases.
Step 2: Calculate the output matrix of the hidden layer nodes H.
Step 3: Calculate the output weight where H † is the Moore-Penrose generalized inverse of the hidden layer output matrix H, Y is the training data target. If matrix H is square and invertible, then β = H −1 Y.

2.
If matrix H is rectangle, then β = H † Y; β is the minimal least square solution; in other words, β = arg min Yβ − T .

3.
If matrix H is singular, then β = H † Y, where H † = H T λI + HH T −1 ; λ is a regularization parameter, which can be set based on a specific norm.
The procedure for training the network weights of the LWNN to solve the European options pricing model is given in Algorithm 2.

Remark 2.
The detailed proof of this theorem can be referred to the related knowledge of the generalized inverse matrix [43].

Numerical Experiment and Discussion
Once the numerical solution is obtained for a specific options pricing model, the main concern of market practitioners becomes its implementation. Whether the options price can be efficiently computed is one of the core criteria to evaluate its characteristics. Therefore, in this section, as far as the validation of the numerical solution's accuracy and the convergence order of the proposed method are concerned, a double barrier knock-out call European option is selected as an experiment. For other kinds of options pricing models, their solutions can be straightforwardly obtained by a simple modification or the parity relationship [75].
In order to validate the feasibility of the proposed method, calculating the solution for α = 1 can be considered as the best way, and comparing it with a benchmark solution based on the implicit difference method (IDM) [22], where all numerical experiments were implemented with MATLAB R2013b on a desktop with i7-6700K CPU, 8GB memory, 1TB HDD, and a Windows 7 operating system. From Figure 2, it can be seen that the two groups of options price coincide perfectly, which can affirm the correctness of the proposed method. Besides, the calculation times for different expirations is presented in Table 1; it is clear that the time consumption of the LWNN-ELM method decreases significantly compared with the IDM, which can save computing resources.  Here, the influence of different α on the options price is investigated. From Equation (6), intuitively, α can be regarded as a metric of the dependence of U(S, t) on the stock price at the same underlying level from the current time to the expiration; in the practical case, T − t 1; it is not difficult to observe that the larger α is, the weaker the dependence becomes. This may explain the phenomena illustrated in Figure 3, wherein the B-S model tends to overprice a double knock-out call when the underlying is close to the lower barrier. On the other hand, it can be observed that from a certain underlying value and onward, this B-S model prefers to underprice the option, which possibly results in the larger impact of the nonzero payoff value on the options price for smaller α. Table 2 offers the calculation time of the LWNN-ELM for different expirations; it can be seen that for a longer expiration, the time consumption has no apparent increase which is beneficial for some complex options pricing models.   For the European put option, the final and boundary conditions are Π(S) = max{K − S, 0}, Q(t) = q = 0, and P(t) = K exp(−r(T − t)). The put options prices for different M are presented in Figure 4; it is clear that the numerical solutions of the LWNN-ELM for different M are quite close; meanwhile, the time consumption is 4.8391 (M = 12) and 6.1726 (M = 16), respectively. In addition, the numerical solutions based on the LWNN-ELM and IDM are provided in Table 3; it is clear that the numerical solutions obtained by the LWNN-ELM agree very well with the benchmark method. Although the numerical solution for M = 12 has a faster computation rate compared to that of M = 16, the numerical solution for M = 16 has a higher precision. Consequently, for a more accurate description of financial dynamics, a fast and efficient scheme to solve the options price is preferred. In the end, the order of convergence is utilized to explore the solution accuracy [76], which is defined as where e (∆τ) = U (∆τ) − U (∆τ/2) and U (∆τ) is the numerical solution at the final time with step size ∆τ. In the following experiments, N s = 460 was selected to saturate the stock price discretization so that the error is mainly determined by the time discretization. The order of convergence of the European option is illustrated in Table 4. It can be observed that the LWNN-ELM has first-order convergence for different α.

Conclusions
Being a generalization of the classic B-S model, the "global" characteristic of the time fractional B-S model makes it more difficult to calculate the analytical or numerical solution than the integer-order model. In this paper, the options price model is governed by a time fractional derivative B-S equation when the underlying price change is considered as a fractal transmission system, then the fractional derivative operational matrix of the two-dimensional Legendre wavelet is derived, and the European options pricing model (including final and boundary conditions) is transformed to a group of algebraic equations. The LWNN-ELM is employed to train the output layer weights owing to its fast learning rate and moderate fitting characteristic. In the end, numerical experiments are provided to price the European option with double barriers based on the LWNN-ELM and the implicit differential method, and the experimental results illustrate that the proposed method has less time consumption compared with the benchmark method; besides, the first-order convergence is demonstrated when α is less than 1. Furthermore, the proposed scheme is simple and easy to implement and can be available to other similar fractional models for pricing different European options.