# High-Accurate Numerical Schemes for Black–Scholes Models with Sensitivity Analysis

Authors: Samir Kumar Bhowmik; Jakobin Alam Khan
Journal: Journal of Mathematics
 (2022)

---

## Abstract

The significance of both the linear and nonlinear Black-Scholes partial differential equation model is huge in the field of financial analysis. In most cases, the exact solution to such a nonlinear problem is very hard to obtain, and in some cases, it is impossible to get an exact solution to such models. In this study, both the linear and the nonlinear Black-Scholes models are investigated. This research mainly focuses on the numerical approximations of the Black-Scholes (BS) model with sensitivity analysis of the parameters. It is to note that most applied researchers use finite difference and finite element-based schemes to approximate the BS model. Thus, an urge for a high accurate numerical scheme that needs fewer grids/nodes is huge. In this study, we aim to approximate and analyze the models using two such higher-order schemes. To be specific, the Chebyshev spectral method and the differential quadrature method are employed to approximate the BS models to see the efficiency of such highly accurate schemes for the option pricing model. First, we approximate the model using the mentioned methods. Then, we move on to use the numerical results to analyze different aspects of stock market through sensitivity analysis. All the numerical schemes have been illustrated through some graphics and relevant discussions. We finish the study with some concluding remarks.


---

## Body

## 1. Introduction

The Black-Scholes parabolic partial differential equation model is considered to be one of the major appliances for option pricing as it is efficient to compute the price of an option more correctly [1–5]. An option is a financial derivative, and it is a contract that gives the owner the right to buy or sell the asset for a specified price within a specified time [4, 6–9]. The fixed price is called the strike price or exercise price, and it is denoted by K; and the fixed time is called the expiration time or the maturity time, and it is denoted by T [3, 5, 7, 8]. There are mainly two types of options available in the derivative market, namely, call option and put option [8]. The call option gives the owner the right to buy the underlying real asset, and the put option gives the owner the right to sell the underlying real asset [3, 4, 7, 10].Though there are various types of options [8], the most commonly applied options are Asian option, European option, and American option [4, 7, 8, 11]. The option popularly known as European option can only be applicable at the expiration time T [7, 9, 10, 12, 13], whereas American option can be exercised at any time until the expiration time T [4, 7, 9, 10]. It is true that there is flexibility in American option. The average value of the underlying asset price is considered in Asian option [11]; it is an option that can be partitioned into two classes [8, 14]; and it contains less volatility [8]. In this paper, we work on European option.Fisher Black, Myron Scholes, and Robert Merton, these three great economists [7] contributed to build up the model in 1970 [3–5, 15]. As the Black-Scholes model is capable to provide a fair value of an option, a large number of researchers and investors have shown their interests in the study of the Black-Scholes model [5, 7].Basically, the BS model represents a PDE [4, 9, 13, 16, 17], and it is also one kind of parabolic equation [1, 7]. The Black-Scholes (BS) equation depends on two independent variables; the first is the time; and the second is the stock price [4, 5]. The BS equation can be derived from Brownian motion with the help of Ito’s lemma [3–5, 9, 18], what we discuss briefly in the next section, where we assume some facts and conditions while we consider the Black-Scholes pricing model [5, 9, 18]. Through a right choice of parameters, one may have a linear and nonlinear Black-Scholes equation where some basic assumptions are almost similar [9, 17]. Usually in most of the cases, simply the BS equation indicates the linear BS model. It is to be noted that when one wants to find the solution to a partial differential equation (PDE), it might be ended up with an infinite number of solutions [9]. That is why to get a unique solution, we impose initial/terminal and boundary conditions in the BS model [4, 9, 13].There are certain parameters that control the nonlinear Black-Scholes model such as transaction costs, volatility, market liquidity, uncertain risks, and some other issues [4, 7, 9, 16]. The nonlinearity is considered as one of the results of transaction costs [4, 9]. The nonlinear BS model is reduced to a nonlinear parabolic type equation [4, 9, 19]. Though the solution to the nonlinear Black-Scholes model sometimes shows asymptotic behavior [9, 20], the solution to the nonlinear BS model can converge even for large market frictions [4, 9, 20].The nonlinear model is considered as a model with transaction cost, whereas the linear model does not consider transaction cost [9, 21]. In the linear model, one restricts himself/herself by considering some assumptions, which can sometimes become unreal [16, 21]. On the other hand, the nonlinear BSmodel is more relaistic because of the parameter consideration [16, 21]. For example, in the nonlinear model, we consider volatility is nonconstant due to transaction cost [4, 9, 16, 21, 22]. That is why considering all the circumstances, we are having the priority to consider the nonlinear BS model over the linearized version.Approximating the nonlinear BS model and analyzing the parameter dynamics of the model are both very challenging. A few articles can be found where researchers approximate the nonlinear Black-Scholes model with a volatility corrector as a variable, and the corrector contains a function whose solution can only be found by solving a relevant initial value problem. Most applied researchers useOΔx2 schemes for the spatial approximation of the model, and thus, the solutions and the relevant sensitivity analysis lack sufficient accuracy. Keeping this fact in mind, this article focuses on to propose high-accurate easily applicable schemes for the model so that applied researchers can use it to analyze the model well. To be specific, in this study, the Chebyshev spectral method and the differential quadrature method (DQM) have been employed to approximate the BS (nonlinear/linear) model as these two schemes are very easy and powerful schemes in financial engineering. Different types of sensitivities (Greeks) of European option are also studied and highlighted in this study.Finance represents an emergent applied field where parameters play an important role. The motivation in studying sensitivity analysis of parameters in finance lies in the need for modeling the vagueness and uncertainty that are mostly implicit in many practical situations. Finance and related fields have developed over the years by adopting mathematical and computational tools and techniques.Estimation of sensitivity analysis of a model with respect to the changes in the parameter values is a crucial part of signifying the validity of the conclusions. Thus, the sensitivities of parameters of the BS model are represented by different Greek alphabets. These parameters play a decisive role in the definition of the shape of the option price.The rest of the article has been organized in the following way:(i)
Section2 contains a short discussion on the Black-Scholes model.(ii)
We propose exponential accurate numerical schemes (the Chebyshev spectral method and the differential quadrature method) for the spatial approximation of the model in Section3.(iii)
A brief discussion about the convergence of the proposed scheme has been included in Section4.(iv)
Then, we move on to present a numerical discussion on the BS model in Section5.(v)
Section6 contains sensitivity analysis for European option.(vi)
We summarize this study in Section7.


## 2. Preliminaries of the Black–Scholes Model

In this study, we aim to develop a fast numerical scheme to approximate both the linear and the nonlinear Black-Scholes models. Before we move on to the main topic, we need to discuss some basic issues to facilitate the reading of the article. LetS be underlying asset price or stock price, V be option price, and T be expiration time (or the time of maturity). It is assumed that the asset price abides by a geometric Brownian motion, which is(1)dS=μSdt+σSdW,where W represents the stochastic variable. W and dW work as uncertainty in stock market where Wt is random. μdt is the expected value, and σ2dt is the variance. Using Ito’s lemma, we have(2)dV=∂V∂t+12σ2S2∂2V∂S2+μS∂V∂Sdt+σS∂V∂SdW.Now using some changes in variables and some simple calculus, one can derive an uncertainty ΔW and W free risk-free model as(3)∂V∂t+12σ2S2∂2V∂S2+rS∂V∂S−rV=0,which is popularly known as the Black-Scholes partial differential equation. From equation (3), one may notice that the derivatives ∂V/∂t, ∂2V/∂S2, and ∂V/∂S must exist.Let call option is big and put option is small. Then, we have(4)VfS,t=VecS,t−VepS,t.Again, from the terminal condition of the Black-Scholes model, we have(5)VfS,t=S−Ke−rT−t.From the last two equations, it is evident that(6)VecS,t−VepS,t=S−Ke−rT−t,which is known as the put-call parity, and the formula is as follows:(7)Call−Put=Asset−Forward.The figure below shows the payoff function for both European call option and European put option for the inputs K=40, 0≤S≤80.Three situations are visible in Figure1, which are out-the-money, at-the-money, and in-the-money. At “out-the-money” and “at-the-money” situations, the holder does not exercise the option because cash flow ST−K is negative and zero, respectively [4]. On the other hand, at “in-the-money” position, the holder exercises the option as the cash flow ST−K is positive [4].Figure 1
Payoff function of European option.
(a)(b)In the nonlinear BS model, we assume that driftμ is constant and volatility σ is nonconstant [4, 16, 23, 24], which is defined as σ˜2=σ˜2t,S,VS,VSS [16, 25, 26]. Like the linear BS model, there are also some assumptions on the nonlinear BS model. Because of transaction costs, the drift μ and volatility σ depend on some other parameters, such as stock price S, time t, derivatives of the option price V, etc. It means that the nonlinearity arises as an output of transaction costs. The nonlinear Black-Scholes equation is as follows [16, 23, 25–27]:(8)Vt+rSVS+12σ˜2t,S,VS,VSSS2VSS−rV=0;0≤S≤S∗,0≤t≤T.The domain for the problem (8) is S∈0,S∗, t∈0,T, where S∗ is sufficiently large.The terminal and boundary conditions of the nonlinear Black-Scholes model for European call option are defined in the following way. The terminal condition is as follows:(9)VS,T=S−K+=maxS−K,0;0≤S≤S∗.The boundary condition is as follows:(10)V0,t=0;0≤t≤T,VS∗,t=S∗−Ke−rT−t;0≤t≤T.As the analytical solution to the nonlinear BS model (transaction cost models) is unknown [16], we focus on to develop numerical schemes of different volatility models. Some transaction cost models with different volatilities have been studied. The following three types of volatility models [4, 16, 27] have been considered, where κ represents round-trip transaction cost for per unit dollar of transaction, σ represents original volatility, and δt is transaction frequency.


### 2.1. Leland’s Model

In Leland’s model, the volatility of the nonlinear Black-Scholes equation is defined as follows [4, 9, 16, 25, 27]:(11)σ˜2=σ21+LesignVSS,where the Leland number is Le=2/πκ/σδt and signx=1,ifx>0,0,ifx=0,−1,ifx<0.


### 2.2. Barles’ and Soner’s Model

The volatility of Barles’ and Soner’s model is defined as follows [4, 9, 16, 23, 25, 27]:(12)σ˜2=σ21+ΨerT−ta2S2VSS,where a=κ/ε and Ψx satisfy the following equation [9, 23, 25]:(13)ddxΨx=Ψx+1−x+2xΨx,x≠0,satisfying Ψ0=0,limx⟶∞Ψx/x=1,limx⟶−∞Ψx=−1.For identity,Ψx=x. Then, the volatility of equation (12) becomes the following [9, 16, 23, 25]:(14)σ˜2=σ21+erT−ta2S2VSS.


### 2.3. Risk-Adjusted Pricing Methodology Model (RAPM Model)

The volatility of this model is as follows [4, 9, 16, 25–27]:(15)σ˜2=σ21+3C2M2πSVSS1/3,where M≥0 denotes the measure of transaction cost, and C≥0 denotes the measure of risk premium [4, 9, 16, 25].Now we move on to write the model (8) to facilitate the numerical integration of the model problem. To that end, we transfer the nonlinear BS equation of the European call option into a convection-diffusion equation through the following transformation of variables [16, 26].(16)x=lnSK,τ=12σ2T−t,VS,t=Kexux,τ.With the changes in variables, it yields(17)∂V∂t=−12σ2Suτ,∂V∂S=ux+u,∂2V∂S2=1Suxx+ux,and thus, the nonlinear BS model (8) transforms into(18)−12σ2Suτ+rSux+u+12σ˜2S21Suxx+ux−ruS=0,which can be written by the following simplified form(19)uτ=σ˜2σ2uxx+ux+Bux.Let B=2r/σ2 in equation (19) can be integrated to find the solution ux,τ for the domain x∈ℝ (i.e., −∞<x<∞) and 0≤τ≤T˜=σ2T/2; then, one may return to the original solution by backward substitution. This process confirms a need of numerical quadrature as it involves the nonlinearity σ˜. After the changes in variables, the volatility for Leland’s model becomes the following [16]:(20)σ˜2=σ21+Lesignuxx+ux.The volatility of Barles’ and Soner’s model is as follows [16]:(21)σ˜2=σ21+Ψe2rτ/σ2a2Kexuxx+ux,and the volatility of the RAPM model is as follows [16, 26]:(22)σ˜2=σ21+3C2M2πuxx+ux1/3,The initial and boundary conditions for European call option are defined in the following way [16, 26]. The initial condition is ux,0=1−e−x+=max1−e−x,0;x∈ℝ. The boundary condition is ux,τ=0;x⟶−∞, and ux,τ=1−e−Bτ−x;x⟶∞. Here, the main difficulty in approximating the nonlinear BS model (19) is to find a suitable scheme for the nonlinearity σ˜. We propose two different schemes for the nonlinearity in the next section.Now to get rid of using the nonlinearityσ˜, there is a strong evidence of using σ˜ as a constant in equation (8), which is the well-established linear Black-Scholes (BS) partial differential (3) [4, 9, 13]. The solution to BS equation (3) is defined on the domain 0≤S≤S∗ and 0≤t≤T, where S∗ is sufficiently large [10, 13, 18]. That is, the terminal condition is given for the time t=T [10] and the boundary conditions are taken for S=0 and S=S∗ [13]. For European call option, we have the following conditions [4, 5, 7, 10, 13, 28]:(23)VS,T=maxS−K,0=S−K+,for0≤S≤S∗,V0,t=0,fort∈0,T,VS∗,t=S∗−Ke−rT−t,fort∈0,T.The analytic solution to the linear BS equation can be obtained with the help of Fourier transformation from heat equation (13). The solution that we get is unique [5]. The theoretical solution to the BS model for European call option is given by(24)VS,t=VcS,t=CS,t=SNd1−Ke−rT−tNd2,and for European put option is given by(25)VS,t=VpS,t=PS,t=Ke−rT−tN−d2−SN−d1,where Nd is the standard normal distribution of the function d [4, 7, 13] with(26)d1=lnS/K+T−tr+σ2/2σT−t,d2=lnS/K+T−tr−σ2/2σT−t=d1−σT−t.In this case, we apply similar transformations in equation (12) of the nonlinear BS model. Then, we get the following convection-diffusion equation for European call option.(27)uτ=uxx+1+Bux,where B=2r/σ2. The initial and boundary conditions for the linear BS model are defined follows [16, 26].The initial condition isux,0=1−e−x+=max1−e−x,0;x∈ℝ, and the boundary conditions are ux,τ=0;x⟶−∞, and ux,τ=1−e−Bτ−x;x⟶∞.


## 2.1. Leland’s Model

In Leland’s model, the volatility of the nonlinear Black-Scholes equation is defined as follows [4, 9, 16, 25, 27]:(11)σ˜2=σ21+LesignVSS,where the Leland number is Le=2/πκ/σδt and signx=1,ifx>0,0,ifx=0,−1,ifx<0.


## 2.2. Barles’ and Soner’s Model

The volatility of Barles’ and Soner’s model is defined as follows [4, 9, 16, 23, 25, 27]:(12)σ˜2=σ21+ΨerT−ta2S2VSS,where a=κ/ε and Ψx satisfy the following equation [9, 23, 25]:(13)ddxΨx=Ψx+1−x+2xΨx,x≠0,satisfying Ψ0=0,limx⟶∞Ψx/x=1,limx⟶−∞Ψx=−1.For identity,Ψx=x. Then, the volatility of equation (12) becomes the following [9, 16, 23, 25]:(14)σ˜2=σ21+erT−ta2S2VSS.


## 2.3. Risk-Adjusted Pricing Methodology Model (RAPM Model)

The volatility of this model is as follows [4, 9, 16, 25–27]:(15)σ˜2=σ21+3C2M2πSVSS1/3,where M≥0 denotes the measure of transaction cost, and C≥0 denotes the measure of risk premium [4, 9, 16, 25].Now we move on to write the model (8) to facilitate the numerical integration of the model problem. To that end, we transfer the nonlinear BS equation of the European call option into a convection-diffusion equation through the following transformation of variables [16, 26].(16)x=lnSK,τ=12σ2T−t,VS,t=Kexux,τ.With the changes in variables, it yields(17)∂V∂t=−12σ2Suτ,∂V∂S=ux+u,∂2V∂S2=1Suxx+ux,and thus, the nonlinear BS model (8) transforms into(18)−12σ2Suτ+rSux+u+12σ˜2S21Suxx+ux−ruS=0,which can be written by the following simplified form(19)uτ=σ˜2σ2uxx+ux+Bux.Let B=2r/σ2 in equation (19) can be integrated to find the solution ux,τ for the domain x∈ℝ (i.e., −∞<x<∞) and 0≤τ≤T˜=σ2T/2; then, one may return to the original solution by backward substitution. This process confirms a need of numerical quadrature as it involves the nonlinearity σ˜. After the changes in variables, the volatility for Leland’s model becomes the following [16]:(20)σ˜2=σ21+Lesignuxx+ux.The volatility of Barles’ and Soner’s model is as follows [16]:(21)σ˜2=σ21+Ψe2rτ/σ2a2Kexuxx+ux,and the volatility of the RAPM model is as follows [16, 26]:(22)σ˜2=σ21+3C2M2πuxx+ux1/3,The initial and boundary conditions for European call option are defined in the following way [16, 26]. The initial condition is ux,0=1−e−x+=max1−e−x,0;x∈ℝ. The boundary condition is ux,τ=0;x⟶−∞, and ux,τ=1−e−Bτ−x;x⟶∞. Here, the main difficulty in approximating the nonlinear BS model (19) is to find a suitable scheme for the nonlinearity σ˜. We propose two different schemes for the nonlinearity in the next section.Now to get rid of using the nonlinearityσ˜, there is a strong evidence of using σ˜ as a constant in equation (8), which is the well-established linear Black-Scholes (BS) partial differential (3) [4, 9, 13]. The solution to BS equation (3) is defined on the domain 0≤S≤S∗ and 0≤t≤T, where S∗ is sufficiently large [10, 13, 18]. That is, the terminal condition is given for the time t=T [10] and the boundary conditions are taken for S=0 and S=S∗ [13]. For European call option, we have the following conditions [4, 5, 7, 10, 13, 28]:(23)VS,T=maxS−K,0=S−K+,for0≤S≤S∗,V0,t=0,fort∈0,T,VS∗,t=S∗−Ke−rT−t,fort∈0,T.The analytic solution to the linear BS equation can be obtained with the help of Fourier transformation from heat equation (13). The solution that we get is unique [5]. The theoretical solution to the BS model for European call option is given by(24)VS,t=VcS,t=CS,t=SNd1−Ke−rT−tNd2,and for European put option is given by(25)VS,t=VpS,t=PS,t=Ke−rT−tN−d2−SN−d1,where Nd is the standard normal distribution of the function d [4, 7, 13] with(26)d1=lnS/K+T−tr+σ2/2σT−t,d2=lnS/K+T−tr−σ2/2σT−t=d1−σT−t.In this case, we apply similar transformations in equation (12) of the nonlinear BS model. Then, we get the following convection-diffusion equation for European call option.(27)uτ=uxx+1+Bux,where B=2r/σ2. The initial and boundary conditions for the linear BS model are defined follows [16, 26].The initial condition isux,0=1−e−x+=max1−e−x,0;x∈ℝ, and the boundary conditions are ux,τ=0;x⟶−∞, and ux,τ=1−e−Bτ−x;x⟶∞.


## 3. Spatial Integration of the Model

An efficient and high-accurate scheme is important for any differential equation. For the BS model, spatial approximation means approximation of spatial derivatives. The finite difference method and the linear finite element method are popular schemes to serve the purpose withOΔx2 accuracy. Thus, one needs sufficient spatial grid points to confirm the desired accuracy of the solutions [22, 29–32]. An efficient alternative to that is to use some function approximations, especially global polynomial approximation can significantly reduce the consideration of grid points [22, 29, 30, 33]. Keeping this basic in mind, in this section, we focus on to approximate the model using global polynomial-based schemes. The main goal here is to use a high-accurate scheme that is easily applicable. By that it is meant a situation where a differential operator of any order can be replaced by a matrix operator so that an applied researcher/finance people can easily apply the scheme and get a better market prediction. There may have several such schemes for equations (19) and (27) subject to some appropriate initial and boundary conditions [30]. Here, Chebyshev spectral scheme and differential quadrature scheme have been employed to serve the purpose. This work focuses on to develop two exponential accurate schemes for the spatial approximation. To that end, we begin with replacing x∈ℝ and the temporal domain by x∈−1,1 and 0,T˜, respectively [10, 16]. In the Chebyshev method and the differential quadrature method, the spatial domain xi is defined for −1 to 1. That is why we are taking the spatial domain x∈−1,1 to apply the spectral methods (Chebyshev and DQM) more effectively.As in the convection-diffusion equation, we have first- and second-order derivatives in space, and we replace these derivatives by the spectral method. That is,∂d/∂xdu is replaced by Ddu, where d is the order of the derivative, and D is a spectral matrix.


### 3.1. Chebyshev Spectral Method

The spectral method is an efficient tool to solve ODE or PDE [22, 32, 33]. If we apply the spectral method in the financial field, then the problem becomes more smooth. A very simple example is presented in this paper to highlight the role of the spectral method in accuracy and computational scale. This method is very simple, more accurate, and computationally inexpensive [22, 29, 32–34].There are several methods resorted to solve the Black-Scholes equation. These methods are successful for their framework. The spectral method is one of them [34]. The spectral method has gained mathematical maturity in financial areas [35]. We get impressive results from the spectral collocation method [22, 29, 33]. We have proposed the Chebyshev spectral method. Chebyshev points are essential in the spectral method [32, 34].We take the following Chebyshev collocation or grid points [22, 29, 33, 34].(28)xi=cosiπN;i=0,1,2,…,N.These points are used to form the Chebyshev differentiation matrix. We get Chebyshev collocation points for the interval −1,1 [33]. We take the Lagrangian form for polynomial px [22, 33]. The basic idea of spectral approximation is to use spectral derivative in the place of differential operator [22, 29, 33]. We denote N+1×N+1-sized Chebyshev matrix by DN [29, 33]. Let v=px be a polynomial and N=2. Then, Chebyshev points are as follows: x0=1, x1=0, and x2=−1. As we get 3 points, so the interpolating polynomial px will be quadratic. We define px as [22, 33](29)px=12x1+xu0+1+x1−xu1+12xx−1u2.Then,(30)p′x=x+12u0−2xu1+x−12u2.1st k=1 term in equation (30) is p1=x+1/2,(31)x0=1⟹p1=1+12=32,x1=0⟹p1=0+12=12,x2=−1⟹p1=−1+12=−12.So, the 1st k=1 column of DN matrix is DN=p1x0⋯⋯p1x1⋯⋯p1x2⋯⋯=3/2⋯⋯1/2⋯⋯−1/2⋯⋯.Similarly, 2nd k=2 term in equation (30) is p2=−2x,(32)x0=1⟹p2=−2×1=−2,x1=0⟹p2=−2×0=0,x2=−1⟹p2=−2−1=2.So, the 2nd k=2 column of DN matrix is DN=⋯p2x0⋯⋯p2x1⋯⋯p2x2⋯=⋯−2⋯⋯0⋯⋯2⋯.Finally, 3rd k=3 term in equation (30) is p3=x−1/2,(33)x0=1⟹p3=1−12=12,x1=0⟹p3=0−12=−12,x2=−1⟹p3=−1−12=−32.Thus, the 3rd k=3 column of DN matrix is DN=⋯⋯p3x1⋯⋯p3x2⋯⋯p3x3=⋯⋯1/2⋯⋯−1/2⋯⋯−3/2.So, DN matrix is [22, 33](34)DN=32−212120−12−122−32.In general, the k th column of DN matrix is computed by putting the Chebyshev points at the k th term of the polynomial p′x [22, 33]. There are formulas for the entries of DN matrix [22, 29, 33, 34]. The explicit methods are helpful for the Chebyshev differentiation matrix [22]. Let N≥1. Then, the rows and columns of the Chebyshev differentiation matrix DN can be obtained for index 0 to N [22, 29]. The entries of DN matrix for the first-order derivative can be written as follows [22, 29, 33, 34]:(35)DN00=2N2+16,DNNN=−2N2+16,DNwk=−xk21−xk2;k=1,2,…,N−1,DNwk=gi−1w+kgjxw−xk;w≠k&w,k=1,2,…,N−1.Here,(36)gw=2;w=0,N,1;otherwise.The spatial derivative ∂/∂x is approximated by the matrix operator D. For example, the second-order derivative ∂2/∂x2 is computed by D2, i.e., we square D, which a is very simple and easily applicable procedure.


### 3.2. Differential Quadrature Method

The differential quadrature method (DQM) is a simple and straightforward technique to implement [36, 37]. The DQM is known for its simplicity, accuracy, and efficiency [36–38]. We need a very little computational endeavor to apply DQM [36–38]. We can use DQM to solve the initial value problem and boundary value problem. This method is significant to solve critical problems in option pricing. So this method can be a suitable choice for option pricing problems. The DQM is similar to the collocation or the pseudospectral method [36]. The DQM was developed by Richard Bellman and his associates Kashef and Casti in 1971 [36, 38]. The DQM is used to discretize spatial derivatives. For our problem, ∂/∂x and ∂2/∂x2 are replaced by differential quadrature matrix D and D2, respectively.Grid points taken in DQM have a swift convergence rate rather than equally spaced mesh points [36, 37]. If we choose nonuniform grid points, we can get stable results [37, 39]. That is why we choose Chebyshev-Gauss-Lobatto (CGL) grid points in spatial discretization to have a stable result [37]. The Chebyshev-Gauss-Lobatto distribution is defined as follows [38, 40, 41]:(37)xi=121−cosi−1πN−1;i=1,2,…,N.We solve the BS equation for the intervalx∈−1,1. For Chebyshev-Gauss-Lobatto points, we get the grid points in the interval 0,1. To transfer the interval from 0,1 to −1,1, we use the linear transformation 2x−1.In the DQM, we need weighting coefficients [36–38, 40, 41]. We use an approximating function to find the weighting coefficients [41]. This approximate function is known as test or trial functions. We approximate the spatial derivatives at the grid point xi in the following manner [36–38, 40]. To that end, it is considered that(38)pkx=Lnxx−xkLn1xk,wherek=11n,and(39)Lnx=∏j=1nx−xj,Ln1xi=∏j=1i≠jnxi−xj,Then, the unknown function can be approximated by(40)ux,τ=∑j=1npkxuxj,τ,where xj is the discrete sampling points in the spatial domain of computation. The derivatives of such a function at x=xi can be presented by(41)∂q∂xquxi,τ=∑j=1NAijquxj,τ;i=1,2,…,N,j=1,2,…,N,q=1,2,where Aijq, q=1,2, is the weighting coefficient at the grid point xi. Lagrange polynomials pkx are used to compute the weighting coefficient Aijq. Now we express the exact values of approximate derivatives and start by presenting entries of the first derivative using the base function pkx, k=11n. There are formulas to express the diagonal entries and off-diagonal entries for the weighting coefficients of the first-order derivative [36], which is given by(42)Aij1=∏k=1k≠iNxi−xkxj−xi∏k=1k≠jNxj−xk;i≠j,∑k=1k≠iN1xi−xk;i=j,In a similar fashion, the formula for weighting coefficients of the second-order derivative is presented below(43)Aij2=2Aij1Aii1−1xi−xj;i≠j,∑k=1k≠iNAik2;i=j,One may use lower-order weighting coefficients defined above to express higher-order weighting coefficients. We can compute the weighting coefficients for higher-order derivatives by using the following recurrence formula(44)Aijq=Aij1Aiiq−1−Aijq−1xi−xj;i≠j,−∑k=1k≠iNAikq;i=j,For further details, one may consult [36, 41]. For this study, it is assumed that ∂/∂x≡D≡A and ∂2/∂x2≡D2≡A2.After spatial integration, equations (19) and (27) can be written as the following semidiscrete system of equations of the form(45)uτ=fA,u,where A is a matrix, and u is a solution vector that depends on τ. The equation given by equation (45) is decomposed by δi‘s and its time derivatives δ˙i‘s by the following Crank-Nicolson formula(46)δi=δin+1+δin2,and usual finite difference approximation(47)δ˙i=δin+1−δinΔτ.This approximation leads to the following recurrence relationship between two time levels n and n+1 relating to δin+1 and δin.That is, we replace the derivatives∂/∂x and ∂2/∂x2 in equation (19) and in volatility equations (20)–(22) by D and D2 matrices, respectively, which are taken from the Chebyshev spectral method and the differential quadrature method. Equation (19) is iteratively executed for each temporal integration.


### 3.3. Computation Algorithms

We investigate the convection-diffusion equation (19) for the transformed mesh/domain −1,1 and 0≤τ≤T˜=σ2T/2 [16, 26]. We consider the grid points as xi,τn, where spatial points xi are taken from the Chebyshev spectral method and the differential quadrature method. For the linear Black-Scholes model, we apply the same procedure that we did with the nonlinear BS model. That is, equation (27) is solved by applying the Chebyshev spectral scheme and the differential quadrature method in space with the explicit FDM in time. The finite difference method (FDM) can be successfully implemented in time integration as this method is very easy to implement and provide the most accurate results [1, 42, 43]. The explicit FDM method has truncation error of OΔτ [1, 21, 42, 43]. We replace the derivatives ∂/∂x and ∂2/∂x2 in equation (27) by D and D2 matrices, respectively, which are taken from the Chebyshev spectral method and the differential quadrature method. For each temporal point, equation (27) is iteratively executed. We present a simple computational algorithm/process to facilitate the readership, and complete MATLAB codes will be available upon request (MATLAB codes for the nonlinear model will be made public (will put on the article/web page) once the article is accepted if the honorable reviewers agree).Algorithm 1.
The nonlinear BS model.Step 1. We computeσ˜2 (equations (20)–(22) for Leland’s model, Barles’-Soner’s model, and RAPM model, respectively, in equation (19)).Step 2. We replaceuxx and ux by D2u and Du, respectively, in(48)uτ=L,whereL=σ˜2/σ2uxx+ux+Bux.Step 3. We replace the time derivativeuτ in equation (48) by the FDM formula uτ=uin+1−uin/Δτ or the Crank-Nicolson method to get(49)un+1=un+ΔτL,Step 4. We iteratively execute equation (49) for each time grid to approximate the solution.Step 5. After finding the solution to the convection-diffusion equation (19), we return back to the original solution to the nonlinear model with the help of the backward substitution equation (16).Algorithm 2.
The linear BS model.Step 1. We replace the derivativesuxx and ux by D2u and Du, respectively, in(50)uτ=L,whereL=uxx+1+Bux.Step 2. We replace the time derivativeuτ in equation (50) by the formula uτ=uin+1−uin/Δτ or the Crank-Nicolson scheme to achieve(51)un+1=un+ΔτL.Step 3. We iterate equation (51) to get a solution for each time stepping.Step 4. We use back substitutions of variables to return back to the original solution to the linear Black-Scholes model by the backward substitution of equation (16).


## 3.1. Chebyshev Spectral Method

The spectral method is an efficient tool to solve ODE or PDE [22, 32, 33]. If we apply the spectral method in the financial field, then the problem becomes more smooth. A very simple example is presented in this paper to highlight the role of the spectral method in accuracy and computational scale. This method is very simple, more accurate, and computationally inexpensive [22, 29, 32–34].There are several methods resorted to solve the Black-Scholes equation. These methods are successful for their framework. The spectral method is one of them [34]. The spectral method has gained mathematical maturity in financial areas [35]. We get impressive results from the spectral collocation method [22, 29, 33]. We have proposed the Chebyshev spectral method. Chebyshev points are essential in the spectral method [32, 34].We take the following Chebyshev collocation or grid points [22, 29, 33, 34].(28)xi=cosiπN;i=0,1,2,…,N.These points are used to form the Chebyshev differentiation matrix. We get Chebyshev collocation points for the interval −1,1 [33]. We take the Lagrangian form for polynomial px [22, 33]. The basic idea of spectral approximation is to use spectral derivative in the place of differential operator [22, 29, 33]. We denote N+1×N+1-sized Chebyshev matrix by DN [29, 33]. Let v=px be a polynomial and N=2. Then, Chebyshev points are as follows: x0=1, x1=0, and x2=−1. As we get 3 points, so the interpolating polynomial px will be quadratic. We define px as [22, 33](29)px=12x1+xu0+1+x1−xu1+12xx−1u2.Then,(30)p′x=x+12u0−2xu1+x−12u2.1st k=1 term in equation (30) is p1=x+1/2,(31)x0=1⟹p1=1+12=32,x1=0⟹p1=0+12=12,x2=−1⟹p1=−1+12=−12.So, the 1st k=1 column of DN matrix is DN=p1x0⋯⋯p1x1⋯⋯p1x2⋯⋯=3/2⋯⋯1/2⋯⋯−1/2⋯⋯.Similarly, 2nd k=2 term in equation (30) is p2=−2x,(32)x0=1⟹p2=−2×1=−2,x1=0⟹p2=−2×0=0,x2=−1⟹p2=−2−1=2.So, the 2nd k=2 column of DN matrix is DN=⋯p2x0⋯⋯p2x1⋯⋯p2x2⋯=⋯−2⋯⋯0⋯⋯2⋯.Finally, 3rd k=3 term in equation (30) is p3=x−1/2,(33)x0=1⟹p3=1−12=12,x1=0⟹p3=0−12=−12,x2=−1⟹p3=−1−12=−32.Thus, the 3rd k=3 column of DN matrix is DN=⋯⋯p3x1⋯⋯p3x2⋯⋯p3x3=⋯⋯1/2⋯⋯−1/2⋯⋯−3/2.So, DN matrix is [22, 33](34)DN=32−212120−12−122−32.In general, the k th column of DN matrix is computed by putting the Chebyshev points at the k th term of the polynomial p′x [22, 33]. There are formulas for the entries of DN matrix [22, 29, 33, 34]. The explicit methods are helpful for the Chebyshev differentiation matrix [22]. Let N≥1. Then, the rows and columns of the Chebyshev differentiation matrix DN can be obtained for index 0 to N [22, 29]. The entries of DN matrix for the first-order derivative can be written as follows [22, 29, 33, 34]:(35)DN00=2N2+16,DNNN=−2N2+16,DNwk=−xk21−xk2;k=1,2,…,N−1,DNwk=gi−1w+kgjxw−xk;w≠k&w,k=1,2,…,N−1.Here,(36)gw=2;w=0,N,1;otherwise.The spatial derivative ∂/∂x is approximated by the matrix operator D. For example, the second-order derivative ∂2/∂x2 is computed by D2, i.e., we square D, which a is very simple and easily applicable procedure.


## 3.2. Differential Quadrature Method

The differential quadrature method (DQM) is a simple and straightforward technique to implement [36, 37]. The DQM is known for its simplicity, accuracy, and efficiency [36–38]. We need a very little computational endeavor to apply DQM [36–38]. We can use DQM to solve the initial value problem and boundary value problem. This method is significant to solve critical problems in option pricing. So this method can be a suitable choice for option pricing problems. The DQM is similar to the collocation or the pseudospectral method [36]. The DQM was developed by Richard Bellman and his associates Kashef and Casti in 1971 [36, 38]. The DQM is used to discretize spatial derivatives. For our problem, ∂/∂x and ∂2/∂x2 are replaced by differential quadrature matrix D and D2, respectively.Grid points taken in DQM have a swift convergence rate rather than equally spaced mesh points [36, 37]. If we choose nonuniform grid points, we can get stable results [37, 39]. That is why we choose Chebyshev-Gauss-Lobatto (CGL) grid points in spatial discretization to have a stable result [37]. The Chebyshev-Gauss-Lobatto distribution is defined as follows [38, 40, 41]:(37)xi=121−cosi−1πN−1;i=1,2,…,N.We solve the BS equation for the intervalx∈−1,1. For Chebyshev-Gauss-Lobatto points, we get the grid points in the interval 0,1. To transfer the interval from 0,1 to −1,1, we use the linear transformation 2x−1.In the DQM, we need weighting coefficients [36–38, 40, 41]. We use an approximating function to find the weighting coefficients [41]. This approximate function is known as test or trial functions. We approximate the spatial derivatives at the grid point xi in the following manner [36–38, 40]. To that end, it is considered that(38)pkx=Lnxx−xkLn1xk,wherek=11n,and(39)Lnx=∏j=1nx−xj,Ln1xi=∏j=1i≠jnxi−xj,Then, the unknown function can be approximated by(40)ux,τ=∑j=1npkxuxj,τ,where xj is the discrete sampling points in the spatial domain of computation. The derivatives of such a function at x=xi can be presented by(41)∂q∂xquxi,τ=∑j=1NAijquxj,τ;i=1,2,…,N,j=1,2,…,N,q=1,2,where Aijq, q=1,2, is the weighting coefficient at the grid point xi. Lagrange polynomials pkx are used to compute the weighting coefficient Aijq. Now we express the exact values of approximate derivatives and start by presenting entries of the first derivative using the base function pkx, k=11n. There are formulas to express the diagonal entries and off-diagonal entries for the weighting coefficients of the first-order derivative [36], which is given by(42)Aij1=∏k=1k≠iNxi−xkxj−xi∏k=1k≠jNxj−xk;i≠j,∑k=1k≠iN1xi−xk;i=j,In a similar fashion, the formula for weighting coefficients of the second-order derivative is presented below(43)Aij2=2Aij1Aii1−1xi−xj;i≠j,∑k=1k≠iNAik2;i=j,One may use lower-order weighting coefficients defined above to express higher-order weighting coefficients. We can compute the weighting coefficients for higher-order derivatives by using the following recurrence formula(44)Aijq=Aij1Aiiq−1−Aijq−1xi−xj;i≠j,−∑k=1k≠iNAikq;i=j,For further details, one may consult [36, 41]. For this study, it is assumed that ∂/∂x≡D≡A and ∂2/∂x2≡D2≡A2.After spatial integration, equations (19) and (27) can be written as the following semidiscrete system of equations of the form(45)uτ=fA,u,where A is a matrix, and u is a solution vector that depends on τ. The equation given by equation (45) is decomposed by δi‘s and its time derivatives δ˙i‘s by the following Crank-Nicolson formula(46)δi=δin+1+δin2,and usual finite difference approximation(47)δ˙i=δin+1−δinΔτ.This approximation leads to the following recurrence relationship between two time levels n and n+1 relating to δin+1 and δin.That is, we replace the derivatives∂/∂x and ∂2/∂x2 in equation (19) and in volatility equations (20)–(22) by D and D2 matrices, respectively, which are taken from the Chebyshev spectral method and the differential quadrature method. Equation (19) is iteratively executed for each temporal integration.


## 3.3. Computation Algorithms

We investigate the convection-diffusion equation (19) for the transformed mesh/domain −1,1 and 0≤τ≤T˜=σ2T/2 [16, 26]. We consider the grid points as xi,τn, where spatial points xi are taken from the Chebyshev spectral method and the differential quadrature method. For the linear Black-Scholes model, we apply the same procedure that we did with the nonlinear BS model. That is, equation (27) is solved by applying the Chebyshev spectral scheme and the differential quadrature method in space with the explicit FDM in time. The finite difference method (FDM) can be successfully implemented in time integration as this method is very easy to implement and provide the most accurate results [1, 42, 43]. The explicit FDM method has truncation error of OΔτ [1, 21, 42, 43]. We replace the derivatives ∂/∂x and ∂2/∂x2 in equation (27) by D and D2 matrices, respectively, which are taken from the Chebyshev spectral method and the differential quadrature method. For each temporal point, equation (27) is iteratively executed. We present a simple computational algorithm/process to facilitate the readership, and complete MATLAB codes will be available upon request (MATLAB codes for the nonlinear model will be made public (will put on the article/web page) once the article is accepted if the honorable reviewers agree).Algorithm 1.
The nonlinear BS model.Step 1. We computeσ˜2 (equations (20)–(22) for Leland’s model, Barles’-Soner’s model, and RAPM model, respectively, in equation (19)).Step 2. We replaceuxx and ux by D2u and Du, respectively, in(48)uτ=L,whereL=σ˜2/σ2uxx+ux+Bux.Step 3. We replace the time derivativeuτ in equation (48) by the FDM formula uτ=uin+1−uin/Δτ or the Crank-Nicolson method to get(49)un+1=un+ΔτL,Step 4. We iteratively execute equation (49) for each time grid to approximate the solution.Step 5. After finding the solution to the convection-diffusion equation (19), we return back to the original solution to the nonlinear model with the help of the backward substitution equation (16).Algorithm 2.
The linear BS model.Step 1. We replace the derivativesuxx and ux by D2u and Du, respectively, in(50)uτ=L,whereL=uxx+1+Bux.Step 2. We replace the time derivativeuτ in equation (50) by the formula uτ=uin+1−uin/Δτ or the Crank-Nicolson scheme to achieve(51)un+1=un+ΔτL.Step 3. We iterate equation (51) to get a solution for each time stepping.Step 4. We use back substitutions of variables to return back to the original solution to the linear Black-Scholes model by the backward substitution of equation (16).


## 4. Convergence of the Spatiotemporal Scheme

In this part of the study, we aim to analyze the error while integrating a generalized BS equation, which is a nonlinear PDE, using the spectral method in space followed by a one-step scheme for the temporal approximation. It is to note that the performance of a numerical scheme depends on its error reduction and computational stability. The global polynomial schemes are very much popular for their priori and a posteriori estimates of the rate of convergence. Most theories depend on a functional analysis framework, which has been discussed in detail in numerous research write-ups. Here, we review some relevant concepts and key outcomes without a formal proof. For clarity, we cite relevant sources. To be specific, we opt to a brief survey of the error estimates of the spatiotemporal scheme we employed here. One may consult [22, 32, 44–46] and many other well-known references for detail. It is to mention that we employ some constants Ci≥0 in this study, which are necessarily not the same for all the cases. The error of such an approximation is of exponential order if the approximating function u is smooth enough. However, if the function is analytic piecewise but discontinuous then the error is of O1. Near the point where it breaks the continuity, the error occurs of O1/N, which is known as Gibbs phenomenon.Usually, global polynomial interpolations are employed to find polynomial approximations of the solutions of ordinary and partial differential equations considering a simple computational domain, and unknown curves are considered to be smooth enough. Polynomials are smooth, which is very important in the approximation theory. It assists researchers to approximate and analyze solutions using the polynomial basis functions [22]. Let u be the approximation and INu be the interpolating polynomial. For the spectral method, the difference between u and INu exponentially goes to 0. Let C be a function of u. Then, upper bounds of u−INu can be written as(52)u−INuLw2≤CNm∑i=0muiLw2,u−INu∞≤CNm−1/2∑i=0muiLw2.If we increase the power ofm, then these upper bounds decrease. That is, if we increase the power of N, then the error becomes less, i.e., the difference u−INu⟶0 if N increases. Then, we can say the solution u converges. In short, if N⟶∞, then there is kind of guaranteed convergence of the spectral method.So this approach is simpler and efficient to compute solutions. Here, it is well known that the typical convergence rate isON−m for every m for functions that are smooth and OcN, 0<c<1, for every analytic function. Thus, if one approximates a smooth function wx by spectral polynomials whx, then it yields(53)wx−whx=ON−m.For time, a forward difference scheme has been used, which is accurate ofOΔτ in L20,T norm for some T>0 [22, 46]. So for the space-time discretization, the error bound is of the form(54)ux,τ−uhx,τ≤C1N−m+C2Δτ,for a suitable C1≥0 and C2≥0.


## 5. Numerical Discussion

In this part of the study, we focus on to discuss the computer-implemented results of the model problem. Both the schemes have been implemented for the nonlinear BSM and the linear BSM. For convenience, we build up a MATLAB routine for the Chebyshev spectral matrix and the differential quadrature matrix to run the calculation. For the temporal integration of the semidiscrete time-dependent model, we use the finite difference scheme (explicit FDM).


### 5.1. Nonlinear Black–Scholes Model

The graphs of the volatility models of European call option are presented here. At timet=0 [10], we plot the volatility models and some numerical values are also presented in a tabulated form. We consider the following parameter values σ=0.2, r=0.1, K=75, T=1, N=Sizeofmatrix=40, number of temporal points, m=2500, a=0.02, Le=0.6, M=2, and C=0.01, and start by illustrating the solutions of the nonlinear BS model by using the DQM. We present solutions of the different volatility models in Figure 2.Figure 2
Different types of volatility models by differential quadrature method. Here, we useK=75, T=1, R=1, N=40, m=2500, σ=0.2, r=0.1, a=0.02, Le=0.6, M=2, and C=0.01.Figure 3
Here, we illustrateVNL−VL considering the volatility models by the differential quadrature method, where NL stands for the nonlinear model and L stands for the linear version of the same.In Figure3 we present the difference between European call option value with the nonlinear model and European call option value with the linear model. This computation confirms the economical deviation of prices between the models. It is observed that the difference is minimum at the expiry values and maximum at about S=80 when t=0.It is noticed from Figure4 that the highest price is presented by Leland’s model, after that Barles’-Soner’s model, next RAPM model, followed by Barles’-Soner’s (identity) model at the time to expiration 1 year (i.e., t=0). To compute the errors for volatility models, we compare the volatility models with the linear Black-Scholes scheme. That is, we take the absolute value of the difference between option value of volatility models and option value of the linear BS model to measure the deviation/error. Table 1 shows absolute errors of different volatility models executed by DQM at time t=0.Figure 4
2D plot of the volatility models by the differential quadrature method att=0.Table 1
Absolute errors of volatility models by differential Quadrature scheme att=0. Here, BS means the Black-Scholes model, and BSI means the Barles’-Soner’s (identity).
Stock price (S)Linear BS (exact values)RAPM valueErrorLeland’s valueErrorBarles’-Soner’s valueErrorBSI valueError30.939450.000090.000170.000080.003260.003170.000250.000160.000120.0000332.210680.000220.000390.000170.005850.005640.000530.000320.000290.0000733.718460.000560.000980.000420.010970.010410.001380.000810.000740.00017Then, we move on to illustrate solutions using the spectral method. Figure5 represents different volatility models executed by the Chebyshev method. In Figure 6, we compare solutions between the volatility models (transaction cost models) and the linear BS model. We plot the absolute difference between the transaction cost model and the linear Black-Scholes model in Figure 6. From Figure 6, we observe that the difference is significant when the time to expiration is 1 year (i.e., t=0) for all the volatility models [16, 27], where the stock price approximately lies between S=50 and S=130. That is, the nonlinear price (model with transaction cost) is significantly bigger than the linear price (model without transaction cost) at the very beginning of the time (i.e., t=0 or time to expiration 1 year) [16, 27]. But we note from this illustration that the difference decreases when the time to expiration decreases (i.e., t increases) [27].Figure 5
Different types of volatility models by Chebyshev spectral method. Here, we user=0.1, K=75, T=1, σ=0.2, R=1, N=40, m=2500, M=2, C=0.01, a=0.02, and Le=0.6.Figure 6
Here, we illustrateVNL−VL considering volatility models by Chebyshev spectral method, where NL stands for nonlinear model and L stands for linear model.Figure 7
2D plot of volatility models by Chebyshev spectral method att=0.From Figure7 we can see that for Chebyshev spectral method, we can arrange the volatility models according to highest to lowest prices at t=0 in the following serial: first, Leland’s price, next Barles’-Soner’s value, after that RAPM value, and finally Barles’-Soner’s (identity) price. In Table 2, we present the error VNL−VL of volatility models by the Chebyshev spectral scheme at t=0.Table 2
The absolute error of volatility models by Chebyshev spectral method att=0. Here, BSI means Barles’-Soner’s (identity).
Stock price (S)Linear BS model (exact values)RAPM valueErrorLeland’s valueErrorBarles’-Soner’s valueErrorBSI valueError37.515860.004360.007070.002710.042740.038370.009810.005440.005540.0011839.846830.012410.019270.006850.086490.074080.024110.01170.015490.0030742.496570.034780.051370.016590.174680.139900.061720.026940.042540.00776


### 5.2. Numerical Results of Linear Black–Scholes Model

Then, we start to apply the scheme for the linear version of the model. For the linear model,σ=0.2, r=0.1, N=Size ofDmatrix=60, m=Number of temporal points=6000, K=75, and T=1 year have been considered along with the European call option in Figure 8. The absolute difference between the exact solutions and the approximate solutions has been presented in Figure 9. It is we noticed from the illustration that the difference between exact option value and numerical option value is significant when S=73.0298 at time to expiration T−t=0.Figure 8
Linear Black-Scholes model by the Chebyshev spectral method (a) and by the differential quadrature method (b).
(a)(b)Figure 9
Here, we plotVExact−VNumeric of the linear Black-Scholes model. (a) One is for the Chebyshev spectral method, and (b) one is for the differential quadrature method.
(a)(b)It is evident from Figure10 that both the schemes provide very close results to the exact values. This makes these two schemes more reliable to approximate the Black-Scholes model. Though the entries in the Chebyshev spectral and differential quadrature methods are similar, their positions are just opposite to each other. Next, we present the approximate values of the linear BS model at t=0 in the following table where absolute error is considered. Here, by error we mean | exact value-approximate value |.Figure 10
2D plot of the linear Black-Scholes model att=0.We have also compared our proposed methods with a benchmark method FDM (finite difference method) to justify the accuracy of spectral methods (Chebyshev and DQM). The discrete solutions of equation (27) are denoted as uxi,τj=uij. In our methods, spatial grid points xi are nonuniform and temporal grid points τj are uniform. We replace the derivatives in convection-diffusion equation (27) to apply FDM [1, 47, 48].(55)uij+1−uijΔτ=2xi+1−xiui−1j−xi+1−xi−1uij+xi−xi−1ui+1jxi−xi−1xi+1−xixi+1−xi−1+1+Bui+1j−ui−1jxi+1−xi−1,i.e.uij+1=Δτ2xi+1−xiui−1j−xi+1−xi−1uij+xi−xi−1ui+1jxi−xi−1xi+1−xixi+1−xi−1+1+Bui+1j−ui−1jxi+1−xi−1+uij.Equation (55) is executed iteratively on time grids τj to approximate the solutions in case of FDM.We notice from Table3 that the approximations are similar for both Chebyshev spectral and differential quadrature schemes, and the errors for these two schemes are significantly small. These two methods approximate the Black-Scholes model way better than FDM. That is why it is well justified to say that spectral methods are reliable to approximate BS model.Table 3
The absolute error of Chebyshev spectral scheme, differential quadrature scheme, and FDM att=0.
Stock price (S)V_ExactV_ChebyshevErrorV_DQMErrorV_FDMError37.332210.003990.004110.000120.004110.000120.005780.0017938.821600.007980.008190.000210.008190.000210.010850.0028840.445830.015890.016270.000380.016270.000380.020410.0045342.211840.031360.032040.000670.032040.000670.038260.0068944.126800.060990.062150.001160.062150.001160.071050.0100646.198030.116210.118140.001930.118140.001930.130140.0139348.432970.215860.218930.003070.218930.003070.233920.0180750.839010.389190.393850.004660.393850.004660.410810.0216253.423420.678610.685330.006720.685330.006720.702020.0234156.193221.141091.150220.009131.150220.009131.163300.02221


## 5.1. Nonlinear Black–Scholes Model

The graphs of the volatility models of European call option are presented here. At timet=0 [10], we plot the volatility models and some numerical values are also presented in a tabulated form. We consider the following parameter values σ=0.2, r=0.1, K=75, T=1, N=Sizeofmatrix=40, number of temporal points, m=2500, a=0.02, Le=0.6, M=2, and C=0.01, and start by illustrating the solutions of the nonlinear BS model by using the DQM. We present solutions of the different volatility models in Figure 2.Figure 2
Different types of volatility models by differential quadrature method. Here, we useK=75, T=1, R=1, N=40, m=2500, σ=0.2, r=0.1, a=0.02, Le=0.6, M=2, and C=0.01.Figure 3
Here, we illustrateVNL−VL considering the volatility models by the differential quadrature method, where NL stands for the nonlinear model and L stands for the linear version of the same.In Figure3 we present the difference between European call option value with the nonlinear model and European call option value with the linear model. This computation confirms the economical deviation of prices between the models. It is observed that the difference is minimum at the expiry values and maximum at about S=80 when t=0.It is noticed from Figure4 that the highest price is presented by Leland’s model, after that Barles’-Soner’s model, next RAPM model, followed by Barles’-Soner’s (identity) model at the time to expiration 1 year (i.e., t=0). To compute the errors for volatility models, we compare the volatility models with the linear Black-Scholes scheme. That is, we take the absolute value of the difference between option value of volatility models and option value of the linear BS model to measure the deviation/error. Table 1 shows absolute errors of different volatility models executed by DQM at time t=0.Figure 4
2D plot of the volatility models by the differential quadrature method att=0.Table 1
Absolute errors of volatility models by differential Quadrature scheme att=0. Here, BS means the Black-Scholes model, and BSI means the Barles’-Soner’s (identity).
Stock price (S)Linear BS (exact values)RAPM valueErrorLeland’s valueErrorBarles’-Soner’s valueErrorBSI valueError30.939450.000090.000170.000080.003260.003170.000250.000160.000120.0000332.210680.000220.000390.000170.005850.005640.000530.000320.000290.0000733.718460.000560.000980.000420.010970.010410.001380.000810.000740.00017Then, we move on to illustrate solutions using the spectral method. Figure5 represents different volatility models executed by the Chebyshev method. In Figure 6, we compare solutions between the volatility models (transaction cost models) and the linear BS model. We plot the absolute difference between the transaction cost model and the linear Black-Scholes model in Figure 6. From Figure 6, we observe that the difference is significant when the time to expiration is 1 year (i.e., t=0) for all the volatility models [16, 27], where the stock price approximately lies between S=50 and S=130. That is, the nonlinear price (model with transaction cost) is significantly bigger than the linear price (model without transaction cost) at the very beginning of the time (i.e., t=0 or time to expiration 1 year) [16, 27]. But we note from this illustration that the difference decreases when the time to expiration decreases (i.e., t increases) [27].Figure 5
Different types of volatility models by Chebyshev spectral method. Here, we user=0.1, K=75, T=1, σ=0.2, R=1, N=40, m=2500, M=2, C=0.01, a=0.02, and Le=0.6.Figure 6
Here, we illustrateVNL−VL considering volatility models by Chebyshev spectral method, where NL stands for nonlinear model and L stands for linear model.Figure 7
2D plot of volatility models by Chebyshev spectral method att=0.From Figure7 we can see that for Chebyshev spectral method, we can arrange the volatility models according to highest to lowest prices at t=0 in the following serial: first, Leland’s price, next Barles’-Soner’s value, after that RAPM value, and finally Barles’-Soner’s (identity) price. In Table 2, we present the error VNL−VL of volatility models by the Chebyshev spectral scheme at t=0.Table 2
The absolute error of volatility models by Chebyshev spectral method att=0. Here, BSI means Barles’-Soner’s (identity).
Stock price (S)Linear BS model (exact values)RAPM valueErrorLeland’s valueErrorBarles’-Soner’s valueErrorBSI valueError37.515860.004360.007070.002710.042740.038370.009810.005440.005540.0011839.846830.012410.019270.006850.086490.074080.024110.01170.015490.0030742.496570.034780.051370.016590.174680.139900.061720.026940.042540.00776


## 5.2. Numerical Results of Linear Black–Scholes Model

Then, we start to apply the scheme for the linear version of the model. For the linear model,σ=0.2, r=0.1, N=Size ofDmatrix=60, m=Number of temporal points=6000, K=75, and T=1 year have been considered along with the European call option in Figure 8. The absolute difference between the exact solutions and the approximate solutions has been presented in Figure 9. It is we noticed from the illustration that the difference between exact option value and numerical option value is significant when S=73.0298 at time to expiration T−t=0.Figure 8
Linear Black-Scholes model by the Chebyshev spectral method (a) and by the differential quadrature method (b).
(a)(b)Figure 9
Here, we plotVExact−VNumeric of the linear Black-Scholes model. (a) One is for the Chebyshev spectral method, and (b) one is for the differential quadrature method.
(a)(b)It is evident from Figure10 that both the schemes provide very close results to the exact values. This makes these two schemes more reliable to approximate the Black-Scholes model. Though the entries in the Chebyshev spectral and differential quadrature methods are similar, their positions are just opposite to each other. Next, we present the approximate values of the linear BS model at t=0 in the following table where absolute error is considered. Here, by error we mean | exact value-approximate value |.Figure 10
2D plot of the linear Black-Scholes model att=0.We have also compared our proposed methods with a benchmark method FDM (finite difference method) to justify the accuracy of spectral methods (Chebyshev and DQM). The discrete solutions of equation (27) are denoted as uxi,τj=uij. In our methods, spatial grid points xi are nonuniform and temporal grid points τj are uniform. We replace the derivatives in convection-diffusion equation (27) to apply FDM [1, 47, 48].(55)uij+1−uijΔτ=2xi+1−xiui−1j−xi+1−xi−1uij+xi−xi−1ui+1jxi−xi−1xi+1−xixi+1−xi−1+1+Bui+1j−ui−1jxi+1−xi−1,i.e.uij+1=Δτ2xi+1−xiui−1j−xi+1−xi−1uij+xi−xi−1ui+1jxi−xi−1xi+1−xixi+1−xi−1+1+Bui+1j−ui−1jxi+1−xi−1+uij.Equation (55) is executed iteratively on time grids τj to approximate the solutions in case of FDM.We notice from Table3 that the approximations are similar for both Chebyshev spectral and differential quadrature schemes, and the errors for these two schemes are significantly small. These two methods approximate the Black-Scholes model way better than FDM. That is why it is well justified to say that spectral methods are reliable to approximate BS model.Table 3
The absolute error of Chebyshev spectral scheme, differential quadrature scheme, and FDM att=0.
Stock price (S)V_ExactV_ChebyshevErrorV_DQMErrorV_FDMError37.332210.003990.004110.000120.004110.000120.005780.0017938.821600.007980.008190.000210.008190.000210.010850.0028840.445830.015890.016270.000380.016270.000380.020410.0045342.211840.031360.032040.000670.032040.000670.038260.0068944.126800.060990.062150.001160.062150.001160.071050.0100646.198030.116210.118140.001930.118140.001930.130140.0139348.432970.215860.218930.003070.218930.003070.233920.0180750.839010.389190.393850.004660.393850.004660.410810.0216253.423420.678610.685330.006720.685330.006720.702020.0234156.193221.141091.150220.009131.150220.009131.163300.02221


## 6. Sensitivity Analysis

Option prices vary with respect to volatility, price of the underlying assets and some other factors [7, 49–51]. These types of changes are measured with some ratios that are expressed by Greek symbols [7, 49–51]. These ratios are called sensitivity, and there are various types of sensitivities such as delta, gamma, vega, theta, and rho [7, 49–51]. All these types are considered as the very significant index in case of option pricing [50]. For market analysis, sensitivity is considered as a powerful tool [50], which is computed through partial derivatives [7, 50]. We study the behaviors of the Greeks. By studying sensitivity of options, we can gather information to predict the situation of stock market [7]. Sensitivity computes the rate of change of option value whenever we change one input considering the other inputs are constant [7]. That is, we measure the partial derivatives of option price with respect to different underlying parameters [49–51]. We present some of the sensitivities (Greeks) for European call option. We have also computed the Greeks numerically with the help of finite difference method and compared this with the exact Greek values. This will surely define the accuracy of our proposed methods for computing the Greeks. For numerical approximation of Greeks, the option values VSi,tj=Vij are taken from Chebyshev and differential quadrature methods.


### 6.1. Delta

Delta,Δ=∂V/∂S=Nd1>0, measures the rate of change of the option value V with respect to the stock price S [7, 49, 51]. Delta varies from 0 to 1 [50]. For nonuniform spatial grid points Si, we apply the following FDM formula [47] to numerically compute delta.(56)ΔSi,tj=Δij=∂V∂Sij=Vi+1j−Vi−1jSi+1−Si−1.To observe the behavior of delta, the parameter values are taken as T=1, r=0.1, σ=0.2, K=75, N=Size of D matrix=40, and m=Number of temporal points=2500.From Figure11, we observe that when the time to maturity is 1 year, the delta value starts to increase when the stock price is S=40 and keeps increasing between S=40 and S=120 and then reaches to its maximum value. At time to expiration 0, delta value suddenly increases when S=70. From this observation, we can conclude that the delta is highly sensitive when the stock price lies between S=70 and S=90 at time to maturity 1 year.Figure 11
Sensitivity analysis (delta) by Chebyshev spectral method (a) and by differential quadrature method (b).
(a)(b)In Figure12, we notice that the closed-form values and numerical values of delta are very close to each other, which shows that our proposed methods are more accurate to study the sensitivity analysis of BS model.Figure 12
Here, we plot exact and numerical values of delta at timet=0. (a) is for the Chebyshev method, and (b) is for the DQM.
(a)(b)


### 6.2. Gamma

Gamma,Γ=∂2V/∂S2=∂Δ/∂S=N′d1/SσT−t=e−d12/2/Sσ2πT−t>0, determines the rate of change (sensitivity) of delta with respect to the underlying stock price [7, 49, 51]. We can predict the errors of delta hedging with the help of gamma [51]. We apply the following formula [48] to compute the second-ordered derivative ∂2V/∂S2 for the nonuniform spatial grid Si.(57)ΓSi,tj=Γij=∂2V∂S2ij=∂∂SΔij=2Si+1−SiVi−1j−Si+1−Si−1Vij+Si−Si−1Vi+1jSi−Si−1Si+1−SiSi+1−Si−1.The parameter values here are T=1, r=0.1, σ=0.2, K=75, N=Size of D matrix=40, and m=Number of temporal points=2500.We note from Figure13 that at the time to expiration is 0, and gamma values start to increase when the stock price is 60. At the same time, gamma reaches to its highest value when S stays at 70; then, gamma slowly keeps decreasing up to S=90; and subsequently, it becomes linear at the time to expiration 0. At maturity time 1 year, gamma initially increases up to S=65 and then decreases if S>65. So, major changes in gamma can be observed between S=60 and S=90. So, gamma is sensitive for S=60 up to S=90.Figure 13
Sensitivity analysis (gamma) by Chebyshev spectral method (a) and by differential quadrature method (b).
(a)(b)It can be observed from Figure14 that the deviation of numerical gamma values from the exact gamma values is very less.Figure 14
2D plot exact and numerical values of gamma at timet=0. (a) is for the Chebyshev method, and (b) is for the DQM.
(a)(b)


### 6.3. Rho

Rho,ρ=∂V/∂r=T−tKe−rT−tNd2>0, determines the rate of change of option price with respect to the risk-free interest rate r [7, 49–51]. We apply the FDM formula ρij=∂V/∂rij=Vijr+Δr−Vijr−Δr/2Δrto numerically evaluate rho values [52], where Δr=0.001. We fix the parameters as T=1, r=0.1, σ=0.2, K=75, N=Size of D matrix=40, and m=Number of temporal points=2500.We observe from Figure15 that initially rho stays linear up to 40 units for almost all the temporal points. But the major change in rho is visible when S is greater than 40 units at the time to expiration 1 year. We compute rho values by changing the interest rate and keeping the other parameters constant in Table 4.Figure 15
Sensitivity analysis (rho) by Chebyshev spectral method (a) and by differential quadrature method (b).
(a)(b)Table 4
Rho values for different interest rates.
Stock price (S)Time to expiration,T−tInterest raterRho numericρ78.08161 year0.0545.36560.149.15550.252.9749From Table4, we can observe that if we increase the interest rate, then there are also changes in rho values. Though the changes in rho for change in interest rate are not huge, we can conclude that rho is sensitive to the changes in interest rate.From Figure16, we can notice that in terms of accuracy, spectral methods (Chebyshev and DQM) are reliable as the exact values and numerical values of rho are almost same.Figure 16
We plot exact and numerical values of rho at timet=0. (a) is for the Chebyshev method, and (b) is for the DQM.
(a)(b)


### 6.4. Theta

Theta,Θ, evaluates the change in option value with respect to time [7, 49–51]. The formula for theta is as follows [49, 51]:(58)Θ=∂V∂t=−∂V∂g=−σSN′d12T−t−rKe−rT−tNd2=−σSe−d12/222πT−t−rKe−rT−tNd2<0,where g=T−t. For numerical computation, we apply the following FDM formula for uniform temporal point tj [1].(59)ΘSi,tj=Θij=∂V∂tij=Vij+1−Vijtj+1−tj.To analyze the sensitivity, we set the parameters as T=1, r=0.1, σ=0.2, K=75, N=Size ofDmatrix=40, and m=Number of temporal points=2500.We notice from Figure17 that theta value is negative, which usually happens for European call option [51].Figure 17
Sensitivity analysis (theta) by Chebyshev spectral method (a) and by differential quadrature method (b).
(a)(b)Figure 18
We plot exact and numerical values of theta forS=99.0580. (a) is for the Chebyshev method, and (b) is for the DQM.
(a)(b)From Figure18 it is evident that closed-form values and numerical values of theta are very close to each other. We also examine theta values by changing time to expiration and all the other remaining parameters as constant. Table 5 shows that when the time to expiration increases, then the theta value also increases, which guarantees the high sensitivity of theta.Table 5
Theta values for different time to expiration.
Stock price (S)Time to expiration,g=T−tTheta numeric,Θ99.05800.4134−7.39290.4158−7.3928


### 6.5. Vega

Vega,ϑ=∂V/∂σ=SN′d1T−t=ST−te−d12/2/2π>0, determines the change in option price with respect to the volatility [7, 49–51]. It is to be noted that call and put have the same vega [7]. We use the formula ϑSi,tj=ϑij=∂V/∂σij=Vijσ+Δσ−Vijσ−Δσ/2Δσ to numerically compute vega, where Δσ=0.001 [52]. To execute vega, T=1, r=0.1, σ=0.2, K=75, N=Size ofDmatrix=40, and m=Number of temporal points=2500 are considered as parameter values to vary vega.From Figure19, we analyze the change in vega values. At the time to expiration 1 year, vega value grows up for stock price S=40 to S=80 units.Figure 19
Sensitivity analysis (vega) by Chebyshev spectral method (a) and by differential quadrature method (b).
(a)(b)It can be noticed from Figure20 that the exact and numerical values of vega are almost same. For different volatilities, we inspect vega keeping the other parameters constant. From Table 6, it can be observed that when volatility increases, vega also increases along with it. This establishes that vega is highly sensitive to the changes in volatility.Figure 20
We plot exact and numerical values of vega att=0. (a) is for the Chebyshev method, and (b) is for the DQM.
(a)(b)Table 6
Vega values for different volatilities.
Stock price (S)Time to expirationg=T−tVolatilityσVega numericϑ174.63151 year0.20.0004640.30.2976


## 6.1. Delta

Delta,Δ=∂V/∂S=Nd1>0, measures the rate of change of the option value V with respect to the stock price S [7, 49, 51]. Delta varies from 0 to 1 [50]. For nonuniform spatial grid points Si, we apply the following FDM formula [47] to numerically compute delta.(56)ΔSi,tj=Δij=∂V∂Sij=Vi+1j−Vi−1jSi+1−Si−1.To observe the behavior of delta, the parameter values are taken as T=1, r=0.1, σ=0.2, K=75, N=Size of D matrix=40, and m=Number of temporal points=2500.From Figure11, we observe that when the time to maturity is 1 year, the delta value starts to increase when the stock price is S=40 and keeps increasing between S=40 and S=120 and then reaches to its maximum value. At time to expiration 0, delta value suddenly increases when S=70. From this observation, we can conclude that the delta is highly sensitive when the stock price lies between S=70 and S=90 at time to maturity 1 year.Figure 11
Sensitivity analysis (delta) by Chebyshev spectral method (a) and by differential quadrature method (b).
(a)(b)In Figure12, we notice that the closed-form values and numerical values of delta are very close to each other, which shows that our proposed methods are more accurate to study the sensitivity analysis of BS model.Figure 12
Here, we plot exact and numerical values of delta at timet=0. (a) is for the Chebyshev method, and (b) is for the DQM.
(a)(b)


## 6.2. Gamma

Gamma,Γ=∂2V/∂S2=∂Δ/∂S=N′d1/SσT−t=e−d12/2/Sσ2πT−t>0, determines the rate of change (sensitivity) of delta with respect to the underlying stock price [7, 49, 51]. We can predict the errors of delta hedging with the help of gamma [51]. We apply the following formula [48] to compute the second-ordered derivative ∂2V/∂S2 for the nonuniform spatial grid Si.(57)ΓSi,tj=Γij=∂2V∂S2ij=∂∂SΔij=2Si+1−SiVi−1j−Si+1−Si−1Vij+Si−Si−1Vi+1jSi−Si−1Si+1−SiSi+1−Si−1.The parameter values here are T=1, r=0.1, σ=0.2, K=75, N=Size of D matrix=40, and m=Number of temporal points=2500.We note from Figure13 that at the time to expiration is 0, and gamma values start to increase when the stock price is 60. At the same time, gamma reaches to its highest value when S stays at 70; then, gamma slowly keeps decreasing up to S=90; and subsequently, it becomes linear at the time to expiration 0. At maturity time 1 year, gamma initially increases up to S=65 and then decreases if S>65. So, major changes in gamma can be observed between S=60 and S=90. So, gamma is sensitive for S=60 up to S=90.Figure 13
Sensitivity analysis (gamma) by Chebyshev spectral method (a) and by differential quadrature method (b).
(a)(b)It can be observed from Figure14 that the deviation of numerical gamma values from the exact gamma values is very less.Figure 14
2D plot exact and numerical values of gamma at timet=0. (a) is for the Chebyshev method, and (b) is for the DQM.
(a)(b)


## 6.3. Rho

Rho,ρ=∂V/∂r=T−tKe−rT−tNd2>0, determines the rate of change of option price with respect to the risk-free interest rate r [7, 49–51]. We apply the FDM formula ρij=∂V/∂rij=Vijr+Δr−Vijr−Δr/2Δrto numerically evaluate rho values [52], where Δr=0.001. We fix the parameters as T=1, r=0.1, σ=0.2, K=75, N=Size of D matrix=40, and m=Number of temporal points=2500.We observe from Figure15 that initially rho stays linear up to 40 units for almost all the temporal points. But the major change in rho is visible when S is greater than 40 units at the time to expiration 1 year. We compute rho values by changing the interest rate and keeping the other parameters constant in Table 4.Figure 15
Sensitivity analysis (rho) by Chebyshev spectral method (a) and by differential quadrature method (b).
(a)(b)Table 4
Rho values for different interest rates.
Stock price (S)Time to expiration,T−tInterest raterRho numericρ78.08161 year0.0545.36560.149.15550.252.9749From Table4, we can observe that if we increase the interest rate, then there are also changes in rho values. Though the changes in rho for change in interest rate are not huge, we can conclude that rho is sensitive to the changes in interest rate.From Figure16, we can notice that in terms of accuracy, spectral methods (Chebyshev and DQM) are reliable as the exact values and numerical values of rho are almost same.Figure 16
We plot exact and numerical values of rho at timet=0. (a) is for the Chebyshev method, and (b) is for the DQM.
(a)(b)


## 6.4. Theta

Theta,Θ, evaluates the change in option value with respect to time [7, 49–51]. The formula for theta is as follows [49, 51]:(58)Θ=∂V∂t=−∂V∂g=−σSN′d12T−t−rKe−rT−tNd2=−σSe−d12/222πT−t−rKe−rT−tNd2<0,where g=T−t. For numerical computation, we apply the following FDM formula for uniform temporal point tj [1].(59)ΘSi,tj=Θij=∂V∂tij=Vij+1−Vijtj+1−tj.To analyze the sensitivity, we set the parameters as T=1, r=0.1, σ=0.2, K=75, N=Size ofDmatrix=40, and m=Number of temporal points=2500.We notice from Figure17 that theta value is negative, which usually happens for European call option [51].Figure 17
Sensitivity analysis (theta) by Chebyshev spectral method (a) and by differential quadrature method (b).
(a)(b)Figure 18
We plot exact and numerical values of theta forS=99.0580. (a) is for the Chebyshev method, and (b) is for the DQM.
(a)(b)From Figure18 it is evident that closed-form values and numerical values of theta are very close to each other. We also examine theta values by changing time to expiration and all the other remaining parameters as constant. Table 5 shows that when the time to expiration increases, then the theta value also increases, which guarantees the high sensitivity of theta.Table 5
Theta values for different time to expiration.
Stock price (S)Time to expiration,g=T−tTheta numeric,Θ99.05800.4134−7.39290.4158−7.3928


## 6.5. Vega

Vega,ϑ=∂V/∂σ=SN′d1T−t=ST−te−d12/2/2π>0, determines the change in option price with respect to the volatility [7, 49–51]. It is to be noted that call and put have the same vega [7]. We use the formula ϑSi,tj=ϑij=∂V/∂σij=Vijσ+Δσ−Vijσ−Δσ/2Δσ to numerically compute vega, where Δσ=0.001 [52]. To execute vega, T=1, r=0.1, σ=0.2, K=75, N=Size ofDmatrix=40, and m=Number of temporal points=2500 are considered as parameter values to vary vega.From Figure19, we analyze the change in vega values. At the time to expiration 1 year, vega value grows up for stock price S=40 to S=80 units.Figure 19
Sensitivity analysis (vega) by Chebyshev spectral method (a) and by differential quadrature method (b).
(a)(b)It can be noticed from Figure20 that the exact and numerical values of vega are almost same. For different volatilities, we inspect vega keeping the other parameters constant. From Table 6, it can be observed that when volatility increases, vega also increases along with it. This establishes that vega is highly sensitive to the changes in volatility.Figure 20
We plot exact and numerical values of vega att=0. (a) is for the Chebyshev method, and (b) is for the DQM.
(a)(b)Table 6
Vega values for different volatilities.
Stock price (S)Time to expirationg=T−tVolatilityσVega numericϑ174.63151 year0.20.0004640.30.2976


## 7. Conclusions

The Black-Scholes model is a mainstream tool for financial mathematics. In this article, we have mainly attempted to solve the widely applied option model (Black-Scholes equation) numerically. We have discussed the solution procedure for both the linear and the nonlinear BS model in detail. The solution procedure has been employed considering European options. We have applied Chebyshev spectral method and differential quadrature method on the nonlinear volatility models and the linear model. It is noticed that both the numerical schemes are very simple to implement for the model problems, and at the same time, we can rely on the Chebyshev method and the DQM to approximate the sensitivities of the Black-Scholes model. Numerical results obtained by these two numerical schemes demonstrate the efficiency of the schemes. We note that for a small number of spatial grid values/points, we get reasonable approximate solutions from both the schemes. For our problem, we obtain similar results for both the Chebyshev spectral scheme and the differential quadrature scheme. At the same time, we have attempted to analyze the sensitivity (Greeks) of the Black-Scholes model in this study. We observe from the study that different parameters have an impact on option pricing. Therefore, we can consider sensitivity analysis as one of the significant elements in stock market. It is noted from the favor of the efficiency of the schemes that a small number of spatial grid values were enough to analyze the sensitivity. This fact will surely encourage the researchers to use the nonlinear BS model along with the numerical schemes presented here to analyze market data. From the computational and theoretical studies, we notice an order one error in temporal integration, which reduces the efficiency of spatial scheme. So spectral scheme for spatial approximation along with a higher-order temporal scheme is recommended for the BS model. A continuation of the same work on American options along with an application of higher-accurate time integration schemes left as future research directions.



---

*Source: 4488082-2022-08-24.xml*