An Analysis of the Fractional-Order Option Pricing Problem for Two Assets by the Generalized Laplace Variational Iteration Approach

Abstract

An option is the right to buy or sell a good at a predetermined price in the future. For customers or financial companies, knowing an option’s pricing is crucial. It is well recognized that the Black–Scholes model is an effective tool for estimating the cost of an option. The Black–Scholes equation has an explicit analytical solution known as the Black–Scholes formula. In some cases, such as the fractional-order Black–Scholes equation, there is no closed form expression for the modified Black–Scholes equation. This article shows how to find the approximate analytic solutions for the two-dimensional fractional-order Black–Scholes equation based on the generalized Riemann–Liouville fractional derivative. The generalized Laplace variational iteration method, which incorporates the generalized Laplace transform with the variational iteration method, is the methodology used to discover the approximate analytic solutions to such an equation. The expression of the two-parameter Mittag–Leffler function represents the problem’s approximate analytical solution. Numerical investigations demonstrate that the proposed scheme is accurate and extremely effective for the two-dimensional fractional-order Black–Scholes Equation in the perspective of the generalized Riemann–Liouville fractional derivative. This guarantees that the generalized Laplace variational iteration method is one of the effective approaches for discovering approximate analytic solutions to fractional-order differential equations.

1. Introduction

The right to purchase or sell a basic product at a specific price in the future is known as an option. Options have a significant presence on marketplaces and exchanges. Determining the prices of an option is important for customers or financial companies. The valuation of options is one of the most important challenges in the field of financial investing. It is well known that the Black–Scholes model [1,2] is an effective instrument for figuring out an option’s cost. There are analytical and numerical approaches used by researchers to solve the Black–Scholes Equations [3,4,5,6,7,8,9,10].

We observe that the Black–Scholes Equations (1)–(3) are partial differential equations with integer-order derivatives. Further study [11,12,13,14] demonstrates that the globalized financial markets are fractal in nature. This illustrates that the traditional Black–Scholes model does not adequately reflect the actual financial market. Studies confirmed the applicability of fractional differential equations many years ago, demonstrating their usefulness for researching aspects linked to fractal geometry and fractal dynamics. Secondly, fractional differential equations provide several benefits in describing significant phenomena in a variety of disciplines, including electromagnetics, fluid flow, acoustics, electrochemistry, as well as material science [15,16,17,18]. Is it reasonable to use the fractional differential equation in the financial market? The answer to the question is “yes”. Fractional derivatives can be used in the financial market because they have a property called “self-similarity”. Further, fractional derivatives respond better to long-term repositories than integer order derivatives. The fractional derivative’s remarkable abilities are employed to solve the fractal complexity in the financial market. At the present, there is an increase in the number of publications that discuss the use of fractional calculus in financial theory [19].

The Black–Scholes equation with two assets of a European call option is defined by:

$$\frac{\partial u}{\partial \tau }+\frac{1}{2}{\sigma }_1^2{S}_1^2\frac{{\partial }^{2}u}{\partial S_1^2}+\frac{1}{2}{\sigma }_1^2{S}_2^2\frac{{\partial }^{2}u}{\partial S_2^2}+\omega {\sigma }_1{\sigma }_2\frac{{\partial }^{2}u}{\partial S_1\partial S_2}+r\left(S_1\frac{\partial u}{\partial S_1}+S_2\frac{\partial u}{\partial S_2}\right)-ru=0,$$

(1)

for ($S_1$,$S_2$,$\tau$ ) ∈ [0,∞ ) × [0,∞ ) × $[0,T]$, with the terminal condition:

$$u(S_1,S_2,T)=max\left\{{\beta }_1{S}_1+{\beta }_2{S}_2-K,0\right\}\mathrm{for}(S_1,S_2)\in [0,\infty )\times [0,\infty ),$$

(2)

where u is the call option depending on the underlying asset prices $S_1$, $S_2$ at time $\tau$;

• $S_1$, $S_2$ are the asset price variables;
• ${\sigma }_1$, ${\sigma }_2$ are the volatility function of underlying assets;
• ${\beta }_1$, ${\beta }_2$ are coefficients so that all risky asset price are at the same level;
• $\omega$ is the volatility of $S_1$ and $S_2$;
• r is the risk-free interest rate;
• T is the expiration date;
• $K=max\{K_1,K_2\}$ where $K_i$ is strike price of the ith underlying asset.

K. Trachoo, W. Sawangtong, and P. Sawangtong [20] researched the two-dimensional Black–Scholes equation with European call option (1) and (2) in 2017. Using the Laplace transform homotopy perturbation approach, they demonstrated that the explicit solution to this issue is represented as a Mellin–Ross function.

P. Sawangtong, K. Trachoo, W. Sawangtong, and B. Wiwattanapataphee [21] investigated the modified Black–Scholes model of (1) and (2) with two assets based on the Liouville–Caputo fractional derivative in 2018. They established, using the Laplace transform homotopy perturbation technique, that the explicit solution to this problem is represented as the Generalized Mittag–Leffler function.

2. The Modified Black–Scholes Equation

The modified Black–Scholes equation in fractional-order derivative form is presented in this section. Let $x=ln\left(S_1\right)-\left(r-\frac{1}{2}{\sigma }_1^2\right)\tau ,y=ln\left(S_2\right)-\left(r-\frac{1}{2}{\sigma }_2^2\right)\tau ,$ $t=T-\tau$ and $u(S_1,S_2,\tau )=e^{-r(T-\tau )}v(x,y,t).$ Readers may find out more information for transformation in [21]. Without loss of generality, we consider the variables x and y by $x\in [0,x_{max}]$ and $y\in [0,y_{max}]$ where $x_{max}$ and $y_{max}$ are positive constants. In the end, the Black–Scholes partial differential equation with two assets of the European call option Equations (1) and (2) is capable of being converted into the following form:

$$v_t=\frac{1}{2}{\sigma }_1^2\frac{{\partial }^{2}v}{\partial x^2}+\frac{1}{2}{\sigma }_2^2\frac{{\partial }^{2}v}{\partial y^2}+\omega {\sigma }_1{\sigma }_2\frac{{\partial }^{2}v}{\partial x\partial y},\mathrm{for}(x,y,t)\in [0,x_{max}]\times [0,y_{max}]\times [0,T],$$

(3)

with the initial condition:

$$v(x,y,0)=\max \left\{c_1{e}^x+c_2{e}^y-K,0\right\},\mathrm{for}(x,y)\in [0,x_{max}]\times [0,y_{max}],$$

(4)

where $c_1$ and $c_2$ are constants defined by

$$c_1={\beta }_1{e}^{(r-\frac{1}{2}{\sigma }_1^2)T}\mathrm{and}{c}_2={\beta }_2{e}^{(r-\frac{1}{2}{\sigma }_2^2)T}.$$

(5)

In this study, we extend the previous work [20] and analyze the general form of the Black–Scholes equation in Equations (3) and (4) by replacing the integer-order time derivative with the fractional-order time derivative. The fractional-order Black–Scholes equation based on the generalized Riemann–Liouville fractional derivative with $\alpha \in (0,1)$ is considered in the form:

$$D_t^{\alpha ,\rho }u(x,y,t;\rho ,\alpha )=\frac{1}{2}{\sigma }_1^2\frac{{\partial }^{2}u}{\partial x^2}+\frac{1}{2}{\sigma }_2^2\frac{{\partial }^{2}u}{\partial y^2}+\omega {\sigma }_1{\sigma }_2\frac{{\partial }^{2}u}{\partial x\partial y},\mathrm{for}(x,y,t)\in [0,x_{max}]\times [0,y_{max}]\times [0,T],$$

(6)

with the fractional integral initial condition:

$$I_t^{1-\alpha ,\rho }u(x,y,0;\rho ,\alpha )=max\left\{c_1{e}^x+c_2{e}^y-K,0\right\},\mathrm{for}(x,y)\in [0,x_{max}]\times [0,y_{max}].$$

(7)

where $\rho >0$, $D_t^{\alpha ,\rho }$ and $I_t^{1-\alpha ,\rho }$ denote the generalized Riemann–Liouville fractional-order derivative with order $\alpha$ and the generalized fractional-order integral with order $1-\alpha$, respectively.

The generalized Laplace variational iteration method is a methodology combining the variational iteration approach with the generalized Laplace transform. Analytical solutions are more complex to obtain than numerical solutions, particularly for fractional partial differential equations. Consequently, the analytical solution offers a valuable instrument for analyzing financial behavior. The generalized Laplace variational iteration approach is used in this research to provide the approximate analytic solution of the time fractional-order Black–Scholes model with two assets for the European call option (6) and (7). In addition, the closed-form analytic solution of the fractional-order Black–Scholes model (6) and (7) is investigated under certain requirements.

The structure of the article is as follows. The definitions of the generalized fractional-order derivative and integral are presented in Section 3. Section 4 discusses the generalized Laplace variational iteration technique’s application and convergence analysis. The explicit solution of the fractional-order Black–Scholes equation is provided in Section 5. In Section 6, numerical results with various parameter values can be seen. This work’s conclusion is provided into Section 7.

3. Basic Definitions

In this section, the generalized Riemann–Liouville fractional integral, the generalized Riemann–Liouville fractional derivative, and the generalized Laplace transform with their some properties have been discussed. For more details, readers can see [22,23]. Throughout this article, we assume that $\alpha$ and $\rho$ are constants with $0<\alpha \le 1$ and $\rho >0,$ and we denote the gamma function by $\mathrm{\Gamma }$.

Definition 1. 

The generalized fractional-order integral α of a continuous function $f:[0,\infty )\to R$ is expressed as

$$I_t^{\alpha ,\rho }f\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}\underset{0}{\overset{t}{\int }}{\left(\frac{t^{\rho }-{\tau }^{\rho }}{\rho }\right)}^{\alpha -1}\frac{f\left(\tau \right)}{{\tau }^{1-\rho }}d\tau .$$

Definition 2. 

The generalized Riemann–Liouville fractional-order derivative of α of a continuous function $f:[0,\infty )\to R$ is given as

$$D_t^{\alpha ,\rho }f\left(t\right)=\frac{1}{\mathrm{\Gamma }(1-\alpha )}\left(t^{1-\rho }\frac{d}{dt}\right)\underset{0}{\overset{t}{\int }}{\left(\frac{t^{\rho }-{\tau }^{\rho }}{\rho }\right)}^{-\alpha }\frac{f\left(\tau \right)}{{\tau }^{1-\rho }}d\tau .$$

We next give some properties that deal with the generalized Riemann–Liouville fractional derivative.

Lemma 1.

Let $f:[0,\infty )\to R$ be a continuous function and c be a constant. Then,

• $I_t^{\alpha ,\rho }{\left(\frac{t^{\rho }}{\rho }\right)}^c=\frac{\mathrm{\Gamma }(c+1)}{\mathrm{\Gamma }(c+\alpha +1)}{\left(\frac{t^{\rho }}{\rho }\right)}^{c+\alpha }$,
• $D_t^{\alpha ,\rho }{I}_t^{\alpha ,\rho }f\left(t\right)=f\left(t\right),$
• $D_t^{\alpha ,\rho }c=\frac{{\rho }^{\alpha -1}}{\mathrm{\Gamma }(1-\alpha )}{t}^{-\alpha \rho }.$

The following part discusses the generalized Laplace transform and some of its properties.

Definition 3. 

The generalized Laplace transform of a continuous function $f:[0,\infty )\to R$ is defined as

$${\mathcal{L}}_{\rho }\left\{f\left(t\right)\right\}\left(s\right)=\underset{0}{\overset{\infty }{\int }}{e}^{-s\frac{t^{\rho }}{\rho }}\frac{f\left(t\right)}{t^{1-\rho }}dt,$$

where s is the Laplace transform parameter.

It is important to note that the generalized Laplace transform can be reduced to the Laplace transform when $\rho =1$.

Lemma 2.

Let $f:[0,\infty )\to R$ be a continuous function and c be a constant. Then,

• ${\mathcal{L}}_{\rho }\left\{{\left(\frac{t^{\rho }}{\rho }\right)}^c\right\}=\frac{\mathrm{\Gamma }\left(c\right)}{s^{c+1}},$
• ${\mathcal{L}}_{\rho }\left\{D_t^{\alpha ,\rho }f\left(t\right)\right\}=s^{\alpha }{\mathcal{L}}_{\rho }\left\{f\left(t\right)\right\}\left(s\right)-I_t^{1-\alpha ,\rho }f\left(0\right)$.

Definition 4. 

Let $f:[0,\infty )\to R$ and $g:[0,\infty )\to R$ be continuous functions. The generalized convolution of f and g is defined by

$$f\left(t\right){\ast }_{\rho }g\left(t\right)={\int }_0^{t}f\left({(t^{\rho }-{\tau }^{\rho })}^{\frac{1}{\rho }}\right)g\left(\tau \right)\frac{d\tau }{{\tau }^{1-\rho }},$$

if the integral exists.

Lemma 3. 

Let $f:[0,\infty )\to R$ and $g:[0,\infty )\to R$ be continuous functions. If ${\mathcal{L}}_{\rho }\left\{f\left(t\right)\right\}\left(s\right)$ and ${\mathcal{L}}_{\rho }\left\{g\left(t\right)\right\}\left(s\right)$ exist, then

$${\mathcal{L}}_{\rho }\left\{f\left(t\right){\ast }_{\rho }g\left(t\right)\right\}\left(s\right)={\mathcal{L}}_{\rho }\left\{f\left(t\right)\right\}\left(s\right){\mathcal{L}}_{\rho }\left\{g\left(t\right)\right\}\left(s\right).$$

In the last part of this section, we will introduce a special function that helps us rewrite complex expressions in a simple form.

Definition 5. 

The two-parameter Mittag–Leffler function is defined as follows:

$$E_{\alpha ,\beta }\left(z\right)=\sum _{k=0}^{\infty }\frac{z^k}{\mathrm{\Gamma }(\alpha k+\beta )},$$

where Γ denotes the Gamma function, $\alpha >0,$ $\beta \in R$ and $z\in C.$

It is important to note that, $E_{1,2}\left(z\right)=\frac{e^z-1}{z}.$

4. The General Methodology of the Generalized Laplace Variational Iteration Method

In this section, we apply the generalized variational iteration method to the nonlinear partial differential equation. Assume that $\mathrm{\Omega }$ is the bounded domain. Let us consider the following general fractional differential equation in the generalized Riemann–Liouville fractional derivative

$$D_t^{\alpha ,\rho }u(x,y,t)+R\left[u(x,y,t)\right]+N\left[u(x,y,t)\right]=f(x,y,t)\mathrm{for}(x,y,t)\in \mathrm{\Omega }\times [0,T],$$

(8)

and the generalized Riemann–Liouville fractional initial condition

$$I_t^{1-\alpha ,\rho }u(x,y,0)=g(x,y)\mathrm{for}(x,y)\in \overline{\mathrm{\Omega }}$$

(9)

where $D_t^{\alpha ,\rho }$ and $I_t^{\alpha ,\rho }$ are the generalized Riemann–Liouville fractional derivative and integral of order $0<\alpha \le 1,$ respectively, $R\left[u\right]$ is a linear term, $N\left[u\right]$ is a nonlinear term and f and g are given functions.

In the first step of the process, we find the suitable Lagrange multiplier $\lambda$ that will be found by using properties of the generalized Riemann–Liouville fractional derivative and integral.

Based on the generalized Riemann–Liouville integration, the correction functional for the nonlinear problem (8) and (9) is defined by

$$u_{n+1}(x,y,t)=u_n(x,y,t)+I_t^{\alpha ,\rho }\lambda \left({\left(t^{\rho }-{\tau }^{\rho }\right)}^{\frac{1}{\rho }}\right)[D_t^{\alpha ,\rho }{u}_n(x,y,\tau )+R\left[u_n\right]+N\left[u_n\right]-f(x,y,t)].$$

(10)

4.1. Lagrange Multipliers

Theorem 1. 

The Lagrange multiplier λ for the fractional-order nonlinear partial differential Equations (8) and (9) can be determined by $\lambda \left({\left(t^{\rho }-{\tau }^{\rho }\right)}^{\frac{1}{\rho }}\right)=-1$.

Proof. 

Let us consider the correction functional (10) for the nonlinear problem (8) and (9):

$$\begin{array}{ccc}\hfill u_{n+1}(x,y,t)& =& u_n(x,y,t)+I_t^{\alpha ,\rho }\lambda \left({\left(t^{\rho }-{\tau }^{\rho }\right)}^{\frac{1}{\rho }}\right)[D_{\tau }^{\alpha ,\rho }{u}_n(x,y,t)+R\left[u_n(x,y,t)\right]\hfill \\ & & -N\left[u_n(x,y,t)\right]-f(x,y,t)]\hfill \\ \hfill u_{n+1}(x,y,t)& =& u_n(x,y,t)+\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_0^t{\left(\frac{t^{\rho }-{\tau }^{\rho }}{\rho }\right)}^{\alpha -1}\lambda \left({\left(t^{\rho }-{\tau }^{\rho }\right)}^{\frac{1}{\rho }}\right)D_{\tau }^{\alpha ,\rho }{u}_n(x,y,\tau )\frac{d\tau }{{\tau }^{1-\rho }}\hfill \\ & & +\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_0^t{\left(\frac{t^{\rho }-{\tau }^{\rho }}{\rho }\right)}^{\alpha -1}\lambda \left({\left(t^{\rho }-{\tau }^{\rho }\right)}^{\frac{1}{\rho }}\right)\left[R\left[u_n\right]-N\left[u_n\right]-f(x,y,\tau )\right]\frac{d\tau }{{\tau }^{1-\rho }}.\hfill \end{array}$$

Let $a\left(t\right)=t^{\rho (\alpha -1)}\lambda \left(t\right)$. Thus,

$$\begin{array}{ccc}\hfill u_{n+1}(x,y,t)& =& u_n(x,y,t)+\frac{{\rho }^{1-\alpha }}{\mathrm{\Gamma }\left(\alpha \right)}a\left(t\right){\ast }_{\rho }{D}_t^{\alpha ,\rho }{u}_n(x,y,t)\hfill \\ & & +\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_0^t{\left(\frac{t^{\rho }-{\tau }^{\rho }}{\rho }\right)}^{\alpha -1}\lambda \left({\left(t^{\rho }-{\tau }^{\rho }\right)}^{\frac{1}{\rho }}\right)[R\left[u_n\right]-N\left[u_n\right]\hfill \\ & & -f(x,y,\tau )]\frac{d\tau }{{\tau }^{1-\rho }},\hfill \end{array}$$

where $a\left(t\right){\ast }_{\rho }{D}_t^{\alpha ,\rho }{u}_n(x,y,t)$ is the generalized convolution of a and $D_t^{\alpha ,\rho }{u}_n$. The generalized Laplace variational iteration correction functional will be defined in the following manner:

$$\begin{array}{ccc}\hfill {\mathcal{L}}_{\rho }\left\{u_{n+1}(x,y,t)\right\}\left(s\right)& =& {\mathcal{L}}_{\rho }\left\{u_n(x,y,t)\right\}\left(s\right)+\frac{{\rho }^{1-\alpha }}{\mathrm{\Gamma }\left(\alpha \right)}{\mathcal{L}}_{\rho }\left\{a\left(t\right){\ast }_{\rho }{D}_{\tau }^{\alpha ,\rho }{u}_n(x,y,t)\right\}\left(s\right)\hfill \\ & & +\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\mathcal{L}}_{\rho }\left\{{\int }_0^t{\left(\frac{t^{\rho }-{\tau }^{\rho }}{\rho }\right)}^{\alpha -1}\lambda \left({\left(t^{\rho }-{\tau }^{\rho }\right)}^{\frac{1}{\rho }}\right)\right[R\left[u_n\right]-N\left[u_n\right]\hfill \\ & & -f(x,y,\tau )\left]\frac{d\tau }{{\tau }^{1-\rho }}\right\}\left(s\right)\hfill \end{array}$$

or equivalenty, upon applying the properties of the Laplace transform, we have

$$\begin{array}{ccc}\hfill {\mathcal{L}}_{\rho }\left\{u_{n+1}(x,y,t)\right\}\left(s\right)& =& {\mathcal{L}}_{\rho }\left\{u_n(x,y,t)\right\}\left(s\right)+\frac{{\rho }^{1-\alpha }}{\mathrm{\Gamma }\left(\alpha \right)}{\mathcal{L}}_{\rho }\left\{a\left(t\right)\right\}\left(s\right)\left\{s^{\alpha }{\mathcal{L}}_{\rho }\left\{u_n(x,y,t)\right\}\left(s\right)-g(x,y)\right\}\hfill \\ & & +\frac{1}{\mathrm{\Gamma }\left(\alpha \right))}{\mathcal{L}}_{\rho }\left\{{\int }_0^t{\left(\frac{t^{\rho }-{\tau }^{\rho }}{\rho }\right)}^{\alpha -1}\lambda \left({\left(t^{\rho }-{\tau }^{\rho }\right)}^{\frac{1}{\rho }}\right)\right[R\left[u_n\right]-N\left[u_n\right]\hfill \\ & & -f(x,y,\tau )\left]\frac{d\tau }{{\tau }^{1-\rho }}\right\}\left(s\right)\hfill \end{array}$$

Taking the variation with respect to $u_n$ of both side of the latter equation, leads to

$$\begin{array}{ccc}\hfill \frac{\delta }{\delta u_n}{\mathcal{L}}_{\rho }\left\{u_{n+1}(x,y,t)\right\}\left(s\right)& =& \frac{\delta }{\delta u_n}{\mathcal{L}}_{\rho }\left\{u_n(x,y,t)\right\}\left(s\right)\hfill \\ & & +\frac{\delta }{\delta u_n}\frac{{\rho }^{1-\alpha }}{\mathrm{\Gamma }\left(\alpha \right)}{\mathcal{L}}_{\rho }\left\{a\left(t\right)\right\}\left(s\right)\left\{s^{\alpha }{\mathcal{L}}_{\rho }\left\{u_n(x,y,t)\right\}\left(s\right)-g(x,y)\right\}\hfill \\ & & +\frac{\delta }{\delta u_n}\frac{1}{\mathrm{\Gamma }\left(\alpha \right))}{\mathcal{L}}_{\rho }\left\{{\int }_0^t{\left(\frac{t^{\rho }-{\tau }^{\rho }}{\rho }\right)}^{\alpha -1}\lambda \left({\left(t^{\rho }-{\tau }^{\rho }\right)}^{\frac{1}{\rho }}\right)\right[R\left[u_n\right]-N\left[u_n\right]\hfill \\ & & -f(x,y,\tau )\left]\frac{d\tau }{{\tau }^{1-\rho }}\right\}\left(s\right)\hfill \end{array}$$

and upopn simplification we obtain

$$\begin{array}{ccc}\hfill {\mathcal{L}}_{\rho }\left\{\frac{\delta }{\delta u_n}{u}_{n+1}(x,y,t)\right\}\left(s\right)& =& {\mathcal{L}}_{\rho }\left\{\frac{\delta }{\delta u_n}{u}_n(x,y,t)\right\}\left(s\right)\hfill \\ & & +\frac{{\rho }^{1-\alpha }}{\mathrm{\Gamma }\left(\alpha \right)}{\mathcal{L}}_{\rho }\left\{a\left(t\right)\right\}\left(s\right)\frac{\delta }{\delta u_n}\left\{s^{\alpha }{\mathcal{L}}_{\rho }\left\{u_n(x,y,t)\right\}\left(s\right)-g(x,y)\right\}\hfill \\ & & +\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\mathcal{L}}_{\rho }\left\{{\int }_0^t{(t-\tau )}^{\alpha -1}\lambda \left({\left(t^{\rho }-{\tau }^{\rho }\right)}^{\frac{1}{\rho }}\right)\frac{\delta }{\delta u_n}\right[R\left[u_n\right]-N\left[u_n\right]\hfill \\ & & -f(x,y,\tau )\left]\frac{d\tau }{{\tau }^{1-\rho }}\right\}\left(s\right)\hfill \end{array}$$

Furthermore, the extra condition of $u_{n+1}$ requires that $\frac{\delta }{\delta u_n}{u}_{n+1}(x,y,t)=0$. Moreover, the terms $R\left[u_n\right]$ and $N\left[u_n\right]$ are considered as restricted variations, which implies $\frac{\delta }{\delta u_n}R\left[u_n\right]=0$ and $\frac{\delta }{\delta u_n}N\left[u_n\right]=0$. We then obtain $1+\frac{{\rho }^{1-\alpha }{s}^{\alpha }}{\mathrm{\Gamma }\left(\alpha \right)}{\mathcal{L}}_{\rho }\left\{a\left(t\right)\right\}\left(s\right)=0$ or

$${\mathcal{L}}_{\rho }\left\{a\left(t\right)\right\}\left(s\right)=-\frac{\mathrm{\Gamma }\left(\alpha \right)}{{\rho }^{1-\alpha }{s}^{\alpha }}.$$

The inverse generalized Laplace transform implies that

$$a\left(t\right)=-\frac{1}{{\rho }^{1-\alpha }}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}=-t^{\rho (\alpha -1)}.$$

By the definition of a, the Lagrange multipliers is $\lambda \left(t\right)=-1.$ □

Note that it follows form Theorem 1 that the correction functional (10) associated with (8) and (9), is formed as:

$$\begin{array}{c}\hfill u_{n+1}(x,y,t)=u_n(x,y,t)-I_t^{\alpha ,\rho }[D_t^{\alpha ,\rho }{u}_n(x,y,\tau )+R\left[u_n\right]+N\left[u_n\right]-f(x,y,t)].\end{array}$$

4.2. Convergence Analysis of the Proposed Method

In this section, we study the convergence of the generalized Laplace variational iteration method, when applied to the nonlinear partial differential Equations (8) and (9). The sufficient conditions for convergence of the method and the error estimate are presented. The main results are proposed in the below theorems.

We next define the operator $A:D\left(A\right)\subseteq H\to H$, where $D\left(A\right)$ is the domain of the operator A and $\left(H,{∥·∥}_H\right)$ is a Banach space, by:

$$A\left[u(x,y,t)\right]=-I_t^{\alpha ,\rho }\left[D_t^{\alpha ,\rho }u(x,y,t)+R\left[u(x,y,t)\right]+N\left[u(x,y,t)\right]-f(x,y,t)\right],$$

(11)

and define the sequence ${\left\{v_n\right\}}_{n=0}^{\infty }$ by:

$$\left.\begin{array}{c}{v}_0(x,y,t)=u_0(x,y,t)\hfill \\ v_1(x,y,t)=A\left[v_0\right]\hfill \\ v_2(x,y,t)=A[v_0(x,y,t)+v_1(x,y,t)]\hfill \\ v_3(x,y,t)=A[v_0(x,y,t)+v_1(x,y,t)+v_2(x,y,t)]\hfill \\ ⋮\hfill \\ v_n(x,y,t)=A[v_0(x,y,t)+v_1(x,y,t)+\dots +v_{n-1}(x,y,t)].\hfill \end{array}\right\}$$

(12)

The relationship between the sequences ${\left\{u_n\right\}}_{n=0}^{\infty }$ and ${\left\{v_n\right\}}_{n=0}^{\infty }$ given by (10) and (12), respectively, is shown in the following lemma.

Lemma 4. 

Let ${\left\{u_n\right\}}_{n=0}^{\infty }$ and ${\left\{v_n\right\}}_{n=0}^{\infty }$ be the sequences constructed by (10) and (12), respectively. Then, $u_n={\displaystyle \sum _{k=0}^n}{v}_k$ for any $n=0,1,2,\dots .$

Proof. 

We can deduce from (10) and (12) that $v_1=A\left[v_0\right]$ and $u_1=u_0+A\left[u_0\right].$ This implies that

$$u_1=u_0+v_1.$$

(13)

We also know that from (10) and (12), $v_2=A[v_0+v_1]$ and $u_2=u_1+A\left[u_1\right]$. By (13), we then get

$$u_2=u_1+v_2.$$

(14)

Once again, we find that $v_3=A[v_0+v_1+v_2]$ and $u_3=u_2+A\left[u_2\right]$. This yields, by (13) and (14), that $u_3=u_2+v_3.$ Throughout this procedure, we finally discover that that $u_n=u_{n-1}+v_n$ for any $n\ge 1.$ This will lead to the desired results. □

It is important to note that, if the limit exists, Lemma 4 enables us to get that${lim}_{n\to \infty }{u}_n(x,y,t)={\displaystyle \sum _{n=0}^{\infty }}{v}_n(x,y,t).$

The next lemma shows the convergence of the infinite series ${\displaystyle \sum _{n=0}^{\infty }}{v}_n(x,y,t).$

Lemma 5. 

Assume that there exists a positive real number γ with $\gamma <1$ such that ${∥v_{n+1}∥}_H\le \gamma {∥v_n∥}_H$ for any $n=0,1,2,3,\dots .$ Then, the infinite series ${\displaystyle \sum _{n=0}^{\infty }}{v}_n(x,y,t)$ given by (12) converges.

Proof. 

Let $S_n$ denote the partial sum of the infinite series ${\displaystyle \sum _{n=0}^{\infty }}{v}_n(x,y,t).$ We would like to show that the sequence ${\left\{S_n\right\}}_{n=0}^{\infty }$ is a Cauchy sequence in the Banach space $H.$ It follows from the assumption of the theorem that we have:

$${∥S_{n+1}-S_n∥}_H={∥v_{n+1}∥}_H\le \gamma {∥v_n∥}_H\le {\gamma }^2{∥v_{n-1}∥}_H\le \dots \le {\gamma }^{n+1}{∥v_0∥}_H.$$

(15)

Let n and m be any natural numbers with $n\ge m$. We consider that by (15):

$$\begin{array}{ccc}\hfill {∥S_n-S_m∥}_H& =& {∥\left(S_n-S_{n-1}\right)+\left(S_{n-1}-S_{n-2}\right)+\dots +\left(S_{m+1}-S_m\right)∥}_H\hfill \\ & \le & {∥S_n-S_{n-1}∥}_H+{∥S_{n-1}-S_{n-2}∥}_H+\dots +{∥S_{m+1}-S_m∥}_H\hfill \\ & \le & {\gamma }^n{∥v_0∥}_H+{\gamma }^{n-1}{∥v_0∥}_H+\dots +{\gamma }^{m+1}{∥v_0∥}_H\hfill \\ & =& \frac{1-{\gamma }^{n-m}}{1-\gamma }{\gamma }^{m+1}{∥v_0∥}_H.\hfill \end{array}$$

We can deduce from the fact that $0<\gamma <1$ that:

$$\underset{m\to \infty }{lim}{∥S_n-S_m∥}_H=0.$$

Therefore, ${\left\{S_n\right\}}_{n=0}^{\infty }$ is a Cauchy sequence in the Banach space $H.$ This information indicates that the infinite series ${\displaystyle \sum _{n=0}^{\infty }}{v}_n(x,y,t)$ determined by (12) converges in the Banach space $H.$ Hence, Lemma 1 is proved completely.

The below theorem demonstrates that the convergent series ${\displaystyle \sum _{n=0}^{\infty }}{v}_n(x,y,t)$ is the solution of the nonlinear Equations (8) and (9). □

Lemma 6. 

Let $\varphi (x,y,t)$ be the function such that the infinite series ${\displaystyle \sum _{n=0}^{\infty }}{v}_n(x,y,t)$, determined by (12), converges to $\varphi (x,y,t).$ Then, the function $\varphi (x,y,t)$ is the solution of the nonlinear partial differential Equations (8) and (9) for any $(x,y,t)\in \overline{\mathrm{\Omega }}\times [0,T]$.

Proof. 

By the property of the convergent series ${\displaystyle \sum _{n=0}^{\infty }}{v}_n(x,y,t)$, we get that ${lim}_{n\to \infty }{v}_n(x,y,t)$$=0.$ Let us consider the following:

$$\sum _{n=0}^k\left[v_{n+1}(x,y,t)-v_n(x,y,t)\right]=v_{k+1}(x,y,t)-v_0(x,y,t)$$

and

$$\begin{array}{ccc}& & \sum _{n=0}^{\infty }\left[v_{n+1}(x,y,t)-v_n(x,y,t)\right]\hfill \\ & =& \underset{k\to \infty }{lim}\sum _{n=0}^k\left[v_{n+1}(x,y,t)-v_n(x,y,t)\right]\hfill \\ & =& \underset{k\to \infty }{lim}{v}_{k+1}(x,y,t)-v_0(x,y,t)\hfill \\ & =& -v_0(x,y,t).\hfill \end{array}$$

(16)

Taking the generalized Riemann–Liouville fractional derivative on both sides of (17), we find that:

$$D_t^{\alpha ,\rho }\sum _{n=0}^{\infty }\left[v_{n+1}(x,y,t)-v_n(x,y,t)\right]=\sum _{n=0}^{\infty }{D}_t^{\alpha ,\rho }\left[v_{n+1}(x,y,t)-v_n(x,y,t)\right]=-D_t^{\alpha ,\rho }{v}_0(x,y,t).$$

(17)

It follows from (11) and (12) and the linearity property of operators that we obtain:

$$\begin{array}{ccc}& & D_t^{\alpha ,\rho }\left[v_{n+1}(x,y,t)-v_n(x,y,t)\right]\hfill \\ & =& D_t^{\alpha ,\rho }A[v_0(x,y,t)+v_1(x,y,t)+\dots +v_n(x,y,t)]\hfill \\ & & -D_t^{\alpha ,\rho }A[v_0(x,y,t)+v_1(x,y,t)+\dots +v_{n-1}(x,y,t)]\hfill \\ & =& -D_t^{\alpha ,\rho }{v}_n(x,y,t)-R\left[v_n\right]-N[v_0+v_1+\dots +v_n]+N[v_0+v_1+\dots +v_{n-1}]\hfill \end{array}$$

(18)

for all $n\ge 1.$ As a result of (19), we now have that:

$$\begin{array}{ccc}\hfill & & \sum _{n=0}^k{D}_t^{\alpha ,\rho }\left[v_{n+1}(x,y,t)-v_n(x,y,t)\right]\hfill \\ & =& D_t^{\alpha ,\rho }\left[v_1(x,y,t)-v_0(x,y,t)\right]+\sum _{n=1}^k{D}_t^{\alpha ,\rho }\left[v_{n+1}(x,y,t)-v_n(x,y,t)\right]\hfill \\ & =& D_t^{\alpha ,\rho }{v}_1(x,y,t)-D_t^{\alpha ,\rho }{v}_0(x,y,t)\hfill \\ & & -D_t^{\alpha ,\rho }{v}_1(x,y,t)-R\left[v_1(x,y,t)\right]-N[v_0(x,y,t)+v_1(x,y,t)]+N\left[v_0(x,y,t)\right]\hfill \\ & & -D_t^{\alpha ,\rho }{v}_2(x,y,t)-R\left[v_2(x,y,t)\right]-N[v_0(x,y,t)+v_1(x,y,t)+v_2(x,y,t)]+N[v_0+v_1]\hfill \\ & & ⋮\hfill \\ & & -D_t^{\alpha ,\rho }{v}_k(x,y,t)-R\left[v_k(x,y,t)\right]-N[v_0+v_1+\dots +v_k]+N[v_0+v_1+\dots +v_{k-1}].\hfill \end{array}$$

(19)

The fact that $D_t^{\alpha ,\rho }{v}_1(x,y,t)=D_t^{\alpha ,\rho }A\left[v_0\right]=-D_t^{\alpha ,\rho }{v}_0(x,y,t)-R\left[v_0(x,y,t)\right]-$$N\left[v_0(x,y,t)\right]+f(x,y,t)$, yields that the Equation (19) becomes:

$$\begin{array}{ccc}\hfill & & \sum _{n=0}^k{D}_t^{\alpha ,\rho }\left[v_{n+1}(x,y,t)-v_n(x,y,t)\right]\hfill \\ & =& -D_t^{\alpha ,\rho }{v}_0(x,y,t)-D_t^{\alpha ,\rho }\left[\sum _{n=0}^k{v}_n(x,y,t)\right]-R\left[\sum _{n=0}^k{v}_n(x,y,t)\right]\hfill \\ & & -N\left[\sum _{n=0}^k{v}_n(x,y,t)\right]+f(x,y,t).\hfill \end{array}$$

(20)

From (17) and (20), we obtain that as $k\to \infty ,$

$$D_t^{\alpha ,\rho }\left[\sum _{n=0}^{\infty }{v}_n(x,y,t)\right]+R\left[\sum _{n=0}^{\infty }{v}_n(x,y,t)\right]+N\left[\sum _{n=0}^{\infty }{v}_n(x,y,t)\right]=f(x,y,t),$$

or

$$D_t^{\alpha ,\rho }\varphi (x,y,t)+R\left[\varphi (x,y,t)\right]+N\left[\varphi (x,y,t)\right]=f(x,y,t).$$

This information implies ${\displaystyle \sum _{n=0}^{\infty }}{v}_n(x,y,t)$ is the infinite series solution of the nonlinear partial differential Equations (8) and (9). □

Lemma 7. 

Assume that the infinite series ${\displaystyle \sum _{n=0}^{\infty }}{v}_n(x,y,t)$, where $v_n$ is defined by (12), converges to the solution u of the nonlinear partial differential Equations (8) and (9). If the approximate analytic solution $u_{\mathit{approx}}$ is the truncated series constructed by $u_{\mathit{approx}}=u_M(x,y,t)={\displaystyle \sum _{n=0}^M}{v}_n(x,y,t)$ for any $(x,y,t)\in \overline{\mathrm{\Omega }}\times [0,T]$, then the maximum error norm can be evaluated as

$${∥u-u_M∥}_H\le \frac{1}{1-\gamma }{\gamma }^{M+1}{∥v_0∥}_H$$

where γ is the real number given in Lemma 5.

Proof. 

Let n and M be any natural numbers with $n\ge M$. As discussed in Theorem 1, we obtain that:

$${∥S_n-S_M∥}_H\le \frac{1-{\gamma }^{n-M}}{1-\gamma }{\gamma }^{M+1}{∥v_0∥}_H$$

where $S_n$ is the partial sum of the infinite series ${\displaystyle \sum _{n=0}^{\infty }}{v}_n(x,y,t).$ As we let$n$ approach to infinity, we obtain that:

$${∥u-\sum _{n=0}^M{v}_n∥}_H\le \frac{1-{\gamma }^{n-M}}{1-\gamma }{\gamma }^{M+1}{∥v_0∥}_H.$$

The definition of real number $\gamma$ implies that $1-{\gamma }^{n-M}<1$ and

$${∥u-\sum _{n=0}^M{v}_n∥}_H\le \frac{1}{1-\gamma }{\gamma }^{M+1}{∥v_0∥}_H.$$

The proof of Lemma 7 is therefore complete. □

The following main theorem is the result from Lemmas 5 and 6.

Theorem 2. 

By using the generalized Laplace variational iteration approach, the approximate analytic solution $u_{\mathit{approx}}$ for the general fractional differential Equation (8) with the integral initial condition (9) can be obtained by the following iteration:

$$\begin{array}{c}{u}_0(x,y,t)\mathit{is}\mathit{an}\mathit{arbitrary}\mathit{function},\hfill \\ u_{n+1}(x,y,t)=u_n(x,y,t)-I_t^{\alpha ,\rho }\left[D_{\tau }^{\alpha ,\rho }{u}_n(x,y,t)+R\left[u_n(x,y,t)\right]+N\left[u_n(x,y,t)\right]-f(x,y,t)\right].\hfill \end{array}$$

for any $(x,y,t)\in \overline{\mathrm{\Omega }}\times [0,T]$ and for all $n\ge 1$. Moreover, if ${lim}_{n\to \infty }{u}_n(x,y,t)$ exists, then the analytic solution u for the general fractional differential Equation (8) with the integral initial condition (9) can be found by

$$u(x,y,t)=\underset{n\to \infty }{lim}{u}_n(x,y,t)\mathit{for}\mathit{any}(x,y,t)\in \overline{\mathrm{\Omega }}\times [0,T].$$

5. An Application of the Generalized Laplace Variational Iteration Method

In this part, we will apply the generalized Laplace variational iteration method to the fractional-order Black–Scholes equation based on the generalized Riemann–Liouville fractional derivative with the fractional integral initial condition (6) and (7). By (8) and (9), we set that $R\left[u_n\right]=-\left(\frac{1}{2}{\sigma }_1^2\frac{{\partial }^2{u}_n}{\partial x^2}+\frac{1}{2}{\sigma }_2^2\frac{{\partial }^2{u}_n}{\partial y^2}+\omega {\sigma }_1{\sigma }_2\frac{{\partial }^2{u}_n}{\partial x\partial y}\right)$, $N\left[u_n\right]=0$, $f(x,y,t)=0$ and $g(x,y)=max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}.$ By Theorem 2 and by choosing

$$u_0(x,y,t)=\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}+\frac{1}{\mathrm{\Gamma }(\alpha +1)}\omega e^{x+y}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha },$$

We obtain that:

$$\begin{array}{c}{u}_0(x,y,t)=\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}+\frac{1}{\mathrm{\Gamma }(\alpha +1)}\omega e^{x+y}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha },\hfill \\ u_{n+1}(x,y,t)=u_n(x,y,t)-\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_0^t{\left(\frac{t^{\rho }-{\tau }^{\rho }}{\rho }\right)}^{\alpha -1}(D_{\tau }^{\alpha ,\rho }{u}_n-\frac{1}{2}{\sigma }_1^2\frac{{\partial }^2{u}_n}{\partial x^2}-\frac{1}{2}{\sigma }_2^2\frac{{\partial }^2{u}_n}{\partial y^2}\hfill \\ -\omega {\sigma }_1{\sigma }_2\frac{{\partial }^2{u}_n}{\partial x\partial y})\frac{d\tau }{{\tau }^{1-\rho }},n\ge 0.\hfill \end{array}$$

Note that $I_t^{1-\alpha ,\rho }{u}_0(x,y,0)=max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}.$ The generalized Laplace transform is then used to aid us discover the terms $u_n$ for $n\ge 1$. Let us consider the following:

$$\begin{array}{ccc}\hfill u_{n+1}(x,y,t)& =& u_n(x,y,t)-\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_0^t{\left(\frac{t^{\rho }-{\tau }^{\rho }}{\rho }\right)}^{\alpha -1}{D}_{\tau }^{\alpha ,\rho }{u}_n(x,y,\tau )\frac{d\tau }{{\tau }^{1-\rho }}\hfill \\ & & +\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\int }_0^t{\left(\frac{t^{\rho }-{\tau }^{\rho }}{\rho }\right)}^{\alpha -1}\left(\frac{1}{2}{\sigma }_1^2\frac{{\partial }^2{u}_n}{\partial x^2}+\frac{1}{2}{\sigma }_2^2\frac{{\partial }^2{u}_n}{\partial y^2}+\omega {\sigma }_1{\sigma }_2\frac{{\partial }^2{u}_n}{\partial x\partial y}\right)\frac{d\tau }{{\tau }^{1-\rho }}\hfill \\ & =& u_n(x,y,t)-\frac{{\rho }^{1-\alpha }}{\mathrm{\Gamma }\left(\alpha \right)}{t}^{\rho (\alpha -1)}{\ast }_{\rho }{D}_t^{\alpha ,\rho }{u}_n(x,y,t)\hfill \\ & & +\frac{{\rho }^{1-\alpha }}{\mathrm{\Gamma }\left(\alpha \right)}{t}^{\rho (\alpha -1)}{\ast }_{\rho }\left(\frac{1}{2}{\sigma }_1^2\frac{{\partial }^2{u}_n}{\partial x^2}+\frac{1}{2}{\sigma }_2^2\frac{{\partial }^2{u}_n}{\partial y^2}+\omega {\sigma }_1{\sigma }_2\frac{{\partial }^2{u}_n}{\partial x\partial y}\right).\hfill \end{array}$$

Taking the generalized Laplace transform on both sides of the above equation:

$$\begin{array}{ccc}\hfill {\mathcal{L}}_{\rho }\left\{u_{n+1}\right\}\left(s\right)& =& {\mathcal{L}}_{\rho }\left\{u_n\right\}\left(s\right)-\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\mathcal{L}}_{\rho }\left\{{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}\right\}\left(s\right){\mathcal{L}}_{\rho }\left\{D_t^{\alpha ,\rho }{u}_n(x,y,t)\right\}\left(s\right)\hfill \\ & & +\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\mathcal{L}}_{\rho }\left\{{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}\right\}\left(s\right){\mathcal{L}}_{\rho }\left\{\frac{1}{2}{\sigma }_1^2\frac{{\partial }^2{u}_n}{\partial x^2}+\frac{1}{2}{\sigma }_2^2\frac{{\partial }^2{u}_n}{\partial y^2}+\omega {\sigma }_1{\sigma }_2\frac{{\partial }^2{u}_n}{\partial x\partial y}\right\}\left(s\right)\hfill \\ & =& {\mathcal{L}}_{\rho }\left\{u_n\right\}\left(s\right)-\frac{1}{s^{\alpha }}\left(s^{\alpha }{\mathcal{L}}_{\rho }\left\{u_n\right\}\left(s\right)-I_t^{1-\alpha ,\rho }{u}_n(x,y,0)\right)\hfill \\ & & +\frac{1}{s^{\alpha }}{\mathcal{L}}_{\rho }\left\{\frac{1}{2}{\sigma }_1^2\frac{{\partial }^2{u}_n}{\partial x^2}+\frac{1}{2}{\sigma }_2^2\frac{{\partial }^2{u}_n}{\partial y^2}+\omega {\sigma }_1{\sigma }_2\frac{{\partial }^2{u}_n}{\partial x\partial y}\right\}\left(s\right).\hfill \end{array}$$

The inverse generalized Laplace transform yields that

$$\begin{array}{ccc}\hfill u_{n+1}(x,y,t)& =& {\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{\alpha }}{I}_t^{1-\alpha ,\rho }{u}_n(x,y,0)\right\}\hfill \\ & & +{\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{\alpha }}{\mathcal{L}}_{\rho }\left\{\frac{1}{2}{\sigma }_1^2\frac{{\partial }^2{u}_n}{\partial x^2}+\frac{1}{2}{\sigma }_2^2\frac{{\partial }^2{u}_n}{\partial y^2}+\omega {\sigma }_1{\sigma }_2\frac{{\partial }^2{u}_n}{\partial x\partial y}\right\}\left(s\right)\right\}\left(t\right)\hfill \end{array}$$

(21)

for $n\ge 1.$ Thus, the generalized Laplace variational iteration method for finding the approximate analytic solution of the fractional-order Black–Scholes Equation (6) with the fractional integral condition (7) is defined by:

$$\begin{array}{c}\hfill \begin{array}{c}{u}_0(x,y,t)=\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}+\frac{1}{\mathrm{\Gamma }(\alpha +1)}\omega e^{x+y}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha }\hfill \\ u_{n+1}(x,y,t)={\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{\alpha }}{I}_t^{1-\alpha ,\rho }{u}_n(x,y,0)\right\}+{\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{\alpha }}{\mathcal{L}}_{\rho }\right\{\frac{1}{2}{\sigma }_1^2\frac{{\partial }^2{u}_n}{\partial x^2}+\frac{1}{2}{\sigma }_2^2\frac{{\partial }^2{u}_n}{\partial y^2}\hfill \\ +\omega {\sigma }_1{\sigma }_2\frac{{\partial }^2{u}_n}{\partial x\partial y}\left\}\left(s\right)\right\}\left(t\right)\mathrm{for}n\ge 0,\hfill \end{array}\end{array}$$

with $I_t^{1-\alpha ,\rho }{u}_n(x,y,0)=max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}$ for all $n\ge 0.$

Theorem 3. 

The approximate analytic solution for the two dimensional fractional-order Black–Scholes Equation (6) with the fractional integral condition (7) can be defined by the following iteration:

$$\begin{array}{ccc}\hfill u_n(x,y,t)& =& \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}\hfill \\ & & +\sum _{k=0}^{n-1}\left(\frac{1}{2^{k+1}}{\sigma }_1^{2(k+1)}max\left\{c_1{e}^x,0\right\}+\frac{1}{2^{k+1}}{\sigma }_2^{2(k+1)}max\left\{c_2{e}^y,0\right\}\right)\hfill \\ & & \times \frac{1}{\mathrm{\Gamma }\left(\right(k+2\left)\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{(k+2)\alpha -1}\hfill \\ & & +\frac{1}{\mathrm{\Gamma }\left(\right(n+1)\alpha +1)}{\left[\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\omega {\sigma }_1{\sigma }_2\right]}^n{\left(\frac{t^{\rho }}{\rho }\right)}^{n\alpha }\hfill \end{array}$$

(22)

for any $(x,y,t)\in [0,x_{max}]\times [0,y_{max}]\times [0,T]$ and for all $n\ge 1.$ Furthermore, if

$$\begin{array}{c}\hfill \frac{1}{\mathrm{\Gamma }\left(\right(n+1)\alpha +1)}{\left[\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\omega {\sigma }_1{\sigma }_2\right]}^n{\left(\frac{t^{\rho }}{\rho }\right)}^{n\alpha }\end{array}$$

approaches zero when n goes to infinity for any fixed $t\in [0,T]$, then the analytic solution u for the fractional-order Black–Scholes Equations (6) and (7) is in the form:

$$\begin{array}{ccc}\hfill u(x,y,t)& =& max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}\hfill \\ & & +\frac{{\sigma }_1^{2}max\left\{c_1{e}^x,0\right\}}{2}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha -1}{E}_{\alpha ,2\alpha }\left(\frac{{\sigma }_1^2}{2}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha }\right)\hfill \\ & & +\frac{{\sigma }_2^{2}max\left\{c_2{e}^y,0\right\}}{2}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha -1}{E}_{\alpha ,2\alpha }\left(\frac{{\sigma }_2^2}{2}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha }\right)\hfill \end{array}$$

(23)

where $E_{\alpha ,\beta }$ denotes the two-parameter Mittag–Leffler function.

Proof. 

Let $u_0(x,y,t)=\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}+\frac{1}{\mathrm{\Gamma }(\alpha +1)}\omega e^{x+y}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha }.$ The generalized Laplace variational iteration procedure is then started by computing the term $u_1.$ By (22), we obtain:

$$\begin{array}{ccc}\hfill u_1& =& {\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{\alpha }}{I}_t^{1-\alpha ,\rho }{u}_0(x,y,0)\right\}\left(t\right)+{\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{\alpha }}{\mathcal{L}}_{\rho }\left\{\frac{1}{2}{\sigma }_1^2\frac{{\partial }^2{u}_0}{\partial x^2}+\frac{1}{2}{\sigma }_2^2\frac{{\partial }^2{u}_0}{\partial y^2}+\omega {\sigma }_1{\sigma }_2\frac{{\partial }^2{u}_0}{\partial x\partial y}\right\}\left(s\right)\right\}\left(t\right)\hfill \\ & =& max\left\{c_1{e}^x+\widetilde{{\beta }_2}{e}^y-K,0\right\}{\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{\alpha }}\right\}\left(t\right)\hfill \\ & & +{\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{\alpha }}{\mathcal{L}}_{\rho }\left\{\frac{1}{2}{\sigma }_1^2\frac{{\partial }^2{u}_0}{\partial x^2}+\frac{1}{2}{\sigma }_2^2\frac{{\partial }^2{u}_0}{\partial y^2}+\omega {\sigma }_1{\sigma }_2\frac{{\partial }^2{u}_0}{\partial x\partial y}\right\}\left(s\right)\right\}\left(t\right)\hfill \\ & =& \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}\hfill \\ & & +\frac{1}{2}{\sigma }_1^{2}max\left\{c_1{e}^x,0\right\}\frac{1}{\mathrm{\Gamma }\left(2\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha -1}+\frac{1}{2}{\sigma }_2^{2}max\left\{c_2{e}^y,0\right\}\frac{1}{\mathrm{\Gamma }\left(2\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha -1}\hfill \\ & & +\frac{1}{2}{\sigma }_1^2\omega e^{x+y}\frac{1}{\mathrm{\Gamma }(2\alpha +1)}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha }+\frac{1}{2}{\sigma }_2^2\omega e^{x+y}\frac{1}{\mathrm{\Gamma }(2\alpha +1)}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha }\hfill \\ & & +\omega {\sigma }_1{\sigma }_2\omega e^{x+y}\frac{1}{\mathrm{\Gamma }(2\alpha +1)}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha }.\hfill \end{array}$$

Then, the function $u_1$ is

$$\begin{array}{ccc}\hfill u_1(x,y,t)& =& max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}\hfill \\ & & +\left[\frac{1}{2}{\sigma }_1^{2}max\left\{c_1{e}^x,0\right\}+\frac{1}{2}{\sigma }_2^{2}max\left\{c_2{e}^y,0\right\}\right]\frac{1}{\mathrm{\Gamma }\left(2\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha -1}\hfill \\ & & +\left[\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\omega {\sigma }_1{\sigma }_2\right]\omega e^{x+y}\frac{1}{\mathrm{\Gamma }(2\alpha +1)}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha }.\hfill \end{array}$$

Note that the function $u_1$ satisfies that $I_t^{1-\alpha ,\rho }{u}_1(x,y,0)=max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}.$ We next find the function $u_2.$ By (22), we get:

$$u_2={\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{\alpha }}{I}_t^{1-\alpha ,\rho }{u}_1(x,y,0)\right\}\left(t\right)+{\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{\alpha }}{\mathcal{L}}_{\rho }\left\{\frac{1}{2}{\sigma }_1^2\frac{{\partial }^2{u}_1}{\partial x^2}+\frac{1}{2}{\sigma }_2^2\frac{{\partial }^2{u}_1}{\partial y^2}+\omega {\sigma }_1{\sigma }_2\frac{{\partial }^2{u}_1}{\partial x\partial y}\right\}\left(s\right)\right\}\left(t\right)$$

Consider ${\mathcal{L}}_{\rho }\left\{\frac{1}{2}{\sigma }_1^2\frac{{\partial }^2{u}_1}{\partial x^2}+\frac{1}{2}{\sigma }_2^2\frac{{\partial }^2{u}_1}{\partial y^2}+\omega {\sigma }_1{\sigma }_2\frac{{\partial }^2{u}_1}{\partial x\partial y}\right\}\left(s\right).$

$$\begin{array}{ccc}& & {\mathcal{L}}_{\rho }\left\{\frac{1}{2}{\sigma }_1^2\frac{{\partial }^2{u}_1}{\partial x^2}+\frac{1}{2}{\sigma }_2^2\frac{{\partial }^2{u}_1}{\partial y^2}+\omega {\sigma }_1{\sigma }_2\frac{{\partial }^2{u}_1}{\partial x\partial y}\right\}\left(s\right)\hfill \\ & =& \frac{1}{2}{\sigma }_1^{2}max\left\{c_1{e}^x,0\right\}{\mathcal{L}}_{\rho }\left\{\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}\right\}\left(s\right)+\frac{1}{2}{\sigma }_2^{2}max\left\{c_2{e}^y,0\right\}{\mathcal{L}}_{\rho }\left\{\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}\right\}\left(s\right)\hfill \\ & & +\frac{1}{2^2}{\sigma }_1^{4}max\left\{c_1{e}^x,0\right\}{\mathcal{L}}_{\rho }\left\{\frac{1}{\mathrm{\Gamma }\left(2\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha -1}\right\}\left(s\right)+\frac{1}{2^2}{\sigma }_2^{4}max\left\{c_2{e}^y,0\right\}{\mathcal{L}}_{\rho }\left\{\frac{1}{\mathrm{\Gamma }\left(2\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha -1}\right\}\left(s\right)\hfill \\ & & +\frac{1}{2}{\sigma }_1^2\left[\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\omega {\sigma }_1{\sigma }_2\right]\omega e^{x+y}{\mathcal{L}}_{\rho }\left\{\frac{1}{\mathrm{\Gamma }(2\alpha +1)}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha }\right\}\left(s\right)\hfill \\ & & +\frac{1}{2}{\sigma }_2^2\left[\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\omega {\sigma }_1{\sigma }_2\right]\omega e^{x+y}{\mathcal{L}}_{\rho }\left\{\frac{1}{\mathrm{\Gamma }(2\alpha +1)}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha }\right\}\left(s\right)\hfill \\ & & +\omega {\sigma }_1{\sigma }_2\left[\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\omega {\sigma }_1{\sigma }_2\right]\omega e^{x+y}{\mathcal{L}}_{\rho }\left\{\frac{1}{\mathrm{\Gamma }(2\alpha +1)}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha }\right\}\left(s\right)\hfill \\ & =& \frac{1}{2}{\sigma }_1^{2}max\left\{c_1{e}^x,0\right\}\frac{1}{s^{\alpha }}+\frac{1}{2}{\sigma }_2^{2}max\left\{c_2{e}^y,0\right\}\frac{1}{s^{\alpha }}\hfill \\ & & +\frac{1}{2^2}{\sigma }_1^{4}max\left\{c_1{e}^x,0\right\}\frac{1}{s^{2\alpha }}+\frac{1}{2^2}{\sigma }_2^{4}max\left\{c_2{e}^y,0\right\}\frac{1}{s^{2\alpha }}\hfill \\ & & +\frac{1}{2}{\sigma }_1^2\left[\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\omega {\sigma }_1{\sigma }_2\right]\omega e^{x+y}\frac{1}{s^{2\alpha +1}}\hfill \\ & & +\frac{1}{2}{\sigma }_2^2\left[\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\omega {\sigma }_1{\sigma }_2\right]\omega e^{x+y}\frac{1}{s^{2\alpha +1}}\hfill \\ & & +\omega {\sigma }_1{\sigma }_2\left[\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\omega {\sigma }_1{\sigma }_2\right]\omega e^{x+y}\frac{1}{s^{2\alpha +1}}.\hfill \end{array}$$

Then, the term $u_2$ is determined by

$$\begin{array}{ccc}\hfill u_2(x,y,t)& =& \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}\hfill \\ & & +\frac{1}{2}{\sigma }_1^{2}max\left\{c_1{e}^x,0\right\}{\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{2\alpha }}\right\}\left(t\right)+\frac{1}{2}{\sigma }_2^{2}max\left\{c_2{e}^y,0\right\}{\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{2\alpha }}\right\}\left(t\right)\hfill \\ & & +\frac{1}{2^2}{\sigma }_1^{4}max\left\{c_1{e}^x,0\right\}{\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{3\alpha }}\right\}\left(t\right)+\frac{1}{2^2}{\sigma }_2^{4}max\left\{c_2{e}^y,0\right\}{\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{3\alpha }}\right\}\left(t\right)\hfill \\ & & +\frac{1}{2}{\sigma }_1^2\omega e^{x+y}{\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{3\alpha +1}}\right\}\left(t\right)+\frac{1}{2}{\sigma }_2^2\omega e^{x+y}{\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{3\alpha +1}}\right\}\left(t\right)\hfill \\ & & +\omega {\sigma }_1{\sigma }_2\omega e^{x+y}{\mathcal{L}}_{\rho }^{-1}\left\{\frac{1}{s^{3\alpha +1}}\right\}\left(t\right)\hfill \\ & =& \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}\hfill \\ & & +\left[\frac{1}{2}{\sigma }_1^{2}max\left\{c_1{e}^x,0\right\}+\frac{1}{2}{\sigma }_2^{2}max\left\{c_2{e}^y,0\right\}\right]\frac{1}{\mathrm{\Gamma }\left(2\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha -1}\hfill \\ & & +\left[\frac{1}{2^2}{\sigma }_1^{4}max\left\{c_1{e}^x,0\right\}+\frac{1}{2^2}{\sigma }_2^{4}max\left\{c_2{e}^y,0\right\}\right]\frac{1}{\mathrm{\Gamma }\left(3\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{3\alpha -1}\hfill \\ & & +{\left[\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\omega {\sigma }_1{\sigma }_2\right]}^2\omega e^{x+y}\frac{1}{\mathrm{\Gamma }(3\alpha +1)}{\left(\frac{t^{\rho }}{\rho }\right)}^{3\alpha }.\hfill \end{array}$$

Note that the function $u_2$ satisfies that $I_t^{1-\alpha ,\rho }{u}_2(x,y,0)=max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}.$ As a result of the preceding explanation, we now have term $u_3$ as follows:

$$\begin{array}{ccc}\hfill u_3(x,y,t)& =& \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}\hfill \\ & & +\left[\frac{1}{2}{\sigma }_1^{2}max\left\{c_1{e}^x,0\right\}+\frac{1}{2}{\sigma }_2^{2}max\left\{c_2{e}^y,0\right\}\right]\frac{1}{\mathrm{\Gamma }\left(2\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha -1}\hfill \\ & & +\left[\frac{1}{2^2}{\sigma }_1^{4}max\left\{c_1{e}^x,0\right\}+\frac{1}{2^2}{\sigma }_2^{4}max\left\{c_2{e}^y,0\right\}\right]\frac{1}{\mathrm{\Gamma }\left(3\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{3\alpha -1}\hfill \\ & & +\left[\frac{1}{2^3}{\sigma }_1^{6}max\left\{c_1{e}^x,0\right\}+\frac{1}{2^3}{\sigma }_2^{6}max\left\{c_2{e}^y,0\right\}\right]\frac{1}{\mathrm{\Gamma }\left(4\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{4\alpha -1}\hfill \\ & & +{\left[\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\omega {\sigma }_1{\sigma }_2\right]}^3\omega e^{x+y}\frac{1}{\mathrm{\Gamma }(4\alpha +1)}{\left(\frac{t^{\rho }}{\rho }\right)}^{4\alpha }.\hfill \end{array}$$

Note that the function $u_3$ satisfies that $I_t^{1-\alpha ,\rho }{u}_3(x,y,0)=max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}.$ Using the same manner, we can find the expression of $u_n$ in the following form:

$$\begin{array}{ccc}\hfill u_n(x,y,t)& =& \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}\hfill \\ & & +\sum _{k=0}^{n-1}\left(\frac{1}{2^{k+1}}{\sigma }_1^{2(k+1)}max\left\{c_1{e}^x,0\right\}+\frac{1}{2^{k+1}}{\sigma }_2^{2(k+1)}max\left\{c_2{e}^y,0\right\}\right)\hfill \\ & & \times \frac{1}{\mathrm{\Gamma }\left(\right(k+2\left)\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{(k+2)\alpha -1}\hfill \\ & & +\frac{1}{\mathrm{\Gamma }\left(\right(n+1)\alpha +1)}{\left[\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\omega {\sigma }_1{\sigma }_2\right]}^n{\left(\frac{t^{\rho }}{\rho }\right)}^{n\alpha },\hfill \end{array}$$

(24)

and $I_t^{1-\alpha ,\rho }{u}_n(x,y,0)=max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}$ for all $n\ge 1.$ As a consequence of the generalized Laplace variational iteration process, we derive that, under the assumption of the theorem, the analytic solution for the fractional Black–Scholes equation is as follows:

$$\begin{array}{ccc}\hfill u(x,y,t)& =& \underset{n\to \infty }{lim}{u}_n(x,y,t)\hfill \\ & =& max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}\hfill \\ & & +\sum _{k=0}^{\infty }\left(\frac{1}{2^{k+1}}{\sigma }_1^{2(k+1)}max\left\{c_1{e}^x,0\right\}+\frac{1}{2^{k+1}}{\sigma }_2^{2(k+1)}max\left\{c_2{e}^y,0\right\}\right)\hfill \\ & & \times \frac{1}{\mathrm{\Gamma }\left(\right(k+2\left)\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{(k+2)\alpha -1}\hfill \\ & =& max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}\hfill \\ & & +\frac{{\sigma }_1^{2}max\left\{c_1{e}^x,0\right\}}{2}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha -1}\sum _{k=0}^{\infty }\frac{1}{\mathrm{\Gamma }(\alpha k+2\alpha )}{\left[\frac{{\sigma }_1^2}{2}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha }\right]}^k\hfill \\ & & +\frac{{\sigma }_2^{2}max\left\{c_2{e}^y,0\right\}}{2}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha -1}\sum _{k=0}^{\infty }\frac{1}{\mathrm{\Gamma }(\alpha k+2\alpha )}{\left[\frac{{\sigma }_2^2}{2}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha }\right]}^k\hfill \\ & =& max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}\hfill \\ & & +\frac{{\sigma }_1^{2}max\left\{c_1{e}^x,0\right\}}{2}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha -1}{E}_{\alpha ,2\alpha }\left(\frac{{\sigma }_1^2}{2}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha }\right)\hfill \\ & & +\frac{{\sigma }_2^{2}max\left\{c_2{e}^y,0\right\}}{2}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha -1}{E}_{\alpha ,2\alpha }\left(\frac{{\sigma }_2^2}{2}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha }\right).\hfill \end{array}$$

Therefore, the Theorem 3 is proved completely. □

Corollary 1. 

The analytic solution for the two-dimensional fractional-order Black–Scholes equation based on the Riemann–Liouville fractional derivative with the Riemann–Liouville fractional integral condition is given in the following:

$$\begin{array}{ccc}\hfill u(x,y,t)& =& max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}\frac{t^{\alpha -1}}{\mathrm{\Gamma }\left(\alpha \right)}+\frac{{\sigma }_1^{2}max\left\{c_1{e}^x,0\right\}}{2}{t}^{2\alpha -1}{E}_{\alpha ,2\alpha }\left(\frac{{\sigma }_1^2}{2}{t}^{\alpha }\right)\hfill \\ & & +\frac{{\sigma }_2^{2}max\left\{c_2{e}^y,0\right\}}{2}{t}^{2\alpha -1}{E}_{\alpha ,2\alpha }\left(\frac{{\sigma }_2^2}{2}{t}^{\alpha }\right),\hfill \end{array}$$

(25)

for any $(x,y,t)\in [0,x_{max}]\times [0,y_{max}]\times [0,T].$

Proof. 

By setting $\rho =1$, this corollary can be obtained immediately from Theorem 3. □

Corollary 2. 

The analytic solution of the classical two-dimensional Black–Scholes equation with the European call option is:

$$\begin{array}{c}\hfill u(x,y,t)=max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}+max\left\{c_1{e}^x,0\right\}\left(e^{\frac{{\sigma }_1^2}{2}t}-1\right)+max\left\{c_2{e}^y,0\right\}\left(e^{\frac{{\sigma }_2^2}{2}t}-1\right),\end{array}$$

(26)

for any $(x,y,t)\in [0,x_{max}]\times [0,y_{max}]\times [0,T].$

Proof. 

By Theorem 3, we have that

$$\begin{array}{ccc}\hfill u(x,y,t)& =& max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha -1}+\frac{{\sigma }_1^{2}max\left\{c_1{e}^x,0\right\}}{2}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha -1}\hfill \\ & & \times E_{\alpha ,2\alpha }\left(\frac{{\sigma }_1^2}{2}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha }\right)+\frac{{\sigma }_2^{2}max\left\{c_2{e}^y,0\right\}}{2}{\left(\frac{t^{\rho }}{\rho }\right)}^{2\alpha -1}{E}_{\alpha ,2\alpha }\left(\frac{{\sigma }_2^2}{2}{\left(\frac{t^{\rho }}{\rho }\right)}^{\alpha }\right).\hfill \end{array}$$

Since the generalized Riemann–Liouville fractional derivative can be reduced to the usual derivative when $\rho =\alpha =1,$ we obtain:

$$\begin{array}{ccc}\hfill u(x,y,t)& =& max\left\{c_1{e}^x+c_2{e}^y-K,0\right\}\hfill \\ & & +\frac{{\sigma }_1^{2}max\left\{c_1{e}^x,0\right\}}{2}t{E}_{1,2}\left(\frac{{\sigma }_1^2}{2}t\right)+\frac{{\sigma }_2^{2}max\left\{c_2{e}^y,0\right\}}{2}t{E}_{1,2}\left(\frac{{\sigma }_2^2}{2}t\right),\hfill \end{array}$$

for any $(x,y,t)\in [0,x_{max}]\times [0,y_{max}]\times [0,T].$ □

6. Numerical Results

We will assume throughout this section that the strike price K is $70.$ The risk-free annual interest rate is $5\%,$ thus $r=0.05;$ the maturity time is $T=1$ year; and the volatilities of the underlying assets $x\in [0,5]$ and $y\in [0,5]$ are ${\sigma }_1=5\%$ and ${\sigma }_2=10\%$, respectively. In the following, we will demonstrate the European option prices computed from the previously described analytic solutions. Figure 1 shows the European call-option pricing with various values of the parameter $\alpha$ obtained by (24) with the asset pricing $y=5$ and $\rho =1$.

Figure 2 demonstrates the European call-option pricing with various values of the parameter $\rho$ given by (24) with the asset pricing $y=5$ and $\alpha =0.5$.

In case $\rho =1$, the graphs of the option pricing by (24) with four different values for fractional order $\alpha =0.3,0.5,0.7,$ and $0.9$ are shown in Figure 3.

It is well known that in the classical Black–Scholes model, the option price depends on five variables. These are volatility, the price of the underlying asset, the strike price of the option, the time until expiration of the option, and the risk-free interest rate. With these variables, sellers of options could, in theory, set prices for the options they sell that make sense. However, in the fractional-order Black–Scholes equation proposed here, there are two parameters added: $\alpha$ and $\rho$. In each of the figures, it is shown that the value of call options will increase according to the values of two parameters. Therefore, if we can figure out the right values for these two factors, the modified Black–Scholes equation will give us option prices that are close to what they are worth on the market.

7. Conclusions

It is well known that in option pricing theory, the Black–Scholes equation is one of the most significant models for option pricing. In this manuscript, we developed the classical Black–Scholes equation in the form of the fractional-order Black–Scholes equation based on the generalized Riemann–Liouville derivative. This article provides the approximate analytic solution to the fractional-order Black–Scholes equation via the generalized Laplace variational iteration method. Moreover, we show that the solution to the classical Black–Scholes equation is achieved as a special case of the proposed approximate analytic solution. This demonstrates that the generalized Laplace variational iteration method is one of the effective approaches for discovering approximate analytic solutions to fractional-order differential equations. The advantage of the modified Black–Scholes equation is that it has two parameters occurring in the definition of the fractional derivative, that is $\alpha$ and $\rho .$ If we can correctly estimate the values of these two factors, the option prices produced from the modified form will be close to the market value of option prices. We may utilize the genetic algorithm and the actual value of the option to determine the proper values for two parameters. Option pricing may be determined by the solution to the modified Black–Scholes equation after the proper values for these two parameters are already known.