# Asymptotic Convergence Analysis and Error Estimate for Black-Scholes Model of Option Pricing

Authors: Juan He; Wei Tu; Aiqing Zhang
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this work, we discuss the numerical method for the solution of the Black-Scholes model. First of all, the asymptotic convergence for the solution of Black-Scholes model is proved. Second, we develop a linear, unconditionally stable, and second-order time-accurate numerical scheme for this model. By using the finite difference method and Legendre-Galerkin spectral method, we construct a time and space discrete scheme. Finally, we prove that the scheme has second-order accuracy and spectral accuracy in time and space, respectively. Several numerical experiments further verify the convergence rate and effectiveness of the developed scheme.


---

## Body

## 1. Introduction

The Black-Scholes (B–S) model is a classical pricing model of European options, which was developed by the economists Black and Scholes [1]. Because of their special contribution to the model, Merton and Scholes were awarded the Nobel Prize in 1997. This model can effectively price financial derivatives such as stocks, currencies, and bonds. Merton [2, 3] constructed the theoretical framework of option pricing and gave corresponding examples to verify the theory. Nowadays, B–S model is playing an increasingly important role in derivative financial instruments; it has been widely used in actual trading, which has inspired many options, trading strategies, and hedging strategies [4, 5, 5–7].Based on the B–S theoretical framework, Duffie [8] proposed a formula of the market value of securities, in which the arbitrage-free value of derivative securities is obtained by solving partial differential equations. Dennis and Antonio [9] presented a new method to approximately assess price under the background of continuous time model. They discussed this method under various conditions, including option pricing with random fluctuation, Greek calculation, and term structure of interest rate. He [10] discussed the convergence of contingent claim price model from discrete time to continuous time. These results show that the contingent claim price and dynamic replication portfolio strategy converge to the continuous time limit.Kim et al. [11] considered the pricing of a European option on a new multiscale mixed structure of underlying assess price fluctuation, and the result indicates that the option price has an ideal modification to Black-Scholes formula. In addition, this kind of correction can bring significant improvement in fitting the surface of implied volatility through calibration exercises. Lai [12] studied the influence of time discretization on European option pricing. The correction and discrete-time rebalancing strategies caused by discrete transactions are reconsidered, and the higher-order terms are expanded by Taylor series, and the corresponding correction source terms are derived. Yousuf et al. [13] proposed a second-order exponential time differencing scheme to solve nonlinear Black-Scholes model with transaction cost.In this paper, we will study the asymptotic property of the solution of B–S model. In addition, we develop a time and space discrete scheme, where the time direction has second-order accuracy and the space direction has spectral accuracy. The stability and convergence of fully discrete scheme are also proved. Finally, numerical examples were shown to verify the accuracy of theoretical analysis.This paper is organized as follows. In Section2, the B–S model will be introduced. In Section 3, we will present second-order numerical method. The error estimates will be presented in Section 4. In Section 5, several numerical examples will be used to demonstrate the effectiveness of the fully discrete method. We will show some conclusions in the last section.


## 2. Black-Scholes Model

We consider the B–S model as follows:(1)∂VS,t∂t+12σ2S2∂2VS,t∂S2+rS∂VS,t∂S−rVS,t=0,S,t∈0,+∞×0,T,(2)V0,t=pt,V+∞,t=qt,(3)VS,T=wS,where V is the price of the option, S is the price of the underlying asset, r is the interest rate, t is the time, and σ is the volatility of the stock price.Remark 1.
The B–S model is actually a stochastic partial differential equation; this model has a closed solution (we should know that most partial differential equations do not have it). That is to say, the option price can be accurately expressed by a function V =f(S, t), and the option price can be directly calculated by substituting the value of the independent variables S,t. In fact, it is difficult to get the desired theoretical and numerical results by directly analyzing the above equations. Therefore, we need to make some equivalent transformations to transform the equation into a more general form.We sett:=T−τ; and VS,τ=uex,T−τ. Thus, we can convert model (1) into the following form:(4)∂τx,τ−α∂x2ux,τ−β∂xux,τ+rux,τ=0,x,τ∈0,+∞×0,T,where α=1/2σ2,β=r−α, and with the following boundary (barrier) and initial conditions,(5)u−∞,τ=pτ,u+∞,τ=qτ,ux,0=u0x,a<x<b.In order to solve the above model better, we consider constructing a numerical method on the finite intervala,b, and then the considered model becomes(6)∂τux,τ−α∂x2ux,τ−β∂xux,τ+rux,τ=0,x,τ∈a,b×0,T,(7)ua,τ=ub,τ=0,∂xua,τ=∂xub,τ=0,(8)ux,0=u0x,a<x<b.For the solution of models (6)–(8), we have the following asymptotic estimation results.Lemma 1.
Supposeuis the solution of (6)–(8), then we find(9)‖u‖2≤e−2rτ‖u0‖2.Proof.
Taking the inner product of equation (6) with u, we find(10)12ddτu2+α∂xu2+ru2=0.
Multiplying the above equation by2e2rτ,(11)e2rτddτu2+2re2rτu2≤0.
LetHτ=e2rτu2, then we find(12)ddτHτ≤0.
Integrating with respect to time from 0 toτ, we have(13)Hτ≤H0.
Multiplying both sides bye−2rτ, we find(14)u2≤e−2rτu02.
This concludes the proof.


## 3. Second-Order Numerical Method

LetK be a positive integer, Δt=T/K be the time step, and τn=nΔt,n=0,1,…,K be mesh point. Using Crank-Nicolson formula to time discretization, we can obtain the following time-discrete scheme:(15)un+1−unΔt−α∂x2un+1/2−β∂xun+1/2+run+1/2=0.Here, un+1/2=un+1+un/2,n≥0.We can prove the following unconditional energy stability theorem for scheme (15).Theorem 1.
The time discretization scheme (15) is unconditionally stable. It satisfies the following properties:(16)un≤u0,n=1,2,…,K.Proof.
Taking theL2 inner produce of the (15) with Δtun+1+un, we arrive at(17)un+12−un2+2Δtα∂xun+1/22+2Δtrun+1/22=0.
Giving up the last two terms of the left hand side, we have(18)un+12−un2≤0.
This yields (16).


## 4. Error Estimate

In this part, we consider Legendre-Galerkin spectral method for the time-discrete scheme (12). We will present some error estimates for full-discretization schemes in L2 norm. First, we denote SN is the Legendre polynomial space, and denote πN:L2Ω⟶SN is the L2-projection operator which satisfies as follows:(19)πNϕ−ϕ,ψ=0,∀ψ∈SN.We have the following estimate [14]:(20)ϕ−πNϕl≤CNl−mϕm,∀ϕ∈HmΩ,m>l≥0.Next, we begin to analyze the error estimates of the full-discrete scheme (23). We denote the truncation error as follows:(21)Rn+1/2:=u·,τn+1−u·,τnΔt−∂τu·,τn+1/2.We know thatRn+1/2 satisfies(22)Rn+1/2≤CΔτ2.We denote error functions:(23)e1n=πNuτn−uNn,e2n=uτn−πNuτn,eun=e1n+e2n=uτn−uNn.Then, we can develop the following full-discrete scheme:(24)uNn+1−uNnΔτ,ϕN−α∂x2uNn+1/2,ϕN−β∂xuNn+1/2,ϕN+ruNn+1/2,ϕN=0,ϕN∈SN.We now state the stability results for full-discrete scheme (23).Theorem 2.
LetuNnn=1M−1be solution of (23), then we derive that(25)uNn+12≤uN02.Proof.
TakingϕN=ΔtuNn+1+uNn, we find that(26)uNn+12−uNn2+2Δtα∂xuNn+1/22+2ΔtruNn+1/22=0.
Giving up two positive terms on the right hand side, we have(27)uNn+12≤uNn2.
Finally, we obtain the desired result (24).Theorem 3.
For the constructed numerical scheme (24), we have the following error estimate:(28)uτk−uNk≤CΔτ2+N1−m,k=0,1,…,K=T/Δτ.Proof.
Forn=0, equation (23) can be written as(29)1ΔτuN1−uN0,ϕN−α∂x2uN1/2,ϕN−β∂xuN1/2,ϕN+ruN1/2,ϕN=0.
Subtracting (28) from (7) at τ1, we note that(30)∂τu⋅,τ1/2−uN1−uN0Δτ=∂τu⋅,τ1/2−u⋅,τ1−u⋅,τ0Δτ+u⋅,τ1−u⋅,τ0Δτ−πNu⋅,τ1−πNu⋅,τ0Δτ+πNu⋅,τ1−πNu⋅,τ0Δτ−uN1−uN0Δτ=r1/2+1ΔτI−πNu⋅,τ1−u⋅,τ0+1Δτe11−e10,∂x2u⋅,τ1/2−∂x2uN1/2=∂x2u⋅,τ1/2−πN∂x2u⋅,τ1/2+πN∂x2u⋅,τ1/2−∂x2uN1/2=∂x2I−πNu⋅,τ1/2+∂x2e11/2,∂xu⋅,τ1/2−∂xuN1/2=∂xu⋅,τ1/2−πN∂xu⋅,τ1/2+πN∂xu⋅,τ1/2−∂xuN1/2=∂xI−πNu⋅,τ1/2+∂xe11/2,u⋅,τ1/2−uN1/2=u⋅,τ1/2−πNu⋅,τ1/2+πNu⋅,τ1/2−uN1/2=I−πNu⋅,τ1/2+e11/2.
Then, we have(31)e11−e10,ϕN+Δτα∂xe11/2,∂xϕN−Δτβ∂xe11/2,ϕN+Δτre11/2,ϕN=−ΔτR1/2,ϕN+πN−Iu·,τ1−u·,τ0,ϕN+ΔταπN−I∂xu·,τ1/2,∂xϕN−ΔτβI−πNu·,τ1/2,∂xϕN+ΔτrπN−Iu·,τ1/2,ϕN,ϕN∈SN.
SettingϕN=2e11, we have(32)2e112+2Δτα∂xe112+2Δτre112≤Δτ1rR1/22+re112+ΔταπN−I∂xu·,τ1/22+∂xe11/22+ΔτββαI−πNu·,τ1/22+αβ∂xe11/22.
Dropping some positive terms, we find(33)e112≤C1Δτ5+C2N2−2m.
Subtracting (23) from a reformulation of (7) at τn+1/2, we find(34)e1n+1−e1n,ϕN+Δτα∂xeun+1/2,∂xϕN−Δτβ∂xeun+1/2,ϕN+Δτre1n+1/2,ϕN=−ΔτRn+1/2,ϕN+πN−Iu·,τn+1−u·,τn,ϕN+ΔταπN−I∂xu·,τn+1/2,∂xϕN−Δτβ(I−πN)u·,τn+1/2,∂xϕN+Δτr(πN−I)u·,τn+1/2,ϕN.
LettingϕN=2e1n+1/2, we have(35)e1n+12−e1n2+2Δτα∂xe1n+1/22+2Δτre1n+1/22≤−2ΔτRn+1/2,e1n+1/2+2ΔταπN−I∂xu·,τn+1/2,∂xe1n+1/2−2ΔτβI−πNu·,τn+1/2,∂xe1n+1/2≤ΔτRn+1/22+e1n+1/22+ΔταπN−I∂xu·,τn+1/22+∂xe1n+1/22+ΔτββαI−πNu·,τn+1/22+αβ∂xe1n+1/22.
Thus, we have(36)e1n+12−e1n2≤C1Δτ5+C2ΔτN2−2m+C3Δτe1n+12.
Summing up forn=1,…,k, we find(37)e1k+12≤e112+CΔτ4+N2−2m+C3Δτ∑n=1ke1n+12.
By using discrete Gronwall lemma, we can get(38)e1k+12≤CΔτ4+N2−2m.
Note that(39)uτk−uk≤e1k+e2k,and using (15), we can get (27).


## 5. Numerical Examples

### 5.1. Verification of the Convergence Order

In this part, we will test the accuracy and validity of the full-discrete scheme (23). Actually, in the B–S model, researches usually choose α=0.5σ2;β=r−α. Set T=1,N=128 and u0=tanx+x,σ=0.1,r=0.2. In Table 1, we list L2 error and convergence order for different time step size. We can see that the time direction is obvious of second-order accuracy, while the space direction has a good convergence property.Table 1
TheL2, numerical errors, and time convergence order for various temporal steps.
ΔτL2−errorOrder0.020.03652.00450.010.00912.00110.0050.00232.00030.0019.005e-052.00000.00052.2580e-052.0000


### 5.2. Effect of Various Parameters

Next, we will test the asymptotic decay property of the solution. We fixN=50,Δτ=0.01,T=1,Λ=0,6,andr=0.1 and leave the initial value unchanged, and at the same time, let β constantly change. In Figure 1– 4, we list the change process of numerical solution with β=0.2,0.6,1.0,5.0, and we find that when β becomes larger, the solution gradually tends to 0.Figure 1
Change process of numerical solution withβ=0.2.Figure 2
Change process of numerical solution withβ=0.6.Figure 3
Change process of numerical solution withβ=1.0.Figure 4
Change process of numerical solution withβ=5.0.Fixu0=|tanx+x|. In the following numerical experiments, we also test the influence of σ and r on the numerical solution. In Figure 5–6, the variation process of numerical solution with σ and r is given. We find that the numerical solution experiences four peaks and valleys; at the same time, with the increase of R, the fluctuation of numerical solution gradually increases.Figure 5
Variation curve of numerical solution withr=0,0.20.4,0.6,0.8,σ=0.2.Figure 6
Variation curve of numerical solution withr=0.1,0.20.3,0.5,0.7,r=0.2.


## 5.1. Verification of the Convergence Order

In this part, we will test the accuracy and validity of the full-discrete scheme (23). Actually, in the B–S model, researches usually choose α=0.5σ2;β=r−α. Set T=1,N=128 and u0=tanx+x,σ=0.1,r=0.2. In Table 1, we list L2 error and convergence order for different time step size. We can see that the time direction is obvious of second-order accuracy, while the space direction has a good convergence property.Table 1
TheL2, numerical errors, and time convergence order for various temporal steps.
ΔτL2−errorOrder0.020.03652.00450.010.00912.00110.0050.00232.00030.0019.005e-052.00000.00052.2580e-052.0000


## 5.2. Effect of Various Parameters

Next, we will test the asymptotic decay property of the solution. We fixN=50,Δτ=0.01,T=1,Λ=0,6,andr=0.1 and leave the initial value unchanged, and at the same time, let β constantly change. In Figure 1– 4, we list the change process of numerical solution with β=0.2,0.6,1.0,5.0, and we find that when β becomes larger, the solution gradually tends to 0.Figure 1
Change process of numerical solution withβ=0.2.Figure 2
Change process of numerical solution withβ=0.6.Figure 3
Change process of numerical solution withβ=1.0.Figure 4
Change process of numerical solution withβ=5.0.Fixu0=|tanx+x|. In the following numerical experiments, we also test the influence of σ and r on the numerical solution. In Figure 5–6, the variation process of numerical solution with σ and r is given. We find that the numerical solution experiences four peaks and valleys; at the same time, with the increase of R, the fluctuation of numerical solution gradually increases.Figure 5
Variation curve of numerical solution withr=0,0.20.4,0.6,0.8,σ=0.2.Figure 6
Variation curve of numerical solution withr=0.1,0.20.3,0.5,0.7,r=0.2.


## 6. Concluding Remarks

We construct an effective fully discrete scheme for the B–S model based on Legendre-Galerkin scheme for spatial discretization. The scheme is also linear, unconditionally stable, and has second-order accuracy in time, where the Crank-Nicolson method is used for time discretization. Through the implementations of several numerical examples, we demonstrate the accuracy and effectiveness of the developed scheme, numerically.



---

*Source: 6563766-2022-03-27.xml*