The Analytical Solution for the Black-Scholes Equation with Two Assets in the Liouville-Caputo Fractional Derivative Sense

: It is well known that the Black-Scholes model is used to establish the behavior of the option pricing in the ﬁnancial market. In this paper, we propose the modiﬁed version of Black-Scholes model with two assets based on the Liouville-Caputo fractional derivative. The analytical solution of the proposed model is investigated by the Laplace transform homotopy perturbation method.


Introduction
A derivative is one of the financial instruments promising payment at a certain time in the future and the payoff amount depends upon the change of some underlying asset. Its value can be derived from the various sorts of the underlying asset such as shares, bonds, interest rate, commodity, currency, etc. It is unquestioned that options are the primary key of a derivative that is commonly used in the financial market. Hence, the idea of option trading has been continuously developed. The earliest research work was proposed by Corzo et al. [1] using the evidence of the Netherlands having active involvement in trading option. Osborne [2] proposed an option pricing formula using an arithmetic Brownian motion with drift. In the year 1973, Fischer Black and Myron Scholes [3] proposed the Black-Scholes model to investigate the behaviour of the option pricing in a market. Several Mathematical models based on the Black-Scholes equation with five-key components of the strike price, the risk-free rate, the underlying security stock price, the volatility and the mature time have been developed [4][5][6][7].
Several numerical and analytical methods have been studied and developed for finding the solution of Black-Scholes model, for instance, the finite different method [8][9][10][11], finite element method [12] for numerical solutions, and the Mellin transform method [13] and homotopy perturbation method [14,15], Laplace homotopy perturbation method [16] for analytical solutions. Moreover, one of the methods to solve the Black-Scholes equations is radial basis function partition of unity method (RBF-PUM) [17,18] which widely uses to approximate the partial differential equation problem. This method is composed of two methods these are partition of unity (PU) method and radial basis function (RBF) approximants [19,20]. RBF-PUM provides extremely accurate results when applied to data with non-homogeneous density [19].
In general, the Black-Scholes model with 2 assets for option pricing can be written as follows: ∂C ∂τ with the terminal condition: where K = max{K 1 , K 2 }, and the boundary conditions: where: C is the call option depending on the underlying stock prices {S 1 , S 2 } at time τ, q i is the dividend yield on the ith underlying stock, ρ ij is the correlation between the ith and jth underlying stock prices, T is the expiration date, r is the risk-free interest rate to expiration, σ i is the volatility of the ith underlying stock, K i is the strike price of the ith underlying stock, β i is a coefficient so that all risky asset prices are at the same level.
It is noted that the most existing models use the strict assumptions, for example, perfect markets, constant values of both risk-free rate and volatility, log-normal distribution of share price dynamics, no dividends, continuous delta hedging. A divisible number of shares has not adequately represented the reality of the market [21][22][23].
More than three hundred years ago, the fractional differential equation was introduced. Nowadays, fractional calculus could be better for explaining the complicated incidents in the real situation than the traditional calculus. Fractional order models have been used to describe in many areas such as sciences, finance, physics and engineering disciplines [11,12,24,25]. Fractional calculus of the financial market was applied to the Black-Scholes model for expanding the theory of finance.
The time-fractional derivative which is considered in this paper is the Liouville-Caputo derivative because the initial condition for the fractional derivative is similar to the traditional derivative [25].
The Liouville-Caputo-type fractional derivative of order α is defined as [25]: where Γ(·) is a gamma function.
There are many researches that studied the fractional Black-Scholes model with one asset [11,[26][27][28][29][30][31]. The fractional Black-Scholes model is the generalized version of the classical model which extend the limitation of the model. Meng et al. [26] studied the fractional option pricing using Black-Scholes model. They applied the fractional Black-Scholes model to call option price for bank foreign exchange in China. Their results show that the fractional Black-Scholes model is better than the classical Black-Sholes model for estimating the effect in the market mechanism [26].
In this paper, we use the application of the Laplace transform Homotopy Perturbation Method (LHPM) to obtain the explicit solution of the time fractional Black-Scholes model. The LHPM is a method that combines with Homotopy Perturbation Method and Laplace transform. The LHPM gives an explicit solution which is a convergent series. We focus on the time fractional Black-Scholes model with two assets for the European call option. The LHPM is utilized for solving the problem. The analytical solution is a valuable tool to study the behaviors of the solution which are difficult to get the numerical solution especially fractional partial differential equations. The analytical solution provides a useful tool for studying financial behavior.
The paper is organized as follows. We present the time fractional Black-Scholes model in Section 2. The fractional derivatives are described in the sense of the Liouville-Caputo fractional derivative. The LHPM is presented in Section 3. In Section 3.1, the explicit solution of this problem is carried out by using LHPM. In addition, numerical examples and discussions are obtained in Section 4. Finally, a conclusion is presented in Section 5.

Mathematical Model
We consider the standard Black-Scholes partial differential equation with two assets for European-style option, efficient markets, perfect liquidity and no dividends during the option's life. Throughout this paper, we assume that σ 1 , σ 2 , ρ and r are constants.
∂c ∂τ with the terminal condition: and boundary conditions: By changing the variables [32]: Equation (1) can be written as: ∂c ∂τ which satisfies the terminal condition: and boundary conditions: By changing again of variables for eliminating the last term on the left hand side of Equation (2), we define v as: and substitute into Equation (2), we have ∂v ∂τ with the terminal condition: v(x, y, T) = max β 1 e x+(r− 1 and the boundary conditions: To be able to solve the initial boundary value problem (3)-(5), a forward time coordinate t = T − τ is introduced and used in Equation (3). Hence, we obtain ∂v ∂t subject to the initial condition: v(x, y, 0) = max β 1 e x +β 2 e y − K, 0 , and the boundary conditions: 2 )T . By replacing the Liouville Caputo fractional derivative, we obtain the following fractional-time order Black-Scholes model with α ∈ (0, 1] equipped with the initial and boundary conditions: subject to the initial condition: v(x, y, 0) = max β 1 e x +β 2 e y − K, 0 , and boundary conditions:

Basic Ideas of Time Fractional Black-Scholes Model with LHPM
In this section, the general form of the time-dependent differential equation can be expressed in the form of: where u(x, y, t) denotes an unknown function, f (x, y, t) denotes a known analytic function and A is determined by: A general equation can be rewritten as: with the initial condition: u(x, y, 0) = h(x, y) for any (x, y) ∈ R × R and the boundary condition: where B is a boundary operator. Firstly, in Equation (11), taking the Laplace transform with respect to time variable t yields L D α t u(x, y, t) + L Ñ (u(x, y, t)) = L f (x, y, t) .
Using the Laplace transform of Liouville-Caputo fractional derivative [33], we then obtain By taking the inverse Laplace transform on Equation (12), we have where the function G(x, y, t) presents the term arising from the source term and the prescribed initial conditions and boundary conditions. We now using the technique of HPM [14,15] to construct the function ν(x, y, t; p). We set + p ν(x, y, t; p) − G(x, y, t) where p ∈ [0, 1] denotes homotopy parameter or an embedding parameter and ν 0 (x, y, t) denotes the initial approximation of Equation (13). Rearranging (13) gives ν(x, y, t; p) = ν 0 (x, y, t) − p ν 0 (x, y, t) − G(x, y, t) For p = 0 and p = 1, we have H(ν(x, y, t; 1), 1) =ν(x, y, t; 1) − G(x, y, t) From Equations (14) and (15), we have By equating the coefficients corresponding to the power of p on both sides of Equation (16), the sequence ν n is carried out as follows: ν 0 (x, y, t) = ν 0 (x, y, t), When p converges to 1, the approximate solution of Equation (10) can be expressed in the form of: which leads to the explicit solution when infinite series converges.

A Solution of Time Fractional Black-Scholes Model by LHPM
In this section, we find the solution of time fractional Black-Scholes model (7) subjecting to the terminal Condition (8) and boundary Condition (9) by LHPM techniques.
Firstly, we let and express Equation (7) in the general form of: , y, t)).
By taking the Laplace transform with respect to time variable t, we get It follows from Laplace transform of the Liouville-Caputo fractional derivative that Equation (19) becomes, The inverse Laplace transform of an Equation (20) is obtained as: v(x, y, t) =max(β 1 e x +β 2 e y − K, 0) By applying techniques of HPM, we can construct the function v(x, y, t; p) which satisfies the following equation or v(x, y, t; p) = v 0 (x, y, t) − p v 0 (x, y, t) + p max(β 1 e x +β 2 e y − K, 0) where p ∈ [0, 1] is an embedded parameter and v 0 (x, y, t) is an initial approximation of Equation (21) which can be freely chosen [34]. For this model, we choose v 0 (x, y, t) as: v 0 (x, y, t) = max(β 1 e x +β 2 e y − K, 0) + e x+y t α .

Numerical Examples
In this section, an series solution of European call option based on Black-Scholes model with two assets as in Equation (25) is computed by using MATLAB programming. The simulations are carried out using the financial parameters given in Table 1. The graphs of the transformed explicit solution v and original explicit solution c in the case of call option are plotted in Figures 1-5. In Figure 1, the solutions v and c are plotted at a day before an expiration date over a range of 0 ≤ S 1 ≤ 200 and 0 ≤ S 2 ≤ 200 surrounding at the strike price with order α = 0.9. The results show that the option values increase significantly when the stock prices increase.
By setting S 2 = 10, the solutions v and c with order α = 0.9 are plotted in Figure 2a,b. With increasing S 1 from 0 to 50, the option price c reaches to zero. It is similar to c, v also reaches to zero when x increases from 2 to 4. The option price c increases linearly when the stock price is greater than 50. The solution v increases exponentially when x is greater than 4. Figure 3 shows the surface plot of call option with S 1 = 10 over a range of stock price 0 ≤ S 2 ≤ 200 and time 0 ≤ t ≤ 1. With increasing S 2 from 0 to 40, the option price c reaches to zero. v also reaches to zero when x increase from 2 to 3.8. After that the option price c increases linearly when stock price is greater than 40. The solution v increases exponentially when x is greater than 3.8.    At a day before an expiration date, the European call option price over the stock price S 1 and S 2 for various order α = 0.5, 0.7 and 0.9 is investigate. Effects of fractional order α on c for different orders when S 2 = 5 and S 1 = 10 are shown in Figure 4a,b, respectively. The comparison indicates that a higher order α gives a lower call option price. In Figure 5, the solution plot of v is similar trend to c. It is noted that a higher order α gives a lower option price v. Moreover, the option prices v increases rapidly after x = 4.2452 and y = 3.4589 as shown in Figure 5a,b, respectively. Consequently, we can conclude that the effect of time derivative order α has a small effect on the option price.  The values of European call option with the correlation varying from −1 to 1 is presented in Figure 6 at a day before an expiration date. The effects of stock price S 2 with some fixed stock prices S 1 and y with some fixed value x are investigated. Three values of S 1 including 50, 80 and 100 and three values of x including 3.8671, 4.3371 and 4.6556 are chosen for S 2 = 5 and y = 1.6093. The result shows that the relationship between the European call option and the correlation is the increasing linear pattern. In addition, the rate of change of European call options v and c with respect to ρ are also similar to each other as shown in Figure 6a,b.

Conclusions
The fractional Black-Scholes equations is a generalized version of the classical model which extend the restriction of using the model for finding the option price. As shown in the result, you can see that if parameters were not changed in the model, we obtain the difference value of option price with the different order. As a result, the fractional Black-Scholes model is more practical than the integer order model. This work introduced the Laplace homotopy perturbation method in order to find the explicit solution in time-fractional Black-Scholes model with two assets. By using the LHPM, the explicit solution can be written in the form of a special function, generalized Mittag-Leffer function. The benefit of the explicit solution form is easy to use for finding the European call option which depends on two stock price. Moreover, the numerical examples of the solution are presented to illustrate the explicit solution.