A New Approach for the Black–Scholes Model with Linear and Nonlinear Volatilities

: Since ﬁnancial engineering problems are of great importance in the academic community, effective methods are still needed to analyze these models. Therefore, this article focuses mainly on capturing the discrete behavior of linear and nonlinear Black–Scholes European option pricing models. To achieve this, this article presents a combined method; a sixth order ﬁnite difference (FD6) scheme in space and a third–order strong stability preserving Runge–Kutta (SSPRK3) over time. The computed results are compared with available literature and the exact solution. The computed results revealed that the current method seems to be quite strong both quantitatively and qualitatively with minimal computational effort. Therefore, this method appears to be a very reliable alternative and ﬂexible to implement in solving the problem while preserving the physical properties of such realistic processes.


Introduction
In the last few decades, the corporations looked for important tools in terms of financial securities. As part of the financial securities, options are mainly used to assure assets in order to cover the risks generated in the stock prices changes [1]. To properly understand what alternatives we have for these options, Lesmana and Wang [2] stated that there are two main types of options: European options can only be exercised on a given date (expiry date) and American options can be exercised on or before the expiry date.
An economic issue, which is important for both the traders and the investors, is to have a proper method for price options to determine the appropriate price or theoretical value for a call or a put option. During the time, price options were not so commonly explained by the traders as a sufficient financial instrument and valuation of an option has always been a difficult and challenging task. According to Rad et al. [3], options are priced using mathematical models that are often challenging to solve. Of course, in the last 50 years, there have been some attempts to introduce different methods for price options. For instance, a major step was taken by Black and Scholes [4][5][6] in 1973 to propose a mathematical model for calculating a fair value for an option. According to them [4], the Black-Scholes model is a pure log-normal diffusion model and it leads to a parabolic partial differential equation under Ito's calculus.
Another important contribution was made by Merton [7], who extended the model equation proposed by Black and Scholes. All three of them demonstrated that their formulae leading to partial differential equations could help to determine a fair value for a call or a put option. One of the main obstacles identified by Lesmana and Wang [2] is that the Black-Scholes option pricing methodology is no longer valid in the presence of transaction costs on trading in the risk-free security or stock. To overcome this drawback, different models were proposed by Leland [8], Boyle and Vorst [9], Kusuoka [10] or Barles and Soner [11]. In addition, Cox and Ross [12] published an article aimed at finding ways to overcome difficulties with payouts and bankruptcies. Another important contribution is mainly related to the Crank-Nicolson scheme used in numerical pricing within the Black-Scholes framework [13]. A good example of this scheme is the work of McCartin and Labadie [14] in the case of pricing vanilla options. In addition, Ankudinova and Ehrhardt [15] used a Crank-Nicolson method and devised a numerical scheme for the linearized Black-Scholes equation using the frozen values of nonlinear volatility. From the point of view of nonlinear volatility, an approach based on the method of lines and the backward Euler scheme was proposed by Company et al. [16]. Other than these, in recent years, Jeong et al. [17], Koleva et al. [18], Manisha [19] and Sari and Bȃlȃcescu [20], Mashayekhi and Hugger [21] and Hout and Valkov [22] used various schemes based on finite differences to analyze the Black-Scholes model. Jeong et al. [17] proposed a finite difference scheme not using the far-field boundary condition. Koleva et al. [18] construct a fourth order compact finite difference scheme in space for European options modeling markets to liquidity shocks. A high-order difference scheme was considered by Manisha [19] for the Black-Scholes equation, which governs four option styles of European-type variable parameters. To analyze the option pricing model, Sari and Balacescu [20] proposed a fourth-order difference and MacCormack schemes. One of the most comprehensive models, the Barles and Soner [11] model was solved with various difference schemes to compare their precision and order of convergence [21]. Hout and Valkov [22] solved the European two-asset options by a finite difference based numerical method considering non-uniform grids. In addition, Hendrick et al. [23] combined a high-order finite difference scheme with the Alternating Direction Implicit scheme for the parabolic partial differential equation representing the European basket option in a sparse grid setting.
Another work by Markolefas [24] demonstrated the ability of a finite element method to formulate discrete approximation models resulting from the Black-Scholes pricing model. Specifically, his/her work demonstrates a version of Galerkin finite element methods to capture accurate and efficient behaviour: the parabolic partial differential equations, with a complex initial, boundary and/or internal conditions, resulting from various option pricing theories. However, Lin and Chen [25] presented a European option pricing model by applying the Model-Order-Reduction method, which is able to reduce the order of the original finite difference method systems.
Likewise, in the case of a nonlinear parabolic partial differential equation that governs the European option pricing in transaction costs, it is worth pointing out the contribution to Lesmana and Wang [2] who proposed an upwind difference scheme for spatial discretization. In addition, there are more than 30 years since the two experts of banking and financial markets developed the first commercially available pricing formula for options linked to the average price of crude oil; they called this option as Asian option [26]. The Asian option, also known as average options, is an option where the return on return is dependent on the arithmetic average of the price of the underlying asset over its lifetime [27][28][29][30]. Kumar et al. [27] described the implementation of the local radial basis function based on a grid-free method for the numerical solution of the Asian option. In addition, Rad et al. [28] described the valuation of the Asian option with a radial basis function approach based on the finite difference method.
In another representative work created by Sin [31], he/she used some numerical methods to compute barrier options considering both Black-Scholes and Heston models. In addition, he/she computed barrier options using the standard Monte Carlo method. The method was improved by Moon [32] in dealing with the problem. The key idea behind his/her approach, the improved one, is to use the exit probability and uniformly distributed random numbers to predict the first instance of hitting the barriers. In addition, Zeng et al. [33] focused on a comparative study of the Monte Carlo simulation and finite difference method for the European option pricing under a regime-switching framework.
In spite of the great efforts made to deal with the problem in the literature, the researchers chose to use relatively low order methods [34][35][36][37][38][39] for the problem.
In recent years, the option market has become increasingly important. This increases the need for higher computational capacity as the complexity of financial models increases in an operating area where real-time market calculations are required. The well-known Black-Scholes model provides a closed form for the values of certain options such as European style. However, for the solutions of more complex models, either there is no closed form solution or, in the absence of tractable closed form, computationally intensive approximate techniques are required. In this case, the finite difference based methods are commonly employed in price options. Therefore, in this article, a sixth order finite difference scheme (FD6) in space and a third-order strong stability preserving Runge-Kutta (SSP-RK3) in time have been combined to obtain effective numerical solutions of the European put option problem that has an exact closed-form solution. This choice provides a direct and accurate estimation of approximation error. High-order approximations reduce the computational effort to achieve the required accuracy. Therefore, they are highly popular in solving real world problems and for all financial markets increased calculation speed is a strategic advantage. Moreover, the choice of time integration is important in many respects such as accuracy, computational effort and stability. The SSP-RK3 scheme, which is a class of high order SSP-time integration techniques, was developed by Gottlieb et al. [40] to investigate hyperbolic conservation laws with a spatial discretization guaranteeing the stability properties expected of the forward Euler method [40]. Therefore, it was considered that the combination of the FD6 and SSP-RK3 could solve the corresponding model with high accuracy. In addition, this method is quite straightforward to write codes in any programming language. Our computations show that the results produced by the current method approximates the exact solution and available solution in the literature [41] very well. To the best knowledge of the authors, this method has not been implemented for the problem represented by the European put option pricing model.
The rest of the paper is organized as follows: in Section 2, the mathematical formulation of the European put option problem is introduced. In Section 3, numerical methods for space and time discretizations are summarized. The results and illustrations are presented in Section 4. Final remarks are reported in Section 5.

Mathematical Formulation
The Black-Scholes model is of great importance in option pricing theory. The model was developed by Black and Scholes [4] and previously by Merton [7]. The value of an option V(S, t) is a function of stock price S and time t can be determined by the following Black-Scholes partial differential equation, where and r ≥ 0 and σ are positive constant parameters, referred to risk-free interest rate and volatility that are both known functions over the life of the option. Here, the stock price S is modeled by a Geometric Brownian motion, i.e., S satisfies the following stochastic differential equation where µ > 0 is the drift rate of the stock and W is a standard Brownian motion. The Black-Scholes equation is given by Equation (1), which assumes that the market is frictionless, is a trading enviroment and there are no transaction costs and restraints on transactions. In addition, Equation (1) does not include arbitrage opportunities.
The Black-Scholes Equation (1) is very effective for modeling pricing options in a complete market without transaction costs. However, in the presence of transaction cost on trading in the riskless security or stock, the Black-Scholes option pricing methodology is no longer valid, since perfect hedging is impossible [2]. Due to transaction cost where the volatility σ can depend on the time t, the stock price S or the derivatives of the option price V itself, the Black-Scholes model (1) is transformed to the following nonlinear equation: where σ is the modified volatility as a function of t, S, V S , V SS . The problem represented by the model equation is very complex and depends on too many parameters, which makes the model more realistic, and, in this case, careless use of those parameters can lead to unrealistic results. Among these parameters, volatility is one of the most considerable factors in the pricing model. Volatility must be considered carefully, as the option value depends on the future stock price. The option pricing model with transaction cost has been studied in the literature by some researchers, such as Boyle and Vorst [10] proposed an option price with the volatility of the form where µ is the proportional transaction cost, ∆t is the transaction period and σ 0 is the original volatility constant.
Leland [8] obtained that the option price is the solution to Equation (2) with the volatility where Le is the Leland number given by where δt and κ represent the transaction frequency and transaction cost measure, respectively. A more complex model was proposed by Barles and Soner [11]. In their model, the nonlinear volatility is given by where a = κ √ γN with κ being the transaction cost parameter, γ being a risk aversion factor and N being the number of options to be sold. The function Ψ is the solution to the following nonlinear initial value problem: Equation (7) implies that lim x→∞ Ψ(z) z = 1 and lim This property accepts the treatment of the function Ψ(.) as the identity for large arguments and therefore the volatility becomes For solving Equations (1) and (3) uniquely, one final condition and two boundary conditions are required. These conditions depend on both types of options (put or call) and whether or not the stock pays cash dividends. The European put option, which gives the holder the right to sell for the strike price K at maturity time T, is studied in this work, and the boundary and final conditions are as follows: The existence and uniqueness of the solution of the problem were shown in the context of stochastic optimal theory in the study of Barles and Soner [11].
The exact solution of the linear European put option problem is determined by with the parameters where N is the standard normal cumulative probability distribution function and δ is continuous dividend yield [42].

The Solution Method
This section is dedicated to a numerical solution of the linear and nonlinear European put option problems. For the sake of brevity, the nonlinear European put option model is discussed with the nonlinear volatility (9) proposed by Barles and Soner [11]. In spatial discretization, a sixth-order finite difference method (FD6) is applied and a third-order strong stability preserving a Runge-Kutta (SSP-RK3) method is considered in temporal discretization. The FD6 scheme based on high order differences achieves the required accuracy. In addition, the SSP-RK3 scheme, which is a class of the SSP-time integration techniques with the spatial discretization, guarantees the stability [40]. Hence, it is expected that the proposed method will solve the problem effectively.

Spatial Discretization
Spatial derivatives are computed by the FD6 scheme based on the Taylor series expansion [43]. In order to use finite difference approximation, it is started by defining a uniform grid consisting of N points satisfying 0 = x 0 < x 1 < ... < x N = x max . The step size h = x i+1 − x i , i = 0, 1, ..., N − 1 is equal to each other at any point i. The first derivative V i at point i can be approximated by where R and L represent the number of points on the right-hand side and left-hand side for taken nodes, respectively. Hence, Equation (14) involves (R + L + 1) constants, a 0 , a 1 , ..., a R+L . R is equal to L for the considered nodes at internal points, but this is not the case for the boundary points. The coefficients a j are calculated by using the Taylor expansion near point i. The FD6 scheme uses seven points at interior and boundary points as follows: The sixth order scheme (15) can be written in a more compact form as follows: The second order derivative terms are obtained by applying the first operator twice: where V = (V 1 , V 2 , .., V N ) T and

Temporal Discretization
The SSP-RK3 scheme for the discretization of Equation (1) is presented in time. A class of high-order SSP time discretization techniques was introduced by Gottlieb et al. [40] for solving hyperbolic conservation laws with stable spatial discretizations. The SSP methods guarantee the stability properties expected of the forward Euler method [40].
The computational domain for time consists of M points satisfying 0 = t 0 < t 1 < ... < t M = T. The uniform time step dt = t n+1 − t n , n = 0, 1, 2, ..., M − 1 is equidistant at any point n. After applying the FD6 method, Equation (1) can be reduced into a set of ordinary differential equations in time. Then, the equation system can be expressed by The SSP-RK3 scheme integrates the semi-discrete Equation (18) from time t 0 to t 0 + dt through the operations where L is the discretization form of L . For the SSP-RK3, the total variation (TV) of the numerical solution does not increase in time [40];

Numerical Results and Analysis
To show efficiency and accuracy of the proposed methods, numerical experiments on the European put option model have been performed. For the computations through the current schemes, computer codes have been produced in MATLAB R 2019. The validity of the scheme has been verified through the produced results.
Semi-infinite domain [0, +∞) in Equation (1) is replaced by the bounded interval [0, S max ) where S max is an artifical limit will be chosen larger than three to four times the exercise price. In order to use the difference approximation, it is started by defining a uniform grid consisting of N points satisfying 0 = S 0 < S 1 ... < S N = S max . The step size dS = S i+1 − S i , i = 0, 1, ..., N − 1 is equal to each other at any point i. Similarly, the interval (0, T) is divided into M subintervals with mesh nodes satisfying 0 = t 0 < t 1 < ... < t M = T. The step size dt = t n+1 − t n , n = 0, 1, 2, ..., M is equal to each other at any point n. In the nonlinear model, the square of volatility approximation σ 2 given by Equation (9) has been discretized at the node S i and time t, being where and and has been used in the discretized form of Equation ( into consideration in the selection of parameters. In terms of the realistically considered parameters, the computed results are believed to be helpful in making a decision for a viable trading strategy.
The numerical results and exact solution for the linear and numerical results for the nonlinear European put option model are presented in Table 1. The computed results revealed that the current method approximates the exact solution very well and is applicable, effective and easy to use. The qualitative behaviour of the option value is plotted in Figures 1-3. Furthermore, the values of the option for the three different values of a are represented in Figure 4. From the figure, it can be seen that qualitative behavior is in agreement with the literature [2].     After applying the methods to the linear Black-Scholes equation and initial and boundary conditions, the resulting system has M N unknowns. Since the present methods are explicit, the stability condition is given by where a = 1 2 σ 2 S 2 is the coefficient in front of the second derivative term in the Black-Scholes equation. In the current computation, this condition is provided at each time and spatial step. In the nonlinear model, the nonlinear coefficient σ i is frozen at each discretized time level, therefore, the condition (25) is provided.
To show accuracy of the method for the linear model, convergence rates are calculated by where V dt ds represents the solution with the spatial mesh size ds and the time mesh size dt, V exact represents the exact solution for the linear model whilst || . || ∞ and || . || 2 are maximum norm and L 2 -norm are given, respectively, by To determine these rates, a sequence of meshes generated by halving the mesh sizes of the previous ones by starting from a given coarse mesh is accepted. As seen in Table 2, the current method is seen to be effective to enhance the accuracy of the numerical solution. When M and N are halved, the error reduced by a factor that is equal to ten. The order of convergence of the combination of FD6 and SSPRK3 is about 6.6 in || . || ∞ and 4.4 in || . || 2 . For the nonlinear model, since the exact solution is unknown, taking fixed value of dt, the difference V(S, t, ds = k) − V(S, t, ds = k/2) is plotted in Figure 5. In addition, Table 3 presents the errors of the accuracy of the linear and nonlinear models that are calculated by the 2 -error given by Reference solution is denoted by v(x i , 0), which has been obtained by the proposed method on a fine grid dt = 6.25 × 10 −5 and ds = 0.1. The numerical solution V 0 i has been calculated with the previously indicated parameters.  To indicate the influence of transaction cost modeled by volatility (9), the difference V nonlinear (S, t) − V linear (S, t) between the price of the European put option with transaction costs and the price of the European put without transaction costs is plotted in Figure 6. The numerical results present an economically significant price deviation between the linear model and the nonlinear model. It is seen that the difference is not symmetric but decreases closer to the expiry date and is maximal close to the stock price S = 10, where the nonlinear price is higher than the linear price. In Table 4, the required times are compared when the proposed method is applied for the linear and nonlinear option problems. For both situations, the required CPU times increase with respect to the number of spatial and temporal discretizations. Although part of the computational time pertains to the calculation of the nonlinear volatility, the CPU times of the nonlinear model for every N and M discretization point are close to the linear model. The present work provides effective estimates in terms of easy programming and relatively low cost.

Conclusions and Recommendations
The complexity and stochastic properties of option pricing problems make it difficult to determine the value of the option. Therefore, accurate approximation methods are needed to understand the behavior of these problems and those are of great importance for scientific developments in financial markets. In this paper, a combination of a sixth-order finite difference scheme in space and a third-order strong stability preserving Runge-Kutta in time has been implemented to obtain effective numerical solutions of the linear and nonlinear European put option models represented by the Black-Scholes equation. The convergence of the solution has been measured by some error norms and it is confirmed that the present method is asymptotically convergent. In addition, the produced results are in good agreement with the literature and the exact solution. Therefore, the proposed method provides a better perspective to describe behavior of the option pricing model represented by the Black-Scholes equation and can be preferred due to reliability and accuracy with minimal computational effort. This study is believed to help researchers who want to uncover challenging financial and stochastic behaviours in modeling. In terms of the realistically considered parameters, it is believed that the computed results are helpful in making decisions for viable trading strategies. The valuation of the European option pricing problem by the proposed method can be compared with realistic option values.