Application of the Laplace Homotopy Perturbation Method to the Black–Scholes Model Based on a European Put Option with Two Assets

Abstract

In this paper, the Laplace homotopy perturbation method (LHPM) is applied to obtain the approximate solution of Black–Scholes partial differential equations for a European put option with two assets. Different from all other approximation methods, LHPM provides a simple way to get the explicit solution which is represented in the form of a Mellin–Ross function. The numerical examples represent that the solution from the proposed method is easy and effective.

1. Introduction

Over the last half century, financial behavior has become an important thing for humans. The price crashes, momentum, and bubbles are investigated for understanding the dynamics of the financial system. Recently, mathematical models are used to describe and predict the financial behavior of the market in the future [1,2,3,4,5,6,7,8,9,10,11,12]. The stock has played a major role in the financial market. The options are the essential key for the derivative which is generally used in the market. A derivative value can be obtained from many sorts of assets—for instance, interest rate, shares, currency, etc. [13,14].

In the financial market, there are many styles of options such as American, European, Asian, Canary, etc. An option is a financial contract between buyers and sellers. The option price is a premium which is obtained from the writer of the option selling options to the buyer. In general, the options can be divided into two types. These are call options that are paid for conferring the right to buy and put options that are paid for conferring the right to sell.

In this paper, we focus only on a put option for European options that may be exercised only on the expiration date of the contract due to a lack of research proposing the analytical solution of this model with two assets. The dominant benefit of the European option allows holders to buy products from sellers with an agreement involving the price and date before exercising. This process provides a benefit to the seller and buyer because an agreement can be considered before making a decision. Moreover, the criteria of this process can be applied to another process such as the purchasing of products.

To investigate the option pricing in the market, the Black–Scholes model is widely used to forecast the behavior of a European option on a stock. The Black–Scholes model was proposed by Black, Scholes, and Merton [15] in 1973. It is set up by a geometric Brownian motion and governed by a partial differential equation with respect to time and stock price. The Black–Scholes model was studied and modified for predicting the option price in the actual market [16,17,18,19]. The main key components of the Black–Scholes model are the risk-free rate, the underlying stock price, strike price, volatility, and expiration date. Recently, the option price in a market depends on other markets, which results in the multidimensional Black–Scholes model becoming more productive than the one-dimensional model. A number of methods are used to find the solution of the multi-dimensional Black–Scholes model—for example, the finite different method [13,14,20,21,22], the finite element method [23], the Mellin transform method [24], and the adaptive wavelet precise integration method [25]. Moreover, the radial basis function (RBF) based methods that are commonly used to solve the Black–Scholes equations problems; as references papers below give RBF theory [26] and applications [27,28].

The Homotopy Perturbation method (HPM), proposed by He [29], has been successfully applied to find analytical solutions for many problems in engineering, physics, biology, finance, and science [29,30,31,32,33,34,35,36]. The solutions to this method are represented as an infinite series. In particular, the solution from HPM converges rapidly to the exact solution. HPM has become a powerful tool to solve mathematical problems [30,37]. Nevertheless, HPM cannot solve some complicated problems. The Laplace homotopy perturbation method (LHPM) is introduced [38,39,40]. LHPM is a method combined with HPM and Laplace transform. For example, the explicit solutions of the Black–Scholes equation for call options are carried out by LHPM [41,42]. In this work, we study the Black–Scholes equation with the basket option based on a European put option with two assets. The Laplace transformation homotopy perturbation method is used to find the explicit approximate solution of the problem.

The remainder of this paper is organized as follows. In Section 2, we present the mathematical model that is the Black–Scholes equation with a basket option based on a European put option with two assets. Then, the concept of Laplace transform homotopy perturbation method is shown in Section 3. The explicit solution of the Black–Scholes equation for a European put option with two assets, solved by using LHPM, is presented in Section 4. The numerical simulations are obtained in Section 5. Finally, we present conclusions in Section 6.

2. The Mathematical Model

We investigated the standard two-dimensional Black–Scholes partial differential equation for European put options including efficient markets, no dividends and perfect liquidity for all of the option’s life. The two-dimensional European put option depends on the prices of two underlying assets $A_1$ and $A_2$. The Black–Scholes model is based on

p 

the value of put option depending on the stock prices $\{A_1,A_2\}$ at time t,

ρ 

the correlation between the two underlying stock prices $A_1$ and $A_2$,

w_i 

the portions of underlying stock $A_i$ for $i=1,2$,

K_i 

the strike price of the ith underlying stock for $i=1,2$,

r 

the risk-free interest rate,

σ_i 

the volatility of the ith underlying stock for $i=1,2$,

T 

the expiration date.

With the assumptions of the Black–Scholes model [15], the two-asset Black–Scholes option price $p(t,A_1,A_2)$ is in the form of the partial differential equation

$$\frac{\partial p}{\partial \tau }+\frac{1}{2}{\sigma }_1^2{A}_1^2\frac{{\partial }^{2}p}{\partial A_1^2}+\frac{1}{2}{\sigma }_2^2{A}_2^2\frac{{\partial }^{2}p}{\partial A_2^2}+\rho {\sigma }_1{\sigma }_2{A}_1{A}_2\frac{{\partial }^{2}p}{\partial A_1\partial A_2}+r{A}_1\frac{\partial p}{\partial A_1}+r{A}_2\frac{\partial p}{\partial A_2}-rc=0,$$

(1)

for $A_1,A_2\in [0,\infty ),\tau \in [0,T],$ with the terminal condition:

$$p(A_1,A_2,T)=\max (K-(w_1{A}_1+w_2{A}_2),0),$$

(2)

in which $K=max\{K_1,K_2\},$ and boundary condition

$$\begin{array}{cc}\hfill p(0,A_2,\tau )& =max(K{e}^{-r(T-\tau )}-w_2{A}_2,0),\hfill \\ \hfill p(A_1,0,\tau )& =max(K{e}^{-r(T-\tau )}-w_1{A}_1,0),\hfill \\ \hfill p(A_1,A_2,\tau )& =max(K{e}^{-r(T-\tau )}-(w_1{A}_1+w_2{A}_2),0)\mathrm{as}{A}_1\to \infty \mathrm{or}{A}_2\to \infty .\hfill \end{array}$$

(3)

To be able to change the problem to the initial boundary value problem, we introduce a forward time $t=T-\tau$. Thus, Equation (1) can be written as

$$\frac{\partial p}{\partial t}-\frac{1}{2}{\sigma }_1^2{A}_1^2\frac{{\partial }^{2}p}{\partial A_1^2}-\frac{1}{2}{\sigma }_2^2{A}_2^2\frac{{\partial }^{2}p}{\partial A_2^2}-\rho {\sigma }_1{\sigma }_2{A}_1{A}_2\frac{{\partial }^{2}c}{\partial A_1\partial A_2}-r{A}_1\frac{\partial p}{\partial A_1}-r{A}_2\frac{\partial p}{\partial A_2}+rp=0.$$

Without loss of generality, we can consider $A_1$ as x, and $A_2$ as y. Hence, this above equation becomes

$$\frac{\partial p}{\partial t}=\frac{1}{2}{\sigma }_1^2{x}^2\frac{{\partial }^{2}p}{\partial x^2}+\frac{1}{2}{\sigma }_2^2{y}^2\frac{{\partial }^{2}p}{\partial y^2}+\rho {\sigma }_1{\sigma }_{2}xy\frac{{\partial }^{2}p}{\partial x\partial y}+rx\frac{\partial p}{\partial x}+ry\frac{\partial p}{\partial y}-rp.$$

(4)

The terminal condition becomes the initial condition as follows:

$$p(x,y,0)=\max (K-(w_1{A}_1+w_2{A}_2),0).$$

3. The Basic Ideas of Homotopy Perturbation Method with Laplace Transform

In general, the basic ideas of LHPM for the time-dependent differential equation are given as follows:

$$\mathcal{D}\left(u\right(x,y,t\left)\right)-f(x,y,t)=0,$$

(5)

where $\mathcal{D}$ denotes a differential operator, $u(x,y,t)$ denotes an unknown function and $f(x,y,t)$ denotes a known analytic function. $\mathcal{D}$ can be determined as

$$\mathcal{D}\left(u(x,y,t)\right)=\frac{\partial }{\partial t}u(x,y,t)+N\left(u(x,y,t)\right)\mathrm{for}(x,y,t)\in \mathbb{R}\times \mathbb{R}\times [0,T],$$

where N is the remaining part of $\mathcal{D}$.

The general equation on the domain $\mathsf{\Omega }=\mathbb{R}\times \mathbb{R}\times [0,T]$ can be written as the following equation:

$$\frac{\partial }{\partial t}u(x,y,t)+N\left(u(x,y,t)\right)=f(x,y,t),(x,y,t)\in \mathsf{\Omega },$$

(6)

with the initial condition

$$u(x,y,0)=h(x,y)\mathrm{for}\mathrm{any}(x,y)\in \mathbb{R}\times \mathbb{R},$$

and the boundary condition:

$$B\left(u,\frac{\partial u}{\partial x},\frac{\partial u}{\partial y},\frac{\partial u}{\partial t}\right)=0,$$

where B denotes the boundary operator.

Firstly, applying the Laplace transform with respect to t on both sides of Equation (6) yields

$$\mathcal{L}\left\{\frac{\partial }{\partial t}u(x,y,t)\right\}+\mathcal{L}\left\{N\left(u\right(x,y,t\left)\right)\right\}=\mathcal{L}\left\{f(x,y,t)\right\}.$$

By the property of the Laplace transform for differentiation, we obtain

$$\begin{array}{c}\hfill \mathcal{L}\left\{u(x,y,t)\right\}=s^{-1}h(x,y)-s^{-1}\mathcal{L}\left\{N\left(u\right(x,y,t\left)\right)\right\}+s^{-1}\mathcal{L}\left\{f(x,y,t)\right\}.\end{array}$$

(7)

Next, we take the inverse Laplace transform to the Equation (7), so that it becomes

$$u(x,y,t)=G(x,y,t)-{\mathcal{L}}^{-1}\left\{s^{-1}\mathcal{L}\left\{N\left(u\right(x,y,t\left)\right)\right\}\right\},$$

where the function $G(x,y,t)$ represents the term resulting from the source term and the determined initial conditions and boundary conditions.

Applying the homotopy perturbation method (HPM) [31,32], the function v can be constructed

$$v(x,y,t;q):\mathsf{\Omega }\times [0,1]\to \mathbb{R}.$$

We define

$$\begin{array}{ccc}\hfill H(v(x,y,t;q),q)& =& (1-p)[v(x,y,t;q)-\widetilde{v_0}(x,y,t)]+q[v(x,y,t;q)-G(x,y,t)\hfill \\ & & +{\mathcal{L}}^{-1}\left\{s^{-1}\mathcal{L}\left\{N\left(v\right(x,y,t;q\left)\right)\right\}\right\}]=0,\hfill \end{array}$$

(8)

where $q\in [0,1]$ is a homotopy parameter and $\widetilde{v_0}(x,y,t)$ is an initial approximation of the Equation (8).

Equation (8) is called the homotopy equation, and it can be written as the following equation:

$$\begin{array}{ccc}v(x,y,t;q)\hfill & =\hfill & \widetilde{v_0}(x,y,t)-q[\widetilde{v_0}(x,y,t)-G(x,y,t)\hfill \\ & & +{\mathcal{L}}^{-1}\left\{s^{-1}\mathcal{L}\left\{N\left(v\right(x,y,t;q\left)\right)\right\}\right\}].\hfill \end{array}$$

(9)

From Equations (8) and (9), it easy to obtain that

$$\begin{array}{ccc}\hfill q=0\to H\left(v\right(x,y,t;0),0)& =& v(x,y,t;0)-\widetilde{v_0}(x,y,t)=0,\hfill \\ \hfill q=1\to H\left(v\right(x,y,t;1),1)& =& v(x,y,t;1)-G(x,y,t)\hfill \\ & & +{\mathcal{L}}^{-1}\left\{s^{-1}\mathcal{L}\left\{N\left(v\right(x,y,t;1\left)\right)\right\}\right\}=0.\hfill \end{array}$$

By the technique of the homotopy perturbation method, the solution $v(x,y,t;q)$ in the Equation (9) can be expressed by the form

$$\begin{array}{ccc}\hfill v(x,y,t;q)& =& \sum _{i=0}^{\infty }{q}^i{v}_i(x,y,t).\hfill \end{array}$$

(10)

By substituting Equation (10) into Equation (9) and using HPM, it can be expressed as follows:

$$\begin{array}{ccc}\hfill \sum _{i=0}^{\infty }{q}^i{v}_i(x,y,t)& =& \widetilde{v_0}(x,y,t)-q[\widetilde{v_0}(x,y,t)-G(x,y,t)\hfill \\ & & +{\mathcal{L}}^{-1}\left\{s^{-1}\mathcal{L}\left\{N\left(\sum _{i=0}^{\infty }{q}^i{v}_i(x,y,t;p)\right)\right\}\right\}].\hfill \end{array}$$

By comparing the coefficients of the corresponding term of power q on both sides of the above equation, the approximate solution $v_i$ can be written as the following recurrence relations:

$$\begin{array}{ccc}\hfill v_0(x,y,t)& =& \widetilde{v_0}(x,y,t)\hfill \\ \hfill v_1(x,y,t)& =& G(x,y,t)-\widetilde{v_0}(x,y,t)\hfill \\ & & -{\mathcal{L}}^{-1}\left\{s^{-1}\mathcal{L}\left\{N\left(\widetilde{v_0}(x,y,t)\right)\right\}\right\}\hfill \\ \hfill v_m(x,y,t)& =& -{\mathcal{L}}^{-1}\left\{s^{-1}\mathcal{L}\left\{N\left(v_{m-1}(x,y,t)\right)\right\}\right\}\mathrm{when}m\ge 2.\hfill \end{array}$$

According to the solution of Equation (10), we obtain

$$v(x,y,t;q)=v_0(x,y,t)+q{v}_1(x,y,t)+q^2{v}_2(x,y,t)+q^3{v}_3(x,y,t)+\dots .$$

By letting q converge to 1, the approximate solution in this problem (5) can be expressed as follows:

$$u(x,y,t)=v(x,y,t;1)=v_0(x,y,t)+v_1(x,y,t)+v_2(x,y,t)+v_3(x,y,t)+\dots .$$

(11)

Particularly, the series (11) approaches the explicit solution when an infinite series converges.

4. An Analytical Solution of the Black–Scholes Model for European Put Options with Two Assets by Using LHPM

In this section, the Black–Scholes model with two assets in Equation (4) is investigated

$$\begin{array}{c}\hfill \frac{\partial p}{\partial t}=\frac{1}{2}{\sigma }_1^2{x}^2\frac{{\partial }^{2}p}{\partial x^2}+\frac{1}{2}{\sigma }_2^2{y}^2\frac{{\partial }^{2}p}{\partial y^2}+\rho {\sigma }_1{\sigma }_{2}xy\frac{{\partial }^{2}p}{\partial x\partial y}+rx\frac{\partial p}{\partial x}+ry\frac{\partial p}{\partial y}-rp\end{array}$$

(12)

with the initial condition

$$\begin{array}{ccc}\hfill p(x,y,0)& =& \max (K-(w_{1}x+w_{2}y),0).\hfill \end{array}$$

By taking the Laplace transform with respect to t to the Equation (12), we obtain

$$\begin{array}{c}\hfill \mathcal{L}\left\{\frac{\partial p}{\partial t}\right\}=\mathcal{L}\left\{\frac{1}{2}{\sigma }_1^2{x}^2\frac{{\partial }^{2}p}{\partial x^2}+\frac{1}{2}{\sigma }_2^2{y}^2\frac{{\partial }^{2}p}{\partial y^2}+\rho {\sigma }_1{\sigma }_{2}xy\frac{{\partial }^{2}p}{\partial x\partial y}+r\left(x\frac{\partial p}{\partial x}+y\frac{\partial p}{\partial y}\right)-rp\right\},\end{array}$$

(13)

setting

$$N\left(p(x,y,t)\right)=\frac{1}{2}{\sigma }_1^2{x}^2\frac{{\partial }^{2}p}{\partial x^2}+\frac{1}{2}{\sigma }_2^2{y}^2\frac{{\partial }^{2}p}{\partial y^2}+\rho {\sigma }_1{\sigma }_{2}xy\frac{{\partial }^{2}p}{\partial x\partial y}+r\left(x\frac{\partial p}{\partial x}+y\frac{\partial p}{\partial y}\right)-rp.$$

Equation (13) becomes

$$\begin{array}{ccc}\hfill \mathcal{L}\left\{p(x,y,t)\right\}& =& \frac{1}{s}p(x,y,0)+\frac{1}{s}\mathcal{L}\left\{N\left(p\right(x,y,t\left)\right)\right\},\hfill \\ \hfill \mathcal{L}\left\{p(x,y,t)\right\}& =& \frac{1}{s}\max (K-w_{1}x-w_{2}y,0)+\frac{1}{s}\mathcal{L}\left\{N\left(p\right(x,y,t\left)\right)\right\}.\hfill \end{array}$$

Then, taking the inverse Laplace transform to the previous equation, we get

$$\begin{array}{c}\hfill p(x,y,t)=\max (K-w_{1}x-w_{2}y,0)+{\mathcal{L}}^{-1}\left\{\frac{1}{s}\mathcal{L}\left\{N\left(p\right(x,y,t\left)\right)\right\}\right\}.\end{array}$$

Applying the Homotopy Perturbation Method, we can construct the following equation:

$$\begin{array}{cc}\hfill (1-q)\left(p\right(x,y,t;q)& -\widetilde{p_0}(x,y,t))+q[p(x,y,t;q)\hfill \\ & -\max (K-w_{1}x-w_{2}y,0)-{\mathcal{L}}^{-1}\left\{\frac{1}{s}\mathcal{L}\left\{N\left(p\right(x,y,t;q\left)\right)\right\}\right]=0\hfill \end{array}$$

or we can write

$$\begin{array}{ccc}\hfill p(x,y,t;q)& =& \widetilde{p_0}(x,y,t)-q\widetilde{p_0}(x,y,t)+q\max (K-w_{1}x-w_{2}y,0)\hfill \\ & & +q{\mathcal{L}}^{-1}\left\{\frac{1}{s}\mathcal{L}\left\{N\left(p\right(x,y,t;q\left)\right)\right\}\right\},\hfill \end{array}$$

(14)

where $\widetilde{p_0}(x,y,t)$ is an initial approximation of Equation (14) that can be freely chosen.

For this model, we choose $\widetilde{p_0}(x,y,t)$ in the form

$$\widetilde{p_0}(x,y,t)=max(K-w_{1}x-w_{2}y,0)+\rho (x^2+y^2)t.$$

By substituting $\widetilde{p_0}(x,y,t)$ to Equation (14), we have

$$\begin{array}{ccc}\hfill p(x,y,t;q)& =& \max (K-w_{1}x-w_{2}y,0)+\rho (x^2+y^2)t+q[-(x^2+y^2)t\hfill \\ & & +{\mathcal{L}}^{-1}\left\{\frac{1}{s}\mathcal{L}\left\{N\left(p\right(x,y,t;q\left)\right)\right\}\right\}].\hfill \end{array}$$

(15)

By the HPM technique, we assume the solution of problem (15) to be written in the form

$$\begin{array}{c}\hfill p(x,y,t;q)=\sum _{n=0}^{\infty }{q}^n{p}_n(x,y,t).\end{array}$$

(16)

Substituting Equation (16) into Equation (15), we get

$$\begin{array}{ccc}\hfill \sum _{n=0}^{\infty }{q}^n{p}_n(x,y,t)& =& \max (K-w_{1}x-w_{2}y,0)+\rho (x^2+y^2)t\hfill \\ & & +q\left(-(x^2+y^2)t+{\mathcal{L}}^{-1}\left\{\frac{1}{s}\mathcal{L}\right\{\frac{1}{2}{\sigma }_1^2{x}^2\sum _{n=0}^{\infty }{q}^n\frac{{\partial }^2{p}_n}{\partial x^2}\right.\hfill \\ & & +\frac{1}{2}{\sigma }_2^2{y}^2\sum _{n=0}^{\infty }{q}^n\frac{{\partial }^2{p}_n}{\partial y^2}+\rho {\sigma }_1{\sigma }_{2}xy\sum _{n=0}^{\infty }{q}^n\frac{{\partial }^2{p}_n}{\partial x\partial y}\hfill \\ & & +r\left(x\sum _{n=0}^{\infty }{q}^n\frac{\partial p_n}{\partial x}+y\sum _{n=0}^{\infty }{q}^n\frac{\partial p_n}{\partial y}\right)\left.-r\sum _{n=0}^{\infty }{q}^n{p}_n\left\}\right\}\right).\hfill \end{array}$$

Equating the corresponding power of p on both sides, we obtain the following:

$$\begin{array}{ccc}\hfill q^0;p_0(x,y,t)& =& \max (K-w_{1}x-w_{2}y,0)+\rho (x^2+y^2)t,\hfill \\ \hfill q^1;p_1(x,y,t)& =& \frac{t^2\rho }{2}\left[({\sigma }_1^2+r)x^2+({\sigma }_2^2+r)y^2\right]-t[{({\sigma }_1^2+r)}^0{x}^2+{({\sigma }_2^2+r)}^0{y}^2\hfill \\ & & -r(\max (w_1,0)x+\max (w_2,0)y)+r\max (K-w_{1}x-w_{2}y,0)],\hfill \\ \hfill q^2;p_2(x,y,t)& =& \frac{t^3\rho }{6}\left[{({\sigma }_1^2+r)}^2{x}^2+{({\sigma }_2^2+r)}^2{y}^2\right]-\frac{t^2}{2}[{({\sigma }_1^2+r)}^1{x}^2+{({\sigma }_2^2+r)}^1{y}^2\hfill \\ & & +r^2(\max (w_1,0)x+\max (w_2,0)y)-r^2\max (K-w_{1}x-w_{2}y,0)],\hfill \\ \hfill q^3;p_3(x,y,t)& =& \frac{t^4\rho }{24}\left[{({\sigma }_1^2+r)}^3{x}^2+{({\sigma }_2^2+r)}^3{y}^2\right]-\frac{t^3}{6}[{({\sigma }_1^2+r)}^2{x}^2+{({\sigma }_2^2+r)}^2{y}^2\hfill \\ & & -r^3(\max (w_1,0)x+\max (w_2,0)y)+r^3\max (K-w_{1}x-w_{2}y,0)],\hfill \end{array}$$

and so on.

Proceeding in the same process, we obtain

$$\begin{array}{ccc}\hfill q^0;p_0(x,y,t)& =& \max (K-w_{1}x-w_{2}y,0)+\rho (x^2+y^2)t,\hfill \\ \hfill q^n;p_n(x,y,t)& =& \frac{t^{n+1}\rho }{\mathsf{\Gamma }(n+2)}\left[{({\sigma }_1^2+r)}^n{x}^2+{({\sigma }_2^2+r)}^n{y}^2\right]-\frac{t^n}{\mathsf{\Gamma }(n+1)}[{({\sigma }_1^2+r)}^{n-1}{x}^2\hfill \\ & & +{({\sigma }_2^2+r)}^{n-1}{y}^2+{(-1)}^n{r}^n(\max (w_1,0)x+\max (w_2,0)y)\hfill \\ & & -{(-1)}^n{r}^n\max (K-w_{1}x-w_{2}y,0)]\mathrm{when}n\ge 1.\hfill \end{array}$$

Thus, the solution $p(x,y,t;p)$ of problem (15) is given by

$$\begin{array}{ccc}\hfill p(x,y,t;q)& =& \sum _{n=0}^{\infty }{q}^n{p}_n(x,y,t)\hfill \\ & =& p_0(x,y,t)+\sum _{n=1}^{\infty }{q}^n\{\frac{t^{n+1}\rho }{\mathsf{\Gamma }(n+2)}\left[{({\sigma }_1^2+r)}^n{x}^2+{({\sigma }_2^2+r)}^n{y}^2\right]\hfill \\ & & -\frac{t^n}{\mathsf{\Gamma }(n+1)}[{({\sigma }_1^2+r)}^{n-1}{x}^2+{({\sigma }_2^2+r)}^{n-1}{y}^2\hfill \\ & & +{(-1)}^n{r}^n(\max (-w_1,0)x+\max (-w_2,0)y)\hfill \\ & & -{(-1)}^n{r}^n\max (K-w_{1}x-w_{2}y,0)\left]\right\}\hfill \\ & =& p_0(x,y,t)+\sum _{n=0}^{\infty }{q}^{n+1}\{\frac{t^{n+2}\rho }{\mathsf{\Gamma }(n+3)}\left[{({\sigma }_1^2+r)}^{n+1}{x}^2+{({\sigma }_2^2+r)}^{n+1}{y}^2\right]\hfill \\ & & -\frac{t^{n+1}}{\mathsf{\Gamma }(n+2)}[{({\sigma }_1^2+r)}^n{x}^2+{({\sigma }_2^2+r)}^n{y}^2\hfill \\ & & +{(-1)}^{n+1}{r}^{n+1}(\max (-w_1,0)x+\max (-w_2,0)y)\hfill \\ & & -{(-1)}^{n+1}{r}^{n+1}\max (K-w_{1}x-w_{2}y,0)\left]\right\}.\hfill \end{array}$$

By setting q converges to 1, we obtain

$$\begin{array}{ccc}\hfill p(x,y,t)& =& p_0(x,y,t)+\sum _{n=0}^{\infty }\{\frac{t^{n+2}\rho }{\mathsf{\Gamma }(n+3)}\left[{({\sigma }_1^2+r)}^{n+1}{x}^2+{({\sigma }_2^2+r)}^{n+1}{y}^2\right]\hfill \\ & & -\frac{t^{n+1}}{\mathsf{\Gamma }(n+2)}[{({\sigma }_1^2+r)}^n{x}^2+{({\sigma }_2^2+r)}^n{y}^2\hfill \\ & & +{(-1)}^{n+1}{r}^{n+1}(max(-w_1,0)x+max(-w_2,0)y)\hfill \\ & & -{(-1)}^{n+1}{r}^{n+1}max(K-w_{1}x-w_{2}y,0)\left]\right\}.\hfill \end{array}$$

Therefore, we obtain the explicit solution of problem (12) as follows:

$$\begin{array}{cc}\hfill p(x,y,t)=& \max (K-w_{1}x-w_{2}y,0)+x^2\left(\rho ({\sigma }_1^2+r)E_{2,{\sigma }_1^2+r}\left(t\right)-E_{1,{\sigma }_1^2+r}\left(t\right)+\rho t\right)\hfill \\ & +y^2\left(\rho ({\sigma }_2^2+r)E_{2,{\sigma }_2^2+r}\left(t\right)-E_{1,{\sigma }_2^2+r}\left(t\right)+\rho t\right)\hfill \\ & -r\left(\max (-w_1,0)x+\max (-w_2,0)y\right)E_{1,-r}\left(t\right)\hfill \\ & +r\max (K-w_{1}x-w_{2}y,0)E_{1,-r}\left(t\right),\hfill \end{array}$$

(17)

where ${\displaystyle E_{a,b}\left(t\right)=t^a\sum _{k=0}^{\infty }\frac{{\left(bt\right)}^k}{\mathsf{\Gamma }(a+k+1)}}$ is the Mellin–Ross function [43], which a and b being constants.

5. Solution Example

In this section, the explicit solution of the problem (12) as shown in Equation (17) is computed to obtain the value of European put option with the condition as in Equation (2) and $w_1=2,w_2=1$,

$$p(A_1,A_2,T)=max(K-(2A_1+A_2),0),$$

with strike price $K=75$. The risk-free interest rate is $5\%$ per year, $(r=0.05)$, the volatilities of the underlying assets $A_1$ and $A_2$ are ${\sigma }_1=5\%$ and ${\sigma }_2=10\%$, respectively. The expiration time is $T=1$ measured in years.

The solution surface plot of a European put option with the correlation, $\rho =0.8$ at the maturity time is showed in Figure 1. The stock prices are considered over the range $0\le A_1\le 200$ and $0\le A_2\le 200$ which surround the strike price. The result shows that the option price linearly decreases when the stock prices, $A_1$ and $A_2$, increase. After that, the option price reaches zero when the stock prices increase.

At a day before an expiration date, the value of put option price with the influence of stock price $A_1$ and fixed stock price $A_2$ is presented in Figure 2. The result indicated that the option price linearly dropped until the value of stock price $A_1$ was around 60; and the option price p reaches zero. In addition, the value of put option price with the influence of stock price $A_2$ and fixed stock price $A_1$ is also represented in Figure 3. The solution p is decreased linearly when the stock price is lower than 40 and the option price reaches zero when the stock price is greater than 40.

The different values of European put options for non-positive correlation case ${\rho }_1=-0.8$ and ${\rho }_2=0$ and non-negative correlation case ${\rho }_2=0$ and ${\rho }_3=0.8$ are plotted over a range of stock prices $0\le A_1\le 200$ and $0\le A_2\le 200$ at the maturity time, as shown in Figure 4. The results show that the correlation has an effect on the European put option. Furthermore, when stock prices both increase, this leads to an increase in the difference between the value of the European put options.

6. Conclusions

The Black–Scholes model is the most famous and useful mathematical model for describing the behavior of the financial market. In this work, the Black–Scholes model for a European put option with two assets is investigated. The LHPM is applied to find the explicit solution of the Black–Scholes model in the form of the infinite series, which is a special function, called the Mellin–Ross function. The advantage of the explicit solution of the Black–Scholes model is easy to implement for simulating the option price that depends on two asset prices. Moreover, the results can be easy to apply in real world problems for the financial system.