# Application of the Generalized Laplace Homotopy Perturbation Method to the Time-Fractional Black–Scholes Equations Based on the Katugampola Fractional Derivative in Caputo Type

^{1}

^{2}

^{3}

^{*}

Next Article in Journal

Next Article in Special Issue

Next Article in Special Issue

Previous Article in Journal

Previous Article in Special Issue

Previous Article in Special Issue

Department of Mathematics, Faculty of Science, Mahidol University, Bangkok 10400, Thailand

Center of Excellence in Mathematics, Commission on Higher Education, Bangkok 10400, Thailand

Department of Mathematics, Faculty of Applyied Science, King Mongkut’s University of Technology North Bangkok, Bangkok 10800, Thailand

Author to whom correspondence should be addressed.

Academic Editors: Yongwimon Lenbury, Ravi P. Agarwal, Philip Broadbridge and Dongwoo Sheen

Received: 28 January 2021
/
Revised: 22 February 2021
/
Accepted: 4 March 2021
/
Published: 12 March 2021

(This article belongs to the Special Issue Proceedings of the International Conference in Mathematics and Applications 2020_Mahidol University)

In the finance market, the Black–Scholes equation is used to model the price change of the underlying fractal transmission system. Moreover, the fractional differential equations recently are accepted by researchers that fractional differential equations are a powerful tool in studying fractal geometry and fractal dynamics. Fractional differential equations are used in modeling the various important situations or phenomena in the real world such as fluid flow, acoustics, electromagnetic, electrochemistry and material science. There is an important question in finance: “Can the fractional differential equation be applied in the financial market?”. The answer is “Yes”. Due to the self-similar property of the fractional derivative, it can reply to the long-range dependence better than the integer-order derivative. Thus, these advantages are beneficial to manage the fractal structure in the financial market. In this article, the classical Black–Scholes equation with two assets for the European call option is modified by replacing the order of ordinary derivative with the fractional derivative order in the Caputo type Katugampola fractional derivative sense. The analytic solution of time-fractional Black–Scholes European call option pricing equation with two assets is derived by using the generalized Laplace homotopy perturbation method. The used method is the combination of the homotopy perturbation method and generalized Laplace transform. The analytic solution of the time-fractional Black–Scholes equation is carried out in the form of a Mittag–Leffler function. Finally, the effects of the fractional-order in the Caputo type Katugampola fractional derivative to change of a European call option price are shown.

In 1973, Black and Scholes [1] constructed the mathematical model for pricing an options contract, known as the Black-Scholes model. In particular, the Black-Scholes model estimates the variation over time of financial instruments. The Black-Scholes equation with two assets of a European call option is defined by:
with the terminal condition:
and boundary conditions:
and
where

$$\frac{\partial u}{\partial \tau}}+{\displaystyle \frac{1}{2}}{\sigma}_{1}^{2}{S}_{1}^{2}{\displaystyle \frac{{\partial}^{2}u}{\partial {S}_{1}^{2}}}+{\displaystyle \frac{1}{2}}{\sigma}_{2}^{2}{S}_{2}^{2}{\displaystyle \frac{{\partial}^{2}u}{{\partial}^{2}{S}_{2}^{2}}}+\omega {\sigma}_{1}{\sigma}_{2}{S}_{1}{S}_{2}{\displaystyle \frac{{\partial}^{2}u}{\partial {S}_{1}\partial {S}_{2}}}+r\left[{S}_{1}{\displaystyle \frac{\partial u}{\partial {S}_{1}}}+{S}_{2}{\displaystyle \frac{\partial u}{\partial {S}_{2}}}\right]-ru=0,$$

$$u(T,{S}_{1},{S}_{2},)=\mathrm{max}\left({\beta}_{1}{S}_{1}+{\beta}_{2}{S}_{2}-K,0\right),\mathrm{for}\phantom{\rule{4.25pt}{0ex}}{S}_{1},{S}_{2}\in [0,\infty ),\tau \in [0,T],$$

$$u(\tau ,{S}_{1},{S}_{2})=0\phantom{\rule{4.25pt}{0ex}}\mathrm{as}\phantom{\rule{4.25pt}{0ex}}({S}_{1},{S}_{2})\to (0,0),$$

$$u(\tau ,{S}_{1},{S}_{2})={\beta}_{1}{S}_{1}+{\beta}_{2}{S}_{2}-K{e}^{-r(T-\tau )}\phantom{\rule{4.25pt}{0ex}}{S}_{1}\to \infty \phantom{\rule{4.25pt}{0ex}}\mathrm{or}\phantom{\rule{4.25pt}{0ex}}{S}_{2}\to \infty ,$$

- 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 is strike price of the underlying asset.

During the past 10–20 years, fractional calculus is extensively used for real-life problems in a variety of fields like chemistry, biology, physics, earth science and engineering [2,3,4,5,6,7,8,9,10]. Most of the fractional derivatives are represented in the integral form. Two fractional derivatives, which are the most popular, are Riemann-Liouville derivative [11] and Caputo derivative [12]. Katugampola [13,14] defined a generalized fractional derivative including Riemann-Liouville and Hadamard fractional integrals and derivatives. The Katugampola fractional derivative in Caputo type is defined by Almeida et al. [15] and Jarad et al. [16] that some researchers call this fractional derivative type the Caputo generalized fractional derivative.

The Black-Scholes equation with integer-order derivatives has been developed and studied by many researchers in [17,18,19,20]. Those studies show that the financial market has fractal characters both at home and abroad which is not enough to reflect the reality of the financial market. Recently, fractional differential equations are proved and accepted that fractional differential equations are a powerful tool in studying fractal geometry and fractal dynamics. The concept of the fractional derivative model helps to understand the data structure for the real-life matters deeper than the integer-order derivative model. This is the reason why researchers attempted to develop the integer-order Black-Scholes equation into the fractional-order Black-Scholes equation.

There are many analytical methods, used to find solutions of fractional Black-Scholes equations such as homotopy analysis method, homotopy perturbation method [21], modified homotopy perturbation method [22], Adomain decomposition method [23], analysis of decomposition method [24], finite difference method [25,26,27], variational iteration method with coupling of the Sumudu transform [28], Laplace homotopy perturbation method [29,30,31,32], and $\rho $-Laplace homotopy perturbation method [33,34].

In 2017, K. Trachoo, W. Sawangtong and P. Sawangtong [32] studied the two dimensional Black-Scholes equation with the integer-order derivative. They used the Laplace homotopy perturbation method to solve the such equation. The analytic solution of this equation is carried out in the form of the special function, which is Mellin-Ross function.

In 2018, P. Sawangtong, K. Trachoo, W. Sawangtong and B. Wiwattanapataphee [31] solved the Black-Scholes equation with two assets in the Liouville-Caputo fractional derivative sense. By the Laplace transform homotopy perturbation method, the analytic solution obtained is in the form of the Generalized Mittag-Leffler function. Moreover, their presented results ensure that the fractional Black-Scholes model is more practical than the integer-order model.

As in previous descriptions, the characteristic of the finance market is fractal geometry and fractal dynamics, and the fractional differential equation is the powerful tool in studying such characteristics. Moreover, the self-similar property of fractional derivative makes can reply to the long-range dependence better than the integer-order derivative. The motivation of this work comes from the two above papers [31,32] and the advantages of fractional calculus. This article studies the fractional-order Black-Scholes equation with two assets described in the sense of the Katugampola Fractional Derivative in Caputo type. We obtain here the analytic solution of the time-fractional Black-Scholes equation with two assets by the technique of the generalized Laplace homotopy perturbation method. Such the method is the combination between homotopy perturbation method and generalized Laplace transform proposed by Jarad and Abdeljawad [35]. Furthermore, the Katugampola fractional derivative in Caputo type can be transformed to both the Caputo derivative and Caputo-Hadamard derivative [36]. In the particular case, analytic solutions of the time-fractional Black-Scholes equation with two assets in the sense of both Caputo derivative and Caputo-Hadamard derivative are shown. We moreover investigate that the fractional-order of the Caputo type Katugampola fractional derivative has an effect on the price of a European call option.

The paper is organized as follows. We introduce the definitions of fractional derivatives, the generalized Laplace transform and the generalized Mittag-Leffler function in Section 2. We recall the constructive partial differential equation associate with the Black-Scholes models in Section 3. In Section 4, the time fractional Black-Scholes equation with two assets of the European call option is determined and carried out by the generalized Laplace homotopy perturbation method. The analytic solutions of fractional Black-Scholes equations and numerical results are presented in Section 5 and Section 6, respectively. Finally, conclusions and a discussion of current work is given in Section 7.

In this section, we give definitions of fractional derivatives and the generalized Laplace transform. Furthermore, some properties of generalized Laplace transform and two-parameters Mittag-Leffler function used throughout this work are shown.

Let $\alpha ,a,$ and $\rho $ be any real number with $0<\alpha \le 1,a\ge 0,\rho >0$ and $\mathsf{\Gamma}$ be a Gamma function.

([37]). The Riemann-Liouville fractional integral of f with order α has the following form:

$$\left({}_{a}{I}_{t}^{\alpha}f\right)\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma}\left(\alpha \right)}}{\int}_{a}^{t}{(t-u)}^{\alpha -1}f\left(u\right)du.$$

([37]). The Riemann-Liouville fractional derivative of f with order α has the following form:

$$\left({}_{a}{D}_{t}^{\alpha}f\right)\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma}(1-\alpha )}}{\displaystyle \frac{d}{dt}}{\int}_{a}^{t}{(t-u)}^{-\alpha}f\left(u\right)du.$$

([37]). The Caputo fractional derivative of f with order α has the following form:

$$\left({}_{a}^{C}{D}_{t}^{\alpha}f\right)\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma}(1-\alpha )}}{\int}_{a}^{t}{(t-u)}^{-\alpha}{f}^{\prime}\left(u\right)du.$$

([38,39,40]). The Caputo-Hadamard fractional derivative of f with order α has the following form:
where the differential operator δ is defined by $\delta =t{\displaystyle \frac{d}{dt}}$.

$$\left({}_{a}^{CH}{D}_{t}^{\alpha}f\right)\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma}(1-\alpha )}}{\int}_{a}^{t}{\left(log{\displaystyle \frac{t}{u}}\right)}^{-\alpha}\delta f\left(u\right)du,$$

([37]). The Katugampola fractional derivative in Caputo type of f with order α has the following form:
where the differential operator γ is given by $\gamma ={t}^{1-\rho}{\displaystyle \frac{d}{dt}}$.

$$\left({}_{a}^{KC}{D}_{t}^{\alpha ,\rho}f\right)\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma}(1-\alpha )}}{\int}_{a}^{t}{\left({\displaystyle \frac{{t}^{\rho}-{u}^{\rho}}{\rho}}\right)}^{-\alpha}\gamma f\left(u\right){\displaystyle \frac{du}{{u}^{1-\rho}}},$$

Note that if $\rho =1,$ then the fractional derivative as Equation (4) becomes the Caputo fractional derivative with order $\alpha $. On the other hand, if $\rho $ approaches $0,$ then the fractional derivative as Equation (4) results to the Caputo-Hadamard fractional derivative with order $\alpha $.

Next, we present the definition of the generalized Laplace transform.

([35]). Let $f,g:[a,\infty )\to \mathbb{R}$ be real valued functions such that $g\left(t\right)$ is continuous and ${g}^{\prime}\left(t\right)>0$ on $[a,\infty )$. If the generalized Laplace transform of f exists, then
where s is the generalized Laplace transform parameter.

$${\mathcal{L}}_{g}\left\{f\left(t\right)\right\}\left(s\right)={\int}_{a}^{\infty}{e}^{-s\left(g\right(t)-g(a\left)\right)}f\left(t\right){g}^{\prime}\left(t\right)dt,$$

Obverse that if we set $g\left(t\right)=t$ and $a=0$ in Equation (5), then the generalized Laplace transform becomes the classical Laplace transform but if we set $g\left(t\right)={\displaystyle \frac{{t}^{\rho}}{\rho}}$ and $a=0,$ then the generalized Laplace transform reduces to the $\rho -$Laplace transform studied by Fall et al. [33].

Throughout this work, we call the generalized Laplace transform with $g\left(t\right)={\displaystyle \frac{{t}^{\rho}}{\rho}}$ and $a=0$ by $\frac{{t}^{\rho}}{\rho}$-Laplace transform. Then, $\frac{{t}^{\rho}}{\rho}$-Laplace transform is defined by:

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

We next give some properties of $\frac{{t}^{\rho}}{\rho}$-Laplace transform.

([37]).

- 1.
- $${\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\left\{{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{\nu}\right\}\left(s\right)={\displaystyle \frac{\mathsf{\Gamma}(1+\nu )}{{s}^{1+\nu}}},\phantom{\rule{4.25pt}{0ex}}\nu \in \mathbb{R},$$
- 2.
- $${\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\left\{\left({}^{KC}{D}_{t}^{\alpha ,\rho}f\right)\left(t\right)\right\}\left(s\right)={s}^{\alpha}{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\left\{f\left(t\right)\right\}\left(s\right)-{s}^{\alpha -1}f\left(0\right).$$

We give some properties of the two-parameters Mittag-Leffler function.

([43]). The asymptotic expansion of ${E}_{\alpha ,\beta}\left({x}^{\alpha}\right)$ with $\alpha \in (0,1]$ and $\beta \in \mathbb{R}$ is:

$${E}_{\alpha ,\beta}\left({x}^{\alpha}\right)\sim {\displaystyle \frac{{x}^{1-\beta}{e}^{x}}{\alpha}}\phantom{\rule{4.25pt}{0ex}}\mathrm{as}\phantom{\rule{4.pt}{0ex}}x\to \infty .$$

In this section, we let:
and

$$x=\mathrm{ln}\left({S}_{1}\right)-\left(r-{\displaystyle \frac{1}{2}}{\sigma}_{1}^{2}\right)\tau ,\phantom{\rule{4.25pt}{0ex}}y=\mathrm{ln}\left({S}_{2}\right)-\left(r-{\displaystyle \frac{1}{2}}{\sigma}_{2}^{2}\right)\tau ,\phantom{\rule{4.25pt}{0ex}}t=T-\tau ,$$

$$u(\tau ,x,y)={e}^{-r(T-\tau )}v(\tau ,x,y).$$

For more details of transformation, readers can be found in [31,44]. Finally, the Black–Scholes partial differential equation with two assets of the European call option Equation (1) can be transformed in the following form:

$${v}_{t}={\displaystyle \frac{1}{2}}{\sigma}_{1}^{2}{\displaystyle \frac{{\partial}^{2}v}{\partial {x}^{2}}}+{\displaystyle \frac{1}{2}}{\sigma}_{2}^{2}{\displaystyle \frac{{\partial}^{2}v}{\partial {y}^{2}}}+\omega {\sigma}_{1}{\sigma}_{2}{\displaystyle \frac{{\partial}^{2}v}{\partial x\partial y}},\phantom{\rule{4.25pt}{0ex}}(t,x,y)\in [0,T]\times \mathbb{R}\times \mathbb{R}.$$

In this work, we consider the Black–Scholes equation in Equation (8) in the general form by replacing the integer-order time derivative with the time-fractional derivative. The time-fractional Black–Scholes equation is studied based on the Katugampola fractional derivative in Caputo type with $\alpha \in (0,1]$ given by:
with the initial condition:
and the boundary conditions:
where ${c}_{1}={\beta}_{1}{e}^{(r-\frac{1}{2}{\sigma}_{1}^{2})T}$ and ${c}_{2}={\beta}_{2}{e}^{(r-\frac{1}{2}{\sigma}_{2}^{2})T}$.

$${}^{KC}{D}_{t}^{\alpha ,\rho}v={\displaystyle \frac{1}{2}}{\sigma}_{1}^{2}{\displaystyle \frac{{\partial}^{2}v}{\partial {x}^{2}}}+{\displaystyle \frac{1}{2}}{\sigma}_{2}^{2}{\displaystyle \frac{{\partial}^{2}v}{\partial {y}^{2}}}+\omega {\sigma}_{1}{\sigma}_{2}{\displaystyle \frac{{\partial}^{2}v}{\partial x\partial y}},$$

$$v(0,x,y;\rho ,\alpha )=\mathrm{max}\left({c}_{1}{e}^{x}+{c}_{2}{e}^{y}-K,0\right),$$

$$v(t,x,y;\rho ,\alpha )=0\phantom{\rule{4.25pt}{0ex}}\mathrm{as}\phantom{\rule{4.25pt}{0ex}}(x,y)\to -\infty ,\phantom{\rule{4.25pt}{0ex}}\mathrm{and}$$

$$v(t,x,y;\rho ,\alpha )={c}_{1}{e}^{x+\frac{1}{2}{\sigma}_{1}^{2}t}+{c}_{2}{e}^{y+\frac{1}{2}{\sigma}_{2}^{2}t}-K\phantom{\rule{4.25pt}{0ex}}\mathrm{as}\phantom{\rule{4.25pt}{0ex}}x\to \infty \phantom{\rule{4.25pt}{0ex}}\mathrm{or}\phantom{\rule{4.25pt}{0ex}}y\to \infty ,$$

In this section, we present the basic idea of the generalized Laplace homotopy perturbation method with $g\left(t\right)={\displaystyle \frac{{t}^{\rho}}{\rho}}$ and $a=0$ for the time-fractional differential equations.

Firstly, we consider the following problem:
with the initial condition:
where R and N are the linear and nonlinear operators and f and h are determined functions.

$${}^{KC}{D}_{t}^{\alpha ,\rho}u(t,x,y;\rho ,\alpha )+Ru(t,x,y;\rho ,\alpha )+Nu(t,x,y;\rho ,\alpha )=f(t,x,y)\phantom{\rule{4.25pt}{0ex}}\mathrm{for}\phantom{\rule{4.25pt}{0ex}}(t,x,y)\in (0,T]\times \mathbb{R}\times \mathbb{R},$$

$$u(0,x,y;\rho ,\alpha )=h(x,y)\phantom{\rule{4.25pt}{0ex}}\mathrm{for}\phantom{\rule{4.25pt}{0ex}}\mathrm{any}\phantom{\rule{4.25pt}{0ex}}(x,y)\in \mathbb{R}\times \mathbb{R},$$

Step 1. By taking the generalized Laplace transform with respect to time variable t in Equation (12), we obtain:

$${\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\left\{{}^{KC}{D}_{t}^{\alpha ,\rho}u(t,x,y;\rho ,\alpha )\right\}\left(s\right)+{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\{R\left(u(t,x,y;\rho ,\alpha )\right)+N\left(u(t,x,y;\rho ,\alpha )\right)\}\left(s\right)={\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\left\{f(t,x,y)\right\}\left(s\right).$$

Step 2. By using the generalized Laplace transform of Katugampola fractional derivative in Caputo type, we have:

$${\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\left\{u\right\}\left(s\right)={s}^{-1}h(x,y)-{s}^{-\alpha}{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\{R\left(u\right)+N\left(u\right)\}\left(s\right)+{s}^{-\alpha}{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\left\{f(t,x,y)\right\}\left(s\right).$$

Step 3. By taking the inverse generalized Laplace transform on the both sides in Equation (14), we get:
where

$$u(t,x,y;\rho ,\alpha )=G(t,x,y)-{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}^{-1}\left\{{s}^{-\alpha}{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\{R\left(u\right)+N\left(u\right)\}\left(s\right)\right\}\left(t\right),$$

$$G(t,x,y)={\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}^{-1}\left\{{s}^{-1}h(x,y)+{s}^{-\alpha}{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\left\{f(t,x,y)\right\}\right\}\left(t\right).$$

Step 4. By applying the homotopy perturbation technique [45,46], the homotopy function $\nu (t,x,y;p)$ is defined by the homotopy equation:
where $p\in [0,1]$ denote homotopy parameter and ${\tilde{\nu}}_{0}(t,x,y)$ is an arbitrary initial approximation satisfying both the initial and boundary conditions of the problem [47].

$$\nu (t,x,y;p)={\tilde{\nu}}_{0}(t,x,y)-p\left[{\tilde{\nu}}_{0}(t,x,y)-G(t,x,y)+{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}^{-1}\{{s}^{-\alpha}{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\{R\left(\nu (t,x,y;p)\right)+N\left(\nu (t,x,y;p)\right)\}\}\right]$$

Note that we can see that when $p=1$ in Equation (15), the homotopy function $\nu (t,x,y;1)$ is the solution of the problem Equation (12). We assume that:
and the nonlinear term can be applied as:

$$\nu (t,x,y;p)=\sum _{n=0}^{\infty}{p}^{n}{\nu}_{n}(t,x,y)$$

$$N\nu (t,x,y;p)=\sum _{n=0}^{\infty}{p}^{n}{H}_{n}({\nu}_{0},{\nu}_{1},\dots ,{\nu}_{n}).$$

The He’s polynomials ${H}_{n}({\nu}_{0},{\nu}_{1},\dots ,{\nu}_{n})$ can be calculated by using the following formula:

$${H}_{n}({\nu}_{0},{\nu}_{1},{\nu}_{2},\dots ,{\nu}_{n})={\displaystyle \frac{1}{n!}}{\displaystyle \frac{{\partial}^{n}}{\partial {p}^{n}}}N\left({\displaystyle \sum _{i=0}^{n}{p}^{i}{\nu}_{i}}\right){|}_{p=0},n=0,1,2,\dots .$$

Step 5. By substituting Equations (16) and (17) into (15), we obtain:

$$\begin{array}{ccc}{\displaystyle \sum _{n=0}^{\infty}{p}^{n}{\nu}_{n}(t,x,y)}\hfill & =& {\tilde{\nu}}_{0}(t,x,y)-p[{\tilde{\nu}}_{0}(t,x,y)-G(t,x,y)\hfill \\ & +& {\displaystyle {\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}^{-1}\{{s}^{-\alpha}{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\{R(\sum _{n=0}^{\infty}{p}^{n}{\nu}_{n}(t,x,y))+\sum _{n=0}^{\infty}{p}^{n}{H}_{n}({\nu}_{0},{\nu}_{1},\dots ,{\nu}_{n})\}\left(s\right)\}\left(t\right)].}\hfill \end{array}$$

Step 6. By comparing the coefficient of like power of p on both sides of Equation (18), the sequence ${\nu}_{n}$ is carried out:

$$\begin{array}{ccc}{p}^{0}\phantom{\rule{8.5pt}{0ex}}:\phantom{\rule{8.5pt}{0ex}}{\nu}_{0}(t,x,y)\hfill & =& {\tilde{\nu}}_{0}(t,x,y),\hfill \\ {p}^{1}\phantom{\rule{8.5pt}{0ex}}:\phantom{\rule{8.5pt}{0ex}}{\nu}_{1}(t,x,y)\hfill & =& G(t,x,y)-{\tilde{\nu}}_{0}(t,x,y)-{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}^{-1}\{{s}^{-\alpha}{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\{R\left({\tilde{\nu}}_{0}(t,x,y)\right)+{H}_{0}(t,x,y))\}\left(s\right)\}\left(t\right),\hfill \\ {p}^{n}\phantom{\rule{8.5pt}{0ex}}:\phantom{\rule{8.5pt}{0ex}}{\nu}_{n}(t,x,y)\hfill & =& -{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}^{-1}\{{s}^{-\alpha}{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\{R\left({\nu}_{n-1}(t,x,y)\right)+{H}_{n-1}(t,x,y))\}\left(s\right)\}\left(t\right)\phantom{\rule{4.25pt}{0ex}}\mathrm{when}\phantom{\rule{4.25pt}{0ex}}n\ge 2.\hfill \end{array}$$

In this section, we solve the time-fractional Black–Scholes equation Equation (9) subjecting to the initial condition Equation (10) and the boundary conditions Equation (11) by the generalized Laplace homotopy perturbation method.

From Equations (9) and (12), we can see that:
and $Nv(t,x,y;\rho ,\alpha )=0,f(t,x,y)=0,$ with the initial condition as:

$$R\left(v(t,x,y;\rho ,\alpha )\right)={\displaystyle \frac{1}{2}}{\sigma}_{1}^{2}{\displaystyle \frac{{\partial}^{2}v}{\partial {x}^{2}}}+{\displaystyle \frac{1}{2}}{\sigma}_{2}^{2}{\displaystyle \frac{{\partial}^{2}v}{\partial {y}^{2}}}+\omega {\sigma}_{1}{\sigma}_{2}{\displaystyle \frac{{\partial}^{2}v}{\partial x\partial y}},$$

$$v(0,x,y;\rho ,\alpha )=h(x,y)=\mathrm{max}({c}_{1}{e}^{x}+{c}_{2}{e}^{y}-K,0).$$

Note that since the nonlinear term is zero, the He’s polynomials ${H}_{n}=0$ for all $n=1,2,3,\dots $.

From the recurrence relation Equation (19), we choose the function ${\tilde{v}}_{0}$ by:

$${\tilde{v}}_{0}(t,x,y)=\mathrm{max}({c}_{1}{e}^{x}+{c}_{2}{e}^{y}-K,0).$$

We next find the function ${v}_{1}$ by using Equation (19). Then,

$$\begin{array}{ccc}{v}_{1}(t,x,y)\hfill & =& G(t,x,y)-{\tilde{v}}_{0}(t,x,y)-{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}^{-1}\{{s}^{-\alpha}{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\{R\left({\tilde{v}}_{0}(t,x,y)\right)+{H}_{0}(t,x,y))\}\left(s\right)\}\left(t\right)\hfill \\ & =& {\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}^{-1}\left\{{\displaystyle \frac{1}{{s}^{\alpha}}}{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\left\{{\displaystyle \frac{1}{2}}{\sigma}_{1}^{2}{\displaystyle \frac{{\partial}^{2}{v}_{0}}{\partial {x}^{2}}}+{\displaystyle \frac{1}{2}}{\sigma}_{2}^{2}{\displaystyle \frac{{\partial}^{2}{v}_{0}}{\partial {y}^{2}}}+\omega {\sigma}_{1}{\sigma}_{2}{\displaystyle \frac{{\partial}^{2}{v}_{0}}{\partial x\partial y}}\right\}\left(s\right)\right\}\left(t\right)\hfill \\ & =& {\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}^{-1}\left\{{\displaystyle \frac{1}{{s}^{\alpha +1}}}\left({\displaystyle \frac{1}{2}}{\sigma}_{1}^{2}\mathrm{max}({c}_{1}{e}^{x},0)+{\displaystyle \frac{1}{2}}{\sigma}_{2}^{2}\mathrm{max}({c}_{2}{e}^{y},0)\right)\right\}\left(t\right)\hfill \\ & =& \left({\displaystyle \frac{1}{2}}{\sigma}_{1}^{2}\mathrm{max}({c}_{1}{e}^{x},0)+{\displaystyle \frac{1}{2}}{\sigma}_{2}^{2}\mathrm{max}({c}_{2}{e}^{y},0)\right){\displaystyle \frac{1}{\mathsf{\Gamma}(1+\alpha )}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{\alpha}.\hfill \end{array}$$

Therefore,

$${v}_{1}(t,x,y)=\left({\displaystyle \frac{1}{2}}{\sigma}_{1}^{2}\mathrm{max}({c}_{1}{e}^{x},0)+{\displaystyle \frac{1}{2}}{\sigma}_{2}^{2}\mathrm{max}({c}_{2}{e}^{y},0)\right){\displaystyle \frac{1}{\mathsf{\Gamma}(1+\alpha )}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{\alpha}.$$

By Equation (19), we obtain that:

$$\begin{array}{ccc}{v}_{2}(t,x,y)\hfill & =& -{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}^{-1}\{{s}^{-\alpha}{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\{R\left({v}_{1}(t,x,y)\right)+{H}_{1}(t,x,y))\}\left(s\right)\}\left(t\right)\hfill \\ & =& {\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}^{-1}\left\{{\displaystyle \frac{1}{{s}^{\alpha}}}{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\left\{{\displaystyle \frac{1}{2}}{\sigma}_{1}^{2}{\displaystyle \frac{{\partial}^{2}{v}_{1}}{\partial {x}^{2}}}+{\displaystyle \frac{1}{2}}{\sigma}_{2}^{2}{\displaystyle \frac{{\partial}^{2}{v}_{1}}{\partial {y}^{2}}}+\omega {\sigma}_{1}{\sigma}_{2}{\displaystyle \frac{{\partial}^{2}{v}_{1}}{\partial x\partial y}}\right\}\right\}\left(t\right)\hfill \\ & =& \left({\displaystyle \frac{1}{{2}^{2}}}{\sigma}_{1}^{4}\mathrm{max}({c}_{1}{e}^{x},0)+{\displaystyle \frac{1}{{2}^{2}}}{\sigma}_{2}^{4}\mathrm{max}({c}_{2}{e}^{y},0)\right){\displaystyle \frac{1}{\mathsf{\Gamma}(1+2\alpha )}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{2\alpha}.\hfill \end{array}$$

Hence,

$${v}_{2}(t,x,y)=\left({\displaystyle \frac{1}{{2}^{2}}}{\sigma}_{1}^{4}\mathrm{max}({c}_{1}{e}^{x},0)+{\displaystyle \frac{1}{{2}^{2}}}{\sigma}_{2}^{4}\mathrm{max}({c}_{2}{e}^{y},0)\right){\displaystyle \frac{1}{\mathsf{\Gamma}(1+2\alpha )}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{2\alpha}.$$

It follows from Equation (19) that:

$$\begin{array}{ccc}{v}_{3}(t,x,y)\hfill & =& -{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}^{-1}\left\{{\displaystyle \frac{1}{{s}^{\alpha}}}{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\left\{{\displaystyle \frac{1}{2}}{\sigma}_{1}^{2}{\displaystyle \frac{{\partial}^{2}{v}_{2}}{\partial {x}^{2}}}+{\displaystyle \frac{1}{2}}{\sigma}_{2}^{2}{\displaystyle \frac{{\partial}^{2}{v}_{2}}{\partial {y}^{2}}}+\omega {\sigma}_{1}{\sigma}_{2}{\displaystyle \frac{{\partial}^{2}{v}_{2}}{\partial x\partial y}}\right\}\left(s\right)\right\}\left(t\right)\hfill \\ & =& \left({\displaystyle \frac{1}{{2}^{3}}}{\sigma}_{1}^{6}\mathrm{max}({c}_{1}{e}^{x},0)+{\displaystyle \frac{1}{{2}^{3}}}{\sigma}_{2}^{6}\mathrm{max}({c}_{2}{e}^{y},0)\right){\displaystyle \frac{1}{\mathsf{\Gamma}(1+3\alpha )}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{3\alpha}.\hfill \end{array}$$

Then,

$${v}_{3}(t,x,y)=\left({\displaystyle \frac{1}{{2}^{3}}}{\sigma}_{1}^{6}\mathrm{max}({c}_{1}{e}^{x},0)+{\displaystyle \frac{1}{{2}^{3}}}{\sigma}_{2}^{6}\mathrm{max}({c}_{2}{e}^{y},0)\right){\displaystyle \frac{1}{\mathsf{\Gamma}(1+3\alpha )}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{3\alpha}.$$

Like the same manner, we can conclude that by Equation (19),
for $n=1,2,3,\dots $. Therefore, by Equation (15), the homotopy perturbation $v(t,x,y;p)$ which corresponds to problem Equations (12) and (13) is in the form:

$$\begin{array}{ccc}{v}_{n}(t,x,y)\hfill & =& -{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}^{-1}\left\{{\displaystyle \frac{1}{{s}^{\alpha}}}{\mathcal{L}}_{\frac{{t}^{\rho}}{\rho}}\left\{{\displaystyle \frac{1}{2}}{\sigma}_{1}^{2}{\displaystyle \frac{{\partial}^{2}{v}_{n-1}}{\partial {x}^{2}}}+{\displaystyle \frac{1}{2}}{\sigma}_{2}^{2}{\displaystyle \frac{{\partial}^{2}{v}_{n-1}}{\partial {y}^{2}}}+\omega {\sigma}_{1}{\sigma}_{2}{\displaystyle \frac{{\partial}^{2}{v}_{n-1}}{\partial x\partial y}}\right\}\left(s\right)\right\}\left(t\right)\hfill \\ \\ & =& \left({\displaystyle \frac{1}{{2}^{n}}}{\sigma}_{1}^{2n}\mathrm{max}({c}_{1}{e}^{x},0)+{\displaystyle \frac{1}{{2}^{n}}}{\sigma}_{2}^{2n}\mathrm{max}({c}_{2}{e}^{y},0)\right){\displaystyle \frac{1}{\mathsf{\Gamma}(1+n\alpha )}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{n\alpha}.\hfill \end{array}$$

$$\begin{array}{ccc}v(t,x,y;p)\hfill & =& \mathrm{max}\left({c}_{1}{e}^{x}+{c}_{2}{e}^{y}-K,0\right)\hfill \\ & +& {\displaystyle \sum _{n=0}^{\infty}{p}^{n+1}\left\{{\displaystyle \frac{1}{\mathsf{\Gamma}(1+(n+1\left)\alpha \right)}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{(n+1)\alpha}\left({\displaystyle \frac{1}{{2}^{(n+1)}}}{\sigma}_{1}^{2(n+1)}\mathrm{max}\left({c}_{1}{e}^{x},0\right)\right)\right\}}\hfill \\ & +& {\displaystyle \sum _{n=0}^{\infty}{p}^{n+1}\left\{{\displaystyle \frac{1}{\mathsf{\Gamma}(1+(n+1\left)\alpha \right)}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{(n+1)\alpha}\left({\displaystyle \frac{1}{{2}^{(n+1)}}}{\sigma}_{2}^{2(n+1)}\mathrm{max}\left({c}_{2}{e}^{y},0\right)\right)\right\}.}\hfill \end{array}$$

By Equation (20), we obtain the analytic solution of problem Equations (9)–(11) defined by:
or
where ${E}_{\alpha ,\alpha +1}$ is the two-parameters Mittag–Leffler function defined by Equation (6).

$$\begin{array}{ccc}v(t,x,y;\rho ,\alpha )=v(t,x,y;1)\hfill & =& \mathrm{max}\left({c}_{1}{e}^{x}+{c}_{2}{e}^{y}-K,0\right)\hfill \\ & +& {\displaystyle \sum _{n=0}^{\infty}\left\{{\displaystyle \frac{1}{\mathsf{\Gamma}(1+(n+1\left)\alpha \right)}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{(n+1)\alpha}\left({\displaystyle \frac{1}{{2}^{(n+1)}}}{\sigma}_{1}^{2(n+1)}\mathrm{max}\left({c}_{1}{e}^{x},0\right)\right)\right\}}\hfill \\ & +& {\displaystyle \sum _{n=0}^{\infty}\left\{{\displaystyle \frac{1}{\mathsf{\Gamma}(1+(n+1\left)\alpha \right)}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{(n+1)\alpha}\left({\displaystyle \frac{1}{{2}^{(n+1)}}}{\sigma}_{2}^{2(n+1)}\mathrm{max}\left({c}_{2}{e}^{y},0\right)\right)\right\}}\hfill \\ & =& \mathrm{max}\left({c}_{1}{e}^{x}+{c}_{2}{e}^{y}-K,0\right)\hfill \\ & +& {\displaystyle \left({\displaystyle \frac{{\sigma}_{1}^{2}}{2}}\right){\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{\alpha}\mathrm{max}({c}_{1}{e}^{x},0)\sum _{n=0}^{\infty}{\displaystyle \frac{1}{\mathsf{\Gamma}(\alpha +1+n\alpha )}}{\left({\displaystyle \frac{{\sigma}_{1}^{2}}{2}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{\alpha}\right)}^{n}}\hfill \\ & +& {\displaystyle \left({\displaystyle \frac{{\sigma}_{2}^{2}}{2}}\right){\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{\alpha}\mathrm{max}({c}_{2}{e}^{y},0)\sum _{n=0}^{\infty}{\displaystyle \frac{1}{\mathsf{\Gamma}(\alpha +1+n\alpha )}}{\left({\displaystyle \frac{{\sigma}_{2}^{2}}{2}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{\alpha}\right)}^{n},}\hfill \end{array}$$

$$\begin{array}{ccc}v(t,x,y;\rho ,\alpha )\hfill & =& \mathrm{max}\left({c}_{1}{e}^{x}+{c}_{2}{e}^{y}-K,0\right)\hfill \\ & +& \left({\displaystyle \frac{{\sigma}_{1}^{2}}{2}}\right){\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{\alpha}\mathrm{max}\left({c}_{1}{e}^{x},0\right){E}_{\alpha ,\alpha +1}\left({\displaystyle \frac{{\sigma}_{1}^{2}}{2}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{\alpha}\right)\hfill \\ & +& \left({\displaystyle \frac{{\sigma}_{2}^{2}}{2}}\right){\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{\alpha}\mathrm{max}\left({c}_{2}{e}^{y},0\right){E}_{\alpha ,\alpha +1}\left({\displaystyle \frac{{\sigma}_{2}^{2}}{2}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{\alpha}\right),\hfill \end{array}$$

Hence, the analytic solution of the time fractional Black-Scholes equations in the sense of Caputo type Katugampola fractional derivative Equations (9)–(11) is given by Equation (21).

It follows from [36] that as $\rho =1$, the time-fractional Black–Scholes equation with two assets for the European call option based on Caputo type Katugampola fractional derivative becomes the time-fractional Black–Scholes equation with two assets for the European call option based on Caputo fractional derivative:
with the initial condition:
and the boundary conditions:
where ${c}_{1}$ and ${c}_{2}$ are constants.

$${}^{C}{D}_{t}^{\alpha}v={\displaystyle \frac{1}{2}}{\sigma}_{1}^{2}{\displaystyle \frac{{\partial}^{2}v}{\partial {x}^{2}}}+{\displaystyle \frac{1}{2}}{\sigma}_{2}^{2}{\displaystyle \frac{{\partial}^{2}v}{\partial {y}^{2}}}+\omega {\sigma}_{1}{\sigma}_{2}{\displaystyle \frac{{\partial}^{2}v}{\partial x\partial y}},$$

$$v(0,x,y;1,\alpha )=\mathrm{max}({c}_{1}{e}^{x}+{c}_{2}{e}^{y}-K,0),$$

$$v(t,x,y;1,\alpha )=0\phantom{\rule{4.25pt}{0ex}}\mathrm{as}\phantom{\rule{4.25pt}{0ex}}(x,y)\to -\infty ,\phantom{\rule{4.25pt}{0ex}}\mathrm{and}$$

$$v(t,x,y;1,\alpha )={c}_{1}{e}^{x+\frac{1}{2}{\sigma}_{1}^{2}t}+{c}_{2}{e}^{y+\frac{1}{2}{\sigma}_{2}^{2}t}-K\phantom{\rule{4.25pt}{0ex}}\mathrm{as}\phantom{\rule{4.25pt}{0ex}}x\to \infty \phantom{\rule{4.25pt}{0ex}}\mathrm{or}\phantom{\rule{4.25pt}{0ex}}y\to \infty ,$$

By (21), the analytic solution of the time-fractional Black–Scholes equation with two assets for the European call option in Caputo fractional derivative sense is defined:
which differs from the solution studied in [31]. The reason is why we obtain another solution because choosing the initial function ${\tilde{\nu}}_{0}(t,x,y)$ in recurrence relation (19) is different. In the particular case, if we let $\rho =1$ and $\alpha =1,$ then the time-fractional Black–Scholes equation reduces to the classical Black–Scholes equation with two assets as Equation (8). It follows from the property of Mittag–Leffler function with two parpmeters, that is ${E}_{1,2}\left(z\right)=\left({\displaystyle \frac{{e}^{z}-1}{z}}\right)$, that the analytic solution of the classical Black–Scholes equation with two assets for the European call option is given by:
which differs from the solution obtained by [32] because choosing the initial function ${\tilde{\nu}}_{0}(t,x,y)$ in recurrence relation (19) is different.

$$\begin{array}{ccc}\hfill v(t,x,y;1,\alpha )& =& \mathrm{max}\left({c}_{1}{e}^{x}+{c}_{2}{e}^{y}-K,0\right)\hfill \\ & +& \left({\displaystyle \frac{{\sigma}_{1}^{2}{t}^{\alpha}}{2}}\right)\mathrm{max}\left({c}_{1}{e}^{x},0\right){E}_{\alpha ,\alpha +1}\left({\displaystyle \frac{{\sigma}_{1}^{2}{t}^{\alpha}}{2}}\right)+\left({\displaystyle \frac{{\sigma}_{2}^{2}{t}^{\alpha}}{2}}\right)\mathrm{max}\left({c}_{2}{e}^{y},0\right){E}_{\alpha ,\alpha +1}\left({\displaystyle \frac{{\sigma}_{2}^{2}{t}^{\alpha}}{2}}\right),\hfill \end{array}$$

$$\begin{array}{ccc}v(t,x,y;1,1)\hfill & =& \mathrm{max}\left({c}_{1}{e}^{x}+{c}_{2}{e}^{y}-K,0\right)\hfill \\ & +& \left({\displaystyle \frac{{\sigma}_{1}^{2}t}{2}}\right)\mathrm{max}\left({c}_{1}{e}^{x},0\right){E}_{1,2}\left({\displaystyle \frac{{\sigma}_{1}^{2}t}{2}}\right)+\left({\displaystyle \frac{{\sigma}_{2}^{2}t}{2}}\right)\mathrm{max}\left({c}_{2}{e}^{y},0\right){E}_{1,2}\left({\displaystyle \frac{{\sigma}_{2}^{2}t}{2}}\right)\hfill \\ & =& \mathrm{max}\left({c}_{1}{e}^{x}+{c}_{2}{e}^{y}-K,0\right)\hfill \\ & +& \mathrm{max}\left({c}_{1}{e}^{x},0\right)\left({e}^{\frac{{\sigma}_{1}^{2}t}{2}}-1\right)+\mathrm{max}\left({c}_{2}{e}^{y},0\right)\left({e}^{\frac{{\sigma}_{2}^{2}t}{2}}-1\right),\hfill \end{array}$$

We next consider the case when $\rho $ approaches to ${0}^{+}.$

More precisely, as $\rho $ approaches to ${0}^{+},$ the time-fractional Black-Scholes equation with two assets for the European call option based on Caputo type Katugampola fractional derivative Equations (9)–(11) becomes the time-fractional Black-Scholes equation with two assets for the European call option in sense of the Caputo-Hadamard fractional derivative.

Consider the terms ${E}_{\alpha ,\alpha +1}{\left({\left({\displaystyle \frac{{\sigma}_{i}^{2}}{2}}\right)}^{\frac{1}{\alpha}}{\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{\alpha}$ for $i=1\phantom{\rule{4.pt}{0ex}}\mathrm{and}\phantom{\rule{4.pt}{0ex}}2$.

By using the property in Lemma 2 as Equation (7) we have that for each $i=1\phantom{\rule{4.pt}{0ex}}\mathrm{and}\phantom{\rule{4.pt}{0ex}}2$,

$${E}_{\alpha ,\alpha +1}{\left({\left(\frac{{\sigma}_{i}^{2}}{2}\right)}^{\frac{1}{\alpha}}\frac{{t}^{\rho}}{\rho}\right)}^{\alpha}\sim {\displaystyle \frac{1}{\alpha}}{\left({\left(\frac{{\sigma}_{i}^{2}}{2}\right)}^{\frac{1}{\alpha}}\frac{{t}^{\rho}}{\rho}\right)}^{-\alpha}{e}^{{\left(\frac{{\sigma}_{i}^{2}}{2}\right)}^{\frac{1}{\alpha}}\frac{{t}^{\rho}}{\rho}}\phantom{\rule{4.pt}{0ex}}\mathrm{as}\phantom{\rule{4.pt}{0ex}}\rho \to {0}^{+}.$$

It follows from Equation (21) that the analytic solution of the time-fractional Black–Scholes equation in the Caputo–Hadamard fractional derivative is in the form:

$$\begin{array}{ccc}\hfill \underset{\rho \to {0}^{+}}{lim}v(t,x,y;\rho ,\alpha )& =& \mathrm{max}\left({c}_{1}{e}^{x}+{c}_{2}{e}^{y}-K,0\right)\hfill \\ & +& \underset{\rho \to {0}^{+}}{lim}\left({\displaystyle \frac{{\sigma}_{1}^{2}}{2}}\right){\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{\alpha}\mathrm{max}\left({c}_{1}{e}^{x},0\right){E}_{\alpha ,\alpha +1}\left({\displaystyle \frac{{\sigma}_{1}^{2}}{2}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{\alpha}\right)\hfill \\ & +& \underset{\rho \to {0}^{+}}{lim}\left({\displaystyle \frac{{\sigma}_{2}^{2}}{2}}\right){\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{\alpha}\mathrm{max}\left({c}_{2}{e}^{y},0\right){E}_{\alpha ,\alpha +1}\left({\displaystyle \frac{{\sigma}_{2}^{2}}{2}}{\left({\displaystyle \frac{{t}^{\rho}}{\rho}}\right)}^{\alpha}\right)\hfill \\ & =& \mathrm{max}\left({c}_{1}{e}^{x}+{c}_{2}{e}^{y}-K,0\right)+\frac{1}{\alpha}\mathrm{max}\left({c}_{1}{e}^{x},0\right){e}^{{\left(\frac{{\sigma}_{1}^{2}}{2}\right)}^{\frac{1}{\alpha}}\frac{{t}^{\rho}}{\rho}}\hfill \\ & +& \frac{1}{\alpha}\mathrm{max}\left({c}_{2}{e}^{y},0\right){e}^{{\left(\frac{{\sigma}_{2}^{2}}{2}\right)}^{\frac{1}{\alpha}}\frac{{t}^{\rho}}{\rho}}.\hfill \end{array}$$

We see that ${lim}_{\rho \to {0}^{+}}v(t,x,y;\rho ,\alpha )\to \infty $ for all $(t,x,y)\in [0,T]\times \mathbf{R}\times \mathbf{R}.$ This means that ${lim}_{\rho \to {0}^{+}}v(t,x,y;\rho ,\alpha )$ is unbounded for any $(t,x,y)\in [0,T]\times \mathbf{R}\times \mathbf{R}$ which contradicts to the initial condition: $v(0,x,y;\rho ,\alpha )=\mathrm{max}\left({c}_{1}{e}^{x}+{c}_{2}{e}^{y}-K,0\right).$ This implies that the time-fractional Black–Scholes equation with two assets in Caputo–Hadamard sense has no solution.

In this section, we illustrate the analytic solution of the time-fractional Black–Scholes equation with two assets for European call option based on Caputo type Katugampola fractional derivative by using Python programming. The financial parameters are given in Table 1.

The classical values of the options $(\alpha =\rho =1)$ defined by (22) are simulated in Figure 1. In Figure 2, Figure 3 and Figure 4, we have depicted the values of the options in Equation (21) for $\alpha =1$ and $\rho =0.85$, $\alpha =0.75$ and $\rho =1$, and $\alpha =1$ and $\rho =1.15,$ respectively. We note that the analytic solution of the time-fractional Black-Scholes equation described in the sense of Caputo type Katugampola fractional derivative and the analytic solution of the classical Black-Scholes equation with both of the two assets are in good agreement. The European call option prices for different values of $\alpha $ and $\rho $ are shown in Table 2.

We see that in Figure 5 and Figure 6, the values of x and y in the case $\rho <1$ are greater than the values of x and y in the case $\rho =1.$ This means that the order $\rho $ of the Caputo type Katugampola fractional derivative has an acceleration effect in the diffusion process when $\rho <1.$ On the other hand, in Figure 7 and Figure 8, the values of x and y in the case $\rho >1$ are less than the values of x and y in the case $\rho =1.$ This implies that the order $\rho $ of the Caputo type Katugampola fractional derivative has a retardation effect in the diffusion process when $\rho >1$.

In this paper, we have discussed the analytic solution of the time-fractional Black–Scholes equation with two assets in sense of the Katugampola fractional derivative in Caputo type. The method, used here to provide analytical solutions, is called the generalized Laplace homotopy perturbation method. The idea of the generalized Laplace homotopy perturbation method comes from combining the technique of both the homotopy perturbation method and generalized Laplace transform. The result in Equation (21) shows that the analytic solution of the time-fractional Black–Scholes equations in the sense of Caputo type Katugampola fractional derivative can be carried out in the form of the two-parameters Mittag–Leffler function. Since the Katugampola fractional derivative in Caputo type generalizes both Caputo and Caputo–Hadamard derivatives, in the case of the time-fractional Black–Scholes equations in the Caputo derivative sense, the analytic solution is also obtained. Unfortunately, in the case of the time-fractional Black–Scholes equations in the Caputo-Hadamard derivative sense, there is no solution. Moreover, we can observe that the orders $\rho $ and $\alpha $ of the Katugampola fractional derivative in Caputo type have an effect on the price of the European option. Moreover, the European option price is overpriced when $\alpha =0.25$.

Conceptualization, P.S.; methodology, P.S. and S.T.; validation, P.S. and S.T.; formal analysis, S.T., W.S. and P.S.; writing—original draft preparation, S.T.; writing—review and editing, W.S. and P.S.; supervision, P.S.; project administration, P.S.; funding acquisition, P.S. All authors have read and agreed to the published version of the manuscript.

This research was funded by King Mongkut’s University of Technology North Bangkok, 133 contract no. KMUTNB-FF-65-58 and this research was partially supported by the Centre of Excellence in 134 Mathematics, PERDO, Commission on Higher Education, Ministry of Education, Thailand.

Not applicable.

Not applicable.

The information used to support the findings of this study are available in the article, in the links there provided or from the corresponding author upon request.

The authors declare no conflict of interest.

- Black, F.; Scholes, M. The pricing of options and corporate liabilities. J. Political Econ.
**1973**, 81, 637–654. [Google Scholar] [CrossRef] - Atangana, A.; Baleanu, D.; Alsaedi, A. Analysis of time-fractional Hunter-Saxton equation: A model of neumatic liquid crystal. Open Phys.
**2016**, 14, 145–149. [Google Scholar] [CrossRef] - Hristov, T. A transient flow of a non-newtonian fluid modelled by a mixed time-space derivative: An improved integral-balance approach. In Mathematical Methods in Engineering; Springer: Cham, Switzerland, 2019; pp. 153–174. [Google Scholar]
- Mainardi, F. Fractional relaxation-oscillation and fractional diffusion-wave phenomena. Chaos Solitons Fractals
**1996**, 9, 1461–1477. [Google Scholar] [CrossRef] - Podlubny, I. Fractional differential equations. In Mathematics in Science and Engineering; Academic Press: San Diego, CA, USA, 1999. [Google Scholar]
- Rekhviashvili, S.; Pskhu, A.; Agarwal, P.; Jain, S. Application of the fractional oscillator model to describe damped vibrations. Turk. J. Phys.
**2019**, 3, 236–242. [Google Scholar] [CrossRef] - Sene, N. Stokes’s first problem for heated flat plate with Atangana-Baleanu fractional derivative. Chaos Solitons Fractals
**2018**, 117, 68–75. [Google Scholar] [CrossRef] - Sofuoglu, Y.; Ozalp, N. Fractional order bilingualism model without conversion from dominant unilingual group to bilingual group. Differ. Equ. Dyn. Syst.
**2017**, 25, 1–9. [Google Scholar] [CrossRef] - Tarasova, V.V.; Vasily, V.E. Elasticity for economic processes with memory: Fractional differential calculus approach. Fract. Diff. Calc.
**2016**, 6, 219–232. [Google Scholar] [CrossRef] - Tariboon, J.; Ntouyas, S.K.; Agarwal, P. New concepts of fractional quantum calculus and applications to impulsive fractional q-difference equations. Adv. Differ. Equ.
**2015**, 18, 1–19. [Google Scholar] [CrossRef] - Kilbas, A.; Srivastava, H.M.; Trujillo, J.J. Theory and Application of Fractional Differential Equations. In North Holland Mathematics Studies; Elsevier: Amsterdam, The Netherlands, 2006; Volume 204. [Google Scholar]
- Samko, S.G.; Kilbas, A.A.; Marichev, O.I. Fractional Integrals and Derivatives: Theory and Applications; Gordon and Breach: Yverdon, Switzerland, 1993. [Google Scholar]
- Katugampola, U.N. New approach to a generalized factional integral. Appl. Math. Comput.
**2011**, 218, 860–865. [Google Scholar] - Katugampola, U.N. A new approach to generalized factional derivatives. Bull. Math. Anal. Appl.
**2014**, 6, 1–15. [Google Scholar] - Almeida, R.; Malinowska, A.B.; Odzijewicz, T. Fractional differential equations with dependence on the Caputo-Katugampola derivative. J. Comput. Nonlinear Dyn.
**2016**, 11, 061017. [Google Scholar] [CrossRef] - Jarad, F.; Abdeljawad, T.; Baleanu, D. On the generalized fractional derivatives and their Caputo modification. J. Nonlinear Sci. Appl.
**2017**, 10, 2607–2619. [Google Scholar] [CrossRef] - Mandelbrot, B. The variation of certain speculative prices. J. Bus.
**1963**, 36, 394–413. [Google Scholar] [CrossRef] - Peters, E.E. Fractal structure in the capital markets. Financ. Anal. J
**1989**, 45, 32–37. [Google Scholar] [CrossRef] - Li, H.Q.; Ma, C.Q. An empirical study of long-term memory of return and volatility in Chinese stock market. J. Financ. Econ.
**2005**, 31, 29–37. [Google Scholar] - Huang, T.F.; Li, B.Y.; Xiong, J.X. Test on the chaotic characteristic of Chinese futures market. Syst. Eng.
**2012**, 30, 43–53. [Google Scholar] - Kumar, S.; Kumar, D.; Singh, J. Numerical computation of fractional Black–Scholes equation arising in financial market. Egypt. J. Basic Appl. Sci.
**2014**, 1, 177–183. [Google Scholar] [CrossRef] - Yavuz, M.; Özdemir, N. A different approach to the European option pricing model with new fractional operator. Math. Model. Nat. Phenom.
**2018**, 13, 12. [Google Scholar] [CrossRef] - Yavuz, M.; Özdemir, N. A quantitative approach to fractional option pricing problems with decomposition series. Konuralp J. Math.
**2014**, 6, 102–109. [Google Scholar] - Ghandehari, M.A.M.; Ranjbar, M. European option pricing of fractional version of the Black–Scholes model: Approach via expansion in series. J. Nonlinear Sci.
**2014**, 17, 105–110. [Google Scholar] - Huang, J.; Cen, Z.; Zhao, J. An adaptive moving mesh method for a time-fractional Black–Scholes equation. Adv. Differ. Equ.
**2019**, 2019, 516. [Google Scholar] [CrossRef] - Koleva, M.N.; Vulkov, L.G. Numerical solution of time-fractional Black–Scholes equation. J. Comput. Appl. Math.
**2017**, 17, 1699–1715. [Google Scholar] [CrossRef] - Song, L.; Wang, W. Solution of the fractional Black-Scholes option pricing model by finite difference method. In Abstract and Applied Analysis; Hindawi: London, UK, 2013. [Google Scholar]
- Ghandehari, M.A.M.; Ranjbar, M. European option pricing of fractional Black-Scholes model with new Lagrange multipliers. Comput. Meth. Differ. Equ.
**2014**, 2, 1–10. [Google Scholar] - Kumar, S.; Yildirim, A.; Khan, Y.; Jafari, H.; Sayevand, K.; Wei, L. Analytical solution of fractional Black-Scholes European option pricing equation by using Laplace transform. Fract. Calc. Appl. Anal.
**2012**, 2, 1–9. [Google Scholar] - Prathumwan, D.; Trachoo, K. On the solution of two-dimensional fractional Black–Scholes equation for European put option. Adv. Differ. Equ.
**2020**, 2020, 146. [Google Scholar] [CrossRef] - Sawangtong, P.; Trachoo, K.; Sawangtong, W.; Wiwattanapataphee, B. The Analytical Solution for the Black-Scholes Equation with Two Assets in the Liouville-Caputo Fractional Derivative Sense. Mathematics
**2018**, 6, 129. [Google Scholar] [CrossRef] - Trachoo, K.; Sawangtong, W.; Sawangtong, P. Laplace transform homotopy perturbation method for the two dimensional Black Scholes model with European call option. Math. Comput. Appl.
**2017**, 22, 23. [Google Scholar] [CrossRef] - Fall, A.N.; Ndiaye, S.N.; Ndolane, S. Black–Scholes option pricing equations described by the Caputo generalized fractional derivative. Chaos Solitons Fractals
**2019**, 125, 108–118. [Google Scholar] [CrossRef] - Ampun, S.; Sawangtong, P. The approximate analytic solution of the time-fractional Black-Scholes equation with a European option based on the Katugampola fractional derivative. Mathematics
**2021**, 9, 214. [Google Scholar] [CrossRef] - Jarad, F.; Abdeljawad, T. Generalized fractional derivatives and Laplace transform. Disc. Cont. Dyn. Syst. S
**2019**, 13. [Google Scholar] [CrossRef] - Oliveira, D.S.; de Oliveira, E.C. On a Caputo-type fractional derivative. Adv. Pure Appl. Math.
**2019**, 10, 81–91. [Google Scholar] [CrossRef] - Fahd, J.; Abdeljawad, T. A modified Laplace transform for certain generalized fractional operators. Res. Nonlinear Anal.
**2018**, 2, 88–98. [Google Scholar] - Fahd, J.; Abdeljawad, T.; Baleanu, D. Caputo-type modification of the Hadamard fractional derivatives. Adv. Differ. Equ.
**2012**, 142, 1–8. [Google Scholar] - Aljoudi, S.; Ahmad, B.; Nieto, J.J.; Alsaedi, A. A coupled system of Hadamard type sequential fractional differential equations with coupled strip conditions. Chaos Solitons Fractal
**2016**, 91, 39–46. [Google Scholar] [CrossRef] - Etemad, S.; Rezapour, S.; Sakar, F.M. On a fractional Caputo–Hadamard problem with boundary value conditions via different orders of the Hadamard fractional operators. Adv. Differ. Equ.
**2020**, 272, 1–20. [Google Scholar] [CrossRef] - Kilbas, A.A.; Saigo, M.; Saxena, R.K. Generalized Mittag-Leffler function and generalized fractional calculus operators. Integral Transforms Spec. Funct.
**2004**, 14, 31–49. [Google Scholar] [CrossRef] - Gorenflo, R.; Kilbas, A.A.; Mainardi, F.; Rogosin, S.V. Mittag-Leffler Functions, Related Topics and Applications; Springer: Berlin, Germany, 2014. [Google Scholar]
- Simon, T. Mittag-Leffler functions and complete monotonicity. Integral Transform. Spec. Funct.
**2015**, 26, 36–50. [Google Scholar] [CrossRef] - Contreras, M.; Llanquihuén, A.; Villena, M. On the Solution of the Multi-Assets Black–Scholes Model: Correlation, Eigenvalues and Geometry. J. Math. Financ.
**2016**, 71, 1772–1783. [Google Scholar] [CrossRef] - He, J.H. Homotopy perturbation technique. J. Comput. Methods Appl. Mech. Eng.
**1999**, 178, 257–262. [Google Scholar] [CrossRef] - He, J.H. Homotopy perturbation method for bifurcation of nonlinear problems. Int. J. Nonlinear Sci. Numer. Simul.
**2005**, 6, 207–208. [Google Scholar] [CrossRef] - Baholian, E.; Azizi, A.; Saeidian, J. Some notes on using the homotopy perturbation method for solving time-dependent differential equations. Math. Comput. Model.
**2009**, 50, 213–224. [Google Scholar] [CrossRef]

Parameters | Value |
---|---|

strike price (dollars), K | 70 |

risk-free interest rate (per year), r | 0.05 |

expiration date (year), T | 1 |

volatility function of underlying first assets (per year), ${\sigma}_{1}$ | 0.1 |

volatility function of underlying second assets (per year), ${\sigma}_{2}$ | 0.2 |

the volatility of ${S}_{1}$ and ${S}_{2}$, $\omega $ | 0.5 |

${\beta}_{1}$ | 2 |

${\beta}_{2}$ | 1 |

$\mathit{\alpha}$ | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
---|---|---|---|---|---|---|---|

$\rho $ | 0.75 | 1 | 0.75 | 1 | 0.95 | 1 | 1.5 |

r | 0.05 | 0.05 | 0.2 | 0.2 | 0.1 | 0.1 | 0.1 |

$\omega $ | 0.5 | 0.5 | 0.5 | 0.5 | 0.5 | 0.5 | 0.5 |

T | 1 | 1 | 1 | 1 | 1 | 1 | 1 |

${\sigma}_{1}$ | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 |

${\sigma}_{2}$ | 0.2 | 0.2 | 0.2 | 0.2 | 0.2 | 0.2 | 0.2 |

${\beta}_{1}$ | 2 | 2 | 2 | 2 | 2 | 2 | 2 |

${\beta}_{2}$ | 1 | 1 | 1 | 1 | 1 | 1 | 1 |

x | 4.1805 | 4.1805 | 4.1805 | 4.1805 | 4.1805 | 4.1805 | 4.1805 |

y | 1.6093 | 1.6093 | 1.6093 | 1.6093 | 1.6093 | 1.6093 | 1.6093 |

v | 72.7786 | 72.4795 | 95.8850 | 95.5376 | 79.8328 | 79.7846 | 79.5009 |

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. |

© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).