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

Abstract

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-fitting 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.

1. 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 fractional differential equation. An implicit method was proposed by Murio to solve the linear time fractional diffusion equation, constructed by Caputo’s fractional derivative on a finite slab [16]. Langlands offered the stability and accuracy of the implicit numerical method for solving the fractional diffusion equation [17]. Chen proposed an explicit closed-form solution of the double-barrier option based on the time fractional B-S model [18], where the solution is obtained by the convolution of some specific function (Mittag–Leffler function) and finite series expansion with integration. As a consequence, some numerical methods were introduced to simplify the calculation procedure. Song adopted the implicit finite difference method to solve the options pricing problem based on Jumaire’s model [19]. Cen [20] transformed the differential equation to an integral form with a weakly singular kernel and discretized it on an adaptive mesh grid, rigorously analyzed the convergence of the proposed scheme, and took account of possible singular behavior. Rezaei [21] presented the time fractional B-S model under the European double barrier option with the constant of elasticity of variance and presented the rigorous solvability, stability, and convergence of the implicit difference method. Zhang deduced an implicit discrete solution to the time fractional B-S model for the European option with double barriers, which has a spatially second-order accuracy and temporally $2-\alpha$-order accuracy [22]. De Staelen improved the spatial accuracy of the implicit differential method to fourth-order and presented the solvability, stability, and convergence of this method based on the Fourier method [23]. Moreover, some analytical methods have been derived, such as homotopy perturbation and the Laplace transform [24], the homotopy analysis method [25,26,27,28,29,30], the fractional variational iteration method [31,32], the generalized differential transform method [33], and the Adomain decomposition method [34,35,36,37]. The essential aspects of these methods are the convolution of some special function or the integration of infinite series expansion, which are complicated to calculate; what is more, the solvability, stability, and convergence need to be proven. As a consequence, making some improvements to the numerical method becomes practical and necessary.

Compared with standard Brownian motion, fractional Brownian motion can describe the dynamic characteristics of financial markets accurately; meanwhile, some new fractional B-S formulae have been derived. The fractional B-S model for the European option was initially proposed by Wyss [38], then Hu and Oksendal derived a fractional B-S model based on fractional Brownian motion [39]. Li developed a fractional stochastic differential equation to demonstrate the existence of trend memory in options pricing when the Hurst index is between 0.5 and 1 [40]. Jumarie proposed two new families of the fractional B-S model based on fractional Taylor expansion series and the modified Riemann–Liouville derivative [41]. Liu and Chang also proposed a time fractional B-S model with transaction cost based on fractional Brownian motion and provided an approximate solution of the nonlinear Hoggard–Whalley–Wilmott equation [42].

Accompanying the evolution of computer technology, machine learning has gained enormous advancement. Being a branch of artificial intelligence, reliable decisions and results can be obtained through iteration and optimization. In the last decade, a single feed-forward neural network named extreme learning machine [43,44] was proposed by Huang, which is different from traditional learning methods based on the gradient descent principle; meanwhile, the single hidden layer structure makes its learning procedure exceptionally faster compared with that of a deep neural network [45,46]; however, the ELM algorithm has been demonstrated to be sensitive to the activation function [47,48,49]. Furthermore, the wavelet can be regarded as a mathematical microscope for signal analysis (constant ratio of bandwidth and central frequency); its multi-resolution analysis and time–frequency ability can provide a powerful tool for signal analysis and approximation, then different wavelets can be constructed by dilation and translation, such as the Legendre wavelet [50], Chebyshev wavelet [51], Taylor wavelet [52], etc. Aiming at the framework of an adaptive multi-scale wavelet, Zheng pointed out that a sampling rate far from 1/2 can lead to low representation efficiency and proposed a signal Shannon-entropy-based adaptive multi-scale wavelet decomposition for signals analysis [53]. The wavelet neural network, suggested by Zhang in 1992 [54], was initially used for signal approximation [55], which can be applied in image analysis and computer vision, such as video segmentation [56], speech recognition, and data compression [57]. Besides, some improved neural-network-based 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.

2. Double Barrier Option under Time Fractional B-S Model

2.1. 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, $\sigma$, 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

$$\begin{array}{c}\hfill \frac{\partial V}{\partial t}+\frac{1}{2}{\sigma }^2{S}^2\frac{{\partial }^{2}V}{\partial S^2}+(r-q)S\frac{\partial V}{\partial S}-rV=0\end{array}$$

(1)

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

$$\begin{array}{c}\hfill {\int }_t^{T}C(S,\tau )d\tau =S^{d_f-1}{\int }_t^{T}F\left(\tau -t\right)\left[V(S,\tau )-V(S,T)\right]d\tau \end{array}$$

(2)

where $F\left(t\right)$ 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\left(t\right)=A_{\alpha }{t}^{-\alpha }/\Gamma \left(1-\alpha \right)$ is selected as the transmission function with $A_{\alpha }$ and $\alpha$ a constant and the transmission exponent, respectively. Making partial differentiation regarding t on both sides of Equation (2) yields

$$\begin{array}{c}\hfill -C(S,t)=S^{d_f-1}\frac{\partial }{\partial t}{\int }_t^{T}F(\tau -t)\left[V(S,\tau )-V(S,T)\right]d\tau \end{array}$$

(3)

On the other hand, because of the B-S model, we have

$$\begin{array}{c}\hfill C(S,t)=\frac{1}{2}{\sigma }^2{S}^2\frac{{\partial }^{2}V}{\partial S^2}+(r-q)S\frac{\partial V}{\partial S}-rV\end{array}$$

(4)

substituting Equation (3) into (4) yields

$$A_{\alpha }{S}^{d_f-1}\frac{{\partial }^{\alpha }V}{\partial t^{\alpha }}+\frac{1}{2}{\sigma }^2{S}^2\frac{{\partial }^{2}V}{\partial S^2}+(r-q)S\frac{\partial V}{\partial S}-rV=0$$

(5)

where

$$\begin{array}{c}\hfill \frac{{\partial }^{\alpha }V}{\partial t^{\alpha }}=\frac{1}{\Gamma (n-\alpha )}\frac{{\partial }^n}{\partial t^n}{\int }_t^T\frac{V(S,\tau )-V(S,T)}{{(\tau -t)}^{\alpha +1-n}}d\tau ,n-1⩽\alpha <n\end{array}$$

(6)

when $\alpha$ approaches 1, this yields

$$\begin{array}{c}\hfill \underset{\alpha \to 1}{lim}\frac{{\partial }^{\alpha }V}{\partial t^{\alpha }}=\frac{1}{\Gamma (2-1)}\frac{\partial }{\partial t}{\int }_t^T\frac{V(S,\tau )-V(S,T)}{{(\tau -t)}^{1+1-2}}d\tau =\frac{\partial V}{\partial t}\end{array}$$

It is clear that the time fractional B-S model is consistent with the classic B-S model when $A_{\alpha }=d_f=1$. In this paper, $A_{\alpha }=d_f=1$ are selected. As a matter of fact, the solution procedure is applicable to other values of $A_{\alpha }$ 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.

2.2. 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

$$\left\{\begin{array}{c}\frac{{\partial }^{\alpha }V}{\partial t^{\alpha }}+\frac{1}{2}{\sigma }^2{S}^2\frac{{\partial }^{2}V}{\partial S^2}+(r-q)S\frac{\partial V}{\partial S}-rV=0\hfill \\ V(S,T)=\overline{\Pi }\left(S\right)=max(S-k,0)\hfill \\ V(A,t)=\overline{P}\left(t\right)\hfill \\ V(B,t)=\overline{Q}\left(t\right)\hfill \end{array}\right.$$

(7)

where $\overline{P}\left(t\right)$ and $\overline{Q}\left(t\right)$ 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: $t=T-2\tau /\left({\sigma }^2\right)$, $S=K{e}^x$, $V(S,t)=KU(x,\tau )$, $\overline{P}\left(t\right)=KP\left(\tau \right)$, and $\overline{Q}\left(t\right)=KQ\left(\tau \right)$; as a consequence, the European Equation (7) can be rewritten as

$$\begin{array}{cc}\hfill \frac{{\partial }^{\alpha }V}{\partial t^{\alpha }}& =\frac{1}{\Gamma (1-\alpha )}\frac{\partial }{\partial t}{\int }_t^T\frac{V(S,\overline{t})-V(S,T)}{{(\overline{t}-t)}^{\alpha }}d\overline{t}\hfill \\ & =\frac{1}{\Gamma (1-\alpha )}\left(-\frac{{\sigma }^2}{2}\right)\frac{\partial }{\partial \tau }{\int }_{T-\frac{2\tau }{{\sigma }^2}}^T\frac{V(S,\overline{t})-V(S,T)}{{\left[\overline{t}-\left(T-\frac{2\tau }{{\sigma }^2}\right)\right]}^{\alpha }}d\overline{t}\hfill \\ & =\frac{1}{\Gamma (1-\alpha )}\left(-\frac{{\sigma }^2}{2}\right)\frac{\partial }{\partial \tau }{\int }_{\tau }^0\frac{V\left(S,T-\frac{2\overline{\tau }}{{\sigma }^2}\right)-V(S,T)}{{\left[\frac{2}{{\sigma }^2}\left(\tau -\overline{\tau }\right)\right]}^{\alpha }}\left(-\frac{2}{{\sigma }^2}\right)d\overline{\tau }\hfill \\ & =\frac{-K}{\Gamma (1-\alpha )}{\left(\frac{{\sigma }^2}{2}\right)}^{\alpha }\frac{\partial }{\partial \tau }{\int }_0^{\tau }\frac{U(x,\overline{\tau })-U(x,0)}{{\left(\tau -\overline{\tau }\right)}^{\alpha }}d\overline{\tau }=-K{\left(\frac{{\sigma }^2}{2}\right)}^{\alpha }\frac{{\partial }^{\alpha }U}{\partial {\tau }^{\alpha }}\hfill \end{array}$$

(8)

where

$$\begin{array}{c}\hfill \frac{{\partial }^{\alpha }U}{\partial t^{\alpha }}=\frac{1}{\Gamma (1-\alpha )}\frac{\partial }{\partial t}{\int }_0^t\frac{U(x,t)-U(x,0)}{{(\tau -t)}^{\alpha }}d\tau ,0<\alpha <1\end{array}$$

is consistent with the modified Riemann–Liouville derivative in [13]. Substituting Equation (8) into Equation (7) yields

$$\left\{\begin{array}{c}\frac{{\partial }^{\alpha }U}{\partial {\tau }^{\alpha }}=c_1\frac{{\partial }^{2}U}{\partial x^2}+c_2\frac{\partial U}{\partial x}+c_{3}U\hfill \\ U(x,0)=\Pi \left(x\right)\hfill \\ U(B_d,\tau )=P\left(\tau \right)\hfill \\ U(B_u,\tau )=Q\left(\tau \right)\hfill \end{array}\right.$$

(9)

where $\gamma =r{(-{\sigma }^2/2)}^{-\alpha }$, $d=q{({\sigma }^2/2)}^{-\alpha }$, $c_1={({\sigma }^2/2)}^{1-\alpha }$, $c_2=r-d-c_1$, $c_3=-\gamma$, $B_d=ln(A/K)$, and $B_u=ln(B/K)$. It can be observed that $\frac{{\partial }^{\alpha }U}{\partial {\tau }^{\alpha }}$ is consistent with Caputo’s derivative in Definition 2.

3. Legendre Wavelet and Its Properties

Legendre polynomials are constructed by the polynomials’, $\{1,x,\cdots ,x^n,\cdots \}$, orthogonalization with respect to weighting function $w\left(x\right)=1$. In 1814, Rodrigul presented the simple expression as

$$\begin{array}{c}\hfill L_0\left(x\right)=1,L_n\left(x\right)=\frac{1}{2^{n}n!}\frac{d^n}{d{x}^n}\left\{{(x^2-1)}^n\right\},x\in \left[-1,1\right],n\in {\mathbb{Z}}^{+}\end{array}$$

(10)

In [67], the recurrence formula of Legendre polynomials was given as

$$\left\{\begin{array}{c}{L}_0\left(x\right)=1,L_1\left(x\right)=x\hfill \\ L_{n+1}\left(x\right)=\frac{2n+1}{n+1}x{L}_n\left(x\right)-\frac{n}{n+1}{L}_{n-1}\left(x\right)\hfill \end{array}\right.$$

(11)

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 $\tilde{x}=[2x-(l_u+l_d)]/(l_u-l_d)$. Hence, the shifted Legendre polynomials $L_n^{*}\left(x\right)$ can be obtained as $L_n^{*}\left(x\right)=L_n\left(\tilde{x}\right)$. Specifically, the orthogonality characteristic of shifted Legendre polynomials on domain $[0,1]$ is

$${\int }_0^1{L}_m^{*}\left(x\right)L_n^{*}\left(x\right)dx=\left\{\begin{array}{c}\frac{1}{2m+1},m=n\hfill \\ 0,\mathrm{otherwise}\hfill \end{array}\right.$$

(12)

Moreover, the analytical form of shifted Legendre polynomials $L_n^{*}\left(x\right)$ can be expressed as [67]

$$\begin{array}{c}\hfill L_n^{*}\left(x\right)=\sum _{i=0}^n{b}_{n,i}{x}^i\end{array}$$

(13)

where

$$\begin{array}{c}\hfill b_{n,i}={(-1)}^{n+i}\frac{(n+i)!}{(n-i)!{(i!)}^2}\end{array}$$

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:

$$\begin{array}{c}\hfill {\psi }_{a,b}\left(x\right)=\frac{1}{\sqrt{a}}\psi \left(\frac{x-b}{a}\right),a,b\in \mathbb{R},a>0\end{array}$$

(14)

where a and b are the dilation and translation parameters, respectively. When they are constricted to only discrete values $a=a_0^k$, $b=n{a}_0^k{b}_0$, $a_0>1$, $b_0>1$, this yields

$$\begin{array}{c}\hfill {\psi }_{k,n}\left(x\right)=\frac{1}{\sqrt{a_0^k}}\psi \left(\frac{x-n{a}_0^k{b}_0}{a_0^k}\right)=a_0^{-k/2}\psi \left(a_0^{-k}x-n{b}_0\right),n,k\in {\mathbb{Z}}^{+}\end{array}$$

(15)

where ${\psi }_{k,n}\left(x\right)$ represents the basis of $L^2\left(R\right)$. In particular, ${\psi }_{k,n}\left(x\right)$ constructs an orthogonal basis when $a_0=2$ and $b_0=1$. Legendre wavelet ${\psi }_{n,m}\left(t\right)=\psi (k,\widehat{n},m,t)$ is composed of four arguments, $\widehat{n}=2n+1$, $n=0,1,2,\cdots ,2^{k-1}$, $m=0,1,\cdots ,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

$${\psi }_{n,m}\left(x\right)=\left\{\begin{array}{c}\sqrt{2m+1}\sqrt{2^k}{L}_m\left(2^{k+1}x-\widehat{n}\right),\frac{n}{2^k}⩽x⩽\frac{n+1}{2^k}\hfill \\ 0,\mathrm{otherwise}\hfill \end{array}\right.$$

(16)

Theorem 1.

([69]) Supposing that $\mathcal{H}$ and $\mathcal{Y}\subset \mathcal{H}$ are the inner product and complete space, respectively, we define $\{e_0,e_1,\cdots ,e_n\}$ as an orthogonal basis for $\mathcal{H}$. For $f\in \mathcal{H}$, the best approximation in $\mathcal{Y}$:

$$\begin{array}{c}\hfill f_0=\sum _{i=0}^n\left(f,e_i\right)e_i\end{array}$$

such that

$$\begin{array}{c}\hfill \forall y\in \mathcal{Y},{∥f-f_0∥}_2⩽{∥f-y∥}_2\end{array}$$

where $(,)$ denotes the inner product, that is ${∥f∥}_2=\sqrt{(f,f)}$.

From Theorem 1, a function $f\left(x\right)$ on domain $[0,1]$ can be expanded by the Legendre wavelet as

$$\begin{array}{c}\hfill f\left(x\right)=\sum _{n=0}^{\infty }\sum _{m=0}^{\infty }{c}_{n,m}{\psi }_{n,m}\left(x\right)\end{array}$$

(17)

where $c_{n,m}=\left(f\left(x\right),{\psi }_{n,m}\left(x\right)\right)$. If Equation (17) is substituted with finite series expansion, then $f\left(x\right)$ can be written as

$$\begin{array}{c}\hfill f\left(x\right)\simeq \sum _{n=0}^{2^k-1}\sum _{m=0}^{m_0-1}{c}_{n,m}{\psi }_{n,m}\left(x\right)=C^T\Psi \left(x\right)\end{array}$$

(18)

where T represents the matrix transposition operator and C and $\Psi \left(x\right)$ are the $2^k{m}_0$-column vector shown as below.

$$\begin{array}{cc}\hfill C=& {\left[c_{0,0},c_{0,1},\cdots ,c_{0,m_0-1},c_{1,0},c_{1,1},\cdots ,c_{1,m_0-1},\cdots ,c_{2^k-1,0},c_{2^k-1,1},\cdots ,c_{2^k-1,m_0-1}\right]}^T\hfill \\ \hfill \Psi \left(x\right)=& \left[{\psi }_{0,0}\left(x\right),{\psi }_{0,1}\left(x\right),\cdots ,{\psi }_{0,m_0-1}\left(x\right),{\psi }_{1,0}\left(x\right),{\psi }_{1,1}\left(x\right),\cdots ,{\psi }_{1,m_0-1}\left(x\right),\cdots ,{\psi }_{2^{k-1},0}\left(x\right),\right.\hfill \\ & {\left.{\psi }_{2^k-1,1}\left(x\right),\cdots ,{\psi }_{2^k-1,m_0-1}\left(x\right)\right]}^T\hfill \end{array}$$

For simplicity, let $c_i=c_{n,m}$, ${\psi }_i={\psi }_{n,m}$, and $M=2^k{m}_0$, where index i is determined by $i=m_{0}n+m+1$, then Equation (18) can be rewritten as

$$\begin{array}{c}\hfill f\left(x\right)\simeq \sum _{i=1}^M{c}_i{\psi }_i\left(x\right)=C^T\Psi \left(x\right)\end{array}$$

(19)

where

$$\begin{array}{cc}\hfill C=& {\left[c_1,\cdots ,c_{m_0},|,c_{m_0+1},\cdots ,c_{2m_0},|,c_{m_0(2^k-1)+1},\cdots ,c_M\right]}^T\hfill \\ \hfill \Psi \left(x\right)=& \left[{\psi }_1\left(x\right),\cdots ,{\psi }_{m_0}\left(x\right),|,{\psi }_{m_0+1}\left(x\right),\cdots ,{\psi }_{2m_0}\left(x\right),|,\cdots ,|,\right.\hfill \\ & {\left.{\psi }_{m_0(2^k-1)+1}\left(x\right),\cdots ,{\psi }_M\left(x\right)\right]}^T\hfill \end{array}$$

Similarly, an arbitrary function with two variables $f(x,y)\in {\mathbb{L}}^2(\mathbb{R}\times \mathbb{R})$ defined on domain $[0,1]\times [0,1]$ can be approximated with the Legendre wavelet:

$$\begin{array}{c}\hfill f(x,y)\simeq \sum _{i=1}^M\sum _{j=1}^M{\psi }_i\left(x\right)f_{ij}{\psi }_j\left(y\right)={\Psi }^T\left(x\right)F\Psi \left(y\right)\end{array}$$

(20)

where $F={\left[f_{ij}\right]}_{M\times M}$ and $f_{ij}=\left({\psi }_i\left(x\right),\left(f(x,y),{\psi }_j\left(y\right)\right)\right)$.

Theorem 2.

(Convergence analysis [70]) Supposing that a continuous function $f\left(x\right)$ on domain $[0,1]$ is bounded with second derivative $\left|f^{\prime \prime }\left(x\right)\right|⩽K$, it can be uniformly approximated by the infinite sum of Legendre wavelet expansion, that is $f\left(x\right)={\sum }_{n=1}^{\infty }{\sum }_{m=0}^{\infty }{c}_{n,m}{\psi }_{n,m}\left(x\right)$; moreover, coefficients $c_{n,m}$ are bounded with $\left|c_{n,m}\right|⩽K\sqrt{12}/{\left(2n\right)}^{5/2}{(2m-3)}^2$.

Theorem 3.

(Error analysis [71]) Supposing that a continuous function $f(x,y)$ on domain $[0,1]\times [0,1]$ is bounded with four mixed partial derivatives $\left|{\partial }^{4}f(x,y)/\partial x^2\partial y^2\right|⩽D$, then the Legendre wavelet expansion of $f(x,y)$ has a uniform convergence characteristic and

$$\begin{array}{c}\hfill \left|f_{ij}\right|<\frac{12D}{{\left(2n\right)}^5{(2m-3)}^4}\end{array}$$

Definition 1.

Supposing that $A={\left[a_{ij}\right]}_{m\times n}$ and $B={\left[b_{ij}\right]}_{p\times q}$ are two arbitrary matrices, the Kronecker product of A and B is defined as follows [72].

$$\begin{array}{c}\hfill A\otimes B={\left[\begin{array}{cccc}{a}_{11}B& a_{12}B& \cdots & a_{1n}B\\ a_{21}B& a_{22}B& \cdots & a_{2n}B\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}B& a_{m2}B& \cdots & a_{mn}B\end{array}\right]}_{mp\times nq}\end{array}$$

For matrix A, operator Vec $\left(A\right)$ is a column vector made of the column of A stacked on top of each other from left to right [73].

$$\begin{array}{c}\hfill \mathit{Vec}\left(A\right)={\left[a_{11},\cdots ,a_{m1},a_{12},\cdots ,a_{m2},\cdots ,a_{1n},\cdots ,a_{mn}\right]}^T\end{array}$$

Definition 2.

The Caputo’s fractional derivative of order α is defined as [74]

$${}_0^C{D}_x^{\alpha }f(x,y)=\left\{\begin{array}{c}\frac{1}{\Gamma (n-\alpha )}{\int }_0^x\frac{f^{\left(n\right)}(\chi ,y)}{{\left(x-\chi \right)}^{1+\alpha -n}}d\chi ,0⩽n-1<\alpha <n\hfill \\ \frac{{\partial }^{\alpha }f(x,y)}{\partial x^{\alpha }},\alpha =n\hfill \end{array}\right.$$

where operator $D_{*}^{\alpha }$ represents Caputo’s derivative. It can be observed that $D_x^{\alpha }C=0$ (C is a constant); additionally,

$$D_x^{\alpha }{x}^n=\left\{\begin{array}{c}0,n<\lceil \alpha \rceil ,n\in {\mathbb{Z}}^{+}\hfill \\ \frac{\Gamma (n+1)}{\Gamma (n+1-\alpha )}{x}^{n-\alpha },n⩾\lceil \alpha \rceil ,n\in {\mathbb{Z}}^{+}\hfill \end{array}\right.$$

where $\Gamma (·)$ denotes the Gamma function; the ceiling function $\lceil \alpha \rceil$ denotes the smallest integer greater than or equal to α.

4. Operational Matrix of Fractional Derivative

From Equations (13) and (16), for $m<\lceil \alpha \rceil$, the analytical form of the shifted Legendre wavelet can be written as

$$\begin{array}{c}\hfill {\psi }_{n,m}\left(x\right)=2^{k/2}\sqrt{2m+1}\sum _{i=0}^m{b}_{m,i}{\left(2^{k}x-n\right)}^i{I}_{n,k}\left(x\right)\end{array}$$

(21)

where $I_{n,k}=\left[n/2^k,(n+1)/2^k\right]$ ($n=0,1,\cdots ,2^k-1$, $m=0,1,\cdots ,m_0-1$); in addition, the indicator function is defined as

$$I_{n,k}\left(x\right)=\left\{\begin{array}{c}1,x\in \left[\frac{n}{2^k},\frac{n+1}{2^k}\right]\hfill \\ 0,\mathrm{otherwise}\hfill \end{array}\right.$$

It is apparent that $I_{n,k}\left(x\right)$ is zero outside the interval $\left[n/2^k,(n+1)/2^k\right]$. Applying $D_x^{\alpha }$ on both sides of Equation (21) yields

$$\begin{array}{cc}\hfill D_x^{\alpha }{\psi }_{n,m}\left(x\right)=& \sum _{i=0}^m{D}^{\alpha }\left({\left(x-\frac{n}{2^k}\right)}^i\right)=\sum _{i=\lceil \alpha \rceil }^m\frac{a_{m,i}\left(i\right)!}{\Gamma (i-\alpha +1)}{\left(x-\frac{n}{2^k}\right)}^{i-\alpha }{I}_{n,k}\left(x\right)\hfill \\ & +\sum _{i=\lceil \alpha \rceil }^m\frac{a_{m,i}\left(i\right)!}{\Gamma \left(\lceil \alpha \rceil -\alpha \right)\Gamma \left(i-\lceil \alpha \rceil +1\right)}{f}_i\left(x\right)\hfill \end{array}$$

(22)

where

$$\begin{array}{c}\hfill a_{m,i}=2^{k(i+1/2)}\sqrt{2m+1}{b}_{m,i},f_i\left(x\right)={\int }_{n/2^k}^{(n+1)/2^k}\frac{{\left(\chi -\frac{n}{2^k}\right)}^{i-\lceil \alpha \rceil }}{{\left(x-\chi \right)}^{\alpha -\lceil \alpha \rceil +1}}d\chi \end{array}$$

Here, Legendre wavelet expansions are employed to approximate ${(x-n/2^k)}^{i-\alpha }$, that is to say

$$\begin{array}{c}\hfill {\left(x-\frac{n}{2^k}\right)}^{i-\alpha }=\sum _{j=0}^{m_0-1}{e}_{i,j}{\psi }_{n,j}\left(x\right),i=\lceil \alpha \rceil ,\lceil \alpha \rceil +1,\cdots ,m\end{array}$$

(23)

where

$$\begin{array}{cc}\hfill e_{i,j}& ={\int }_{n/2^k}^{(n+1)/2^k}{\left(x-\frac{n}{2^k}\right)}^{i-\alpha }{\psi }_{n,j}\left(x\right)dx\hfill \\ & =\sum _{r=0}^j{c}_{j,r}{\int }_{n/2^k}^{(n+1)/2^k}{\left(2^{k}x-n\right)}^{r+i-\alpha }dx=\sum _{r=0}^j\frac{c_{j,r}}{2^k(r+i-\alpha +1)}\hfill \end{array}$$

and

$$\begin{array}{c}\hfill c_{j,r}=\frac{2^{k/2}\sqrt{2m+1}}{2^{k(i-\alpha )}}{b}_{j,r}\end{array}$$

Furthermore, $f_i\left(x\right)$ in Equation (22) can be approximated by the Legendre wavelet in every interval $I_{l,k}\left(x\right)$, $l=n+1,\cdots ,2^k-1$, yielding

$$\begin{array}{c}\hfill f_i\left(x\right)=\sum _{j=0}^{m_0-1}{d}_{i,j}{\psi }_{i,j}\left(x\right)\end{array}$$

(24)

where

$$\begin{array}{cc}\hfill d_{i,j}& ={\int }_{l/2^k}^{(l+1)/2^k}{f}_i\left(x\right){\psi }_{l,j}\left(x\right)dx\hfill \\ & =\sum _{r=0}^j{2}^{k/2}\sqrt{2m+1}{b}_{j,r}\times {\int }_{l/2^k}^{(l+1)/2^k}{f}_i\left(x\right){\left(2^{k}x-l\right)}^{r}dx\hfill \end{array}$$

Substituting Equations (23) and (24) into Equation (22), we can obtain

$$\begin{array}{c}\hfill D_x^{\alpha }{\psi }_{n,m}\left(x\right)=\sum _{j=0}^{m_0-1}{\Omega }_{\alpha }^{\left(n\right)}(m,j){\psi }_{n,j}\left(x\right)+\sum _{l=n+1}^{2^k-1}\sum _{j=0}^{m_0-1}{\widehat{\Omega }}_{\alpha }^{\left(n\right)}(m,j){\psi }_{l,j}\left(x\right)\end{array}$$

(25)

where

$$\left\{\begin{array}{c}{\Omega }_{\alpha }^{\left(n\right)}(m,j)=\sum _{i=\lceil \alpha \rceil }^m{\Theta }_{m,j,i}\hfill \\ {\Theta }_{m,j,i}=\frac{a_{m,i}\left(i\right)!}{\Gamma (i-\alpha +1)}\times \sum _{r=0}^j\frac{c_{j,r}}{2^k\left(r+i-\alpha +1\right)}\hfill \end{array}\right.$$

and

$$\left\{\begin{array}{c}{\widehat{\Omega }}_{\alpha }^{\left(n\right)}(m,j)=\sum _{i=\lceil \alpha \rceil }^m{\widehat{\Theta }}_{m,j,i}\hfill \\ {\widehat{\Theta }}_{m,j,i}=\frac{a_{m,i}\left(i\right)!}{\Gamma (\lceil \alpha \rceil -\alpha )\Gamma (i-\lceil \alpha \rceil +1)}\times \sum _{r=0}^j{2}^{k/2}\sqrt{2m+1}{b}_{j,r}\times \hfill \\ {\int }_{l/2^k}^{(l+1)/2^k}{f}_i\left(x\right){\left(2^{k}x-l\right)}^{r}dx\hfill \end{array}\right.$$

After simplification, ${\Theta }_{m,j,i}$ and ${\widehat{\Theta }}_{m,j,i}$ can be written as:

$$\begin{array}{c}\hfill {\Theta }_{m,j,i}=\frac{{(-1)}^{m+i+j}{2}^{k\alpha }(2m+1)(m+i)!\Gamma (i-\alpha +2)\Gamma (j-i+\alpha )}{(i-\alpha +1)(m-i)!\left(i\right)!\Gamma (-i+\alpha )\Gamma (i-\alpha +1)\Gamma (j+i-\alpha +2)}\end{array}$$

(26)

and

$$\begin{array}{c}\hfill {\widehat{\Theta }}_{m,j,i}=\frac{2^{k(i+1)}(2m+1)b_{m,i}\left(i\right)!}{\Gamma (\lceil \alpha \rceil -\alpha \left)\Gamma \right(i-\lceil \alpha \rceil +1)}\sum _{r=0}^j{b}_{j,r}{\int }_{l/2^k}^{(l+1)/2^k}{f}_i\left(x\right){(2^k-l)}^{r}dx\end{array}$$

(27)

Consequently, Equation (25) can be written as

$$\begin{array}{cc}\hfill D_x^{\alpha }{\psi }_{n,m}\left(x\right)& =\left[{\Omega }_{\alpha }^{\left(n\right)}(m,0),{\Omega }_{\alpha }^{\left(n\right)}(m,1),\cdots ,{\Omega }_{\alpha }^{\left(n\right)}(m,m_0-1)\right]{\Psi }_n\left(x\right)+\hfill \\ & \sum _{l=n+1}^{2^k-1}\left[{\widehat{\Omega }}_{\alpha }^{\left(n\right)}(m,0),{\widehat{\Omega }}_{\alpha }^{\left(n\right)}(m,1),\cdots ,{\widehat{\Omega }}_{\alpha }^{\left(n\right)}(m,m_0-1)\right]{\Psi }_l\left(x\right)\hfill \end{array}$$

(28)

where

$$\begin{array}{c}\hfill {\Psi }_s\left(x\right)=\left[{\psi }_{s,0}\left(x\right),{\psi }_{s,1}\left(x\right),\cdots ,{\psi }_{s,(m_0-1)}\left(x\right)\right],s=0,1,2\cdots ,2^k-1\end{array}$$

Remark 1.

For $\alpha =\lceil \alpha \rceil$, from Caputo’s derivative, we can obtain

$$D_x^{\alpha }{\psi }_{n,m}\left(x\right)=\left\{\begin{array}{c}0,0⩽m<\lceil \alpha \rceil \hfill \\ \sum _{j=0}^{m_0-1}{\Omega }_{\alpha }^{\left(n\right)}(m,j){\psi }_{n,j}\left(x\right),m⩾\lceil \alpha \rceil \hfill \end{array}\right.$$

(29)

Let $\Psi \left(x\right)$ be the Legendre wavelet vector in Equation (18), $D^{\alpha }$ ($\alpha >0$, and ($\lfloor \alpha \rfloor -1<\alpha <\lfloor \alpha \rfloor$)) be the $\alpha$-order differential operational matrix of the Legendre wavelet, then

$$\begin{array}{c}\hfill D_x^{\alpha }\Psi \left(x\right)\simeq D^{\alpha }\Psi \left(x\right)\end{array}$$

(30)

where

$$\begin{array}{c}\hfill D^{\alpha }={\left[\begin{array}{ccccc}{B}^{\alpha }& F^{\alpha }& \cdots & F^{\alpha }& F^{\alpha }\\ 0& B^{\alpha }& \cdots & F^{\alpha }& F^{\alpha }\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & B^{\alpha }& F^{\alpha }\\ 0& 0& \cdots & 0& B^{\alpha }\end{array}\right]}_{M\times M}\end{array}$$

where matrices $B^{\alpha }$ and $F^{\alpha }$ are given as

$$\begin{array}{c}\hfill B^{\alpha }={\left[\begin{array}{cccc}0& 0& \cdots & 0\\ ⋮& ⋮& \cdots & ⋮\\ 0& 0& \cdots & 0\\ {\Omega }_{\alpha }^{\left(n\right)}(\lceil \alpha \rceil ,0)& {\Omega }_{\alpha }^{\left(n\right)}(\lceil \alpha \rceil ,1)& \cdots & {\Omega }_{\alpha }^{\left(n\right)}(\lceil \alpha \rceil ,m_0-1)\\ {\Omega }_{\alpha }^{\left(n\right)}(\lceil \alpha \rceil +1,0)& {\Omega }_{\alpha }^{\left(n\right)}(\lceil \alpha \rceil +1,1)& \cdots & {\Omega }_{\alpha }^{\left(n\right)}(\lceil \alpha \rceil +1,m_0-1)\\ ⋮& ⋮& \cdots & ⋮\\ {\Omega }_{\alpha }^{\left(n\right)}(m_0-1,0)& {\Omega }_{\alpha }^{\left(n\right)}(m_0-1,1)& \cdots & {\Omega }_{\alpha }^{\left(n\right)}(m_0-1,m_0-1)\end{array}\right]}_{m_0\times m_0}\end{array}$$

and

$$\begin{array}{c}\hfill F^{\alpha }={\left[\begin{array}{cccc}0& 0& \cdots & 0\\ ⋮& ⋮& \cdots & ⋮\\ 0& 0& \cdots & 0\\ {\widehat{\Omega }}_{\alpha }^{\left(n\right)}(\lceil \alpha \rceil ,0)& {\widehat{\Omega }}_{\alpha }^{\left(n\right)}(\lceil \alpha \rceil ,1)& \cdots & {\widehat{\Omega }}_{\alpha }^{\left(n\right)}(\lceil \alpha \rceil ,m_0-1)\\ {\widehat{\Omega }}_{\alpha }^{\left(n\right)}(\lceil \alpha \rceil +1,0)& {\widehat{\Omega }}_{\alpha }^{\left(n\right)}(\lceil \alpha \rceil +1,1)& \cdots & {\widehat{\Omega }}_{\alpha }^{\left(n\right)}(\lceil \alpha \rceil +1,m_0-1)\\ ⋮& ⋮& \cdots & ⋮\\ {\widehat{\Omega }}_{\alpha }^{\left(n\right)}(m_0-1,0)& {\widehat{\Omega }}_{\alpha }^{\left(n\right)}(m_0-1,1)& \cdots & {\widehat{\Omega }}_{\alpha }^{\left(n\right)}(m_0-1,m_0-1)\end{array}\right]}_{m_0\times m_0}\end{array}$$

Here, for a better understanding of the calculation procedure, $B^{\left(1.5\right)}$ and $F^{\left(1.5\right)}$ are provided as below, where $m_0=4$, $k=2$, and $\alpha =1.5$ are selected.

$$\begin{array}{c}\hfill B^{\left(1.5\right)}=\frac{1}{\sqrt{\pi }}\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 640& 128& -\frac{128}{7}& \frac{128}{21}\\ -896& 640& \frac{896}{3}& -\frac{640}{11}\end{array}\right]\end{array}$$

$$\begin{array}{c}\hfill F^{\left(1.5\right)}=\frac{1}{\sqrt{\pi }}\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 1280*(1-\sqrt{2})& 2^8*(3\sqrt{2}-4)& \frac{2^7*(166-118\sqrt{2})}{7}& -\frac{2^8*(538\sqrt{2}-824)}{21}\\ 2^8*(91\sqrt{2}-28)& -2^8*(29\sqrt{2}-51)& -\frac{2^8*(400-281\sqrt{2})}{3}& \frac{2^8*(10245-9557\sqrt{2})}{33}\end{array}\right]\end{array}$$

In order to solve the time fractional B-S model of the European option, $U(x,\tau )$ in Equation (9) can be approximated with

$$\begin{array}{c}\hfill U(x,\tau )\simeq \Psi {\left(x\right)}^{T}U\Psi \left(\tau \right)\end{array}$$

(31)

where $U={\left[u_{i,j}\right]}_{M\times M}$ is an unknown matrix, which needs to be determined; from Equation (30), we can obtain

$$\left\{\begin{array}{cc}\hfill \frac{{\partial }^{\alpha }U(x,\tau )}{\partial {\tau }^{\alpha }}& =\Psi {\left(x\right)}^{T}U{D}^{\alpha }\Psi \left(\tau \right)\hfill \\ \hfill \frac{\partial U(x,\tau )}{\partial x}& =\Psi {\left(x\right)}^T{D}^{T}U\Psi \left(\tau \right)\hfill \\ \hfill \frac{{\partial }^{2}U(x,\tau )}{\partial x^2}& =\Psi {\left(x\right)}^T{\left(D^T\right)}^{2}U\Psi \left(\tau \right)\hfill \end{array}\right.$$

(32)

From Definition 1, Equation (32) can be rewritten as

$$\left\{\begin{array}{cc}\hfill \frac{{\partial }^{\alpha }U(x,\tau )}{\partial {\tau }^{\alpha }}& =Vec\left(\Psi {\left(\tau \right)}^T\otimes \Psi {\left(x\right)}^T\right)\left[{\left(D^{\alpha }\right)}^T\otimes I_M\right]Vec\left(U\right)\hfill \\ \hfill \frac{\partial U(x,\tau )}{\partial x}& =Vec\left(\Psi {\left(\tau \right)}^T\otimes \Psi {\left(x\right)}^T\right)\left[I_M\otimes D^T\right]Vec\left(U\right)\hfill \\ \hfill \frac{{\partial }^{2}U(x,\tau )}{\partial x^2}& =Vec\left(\Psi {\left(\tau \right)}^T\otimes \Psi {\left(x\right)}^T\right)\left[I_M\otimes {\left(D^T\right)}^2\right]Vec\left(U\right)\hfill \end{array}\right.$$

where

$$\begin{array}{c}\hfill Vec\left(\Psi {\left(\tau \right)}^T\otimes \Psi {\left(x\right)}^T\right)={\left[{\psi }_1\left(\tau \right)\Psi \left(x\right),{\psi }_2\left(\tau \right)\Psi \left(x\right),\cdots ,{\psi }_M\left(\tau \right)\Psi \left(x\right)\right]}_{1\times M^2}\end{array}$$

and $I_M$ is an M-dimensional identity matrix. Then, Equation (9) can be converted to

$$\left\{\begin{array}{cc}\hfill {\Psi }^T\left(x\right)U{D}^{\alpha }\Psi \left(\tau \right)& =c_1\Psi {\left(x\right)}^T{\left(D^2\right)}^{T}U\Psi \left(\tau \right)+c_2\Psi {\left(x\right)}^T{D}^{T}U\Psi \left(\tau \right)+\hfill \\ & c_3\Psi {\left(x\right)}^{T}U\Psi \left(\tau \right)\hfill \\ \hfill {\Psi }^T\left(x\right)U\Psi \left(0\right)& =\Pi \left(x\right)\hfill \\ \hfill {\Psi }^T\left(B_d\right)U\Psi \left(\tau \right)& =P\left(\tau \right)\hfill \\ \hfill {\Psi }^T\left(B_u\right)U\Psi \left(\tau \right)& =Q\left(\tau \right)\hfill \end{array}\right.$$

(33)

Define

$$\begin{array}{c}\hfill Z={\left(D^{\alpha }\right)}^T\otimes I_M-c_1\left(I_M\otimes D^T\right)-c_2\left[I_M\otimes {\left(D^T\right)}^2\right]-c_3\left(I_M\otimes I_M\right)\end{array}$$

and let

$$\begin{array}{cc}\hfill \mathsf{\Delta }x& =\frac{B_u-B_d}{N_s},x_i=B_d+i\mathsf{\Delta }x,i=0,1,2,\cdots ,N_s\hfill \\ \hfill \mathsf{\Delta }\tau & =T/N_t,{\tau }_j=j\mathsf{\Delta }\tau ,j=0,1,2,\cdots ,N_t\hfill \end{array}$$

be the step size of stock price and time, respectively. Denoting $U_i^j=U(x_i,{\tau }_j)$ yields

$$\left\{\begin{array}{c}\left[Vec\left(\Psi {\left(x_i\right)}^T\otimes \Psi {\left({\tau }_j\right)}^T\right)Z\right]Vec\left(U\right)=0\hfill \\ \left[Vec\left(\Psi {\left(x_i\right)}^T\otimes \Psi {\left({\tau }_0\right)}^T\right)\right]Vec\left(U\right)=U_i^0\hfill \\ \left[Vec\left(\Psi {\left(x_0\right)}^T\otimes \Psi {\left({\tau }_j\right)}^T\right)\right]Vec\left(U\right)=U_0^j\hfill \\ \left[Vec\left(\Psi {\left(x_{N_s}\right)}^T\otimes \Psi {\left({\tau }_j\right)}^T\right)\right]Vec\left(U\right)=U_{N_s}^j\hfill \end{array}\right.$$

(34)

which can be expressed as

$$\begin{array}{c}\hfill HVec\left(U\right)=Y\end{array}$$

(35)

Let

$$\begin{array}{c}\hfill L={\left[\begin{array}{cccc}{\psi }_1\left(x_0\right)& {\psi }_2\left(x_0\right)& \cdots & {\psi }_M\left(x_0\right)\\ {\psi }_1\left(x_1\right)& {\psi }_2\left(x_1\right)& \cdots & {\psi }_M\left(x_1\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\psi }_1\left(x_{N_s}\right)& {\psi }_2\left(x_{N_s}\right)& \cdots & {\psi }_M\left(x_{N_s}\right)\end{array}\right]}_{(N_s+1)\times M}\end{array}$$

and

$$\begin{array}{c}\hfill K={\left[\begin{array}{cccc}{\psi }_1\left({\tau }_0\right)& {\psi }_2\left({\tau }_0\right)& \cdots & {\psi }_M\left({\tau }_0\right)\\ {\psi }_1\left({\tau }_1\right)& {\psi }_2\left({\tau }_1\right)& \cdots & {\psi }_M\left({\tau }_1\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\psi }_1\left({\tau }_{N_t}\right)& {\psi }_2\left({\tau }_{N_t}\right)& \cdots & {\psi }_M\left({\tau }_{N_t}\right)\end{array}\right]}_{(N_t+1)\times M}\end{array}$$

then

$$\begin{array}{c}\hfill H=\left[\begin{array}{c}\left(L\otimes K\right)Z\\ L\otimes K(1,:)\\ L(1,:)\otimes K\\ L(N_s+1,:)\otimes K\end{array}\right]\end{array}$$

and

$$\begin{array}{c}\hfill Y=\left[\stackrel{(N_s+1)(N_t+1)}{\overbrace{0,\cdots ,0}},\underset{N_s+1}{\underbrace{U_0^0,\cdots ,U_{N_s}^0}},\stackrel{N_t+1}{\overbrace{U_0^0,\cdots ,U_0^{N_t}}},\underset{N_t+1}{\underbrace{U_{N_s}^0,\cdots ,U_{N_s}^{N_t}}}\right]\end{array}$$

5. 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.

Algorithm 1 Extreme learning machine algorithm.

Input: $A={\left\{a_i\right\}}_{i=1}^n$, $B={\left\{b_i\right\}}_{i=1}^n$, $X={\left\{x_i\right\}}_{i=1}^n$ $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 $\beta =H^{†}Y$. where $H^{†}$ is the Moore–Penrose generalized inverse of the hidden layer output matrix H, Y is the training data target.

Theorem 4.

The equation $H\beta =Y$ can be solved in three cases:

1.

If matrix H is square and invertible, then $\beta =H^{-1}Y$.

2.

If matrix H is rectangle, then $\beta =H^{†}Y$; β is the minimal least square solution; in other words, $\beta =argmin∥Y\beta -T∥$.

3.

If matrix H is singular, then $\beta =H^{†}Y$, where $H^{†}=H^T{\left(\lambda I+H{H}^T\right)}^{-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.

Algorithm 2 The procedure for solving the European options pricing model based on the LWNN-ELM.

Input: $x_i=B_d+i\mathsf{\Delta }x$, $\mathsf{\Delta }x=\left(B_u-B_d\right)/N_s$, $i=0,1,2,\cdots ,N_s$ ${\tau }_j=j\mathsf{\Delta }\tau$, $\mathsf{\Delta }\tau =T/N_t$, $j=0,1,2,\cdots ,N_t$ ${\psi }_1\left(x\right)$, ⋯, ${\psi }_M\left(x\right)$, ${\psi }_1\left(\tau \right)$, ⋯, ${\psi }_M\left(\tau \right)$ Step 1: Constructing a numerical solution by the Legendre wavelet as the activation function, that is $U^{*}(x,\tau )={\sum }_{i=0}^n{\sum }_{j=0}^n{\beta }_{i,j}{\Psi }_i\left(x\right){\Psi }_j\left(\tau \right)=\Psi (x,\tau )\beta$. Step 2: Substituting the numerical solution $U^{*}(x,\tau )$ and its derivative into the PDEs, the initial condition, and its boundary. Step 3: Solving $H\beta =Y$ by ELM and obtaining the network weights $\beta =H^{†}Y$, $\beta =argmin∥H\beta -Y∥$. Step 4: To obtain the numerical solution $U^{*}(x,\tau )={\sum }_{i=0}^n{\sum }_{j=0}^n{\beta }_{i,j}{\Psi }_i\left(x\right){\Psi }_j\left(\tau \right)=\Psi (x,\tau )\beta$.

Remark 2.

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

6. 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 $\alpha =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. Time consumption of the LWNN-ELM and IDM for different expirations.

		Method	LWNN-ELM	IDM
	Time Consumption (seconds)
Expiration (years)
1			5.2853	75.2167
0.5			4.3862	62.4791
2			6.7541	83.3754

; 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 $\alpha$ on the options price is investigated. From Equation (6), intuitively, $\alpha$ 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 $\alpha$ 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 $\alpha$.

Table 2. Time consumption of the LWNN-ELM for different expirations.

		Method	LWNN-ELM
	Time Consumption (seconds)
Expiration (years)
1			27.6381
0.5			24.9635
2			31.7652

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 $\overline{\Pi }\left(S\right)=max\{K-S,0\}$, $\overline{Q}\left(t\right)=q=0$, and $\overline{P}\left(t\right)=Kexp(-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. The numerical solutions of the European put option, where the corresponding parameters are $K=50$, $T=1$, $r=0.01$, $\sigma =0.1$, $A=10$, $B=100$, and $M=16$.

(${\mathit{N}}_{\mathit{s}}$,${\mathit{N}}_{\mathit{t}}$)	Method	Stock Price
		$\mathit{S}=\mathbf{30}$	$\mathit{S}=\mathbf{40}$	$\mathit{S}=\mathbf{50}$	$\mathit{S}=\mathbf{60}$	$\mathit{S}=\mathbf{70}$	$\mathit{S}=\mathbf{80}$
(230,500)	LWNN-ELM	17.6089	9.0658	3.5568	1.1969	0.4472	0.1727
	IDM	17.6088	9.0656	3.5571	1.1968	0.4473	0.1726
(460,500)	LWNN-ELM	17.5861	9.2083	3.5714	1.2866	0.4707	0.1757
	IDM	17.5862	9.2083	3.5715	1.2867	0.4706	0.1755

; 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

$$\mathrm{Order}={log}_2\frac{{∥e_{\mathsf{\Delta }\tau }∥}_{\infty }}{{∥e_{(\mathsf{\Delta }\tau /2)}∥}_{\infty }}$$

(36)

where $e_{\left(\mathsf{\Delta }\tau \right)}=U_{\left(\mathsf{\Delta }\tau \right)}-U_{(\mathsf{\Delta }\tau /2)}$ and $U_{\left(\mathsf{\Delta }\tau \right)}$ is the numerical solution at the final time with step size $\mathsf{\Delta }\tau$. 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. The convergence order of the European put options price. $N_t$ is the number of discretization points regarding time; the number of discretization points regarding the stock price is fixed as $N_s=460$.

	$\mathit{\alpha }$	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
${\mathit{N}}_{\mathit{t}}$
250		1.02	1.05	1.02	1.04	1.02	1.04	1.02	1.02	0.98
500		1.01	1.01	1.03	1.03	1.01	1.00	1.04	0.99	1.01
1000		1.01	1.01	1.00	1.02	1.00	1.03	1.02	1.02	0.99
2000		1.00	1.00	1.02	1.01	1.01	1.00	0.99	1.01	0.99

. It can be observed that the LWNN-ELM has first-order convergence for different $\alpha$.

7. 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 $\alpha$ 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.