Solution of the Fractional Black-Scholes Option Pricing Model by Finite Difference Method

Abstract

This work deals with the put option pricing problems based on the time-fractional Black-Scholes equation, where the fractional derivative is a so-called modified Riemann-Liouville fractional derivative. With the aid of symbolic calculation software, European and American put option pricing models that combine the time-fractional Black-Scholes equation with the conditions satisfied by the standard put options are numerically solved using the implicit scheme of the finite difference method.

1. Introduction

The theory and methodology of partial differential equation started to become popular to study option pricing problems, after the classical Black-Scholes equation was proposed. The available research results include mainly two aspects: one is to give values of options using more powerful numerical and analytic methods; the other is to derive new pricing models that reflect the actual financial market more closely.

The Black-Scholes equation has been increasingly attracting interest over the last two decades since it provides effectively the values of options. But the classical Black-Scholes equation was established under some strict assumptions. Therefore, some improved models have been proposed to weaken these assumptions, such as stochastic interest model [1], Jump-diffusion model [2], stochastic volatility model [3], and models with transactions costs [4, 5]. With the discovery of the fractal structure for financial market, the fractional Black-Scholes models [6–9] are derived by replacing the standard Brownian motion involved in the classical model with fractional Brownian motion. These fractional Black-Scholes models are still partial differential equations with integer order derivative, which are reduced to the classical Black-Scholes equation if we let the Hurst exponent 𝐻=1/2. The differential equation involving derivatives of fractional order is a powerful tool for studying fractal geometry and fractal dynamics. As a generalization of the integer-order differential equation, fractional differential equation is used to model important phenomena in various fields such as fluid flow, electromagnetic, acoustics, electrochemistry, cosmology, and material science. Recently, fractional partial differential equation was introduced more and more into financial theory. Wyss [10] gave the fractional Black-Scholes equation with a time-fractional derivative to price European call option. Cartea and del-Castillo-Negrete [11] gave several fractional diffusion models of option prices in markets with jumps and priced barrier option using fractional partial differential equation. Jumarie [12, 13] derived the time- and space-fractional Black-Scholes equations and gave optimal fractional Merton's portfolio. The aim of this paper is to try to combine Jumarie's time-fractional Black-Scholes equation with the terminal and boundary conditions satisfied by the standard put options to study the pricing problems for European and American put options.

In the present work, the option price 𝑉=𝑉(𝑆,𝑡) is suggested to be subject to the time-fractional Black-Scholes equation [12, 13] with the following form: 𝜕^𝛼𝑉𝜕𝑡^𝛼=𝑟𝑉−𝑟𝑆𝜕𝑉𝑡𝜕𝑆Γ(^1−𝛼)−Γ()2−𝛼1+𝛼2𝜎²𝑆²𝜕²𝑉𝜕𝑆²,𝑡>0,0<𝛼≤1, (1) where 𝑟,𝜎 denote risk-free interest rate and volatility, respectively. 𝜕^𝛼𝑉/𝜕𝑡^𝛼 involved in (1) is a modified Riemann-Liouville fractional derivative [12–15], which is defined using the equality 𝐷^𝛼_𝑡𝑓(𝑡)=lim_ℎ→0Δ^𝛼𝑓(𝑡)ℎ^𝛼=lim^∞_{ℎ→0𝑘=0}(−1)^𝑘𝛼𝑘𝑓(𝑡+(𝛼−𝑘)ℎ),0<𝛼≤1 (2) and can be represented as 𝐷^𝛼_𝑡𝑓(𝑡)=1Γ()∫−𝛼^𝑡0(𝑡−𝜉)^−𝛼−1𝑓(𝜉)(0)1−𝑓𝑑𝜉,if𝛼<0,Γ()𝑑1−𝛼∫𝑑𝑡^𝑡0(𝑡−𝜉)^−𝛼𝑓(𝜉)(0)𝑓−𝑓𝑑𝜉,if0<𝛼<1,_{(𝛼−𝑛)}(𝑡)^(𝑛),if𝑛≤𝛼<𝑛+1,𝑛≥1. (3)

The value of European put option is taken as a solution of (1) with the following terminal and boundary conditions: 𝑉()(),𝑉()((,𝑉()𝑆,𝑇=max𝐾−𝑆,00,𝑡=𝐾exp−𝑟𝑇−𝑡))𝑆,𝑡⟶0,𝑆⟶+∞. (4)

Pricing European put option based on (1) and (4) can be implemented by finite differences approximation. For American put option, this idea is operable. The main problem for pricing an American put option is to consider the possibility of early exercise. To avoid arbitrage, the option value at each point in the (𝑆,𝑡) space cannot be less than the intrinsic value [16]. For that, the principle of dynamic programming (at a given time, the optimal strategy corresponds to the maximum of either the exercise value or the value associated with selecting an optimal strategy an instant later) [17] expressed by the following equation is applied, which can be implemented using finite difference method and successive overrelaxation (SOR) method 𝑉(𝑆,𝑡)=max(𝐾−𝑆,𝑉(𝑆+𝑑𝑆,𝑡+𝑑𝑡)). (5)

Fractional Black-Scholes equation (1) is actually a special kind of fractional advection-dispersion equation, which has time-varying coefficient and time-fractional derivative. A general fractional advection-dispersion equation has been studied in many works in the literature. Sousa [18] employed a second-order explicit finite difference method derived by a Lax-Wendroff-type time discretization procedure to solve space-fractional advection-diffusion equation. Mohebbi and Abbaszadeh [19] applied compact finite difference scheme to deal with time-fractional advection-dispersion equation with constant coefficient. Meerschaert and Tadjeran [20] developed explicit and implicit Euler methods for the space-fractional advection-dispersion equation. Liu et al. [21] developed effective numerical methods for solving several space-time fractional advection-dispersion equations. Liu et al. [22] used explicit and implicit difference methods to solve space-time fractional advection-dispersion equation and proved the stability and convergence of the methods. More powerful methods, such as implicit MLS meshless method [23], Adomian decomposition method [24], homotopy perturbation method [25], and spectral regularization method [26], have been used to solve numerical and analytical solutions for various types of fractional advection-dispersion equation. In this paper, we combine time-fractional Black-Scholes equation (1) with conditions satisfied by European and American put options to construct put option pricing models. It is obvious that this fractional derivative model can be reduced to the classical Black-Scholes model if the derivative order 𝛼=1. As a generalization, this time-fractional pricing model of European and American put options is solved numerically using the implicit finite difference technique.

2. Numerical Scheme

Numerical and analytical methods [16–43] are used to study a great quantity of differential equations. Among them, the finite difference method is a direct and effective numerical algorithm. With this technique, the differential equation is transformed into a difference equation by the discretization of derivatives, and numerical solutions are finally obtained. The approach has been extended successfully to deal with various fractional differential equation [18–22, 38–43]. In this section, the implicit finite difference mode is given.

Take the change of variable 𝜏=𝑇−𝑡. (6) Equation (1) can be transformed into the following form: 𝜏^(𝛼−1)(𝑇−𝜏)^(1−𝛼)𝜕^𝛼𝑉𝜕𝜏^𝛼+𝑟𝑉−𝑟𝑆𝜕𝑉()𝜕𝑆𝑇−𝜏^1−𝛼Γ()−Γ()2−𝛼1+𝛼2𝜎²𝑆²𝜕²𝑉𝜕𝑆²=0. (7)

In order to use finite difference approximation, we start by 𝑆∈[0,𝑆_max] and 𝜏∈[0,𝑇). Let ℎ=𝑆_max/𝑀, 𝑘=𝑇/𝑁 (𝑀,𝑁∈𝑍₊) be the grid sizes in space and time; the computational domain is discretized by a uniform grid (𝑆_𝑚,𝜏_𝑛) with 𝑆_𝑚=𝑚ℎ (𝑚=0,1,…,𝑀) and 𝜏_𝑛=𝑛𝑘 (𝑛=0,1,…,𝑁), where 𝑆_𝑀=𝑆_max is a realistic and practical approximation to infinity. 𝑉^𝑛_𝑚 denotes an approximate solution of (7) in 𝑆_𝑚 at the time 𝜏_𝑛.

As pointed out in [12, 13], Jumarie's definition and the so-called Caputo's definition yield the same result when the function is differentiable. So, we take Caputo finite difference approximation [38] for the modified Riemann-Liouville fractional derivative involved in (7), namely, 𝜕^𝛼𝑉𝜕𝜏^𝛼=𝜚^𝑛_{𝛼,𝑘𝑗=1}𝜔_𝑗^(𝛼)𝑉_𝑚^{𝑛−𝑗+1}−𝑉_𝑚^𝑛−𝑗+𝑂(𝑘), (8) where 𝜚_𝛼,𝑘=1/Γ(2−𝛼)𝑘^𝛼 and 𝜔_𝑗^(𝛼)=𝑗^1−𝛼−(𝑗−1)^1−𝛼 that satisfies (i) 𝜔_𝑗^(𝛼)>0, (ii) 𝜔_𝑗^(𝛼)>𝜔^(𝛼)_𝑗+1 [22, 38].

For spatial derivative, we use the following difference approximation: 𝜕𝑉=𝑉𝜕𝑆^𝑛𝑚+1−𝑉^𝑛_𝑚−1(ℎ),𝜕2ℎ+𝑂₂𝑉𝜕𝑆²=𝑉^𝑛_𝑚+1−2𝑉^𝑛_𝑚+𝑉^𝑛_𝑚−1ℎ²+𝑂ℎ². (9)

Substituting (8) and (9) to (7) can derive 𝜏^(𝛼−1)(𝑇−𝜏)^(1−𝛼)𝜚^𝑛_{𝛼,𝑘𝑗=1}𝜔_𝑗^(𝛼)𝑉_𝑚^{𝑛−𝑗+1}−𝑉_𝑚^𝑛−𝑗+𝑟𝑉^𝑛_𝑚𝑉−𝑟𝑚ℎ^𝑛𝑚+1−𝑉^𝑛_𝑚−1×()2ℎ𝑇−𝑛𝑘^1−𝛼Γ()−Γ()2−𝛼1+𝛼2×𝜎²(𝑚ℎ)²𝑉^𝑛_𝑚+1−2𝑉^𝑛_𝑚+𝑉^𝑛_𝑚−1ℎ²=0 (10) or 𝜏^(𝛼−1)(𝑇−𝜏)^(1−𝛼)𝜚^𝑛_{𝛼,𝑘𝑗=1}𝜔_𝑗^(𝛼)𝑉_𝑚^{𝑛−𝑗+1}−𝑉_𝑚^𝑛−𝑗+𝑟𝑉^𝑛_𝑚𝑉−𝑟𝑚^𝑛𝑚+1−𝑉^𝑛_𝑚−12×(𝑇−𝑛𝑘)^1−𝛼Γ()−Γ()2−𝛼1+𝛼2×𝜎²𝑚²𝑉^𝑛_𝑚+1−2𝑉^𝑛_𝑚+𝑉^𝑛_𝑚−1=0, (11) where the local truncation error is 𝑂(𝑘+ℎ²).

From (11), we get the following equality when 𝑛=1: 𝑘^(𝛼−1)(𝑇−𝑘)^(1−𝛼)𝜚_𝛼,𝑘𝑉¹_𝑚−𝑉⁰_𝑚+𝑟𝑉¹_𝑚𝑉−𝑟𝑚^1𝑚+1−𝑉¹_𝑚−12×(𝑇−𝑘)^1−𝛼Γ()−Γ()2−𝛼1+𝛼2×𝜎²𝑚²𝑉¹_𝑚+1−2𝑉¹_𝑚+𝑉¹_𝑚−1=0. (12) After combining like terms, we can derive the following equation: 𝑟𝑚(𝑇−𝑘)^1−𝛼()−Γ()2Γ2−𝛼1+𝛼2𝜎²𝑚²𝑉¹_𝑚−1+𝑘^(𝛼−1)(𝑇−𝑘)^(1−𝛼)𝜚_𝛼,𝑘+𝑟(𝑇−𝑘)^1−𝛼Γ()()𝜎2−𝛼+Γ1+𝛼²𝑚²𝑉¹_𝑚−Γ(1+𝛼)2𝜎²𝑚²+𝑟𝑚(𝑇−𝑘)^1−𝛼()𝑉2Γ2−𝛼^1𝑚+1=𝑘^(𝛼−1)(𝑇−𝑘)^(1−𝛼)𝜚_𝛼,𝑘𝑉⁰_𝑚,𝑚=1,2,…,𝑀−1. (13) For 𝑛≥2, we obtain 𝑟𝑚(𝑇−𝑛𝑘)^1−𝛼()−Γ()2Γ2−𝛼1+𝛼2𝜎²𝑚²𝑉^𝑛_𝑚−1+(𝑛𝑘)^(𝛼−1)(𝑇−𝑛𝑘)^(1−𝛼)𝜚_𝛼,𝑘+𝑟(𝑇−𝑛𝑘)^1−𝛼Γ()()𝜎2−𝛼+Γ1+𝛼²𝑚²𝑉^𝑛_𝑚−Γ(1+𝛼)2𝜎²𝑚²+𝑟𝑚(𝑇−𝑛𝑘)^1−𝛼()𝑉2Γ2−𝛼^𝑛𝑚+1=(𝑛𝑘)^(𝛼−1)(𝑇−𝑛𝑘)^(1−𝛼)𝜚_𝛼,𝑘×𝑉_𝑚^𝑛−1−^𝑛_𝑗=2𝜔_𝑗^(𝛼)𝑉_𝑚^{𝑛−𝑗+1}−𝑉_𝑚^𝑛−𝑗,𝑚=1,2,…,𝑀−1. (14)

If the minimum limit of variable 𝑆 is 𝑆_min, the corresponding implicit difference approximation can be derived by replacing 𝑚 of the above results (13) and (14) with (𝑆_min/ℎ)+𝑚.

From the terminal and boundary conditions of the European put option, we can get 𝑉⁰_𝑚()𝑉=max𝐾−𝑚ℎ,0,𝑚=0,…,𝑀,_𝑛0()𝑉=𝐾exp−𝑟𝑛𝑘,𝑛=1,…,𝑁,_𝑛𝑀=0,𝑛=1,…,𝑁. (15)

In the case of American put option, we perform the above procedure. At the same time, we should check for the possibility of early exercise after computing 𝑉^𝑛_𝑚 and set 𝑉^𝑛_𝑚=max𝑉^𝑛_𝑚,𝐾−𝑚ℎ,for𝑛=1,…,𝑁. (16) Under this implicit scheme, we do not directly apply (16) at each step but use SOR method to complete this procedure, which was suggested in [16].

3. Stability and Convergence

In this section, we will analyze the stability and convergence of implicit finite difference scheme (13) and (14) using Fourier analysis involved in [39–41]. For that, (13) and (14) are rewritten as 𝑎_𝑚,1𝑉¹_𝑚−1+𝑏_𝑚,1𝑉¹_𝑚+𝑐_𝑚,1𝑉¹_𝑚+1=𝑉⁰_𝑚,𝑎_𝑚,𝑛𝑉^𝑛_𝑚−1+𝑏_𝑚,𝑛𝑉^𝑛_𝑚+𝑐_𝑚,𝑛𝑉^𝑛_𝑚+1=1−𝜔₂^(𝛼)𝑉_𝑚^𝑛−1−^𝑛−1_𝑗=2𝜔^(𝛼)_𝑗+1−𝜔_𝑗^(𝛼)𝑉_𝑚^𝑛−𝑗+𝜔_𝑛^(𝛼)𝑉⁰_𝑚, (17) where 𝑎_𝑚,𝑛=(𝑛𝑘)^(1−𝛼)(𝑇−𝑛𝑘)^(𝛼−1)𝜚⁻¹_𝛼,𝑘×𝑟𝑚(𝑇−𝑛𝑘)^1−𝛼()−Γ()2Γ2−𝛼1+𝛼2𝜎²𝑚²,𝑏_𝑚,𝑛=1+(𝑛𝑘)^(1−𝛼)(𝑇−𝑛𝑘)^(𝛼−1)𝜚⁻¹_𝛼,𝑘×𝑟(𝑇−𝑛𝑘)^1−𝛼Γ()()𝜎2−𝛼+Γ1+𝛼²𝑚²,𝑐_𝑚,𝑛=−(𝑛𝑘)^(1−𝛼)(𝑇−𝑛𝑘)^(𝛼−1)𝜚⁻¹_𝛼,𝑘×Γ(1+𝛼)2𝜎²𝑚²+𝑟𝑚(𝑇−𝑛𝑘)^1−𝛼().2Γ2−𝛼 (18)

3.1. Stability Analysis

If 𝑈^𝑛_𝑚 is another approximate solution of (13) and (14), we define the round-off error 𝜀^𝑛_𝑚=𝑉^𝑛_𝑚−𝑈^𝑛_𝑚 which satisfies the following equations according to (17): 𝑎_𝑚,1𝜀¹_𝑚−1+𝑏_𝑚,1𝜀¹_𝑚+𝑐_𝑚,1𝜀¹_𝑚+1=𝜀⁰_𝑚,𝑎_𝑚,𝑛𝜀^𝑛_𝑚−1+𝑏_𝑚,𝑛𝜀^𝑛_𝑚+𝑐_𝑚,𝑛𝜀^𝑛_𝑚+1=1−𝜔₂^(𝛼)𝜀_𝑚^𝑛−1−^𝑛−1_𝑗=2𝜔^(𝛼)_𝑗+1−𝜔_𝑗^(𝛼)𝜀_𝑚^𝑛−𝑗+𝜔_𝑛^(𝛼)𝜀⁰_𝑚,𝜀𝑛=2,…,𝑁,𝑚=1,2,…,𝑀−1,_𝑛0=𝜀^𝑛_𝑀=0,𝑛=0,…,𝑁. (19)(20)

Now, we introduce the grid function in [39]: 𝜀^𝑛(𝑆)=𝜀^𝑛_𝑚,when𝑆_𝑚−ℎ2<𝑆≤𝑆_𝑚+ℎ2ℎ,𝑚=1,…,𝑀−1,0,when0≤𝑆≤2or𝑆_max−ℎ2<𝑆≤𝑆_max. (21)

Recall the result of [39]; defining 𝜀^𝑛=𝜀^𝑛₁,𝜀^𝑛₂,…,𝜀^𝑛_𝑀−1^𝑇,𝜀^𝑛₂=^𝑀−1_𝑚=1ℎ𝜀^𝑛_𝑚^21/2=^𝑆max₀𝜀^𝑛(𝑆)²𝑑𝑆^1/2 (22) can get 𝜀^𝑛2₂=^+∞_𝑙=−∞𝜉_𝑛(𝑙)², (23) where 𝜉_𝑛(𝑙)=(1/𝑆_max)∫^𝑆max₀𝜀^𝑛(𝑆)exp(−𝐼2𝜋𝑙𝑆/𝑆_max)𝑑𝑆.

The solution of (19) and (20) is supposed to have the following form: 𝜀^𝑛_𝑚=𝜉_𝑛exp𝐼𝜇𝑚ℎ, (24) where 𝜇=2𝜋𝑙/𝑆_max, 𝐼=√−1.

Substituting the above expression into (19), we get 𝑎_𝑚,1𝜉₁exp𝐼𝜇(𝑚−1)ℎ+𝑏_𝑚,1𝜉₁exp𝐼𝜇𝑚ℎ+𝑐_𝑚,1𝜉₁exp𝐼𝜇(𝑚+1)ℎ=𝜉₀,𝑎exp𝐼𝜇𝑚ℎ_𝑚,𝑛𝜉_𝑛exp𝐼𝜇(𝑚−1)ℎ+𝑏_𝑚,𝑛𝜉_𝑛exp𝐼𝜇𝑚ℎ+𝑐_𝑚,𝑛𝜉_𝑛()ℎ=exp𝐼𝜇𝑚+11−𝜔₂^(𝛼)𝜉_𝑛−1−exp𝐼𝜇𝑚ℎ_{𝑛−1𝑗=2}𝜔^(𝛼)_𝑗+1−𝜔_𝑗^(𝛼)𝜉_𝑛−𝑗exp𝐼𝜇𝑚ℎ+𝜔_𝑛^(𝛼)𝜉₀exp𝐼𝜇𝑚ℎ. (25) Simplifying (25), we obtain 𝜗_𝑚,1𝜉₁=𝜉₀,𝜗_𝑚,𝑛𝜉_𝑛=1−𝜔₂^(𝛼)𝜉_𝑛−1−^𝑛−1_𝑗=2𝜔^(𝛼)_𝑗+1−𝜔_𝑗^(𝛼)𝜉_𝑛−𝑗+𝜔_𝑛^(𝛼)𝜉₀, (26) where 𝜗_𝑚,𝑛=1+(𝑛𝑘)^(1−𝛼)(𝑇−𝑛𝑘)^(𝛼−1)𝜚⁻¹_𝛼,𝑘[(𝑟(𝑇−𝑛𝑘)^1−𝛼/Γ(2−𝛼))(1−𝐼𝑚sin(𝜇ℎ))+Γ(1+𝛼)𝜎²𝑚²(1−  cos(𝜇ℎ))].

Proposition 1. If 𝜉_𝑛 is a solution of (26), then |𝜉_𝑛|≤|𝜉₀|.

Proof. For 𝑛=1, the first equality of (26) gives |𝜉₁|=(1/|𝜗_𝑚,1|)|𝜉₀|≤|𝜉₀|. If 𝜉_𝑛−1≤𝜉₀, (27) then using the second equality of (26), we obtain 𝜉_𝑛≤1−𝜔₂^(𝛼)−𝜔_𝑛^(𝛼)−𝜔₂^(𝛼)+𝜔_𝑛^(𝛼)𝜗_𝑚,𝑛𝜉₀=1𝜗_𝑚,𝑛𝜉₀≤𝜉₀. (28)

According to (13) and (14), we can obtain the following conclusion using Proposition 1 and equality (23).

Theorem 2. The difference scheme (13) and (14) is unconditionally stable.

3.2. Convergence Analysis

Suppose that 𝑉(𝑆_𝑖,𝑡_𝑘) is an exact solution of (7) at grid point (𝑆_𝑖,𝑡_𝑘) and 𝑉^𝑖_𝑘 is the difference solution of (13) and (14); we define the error 𝜖^𝑘_𝑗=𝑉(𝑆_𝑖,𝑡_𝑘)−𝑉^𝑖_𝑘 which satisfies the following equation according to (17): 𝑎_𝑚,1𝜖¹_𝑚−1+𝑏_𝑚,1𝜖¹_𝑚+𝑐_𝑚,1𝜖¹_𝑚+1=𝑘^𝛼𝑅¹_𝑚,𝑎_𝑚,𝑛𝜖^𝑛_𝑚−1+𝑏_𝑚,𝑛𝜖^𝑛_𝑚+𝑐_𝑚,𝑛𝜖^𝑛_𝑚+1=1−𝜔₂^(𝛼)𝜖_𝑚^𝑛−1−^𝑛−1_𝑗=2𝜔^(𝛼)_𝑗+1−𝜔_𝑗^(𝛼)𝜖_𝑚^𝑛−𝑗+𝑘^𝛼𝑅^𝑛_𝑚,𝜖𝑛=2,…,𝑁,𝑚=1,2,…,𝑀−1,_𝑛0=𝜖^𝑛_𝑀𝜖=0,𝑛=1,…,𝑁,_0𝑚=0,𝑚=0,…,𝑀, (29)(30) where |𝑅^𝑛_𝑚|≤𝐶₁(𝑘+ℎ²) (𝑛=1,2,…,𝑁) and 𝐶₁ is a positive constant.

For completing the proof of convergence, we recall several results below which came from [39].

Similar to the stability analysis, [39] constructed the following grid function: 𝜖^𝑛(𝑆)=𝜖^𝑛_𝑚,when𝑆_𝑚−ℎ2<𝑆≤𝑆_𝑚+ℎ2ℎ,𝑚=1,…,𝑀−1,0,when0≤𝑆≤2or𝑆_max−ℎ2<𝑆≤𝑆_max,𝑅^𝑛(𝑆)=𝑅^𝑛_𝑚,when𝑆_𝑚−ℎ2<𝑆≤𝑆_𝑚+ℎ2ℎ,𝑚=1,…,𝑀−1,0,when0≤𝑆≤2or𝑆_max−ℎ2<𝑆≤𝑆_max. (31)

Defining 𝜖^𝑛=𝜖^𝑛₁,𝜖^𝑛₂,…,𝜖^𝑛_𝑀−1^𝑇,𝑅^𝑛=𝑅^𝑛₁,𝑅^𝑛₂,…,𝑅^𝑛_𝑀−1^𝑇 (32) with the norms 𝜖^𝑛₂=^𝑀−1_𝑚=1ℎ𝜖^𝑛_𝑚^21/2=^𝑆max₀𝜖^𝑛(𝑆)²𝑑𝑆^1/2,𝑅^𝑛₂=^𝑀−1_𝑚=1ℎ𝑅^𝑛_𝑚^21/2=^𝑆max₀𝑅^𝑛(𝑆)²𝑑𝑆^1/2 (33) can obtain 𝜖^𝑛2₂=^+∞_𝑙=−∞𝜁_𝑛(𝑙)²,𝑅^𝑛2₂=^+∞_𝑙=−∞𝜂_𝑛(𝑙)², (34) where 𝜁_𝑛(𝑙)=1𝑆_max^𝑆max₀𝜖^𝑛(𝑆)exp−𝐼2𝜋𝑙𝑆𝑆_max𝜂𝑑𝑆,_𝑛(𝑙)=1𝑆_max^𝑆max₀𝑅^𝑛(𝑆)exp−𝐼2𝜋𝑙𝑆𝑆_max𝑑𝑆. (35)

We can let 𝜖^𝑛_𝑚=𝜁_𝑛,𝑅exp𝐼𝜇𝑚ℎ_𝑛𝑚=𝜂_𝑛exp𝐼𝜇𝑚ℎ, (36) where 𝜇=2𝜋𝑙/𝑆_max, 𝐼=√−1.

Substituting the above expressions into (29) and simplifying them can get 𝜗_𝑚,1𝜁₁=𝑘^𝛼𝜂₁,𝜗_𝑚,𝑛𝜁_𝑛=1−𝜔₂^(𝛼)𝜁_𝑛−1−^𝑛−1_𝑗=2𝜔^(𝛼)_𝑗+1−𝜔_𝑗^(𝛼)𝜁_𝑛−𝑗+𝑘^𝛼𝜂_𝑛, (37) where 𝜗_𝑚,𝑛=1+(𝑛𝑘)^(1−𝛼)(𝑇−𝑛𝑘)^(𝛼−1)𝜚⁻¹_𝛼,𝑘([(𝑟(𝑇−𝑛𝑘)^1−𝛼/Γ(2−𝛼))(1−𝐼𝑚sin(𝜇ℎ))+Γ(1+𝛼)𝜎²𝑚²(1−  cos(𝜇ℎ))].

Proposition 3. There exists a positive constant 𝐶₂, so that |𝜁_𝑛|≤𝐶₂𝑛𝑘^𝛼|𝜂₁|.

Proof. From [39], we know 𝑅^𝑛₂≤𝐶₁𝑆_max𝑘+ℎ²()𝜂𝑛=1,2,…,𝑁,_𝑛≤𝐶₂𝜂₁(𝑛=1,2,…,𝑁). (38)(39) Then the first equality of (37) tells us that 𝜁₁≤1𝜗_𝑚,𝑛𝑘^𝛼𝜂₁≤𝐶₂𝑘^𝛼𝜂₁. (40) If we let 𝜁_𝑛−1≤𝐶₂(𝑛−1)𝑘^𝛼𝜂₁, (41) using the second equality of (37), we can obtain 𝜁_𝑛≤1𝜗_𝑚,𝑛1−𝜔₂^(𝛼))+(𝑛−1_{𝑛−1𝑗=2}𝜔_𝑗^(𝛼)−𝜔^(𝛼)_𝑗+1𝑛−𝑗+1𝐶₂𝑘^𝛼𝜂₁≤1𝜗_𝑚,𝑛1−𝜔₂^(𝛼)()+𝑛−1_{𝑛−1𝑗=2}𝜔_𝑗^(𝛼)−𝜔^(𝛼)_𝑗+1(𝑛−1)+1𝐶₂𝑘^𝛼𝜂₁≤1𝜗_𝑚,𝑛1−𝜔_𝑛^(𝛼)(𝑛−1)+1𝐶₂𝑘^𝛼𝜂₁≤𝑛𝐶₂𝑘^𝛼𝜂₁. (42)

According to equalities (34) and inequality (38), we can obtain the following conclusion using Proposition 3.

Theorem 4. The implicit difference scheme (13) and (14) is 𝐿₂-convergent.

4. Computational Examples

Option pricing model based on the time-fractional differential equation (1) is studied. The implicit difference scheme of (7) is given in Section 2. From the finite difference forms, it is clear that (13) and (14) are just the implicit difference scheme of the classical Black-Scholes equation if we let 𝛼=1. This known implicit difference results are given in [16, 17, 36]. In this section, we focus mainly on investigating the results with 0<𝛼<1. To explain the stability and convergence of the implicit numerical schemes, we firstly take Example 1 with the given terminal and boundary condition as an example.

Example 1. Consider (1) that is subject to the following conditions: 𝑉(𝑆,𝑇)=𝑉(𝑆,𝑇)=𝜎²2𝑇−𝑇^𝛼(𝑆),𝑉𝑆−2ln_min,𝑡=exp(−𝑟(𝑇−𝑡))𝜎²2𝑡−𝑡^𝛼−2ln𝑆_min,𝑉𝑆_max,𝑡=exp(−𝑟(𝑇−𝑡))𝜎²2𝑡−𝑡^𝛼−2ln𝑆_max. (43)
Under the condition 𝑟=𝜎², (1) has an analytical solution 𝑉(𝑆,𝑡)=exp(−𝑟(𝑇−𝑡))𝜎²2𝑡−𝑡^𝛼−2ln(𝑆). (44)

Figure 1 shows analytical solution and numerical solution obtained by the implicit difference method at time 𝑡=0.5 when 𝑟=0.01, 𝜎=0.1, 𝑇=1, 𝑆_min=0.1, 𝑆_max=1, 𝑀=19, and 𝑁=10. Numerical solution compares well with analytical solution, which proves that the implicit scheme is stable. Under the same parameters, Figure 2 gives the absolute error between numerical solutions and analytical solutions, which illustrates that the numerical results are convergent.

(a)

(b)

Figure 1 

Solution curve with 𝑡=0.5. Solid line: numerical solution; Soft dot: analytical solution.

(a)

(b)

Figure 2 

Absolute error between numerical solution and analytical solution. dash-dotted line: 𝑡=0.3; dotted line: 𝑡=0.5; solid line: 𝑡=0.8.

Example 2. European put option pricing model is based on (1) and the condition (4) under the following parameters: 𝑆𝐾=50,𝑟=0.01,𝜎=0.3,𝑇=1,_max=100,𝑀=20,𝑁=10. (45) Through coding, the European put option values with 𝛼=1,1/2,5/7,9/10 are plotted in Figure 3 at 𝑡=0.1/0.5/1.

(a)

(b)

(c)

(d)

Figure 3 

European put option. dash-dotted line: 𝑡=0.1; dotted line: 𝑡=0.5; solid line: 𝑡=1.

Example 3. Consider American put option pricing model. The following parameters are selected for the present study: 𝑆𝐾=60,𝑟=0.01,𝜎=0.4,𝑇=1,_max=120,𝑀=30,𝑁=10. (46) Figure 4 indicates a price comparison of the American put option at 𝑡=0.2/0.6/1 (relaxation parameter and tolerance parameter are 1.2 and 0.001 in the applied process of SOR method).

(a)

(b)

(c)

(d)

Figure 4 

American put option. dash-dotted line: 𝑡=0.2; dotted line: 𝑡=0.6; solid line: 𝑡=1.

Using Examples 1, 2, and 3, we examine the implementation of the implicit finite difference method for the fractional partial differential equation system. According to Figures 1 and 2 in Example 1, we can confirm that the implicit numerical scheme is stable and convergent. From Figures 3 and 4, we can see that the numerical scheme is very effective. Figures 3 and 4 present numerical simulation of the price of the European and American put options when the order of the time-fractional derivative takes different values. Their visible shapes and development trend are similar to the classical put option pricing model based on the standard Black-Scholes equation, which illustrates the essential characteristics of the European and American put options. For making the figures clear, we choose the small values of 𝑀 and 𝑁 and get the similar conclusions if largening the number of steps in time and space. As a generalization of the standard models, these fractional Black-Scholes models are powerful and will be of great interest to researchers in further work.

5. Conclusions

In this work, the finite difference method is employed to solve the time-fractional Black-Scholes equation together with the conditions satisfied by the standard put options. Application of the fractional differential equation to the pricing theory of option is in its beginning stage and needs more further work. This fractional model mentioned in this paper can model the price of other financial derivatives like warrant, swaps, and so on. The successful application of the finite difference method proves that this technique is effective and requires less computational work to solve fractional partial differential equation.

Acknowledgments

The authors would like to express their gratitude to reviewers for the careful reading of the paper and for their constructive comments which greatly improve the quality of this paper. The work is partially supported by the National Natural Science Foundation of China (no. 71171035, no. 71273044, and no. 71271045).