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

e signicance of both the linear and nonlinear Black-Scholes partial dierential equation model is huge in the eld of nancial 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. is 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 nite dierence and nite element-based schemes to approximate the BS model. us, 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 specic, the Chebyshev spectral method and the dierential quadrature method are employed to approximate the BS models to see the eciency of such highly accurate schemes for the option pricing model. First, we approximate the model using the mentioned methods. en, we move on to use the numerical results to analyze dierent aspects of stock market through sensitivity analysis. All the numerical schemes have been illustrated through some graphics and relevant discussions. We nish the study with some concluding remarks.


Introduction
e Black-Scholes parabolic partial di erential equation model is considered to be one of the major appliances for option pricing as it is e cient to compute the price of an option more correctly [1][2][3][4][5]. An option is a nancial derivative, and it is a contract that gives the owner the right to buy or sell the asset for a speci ed price within a speci ed time [4,[6][7][8][9]. e xed price is called the strike price or exercise price, and it is denoted by K; and the xed time is called the expiration time or the maturity time, and it is denoted by T [3,5,7,8]. ere are mainly two types of options available in the derivative market, namely, call option and put option [8]. e 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].
ough there are various types of options [8], the most commonly applied options are Asian option, European option, and American option [4,7,8,11]. e 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 exibility in American option. e 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][4][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]. e Black-Scholes (BS) equation depends on two independent variables; the rst is the time; and the second is the stock price [4,5]. e 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]. rough 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]. at is why to get a unique solution, we impose initial/terminal and boundary conditions in the BS model [4,9,13].
ere 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]. e nonlinearity is considered as one of the results of transaction costs [4,9]. e nonlinear BS model is reduced to a nonlinear parabolic type equation [4,9,19].
ough 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]. e 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]. at 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 use O(Δx 2 ) 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. e 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. us, the sensitivities of parameters of the BS model are represented by different Greek alphabets. ese parameters play a decisive role in the definition of the shape of the option price. e rest of the article has been organized in the following way: (i) Section 2 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 Section 3. (iii) A brief discussion about the convergence of the proposed scheme has been included in Section 4. (iv) en, we move on to present a numerical discussion on the BS model in Section 5. (v) Section 6 contains sensitivity analysis for European option. (vi) We summarize this study in Section 7.

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. Let S 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 where W represents the stochastic variable. W and dW work as uncertainty in stock market where W(t) is random. μdt is the expected value, and σ 2 dt is the variance. Using Ito's lemma, we have Now using some changes in variables and some simple calculus, one can derive an uncertainty ΔW and W free riskfree model as which is popularly known as the Black-Scholes partial differential equation. From equation (3), one may notice that the derivatives (zV/zt), (z 2 V/zS 2 ), and (zV/zS) must exist. Let call option is big and put option is small. en, we have 2 Journal of Mathematics Again, from the terminal condition of the Black-Scholes model, we have From the last two equations, it is evident that which is known as the put-call parity, and the formula is as follows: e figure below shows the payoff function for both European call option and European put option for the inputs K � 40, 0 ≤ S ≤ 80. ree situations are visible in Figure 1, which are outthe-money, at-the-money, and in-the-money. At "out-themoney" and "at-the-money" situations, the holder does not exercise the option because cash flow S(T) − K is negative and zero, respectively [4]. On the other hand, at "in-themoney" position, the holder exercises the option as the cash flow S(T) − K is positive [4].
In the nonlinear BS model, we assume that drift μ is constant and volatility σ is nonconstant [4,16,23,24], which is defined as σ 2 � σ 2 (t, S, V S , V SS ) [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. e nonlinear Black-Scholes equation is as follows [16,23,[25][26][27]: e domain for the problem (8) e terminal and boundary conditions of the nonlinear Black-Scholes model for European call option are defined in the following way. e terminal condition is as follows: e boundary condition is as follows: 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. e 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.

Risk-Adjusted Pricing Methodology Model (RAPM Model).
e volatility of this model is as follows [4,9,16,[25][26][27]: 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].
With the changes in variables, it yields Journal of Mathematics and thus, the nonlinear BS model (8) transforms into which can be written by the following simplified form Let B � 2r/σ 2 in equation (19) can be integrated to find the solution u(x, τ) for the domain x ∈ R (i.e., − ∞ < x < ∞) and 0 ≤ τ ≤ T � (σ 2 T/2); then, one may return to the original solution by backward substitution. is 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]: e volatility of Barles' and Soner's model is as follows [16]: and the volatility of the RAPM model is as follows [16,26]: e initial and boundary conditions for European call option are defined in the following way [16,26]. e initial condition is u( 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]. e 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].
at 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]: e analytic solution to the linear BS equation can be obtained with the help of Fourier transformation from heat equation (13). e solution that we get is unique [5]. e theoretical solution to the BS model for European call option is given by and for European put option is given by where N(d) is the standard normal distribution of the function d [4,7,13] with In this case, we apply similar transformations in equation (12) of the nonlinear BS model. en, we get the following convection-diffusion equation for European call option.
where B � 2r/σ 2 . e initial and boundary conditions for the linear BS model are defined follows [16,26]. e initial condition is u(

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. e finite difference method and the linear finite element method are popular schemes to serve the purpose with O(Δx 2 ) accuracy.
us, one needs sufficient spatial grid points to confirm the desired accuracy of the solutions [22,[29][30][31][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. e 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. ere 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. is work focuses on to develop two exponential accurate schemes for the spatial approximation. To that end, we begin with replacing x ∈ R 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 x i is defined for − 1 to 1. at 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 firstand second-order derivatives in space, and we replace these derivatives by the spectral method.
at is, (z d /zx d )u is replaced by D d u, where d is the order of the derivative, and D is a spectral matrix.

Chebyshev Spectral Method.
e 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.
ere are several methods resorted to solve the Black-Scholes equation.
ese methods are successful for their framework. e spectral method is one of them [34]. e 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].

Journal of Mathematics
So, the 2nd (k � 2) column of D N matrix is In general, the k th column of D N matrix is computed by putting the Chebyshev points at the k th term of the polynomial p ′ (x) [22,33]. ere are formulas for the entries of D N matrix [22,29,33,34]. e explicit methods are helpful for the Chebyshev differentiation matrix [22]. Let N ≥ 1. en, the rows and columns of the Chebyshev differentiation matrix D N can be obtained for index 0 to N [22,29]. e entries of D N matrix for the first-order derivative can be written as follows [22,29,33,34]: Here, e spatial derivative z/zx is approximated by the matrix operator D. For example, the second-order derivative (z 2 /zx 2 ) is computed by D 2 , i.e., we square D, which a is very simple and easily applicable procedure.

Differential Quadrature Method.
e differential quadrature method (DQM) is a simple and straightforward technique to implement [36,37]. e DQM is known for its simplicity, accuracy, and efficiency [36][37][38]. We need a very little computational endeavor to apply DQM [36][37][38]. We can use DQM to solve the initial value problem and boundary value problem. is method is significant to solve critical problems in option pricing. So this method can be a suitable choice for option pricing problems. e DQM is similar to the collocation or the pseudospectral method [36]. e DQM was developed by Richard Bellman and his associates Kashef and Casti in 1971 [36,38]. e DQM is used to discretize spatial derivatives. For our problem, z/zx and z 2 /zx 2 are replaced by differential quadrature matrix D and D 2 , 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].
In the DQM, we need weighting coefficients [36-38, 40, 41]. We use an approximating function to find the weighting coefficients [41]. is approximate function is known as test or trial functions. We approximate the spatial derivatives at the grid point x i in the following manner [36][37][38]40]. To that end, it is considered that and x − x j , en, the unknown function can be approximated by where x j is the discrete sampling points in the spatial domain of computation. e derivatives of such a function at x � x i can be presented by 6 Journal of Mathematics where A ij . Now we express the exact values of approximate derivatives and start by presenting entries of the first derivative using the base function p k (x), k � 1(1)n. ere 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 In a similar fashion, the formula for weighting coefficients of the second-order derivative is presented below 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 For further details, one may consult [36,41]. For this study, it After spatial integration, equations (19) and (27) can be written as the following semidiscrete system of equations of the form where A is a matrix, and u is a solution vector that depends on τ. e equation given by equation (45) is decomposed by δ i 's and its time derivatives _ δ i 's by the following Crank-Nicolson formula and usual finite difference approximation is approximation leads to the following recurrence relationship between two time levels n and n + 1 relating to δ n+1 i and δ n i . at is, we replace the derivatives z/zx and z 2 /zx 2 in equation (19) and in volatility equations (20)-(22) by D and D 2 matrices, respectively, which are taken from the Chebyshev spectral method and the differential quadrature method. Equation (19) is iteratively executed for each temporal integration.

Computation Algorithms.
We investigate the convection-diffusion equation (19) for the transformed mesh/domain [− 1, 1] and 0 ≤ τ ≤ T � (σ 2 T/2) [16,26]. We consider the grid points as (x i , τ n ), where spatial points x i 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.
at is, equation (27) is solved by applying the Chebyshev spectral scheme and the differential quadrature method in space with the explicit FDM in time. e 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]. e explicit FDM method has truncation error of O(Δτ) [1,21,42,43]. We replace the derivatives z/zx and z 2 /zx 2 in equation (27) by D and D 2 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.
e nonlinear BS model.
Step 2. We replace u xx and u x by D 2 u and Du, respectively, in where L � (σ 2 /σ 2 )(u xx + u x ) + Bu x .
Step 3. We replace the time derivative u τ in equation (48) by the FDM formula u τ � (u n+1 i − u n i )/(Δτ) or the Crank-Nicolson method to get Step 4. We iteratively execute equation (49) for each time grid to approximate the solution.
Step 5. After finding the solution to the convectiondiffusion equation (19), we return back to the original Journal of Mathematics 7 solution to the nonlinear model with the help of the backward substitution equation (16).
Algorithm 2. e linear BS model.
Step 1. We replace the derivatives u xx and u x by D 2 u and Du, respectively, in where L � u xx + (1 + B)u x .
Step 2. We replace the time derivative u τ in equation (50) by the formula u τ � (u n+1 i − u n i )/(Δτ) or the Crank-Nicolson scheme to achieve 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).

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 onestep scheme for the temporal approximation. It is to note that the performance of a numerical scheme depends on its error reduction and computational stability. e 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][45][46] and many other well-known references for detail. It is to mention that we employ some constants C i ≥ 0 in this study, which are necessarily not the same for all the cases. e 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 O(1). Near the point where it breaks the continuity, the error occurs of O(1/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 I N u be the interpolating polynomial. For the spectral method, the difference between u and I N u exponentially goes to 0. Let C be a function of u. en, upper bounds of ‖u − I N u‖ can be written as If we increase the power of m, then these upper bounds decrease. at is, if we increase the power of N, then the error becomes less, i.e., the difference |u − I N u| ⟶ 0 if N increases.
en, 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 is O(N − m ) for every m for functions that are smooth and O(c N ), 0 < c < 1, for every analytic function. us, if one approximates a smooth function w(x) by spectral polynomials (w h (x)), then it yields For time, a forward difference scheme has been used, which is accurate of O(Δτ) in L 2 ([0, T]) norm for some T > 0 [22,46]. So for the space-time discretization, the error bound is of the form for a suitable C 1 ≥ 0 and C 2 ≥ 0.

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).

Nonlinear Black-Scholes
Model. e graphs of the volatility models of European call option are presented here. At time t � 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.
In Figure 3 we present the difference between European call option value with the nonlinear model and European call option value with the linear model. is 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 Figure 4 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. at 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. en, we move on to illustrate solutions using the spectral method. Figure 5 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. at 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].
From Figure 7 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 |V NL − V L | of volatility models by the Chebyshev spectral scheme at t � 0.   � 6000, K � 75, and T � 1 year have been considered along with the European call option in Figure 8. e 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.

Numerical Results of Linear
It is evident from Figure 10 that both the schemes provide very close results to the exact values. is makes these two schemes more reliable to approximate the Black-Scholes model. ough 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 |.
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). e discrete solutions of equation (27) are denoted as u(x i , τ j ) � u j i . In our methods, spatial grid points x i are nonuniform and temporal grid points τ j are uniform. We    replace the derivatives in convection-diffusion equation (27) to apply FDM [1,47,48].    Figure 9: Here, we plot |V Exact − V Numeric | of the linear Black-Scholes model. (a) One is for the Chebyshev spectral method, and (b) one is for the differential quadrature method.
Equation (55) is executed iteratively on time grids τ j to approximate the solutions in case of FDM. We notice from Table 3 that the approximations are similar for both Chebyshev spectral and differential quadrature schemes, and the errors for these two schemes are significantly small. ese two methods approximate the Black-Scholes model way better than FDM. at is why it is well justified to say that spectral methods are reliable to approximate BS model.

Sensitivity Analysis
Option prices vary with respect to volatility, price of the underlying assets and some other factors [7,[49][50][51]. ese types of changes are measured with some ratios that are expressed by Greek symbols [7,[49][50][51].
ese ratios are called sensitivity, and there are various types of sensitivities such as delta, gamma, vega, theta, and rho [7,[49][50][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]. at is, we measure the partial derivatives of option price with respect to different underlying parameters [49][50][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. is will surely define the accuracy of our proposed methods for computing the Greeks. For numerical approximation of Greeks, the option values V(S i , t j ) � V j i are taken from Chebyshev and differential quadrature methods.  [50]. For nonuniform spatial grid points S i , we apply the following FDM formula [47] to numerically compute delta.
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 Figure 11, 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.
In Figure 12, 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.
) > 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 z 2 V/zS 2 for the nonuniform spatial grid S i .
We note from Figure 13 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.
It can be observed from Figure 14 that the deviation of numerical gamma values from the exact gamma values is very less.    We observe from Figure 15 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.
From Table 4, we can observe that if we increase the interest rate, then there are also changes in rho values. ough 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 Figure 16, 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.

eta.
eta, Θ, evaluates the change in option value with respect to time [7,[49][50][51]. e formula for theta is as follows [49,51]: where g � T − t. For numerical computation, we apply the following FDM formula for uniform temporal point t j [1].
To analyze the sensitivity, we set 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 notice from Figure 17 that theta value is negative, which usually happens for European call option [51].
From Figure 18 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.
) > 0, determines the change in option price with respect to the volatility [7,[49][50][51]. It is to be noted that call and put have the same vega [7]. We use the formula ϑ(S i , t j ) � ϑ to numerically compute vega, where Δσ � 0.001 [52]. To execute vega, T � 1, r � 0.1, σ � 0.2, K � 75, N � Size of D matrix � 40, and m � Number of temporal points � 2500 are considered as parameter values to vary vega.      Figure 19: Sensitivity analysis (vega) by Chebyshev spectral method (a) and by differential quadrature method (b).  Figure 17: Sensitivity analysis (theta) by Chebyshev spectral method (a) and by differential quadrature method (b). 16 Journal of Mathematics From Figure 19, 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.
It can be noticed from Figure 20 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. is establishes that vega is highly sensitive to the changes in volatility.

Conclusions
e 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. e 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. erefore, 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. is 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.

Data Availability
No data were used to support the findings of the study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.