Homotopy Perturbation Method for Fractional Black-Scholes European Option Pricing Equations Using Sumudu Transform

Abstract

The homotopy perturbation method, Sumudu transform, and He’s polynomials are combined to obtain the solution of fractional Black-Scholes equation. The fractional derivative is considered in Caputo sense. Further, the same equation is solved by homotopy Laplace transform perturbation method. The results obtained by the two methods are in agreement. The approximate analytical solution of Black-Scholes is calculated in the form of a convergence power series with easily computable components. Some illustrative examples are presented to explain the efficiency and simplicity of the proposed method.

1. Introduction

Fractional differential equations have attracted much attention, recently, see, for instance [1–5]. This is mostly due to the fact that fractional calculus provides an efficient and excellent instrument for the description of many practical dynamical phenomena arising in engineering and scientific disciplines such as, physics, chemistry, biology, economy, viscoelasticity, electrochemistry, electromagnetic, control, porous media, and many more, see, for example, [6–9].

Many partial differential equations of fractional order have been studied and solved. For example many researchers studied the existence of solutions of the Black-Scholes model using many methods, see [10–14].

The homotopy perturbation method was first introduced and applied by He [15–17]. This method has been applied by many authors in many fields, for example, it is applied to nonlinear oscillator [18], nonlinear wave equation [19], nonlinear partial differential equations [20], integro-differential equation of fractional order [21], fuzzy differential equation [22], and other fields [23, 24]. Further homotopy perturbation methods are combined with Laplace transform to solve many problems such as one dimensional nonhomogeneous partial differential equations with a variable coefficient [25], Black-Scholes of fractional order [26], and parabolic partial differential equations [27]. The homotopy perturbation method coupled with Sumudu transform basically illustrates how Sumudu transform can be used to approximate the solutions of the linear and nonlinear differential equations by manipulating the homotopy perturbation method. In [28] Singh et al. studied the solution of linear and nonlinear partial differential equations by using the homotopy perturbation method coupled with Sumudu transform. Further, in [29] the authors proposed the homotopy perturbation method coupled with Sumudu transform to solve nonlinear fractional gas dynamics equation.

The Black-Scholes equation is one of the most significant mathematical models for a financial market. It is a second-order parabolic partial differential equation that governs the value of financial derivatives. This Black-Scholes model for the value of an option is described by the following equation: 𝜕𝑣+𝜕𝑡𝜎𝑥²2𝜕²𝑣𝜕𝑥²+𝑟(𝑡)𝑥𝜕𝑣(𝑡)()𝜕𝑥−𝑟𝑣=0,𝑥,𝑡∈𝑅⁺×(0,𝑇),0<𝛼≤1, (1) where 𝑣(𝑥,𝑡) is the European call option price at asset price 𝑥 and at time 𝑡,𝑇 is the maturity,𝑟(𝑡) is the risk free interest rate, and 𝜎(𝑥,𝑡) represents the volatility function of underlying asset. The payoff functions are 𝑣_𝑐(𝑥,𝑡)=max(𝑥−𝐸,0);𝑣_𝑝(𝑥,𝑡)=max(𝐸−𝑥,0), (2) where 𝑣_𝑐(𝑥,𝑡) and 𝑣_𝑝(𝑥,𝑡) are the value of the European call and put options, respectively, 𝐸 denotes the expiration price for the option, and the function max(𝑥,0) gives the large value between 𝑥 and 0.

In this paper, we consider the following fractional Black-Scholes of the form 𝜕^𝛼𝑣𝜕𝑡^𝛼+𝜎𝑥²2𝜕²𝑣𝜕𝑥²+𝑟(𝑡)𝑥𝜕𝑣(𝑡)()𝜕𝑥−𝑟𝑣=0,𝑥,𝑡∈𝑅⁺×(0,𝑇),0<𝛼≤1. (3)

In [29] Singh et al. used homotopy perturbation method coupled with Sumudu transform to solve fractional gas dynamics equation. The aim of this paper is to applied the homotopy perturbation method for fractional Black-Scholes equation by using He's polynomials and Sumudu transform.

2. Sumudu Transform

The Sumudu transform was first introduced and applied by Watugala [30] in (1998). For further details and properties of Sumudu transform see [31–34]. The Sumudu transform is defined over the set of functions: 𝐴=𝑓(𝑡)∶∃𝑀,𝜏₁,𝜏₂>0,𝑓(𝑡)<𝑀𝑒^{𝑡/𝜏𝑗},if𝑡∈(−1)^𝑗×[0,∞) (4) by the following formula 𝑓(𝑢)=𝑆𝑓(𝑡);𝑢=∶^∞₀𝑓(𝑢𝑡)𝑒⁻𝑡𝑑𝑡,𝑢∈(−𝜏,𝜏). (5) The existence and uniqueness was discussed in [35]. For further properties of Sumudu transform and its derivatives, see [36]. Some fundamental further established properties of Sumudu transform can be found in [31].

Similarly, this new transform was applied to one-dimensional neutron transport equation [37]. In [34] Kılıçman et al. show that there is a strong relationship between Sumudu and other integral transforms. Further in [33] the Sumudu transform was extended to the distributions, and some of their properties were also studied in [38]. Recently Kılıçman et al. applied this transform to solve system of differentials equations, for more details see [34, 35, 37–39].

3. Basic Definitions of Fractional Calculus

In this section, we give some basic definitions and properties of fractional calculus theory which will be used in this paper.

Definition 1. The Riemann-Liouville fractional integral operator of order 𝛼≥0 of a function 𝑓∈𝐶𝜇, 𝜇≥−1 is defined as follows: 𝐽^𝛼𝑓(𝑥)=1Γ(𝛼)^𝑥₀(𝑥−𝑡)^𝛼−1𝑓(𝑡)𝑑𝑡,𝛼>0,𝑡>0 (6) in particular 𝐽⁰𝑓(𝑥)=𝑓(𝑥).

For Riemann-Liouville fractional integral, one has 𝐽^𝛼𝑥^𝛾=Γ𝛾+1Γ𝑥𝛼+𝛾+1^𝛼+𝛾. (7)

Definition 2. The Caputo fractional derivative of 𝑓∈𝐶^𝑚₋₁, 𝑚∈𝑁 is defined as follows: 𝐷^𝛼𝑓(𝑥)=1Γ(𝑚−𝛼)^𝑥₀(𝑥−𝑡)^{𝑚−𝛼−1}𝑓^𝑚(𝑡)𝑑𝑡,𝑚−1<𝛼≤𝑚. (8)

Lemma 3. If 𝑚−1<𝛼≤𝑚, 𝑚∈𝑁, 𝑓∈𝐶^𝑚_𝜇, 𝜇>−1 then the following two properties hold: (1)𝐷^𝛼[𝐽^𝛼𝑓(𝑥)]=𝑓(𝑥), (2)𝐽^𝛼[𝐷^𝛼𝑓(𝑥)]=𝑓(𝑥)−∑^𝑚−1_𝑘=1𝑓^𝑘(0)(𝑥^𝑘/𝑘!).

Definition 4. The Mittag-Leffler function 𝐸_𝛼(𝑧) with 𝛼>0 is defined by the following series representation, valid in the whole complex plane: 𝐸_𝛼(𝑧)=^∞₀𝑧^𝑛Γ().𝛼𝑛+1 (9)

Definition 5. The Sumudu transform of the Caputo fractional derivative is defined as follows [40]: 𝑆𝐷^𝛼_𝑡𝑓(𝑡)=𝑢^−𝛼𝑆𝑓(𝑡)−^𝑚−1_𝑘=0𝑢^−𝛼+𝑘𝑓^(𝑘)0⁺,(𝑚−1<𝛼≤𝑚). (10)

4. Homotopy Perturbation Method

To illustrate the basic idea of this method, we consider the following nonlinear differential equation: 𝐴(𝑢)−𝑓(𝑟)=0,𝑟∈Ω (11) with boundary conditions 𝐵𝑢,𝜕𝑢𝜕𝑛=0,𝑟∈Γ, (12) where 𝐴 is a general differential operator, 𝐵 is a boundary operator, 𝑓(𝑟) is a known analytic function, and Γ is the boundary of the domain Ω.

In general, the operator 𝐴 can be divided into two parts 𝐿 and 𝑁, where 𝐿 is linear, while 𝑁 is nonlinear. Equation (11) therefor can be rewritten as follows: 𝐿(𝑢)+𝑁(𝑢)−𝑓(𝑟)=0. (13) By the homotopy technique [41, 42] we construct a homotopy 𝑣(𝑟,𝑝)∶Ω×[0,1]→𝑅 which satisfies 𝐻𝑣,𝑝=1−𝑝𝐿(𝑣)−𝐿𝑢₀𝐴(𝑣)(𝑟[]+𝑝−𝑓)=0𝑝∈0,1,𝑟∈Ω (14) or 𝐻𝑣,𝑝=𝐿(𝑣)−𝐿𝑢₀+𝑝𝐿𝑢₀+𝑝𝑁(𝑣)−𝑓(𝑟)=0, (15) where 𝑝∈[0,1] is an embedding parameter, and 𝑢₀ is an initial approximation of (11) which satisfies the boundary conditions.

From (14) and (15) we have 𝐻(𝑣,0)=𝐿(𝑣)−𝐿𝑢₀𝐻()(𝑣)(𝑟)=0,𝑣,1=𝐴−𝑓=0. (16) The changing in the process of 𝑝 from zero to unity is just that of 𝑣(𝑟,𝑝) from 𝑢₀(𝑟) to 𝑢(𝑟). In topology this is called deformation, and 𝐿(𝑣)−𝐿(𝑢₀) and 𝐴(𝑣)−𝑓(𝑟) are called homotopic.

Now, assume that the solution of (14), (15) can be expressed as 𝑣=𝑣₀+𝑝𝑣₁+𝑝²𝑣₂+⋯. (17) Setting 𝑝=1 results in the approximate solution of (11): 𝑢=lim_𝑝→1𝑣=𝑣₀+𝑣₁+𝑣₂+⋯. (18)

5. Homotopy Perturbation Method Coupled with Sumudu Transform

To illustrate the basic idea of this method, we consider the following nonlinear fractional differential equation: 𝐷^𝛼_𝑡𝑢()[𝑥]𝑢()[𝑥]𝑢()()𝑥,𝑡+𝐿𝑥,𝑡+𝑁𝑥,𝑡=𝑞𝑥,𝑡,𝑡>0,𝑚−1<𝛼≤𝑚, (19) where 𝐷^𝛼_𝑡=𝜕^𝛼/𝜕𝑡^𝛼 is the fractional Caputo derivative of the function 𝑢(𝑥,𝑡), 𝐿 is the linear differential operator, 𝑁 is the nonlinear differential operator, and 𝑞(𝑥,𝑡) is the source term.

Now, applying the Sumudu transform on both sides of (19), we have 𝑆𝐷^𝛼_𝑡𝑢([𝐿[𝑥]𝑢()[𝑥]𝑢(𝑞(.𝑥,𝑡)+𝑆𝑥,𝑡+𝑁𝑥,𝑡)]=𝑆𝑥,𝑡) (20) Using the differential property of Sumudu transform, we have 𝑆[𝑢((𝑥)𝑥,𝑡)]=𝑓−𝑢^𝛼𝑆[𝐿[𝑥]𝑢(𝑥,𝑡)+𝑁[𝑥]𝑢(𝑥,𝑡)]+𝑢^𝛼𝑆𝑞(𝑥,𝑡). (21) Operating with Sumudu inverse on both sides of (21) 𝑢(𝑥,𝑡)=𝑄(𝑥,𝑡)−𝑆⁻¹𝑢^𝛼𝑆(𝐿[𝑥]𝑢(𝑥,𝑡)+𝑁[𝑥]𝑢(𝑥,𝑡)), (22) where 𝑄(𝑥,𝑡) represents the term arising from the source term and the prescribed initial conditions.

Now, applying the classical homotopy perturbation technique, the solution can be expressed as a power series in 𝑝 as given below: 𝑢(𝑥,𝑡)=^∞_𝑛=0𝑝^𝑛𝑢_𝑛(𝑥,𝑡), (23) where the homotopy parameter 𝑝 is considered as a small parameter (𝑝∈[0,1]).

We can decompose the nonlinear term as 𝑁𝑢(𝑥,𝑡)=^∞_𝑛=0𝑝^𝑛𝐻_𝑛(𝑢), (24) where 𝐻_𝑛 are He’s polynomials of 𝑢₀,𝑢₁,𝑢₂,…,𝑢_𝑛 [43–45], and it can be calculated by the following formula: 𝐻𝑢₀,𝑢₁,𝑢₂,…,𝑢_𝑛=1𝜕𝑛!𝜕𝑝^𝑛𝑛𝑁^∞_𝑛=0𝑝^𝑖𝑢_𝑖𝑝=0,𝑛=0,1,2,…. (25) By substituting (23) and (24) and using HPM [15] we get ^∞_𝑛=1𝑝^𝑛𝑢_𝑛()()𝑆𝑥,𝑡=𝑄𝑥,𝑡−𝑝⁻¹𝑢^𝛼𝑆𝐿^∞_𝑛=0𝑝^𝑛𝑢_𝑛(𝑥,𝑡)+^∞_𝑛=0𝑝^𝑛𝐻_𝑛(𝑢). (26) This is coupling of Sumudu transform and homotopy perturbation method using He’s polynomials. By equating the coefficient of corresponding power of 𝑝 on both sides, the following approximations are obtained as 𝑝⁰∶𝑢₀()(),𝑝𝑥,𝑡=𝑄𝑥,𝑡₁∶𝑢₁(𝑥,𝑡)=−𝑆⁻¹𝑢^𝛼𝑆𝐿[𝑥]𝑢₀(𝑥,𝑡)+𝐻₀(𝑢,𝑝₂)∶𝑢₂(𝑥,𝑡)=−𝑆⁻¹𝑢^𝛼𝑆𝐿[𝑥]𝑢₁(𝑥,𝑡)+𝐻₁(𝑢,𝑝₃)∶𝑢₃(𝑥,𝑡)=−𝑆⁻¹𝑢^𝛼𝑆𝐿[𝑥]𝑢₂(𝑥,𝑡)+𝐻₂(𝑢,⋮) (27) Proceeding in the same manner, the rest of the components 𝑢_𝑛(𝑥,𝑡) can be completely obtained, and the series solution is thus entirely determined. Finlay we approximate the solution 𝑢(𝑥,𝑡) by truncated series 𝑢(𝑥,𝑡)=lim^𝑁_{𝑁→∞𝑛=0}𝑢_𝑛(𝑥,𝑡). (28) This series solutions generally converge very rapidly.

6. Examples

In this section, we discuss the implementation of the proposed method.

Example 6. We consider the following fractional Black-Scholes option pricing equation as follows: 𝜕^𝛼𝑣𝜕𝑡^𝛼=𝜕²𝑣𝜕𝑥²+(𝑘−1)𝜕𝑣𝜕𝑥−𝑘𝑣,0<𝛼≤1 (29) subject to initial condition 𝑣(𝑥,0)=max𝑒^𝑥−1,0. (30) Applying Sumudu transform on both sides of (29) subject to initial condition (30), we get 𝑆[𝑣(𝑒𝑥,𝑡)]=max^𝑥−1,0+𝑢^𝛼𝑆𝑣_𝑥𝑥+(𝑘−1)𝑣_𝑥−𝑘𝑣. (31) Operating the inverse Sumudu transform on both sides in (31), we have 𝑣(𝑥,𝑡)=max𝑒^𝑥−1,0−𝑆⁻¹𝑢^𝛼𝑆𝑣_𝑥𝑥+(𝑘−1)𝑣_𝑥−𝑘𝑣. (32) Now, applying homotopy perturbation method ^∞_𝑛=0𝑝^𝑛𝑣_𝑛()𝑒𝑥,𝑡=max^𝑥−1,0−𝑝𝑆⁻¹𝑢^𝛼𝑆^∞_𝑛=0𝑝^𝑛𝐻_𝑛(𝑣), (33) where 𝐻_𝑛=𝑣_𝑛𝑥𝑥+(𝑘−1)𝑣_𝑛𝑥+𝑘𝑣_𝑛,𝑛∈𝑁. (34)
Equating the corresponding power of 𝑝 on both sides in (38), we have 𝑝⁰∶𝑣₀(𝑥,𝑡)=max𝑒^𝑥,𝑝−1,0₁∶𝑣₁(𝑥,𝑡)=𝑆⁻¹𝑢^𝛼𝑆𝐻₀(𝑣𝑒)=−max^𝑥,0−𝑘𝑡^𝛼Γ()𝑒𝛼+1+max^𝑥−1,0−𝑘𝑡^𝛼Γ(),𝑝𝛼+1₂∶𝑣₂(𝑥,𝑡)=𝑆⁻¹𝑢^𝛼𝑆𝐻₁(𝑣𝑒)=max^𝑥,0−𝑘𝑡^𝛼2Γ()𝑒2𝛼+1+max^𝑥−1,0−𝑘𝑡^𝛼2Γ(),⋮𝑝2𝛼+1_𝑛∶𝑣_𝑛(𝑥,𝑡)=𝑆⁻¹𝑢^𝛼𝑆𝐻_𝑛(𝑣𝑒)=max^𝑥,0−𝑘𝑡^𝛼𝑛Γ()𝑒𝑛𝛼+1+max^𝑥−1,0−𝑘𝑡^𝛼𝑛Γ().𝑛𝛼+1 (35) So that the solution 𝑣(𝑥,𝑡) of the problem is given by 𝑣(𝑥,𝑡)=lim^∞_{𝑝→1𝑖=0}𝑝^𝑖𝑢_𝑖()𝑒𝑥,𝑡=max^𝑥−1,0𝐸_𝛼−𝑘𝑡^𝛼+max𝑒^𝑥,01−𝐸_𝛼−𝑘𝑡^𝛼, (36) where 𝐸_𝛼(𝑧) is Mittag-Leffler function in one parameter. For special case 𝛼=1, we get 𝑣(𝑥,𝑡)=max𝑒^𝑥−1,0𝑒^−𝑘𝑡+max𝑒^𝑥,01−𝑒^−𝑘𝑡, (37) which is an exact solution of the given Black-Scholes equation (29) for 𝛼=1.
The behaviour of 𝑣(𝑥,𝑡) with respect to 𝑥 and 𝑡 when 𝛼=1 is given in Figure 1.

Figure 1 

The surface shows the 𝑣(𝑥,𝑡) for (29) with respect to 𝑥 and 𝑡 for 𝛼=1.

Example 7. We consider the following fractional Black-Scholes option pricing equation as follows: 𝜕^𝛼𝑣𝜕𝑡^𝛼+0.08(2+sin𝑥)²𝑥²𝜕²𝑣𝜕𝑥²+0.06𝜕𝑣𝜕𝑥−0.06𝑣=0,0<𝛼≤1 (38) subject to initial condition 𝑣(𝑥,0)=max𝑥−25𝑒^−0.06,0. (39) Firstly, applying Sumudu transform on both sides of (38) subject to initial condition (39), we get 𝑆[𝑣(𝑥,𝑡)]=max𝑥−25𝑒^−0.06,0−𝑢^𝛼𝑆×0.08(2+sin𝑥)²𝑥²𝑣_𝑥𝑥+0.06𝑣_𝑥−0.06𝑣. (40) Operating the inverse Sumudu transform on both sides in (40), we have 𝑣(𝑥,𝑡)=max𝑥−25𝑒^−0.06,0−𝑆⁻¹×𝑢^𝛼𝑆0.08(2+sin𝑥)²𝑥²𝑣_𝑥𝑥+0.06𝑣_𝑥−0.06𝑣. (41) Now, applying the homotopy perturbation method we have ^∞_𝑛=0𝑝^𝑛𝑣_𝑛(𝑥,𝑡)=max𝑥−25𝑒^−0.06𝑆−1,0−𝑝⁻¹𝑢^𝛼𝑆^∞_𝑛=0𝑝^𝑛𝐻_𝑛(𝑣), (42) where 𝐻_𝑛=0.08(2+sin𝑥)²𝑥²𝑣_𝑛𝑥𝑥+0.06𝑣_𝑛𝑥−0.06𝑣_𝑛,𝑛∈𝑁. (43) Equating the corresponding power of 𝑝 on both sides in (42), we have 𝑝⁰∶𝑣₀(𝑥,𝑡)=max𝑥−25𝑒^−0.06,𝑝,0₁∶𝑣₁(𝑥,𝑡)=𝑆⁻¹𝑢^𝛼𝑆𝐻₀(𝑣)=−𝑥−0.06𝑡^𝛼Γ(𝛼+1)+max𝑥−25𝑒^−0.06,0−0.06𝑡^𝛼Γ(),𝑝𝛼+1₂∶𝑣₂(𝑥,𝑡)=𝑆⁻¹𝑢^𝛼𝑆𝐻₁(𝑣)=−𝑥−0.06𝑡^𝛼2Γ(2𝛼+1)+max𝑥−25𝑒^0.06,0−0.06𝑡^𝛼2Γ(),⋮𝑝2𝛼+1_𝑛∶𝑣_𝑛(𝑥,𝑡)=𝑆⁻¹𝑢^𝛼𝑆𝐻_𝑛(𝑣)=−𝑥−0.06𝑡^𝛼𝑛Γ(𝑛𝛼+1)+max𝑥−25𝑒^𝑥,0−0.06−𝑡^𝛼𝑛Γ().𝑛𝛼+1 (44) So that the solution 𝑣(𝑥,𝑡) of the problem is given by 𝑣(𝑥,𝑡)=lim^∞_{𝑝→1𝑖=0}𝑝^𝑖𝑢_𝑖(𝑥,𝑡)=𝑥1−𝐸_𝛼−0.06𝑡^𝛼+max𝑥−25𝑒^−0.06,0𝐸_𝛼−0.06𝑡^𝛼. (45) This is the exact solution of the given option pricing equation (38). The solution of (38) at the special case 𝛼=1 is 𝑣(𝑥,𝑡)=𝑥1−𝑒^0.06𝑡−1,0+max𝑥−25𝑒^−0.06,0𝑒^−0.06𝑡. (46) The behaviour of 𝑣(𝑥,𝑡) with respect to 𝑥 and 𝑡 when 𝛼=1 is given in Figure 2.

Figure 2 

The surface shows the 𝑣(𝑥,𝑡) for (38) with respect to 𝑥 and 𝑡 for 𝛼=1.

7. Conclusion

In this paper, the homotopy perturbation Sumudu transform method (HPSTM) is successfully applied for getting the analytical solution of the fractional Black-Scholes option pricing equation. Two examples from the literature [26] are presented. The results of the illustrated examples are in agreement with the results of the method presented in [26]. In conclusion, HPSTM is a very powerful and efficient method to find approximate solutions as well as numerical solutions.

Acknowledgment

The authors express their sincere thanks to the referees for the careful and noteworthy reading of the paper, and the very helpful suggestions that improved the paper substantially.