**Journal of Mathematical Finance**

Vol.08 No.02(2018), Article ID:84599,10 pages

10.4236/jmf.2018.82024

A Study on Numerical Solution of Black-Scholes Model

Md. Nurul Anwar^{1,2*}, Laek Sazzad Andallah^{1}^{ }

^{1}Department of Mathematics, Jahangirnagar University, Savar, Dhaka, Bangladesh

^{2}Department of Natural Science, BGMEA University of Fashion & Technology, Dhaka, Bangladesh

Copyright © 2018 by authors and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

Received: March 19, 2018; Accepted: May 15, 2018; Published: May 18, 2018

ABSTRACT

In the history of option pricing, Black-Scholes model is one of the most significant models. In this article, the main concern is the numerical solution of the Black-Scholes model (a.k.a. Black/Scholes/Merton) for the European call option in a different way. The model is described and an explicit difference scheme was used for the numerical approximation. The stability condition of the scheme is established through convex combination. A different way was used to obtain the numerical value of the model. Estimation of the relative error was calculated in L_{1}-norm in order to test the accuracy of the scheme. Finally, a comparison of the numerical outcomes with the value obtained by another scheme is given.

**Keywords:**

Black-Scholes Model, Numerical Solution, Approximation, Call Option, Error Estimation

1. Introduction

In the historical backdrop of option pricing model, the Black-Scholes or Black-Scholes-Merton model [1] [2] is a standout amongst the most generous model. This model showed that the significance that mathematics plays an important role in the field of finance.

The Black-Scholes model was first published by Fischer Black and Myron Scholes in their 1973 seminal paper [1] , “The Pricing of Options and Corporate Liabilities”, published in the Journal of Political Economy. In the same year, they derived a partial differential equation, now called the Black-Scholes equation, which estimates the price of the option over time. Robert C. Merton was the first to publish a paper escalating the mathematical understanding of the options pricing model, and created the term “Black-Scholes options pricing model”.

Let us consider S as the price of the stock, which we consider as a random variable. V(S, t) be the value of an option as a function of time and stock price, K be the strike price, r be the risk-free interest rate, $\sigma $ be the Volatility/the standard deviation of the stock return, and t be the time in years. Then the famous Black-Scholes equation that was developed by Fisher Black and Myron Scholes is

$\frac{\partial V}{\partial t}+rS\frac{\partial V}{\partial S}+\frac{1}{2}{\sigma}^{2}{S}^{2}\frac{{\partial}^{2}V}{\partial {S}^{2}}-rV=0$ (1)

The above equation is a second-order parabolic partial differential equation known as Black-Scholes equation is actually a variation of a famous equation in physics that models the transfer of heat.

This equation can be solved by both analytically and numerically. Black and Scholes (1973) of M.I.T. first found the solution by taking advantage of previous research on option pricing that gave an idea of what the solution would look like. In [3] Mellin transformation was utilized to explain this model. They did not require variable change or explaining dispersion condition. R. Company, A.L. Gonzalez, L. Jodar [4] solved the Black-Scoles model which was modified with discrete dividend. They utilized a delta-characterizing grouping of generalized Dirac-Delta function and connected the Mellin transformation to acquire an integral formula. At last, the solution was approximated they approximated by utilizing numerical quadrature estimation.

In this article, a finite difference scheme is used to approximate the solution. This article is organized in the following way. Section 2 consists of the discussion of the dynamics of the model problem, and then the formulation of an explicit difference method to approximate the model and establishment of the stability of the scheme by convex combination in Section 3. In Section 4, we progress a programming code in MATLAB for the implementation of the numerical results for explicit finite difference scheme for Black-Scholes model. We compare the numerical solution with the analytical solution. We estimate the relative error of the numerical scheme in L_{1}-norm that leads to show the rate of convergence graphically. We also compare the explicit method with the result obtained by another work using semi-implicit method [5] .

2. The Black-Scholes Model

The Black-Scholes model is a powerful tool for valuation of equity options. In the early 1970’s, Myron Scholes, Robert Merton, and Fisher Black made an important breakthrough in the pricing of complex financial instruments by developing what has become known as the Black-Scholes model. This model displayed the importance that mathematics plays in the field of finance.

The Black-Scholes model is used to calculate the theoretical price of European put and call options, ignoring any dividends paid during the option’s lifetime. While the original Black-Scholes model did not take into consideration the effects of dividends paid during the life of the option, the model can be adapted to account for dividends by determining the ex-dividend date value of the underlying stock.

This section consists of the basic assumptions in the Black-Scholes model [6] . The hugest is that volatility, a measure of how much a stock can be relied upon to move in the close term, is a constant over time. The Black-Scholes model additionally assumes stock move in a way referred to as an arbitrary stroll; at any given time, they are as liable to climb as they are to move down. These suppositions are blended with the rule that option pricing ought to give no instant gain up to either merchant or purchaser. The correct six suppositions of the Black-Scholes model are as per the following:

1) No dividends are paid, 2) option must be exercised upon expiration, 3) market movement course cannot be anticipated, 4) no commission is charged in the exchange, 5) interest rate stays steady and 6) stock returns are ordinarily distributed, hence volatility is consistent.

The equation of Black-Scholes is

$\frac{\partial V}{\partial t}+rS\frac{\partial V}{\partial S}+\frac{1}{2}{\sigma}^{2}{S}^{2}\frac{{\partial}^{2}V}{\partial {S}^{2}}-rV=0$ (2)

It is a 2^{nd} order partial differential equation (PDE) in S-space and 1^{st} order in time. We generally use: at present t = 0, at expiry t = T.

For European call option [5] , let us denote the function value V by C. Now we define the initial and boundary conditions. In particular, we must consider final and boundary condition in the case of European Call option. The value at final time t = T can be calculated from the definition of call option. If at final time stock price is greater than strike price $\left(S>K\right)$ the call option will be worth $\left(S-K\right)$ because the buyer can buy the stock for K and sell it instantly for S and gain profit. But if $S<K$ then the buyer will not exercise the option and it will be worthless. So at t = T the option value is known is called the final condition and it is expressed as

$C\left(S,T\right)=\mathrm{max}\left(S-K,0\right)$ (3)

For boundary conditions, we consider the value of C at S = 0 and as $S\to \infty $ . If S = 0 then from the Equation (1) the payoff must also be 0. So the consequent boundary condition when S = 0:

$C\left(0,t\right)=0$ (4)

when $S\to \infty $ it is more likely that the option will be exercised and the value will be $S-K$ . As $S\to \infty $ the exercise price has not any impact on the option value, so the option value is equivalent to

$C\left(S,t\right)=S-K{\text{e}}^{-r\left(T-t\right)}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{as}\text{\hspace{0.17em}}S\to \infty $ (5)

This is the right boundary condition. Finally, the Black-Scholes initial [final] boundary value problem for European call option is

$\begin{array}{l}\frac{\partial V}{\partial t}+rS\frac{\partial V}{\partial S}+\frac{1}{2}{\sigma}^{2}{S}^{2}\frac{{\partial}^{2}V}{\partial {S}^{2}}-rV=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}0\le S<\infty ;0\le t\le T\\ \text{with}\text{\hspace{0.17em}}\text{\hspace{0.17em}}V\left(S,T\right)=\mathrm{max}\left(S-k,0\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}0\le S<\infty \\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}V\left(0,t\right)=0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}0\le t\le T\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}V\left(S,t\right)=S-K{\text{e}}^{-r\left(T-t\right)}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{as}\text{\hspace{0.17em}}S\to \infty \end{array}$ (6)

3. Objective and Methodology

The main objective of this article is to study the Black-Scholes partial differential equation and find the solution numerically for a European call option so that the dynamics of this model can be understood. A finite difference scheme will be introduced to approximate the numerical solutions. The boundary conditions will be used for the European call option. The stability of this numerical scheme is also studied. The study will be conducted theoretically. For this purpose, some computer programs such as MATLAB, Mathematica was used.

4. Numerical Solution

In most cases, it is very difficult to obtain the analytical solution of partial differential equation, even we obtain an analytical solution but in a very complex form. Therefore one needs to solve the PDE numerically. Finite difference method [7] [8] is one of the popular methods that have been used to solve partial differential equations. In this section, a finite difference scheme is developed in order to obtain to solve the Black-Scholes model numerically.

4.1. Explicit Difference Scheme

Time and Space discretization:

The interval $\left[0,T\right]$ is divided into N equally length subintervals having length ∆t. The underlying asset/stock will take prices in the interval $\left[0,\infty \right]$ . However, an artificial limit ${S}_{\mathrm{max}}$ is introduced which is normally three or four times bigger than the strike price K. This interval $\left[0,{S}_{\mathrm{max}}\right]$ is also divided into M equally sized subintervals having length ∆S. Therefore the space $\left[0,T\right]\times \left[0,{S}_{\mathrm{max}}\right]$ is approximated by a grid

$\left(n\Delta t,i\Delta S\right)\in \left[0,N\Delta t\right]\times \left[0,M\Delta S\right]$ (7)

where n = 0,.., N and i = 0,...,M . Whenever it is written C^{n} it’s actually
mean the numerical value of C(nΔt,iΔS )

Now to formulate an explicit difference scheme, the Black-Scholes partial differential equation

$\frac{\partial C}{\partial t}+rS\frac{\partial C}{\partial S}+\frac{1}{2}{\sigma}^{2}{S}^{2}\frac{{\partial}^{2}C}{\partial {S}^{2}}-rC=0$ (8)

is discretized using the following formulae:

The discretization of the time derivative $\frac{\partial C}{\partial t}$ is obtained by backward difference formula as below:

[Note: As many Partial Differential Equations have an initial condition, the new time step solution is approximate by discretizing the time derivative by forward difference formula. But in this case, the value of the option at the expiration time or at the final time is known which can be explained by the definition of call option. That is why backward difference formula is used for time derivative to calculate the previous time step solution, in this case, the option value at present-time.]

$\frac{\partial C}{\partial t}\left(n\Delta t,i\Delta S\right)\approx \frac{{C}_{i}^{n}-{C}_{i}^{n-1}}{\Delta t}$ (9)

Now the discretization of the first order derivative of S is obtained by central difference approximation and the second order derivative of S is obtained by a symmetric central difference approximation.

$\frac{\partial C}{\partial S}\left(n\Delta t,i\Delta S\right)\approx \frac{{C}_{i+1}^{n}-{C}_{i-1}^{n}}{2\Delta S}$ (10)

$\frac{{\partial}^{2}C}{{\partial}^{2}S}\left(n\Delta t,i\Delta S\right)\approx \frac{{C}_{i+1}^{n}-2{C}_{i}^{n}+{C}_{i-1}^{n}}{{\left(\Delta S\right)}^{2}}$ (11)

Now using Equations (9)-(11) in Equation (8) we get,

$\frac{{C}_{i}^{n}-{C}_{i}^{n-1}}{\Delta t}+r{S}_{i}\frac{{C}_{i+1}^{n}-{C}_{i-1}^{n}}{2\Delta S}+\frac{1}{2}{\sigma}^{2}{S}_{i}{}^{2}\frac{{C}_{i+1}^{n}-2{C}_{i}^{n}+{C}_{i-1}^{n}}{{\left(\Delta S\right)}^{2}}=r{C}_{i}^{n}$

Therefore

$\begin{array}{c}{C}_{i}^{n-1}={C}_{i}^{n}\left(1-r\Delta t-{\sigma}^{2}\frac{\Delta t}{{\left(\Delta S\right)}^{2}}{S}_{i}^{2}\right)+{C}_{i+1}^{n}\left(r\frac{\Delta t}{2\Delta S}{S}_{i}+\frac{1}{2}{\sigma}^{2}\frac{\Delta t}{{\left(\Delta S\right)}^{2}}{S}_{i}^{2}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{C}_{i-1}^{n}\left(\frac{1}{2}{\sigma}^{2}\frac{\Delta t}{{\left(\Delta S\right)}^{2}}{S}_{i}^{2}-r\frac{\Delta t}{2\Delta S}{S}_{i}\right)\end{array}$ (12)

This is the explicit difference scheme for Black-Scholes equation.

4.2. Stability Condition

Here the establishment of the stability condition for numerical scheme of Black-Scholes equation is discussed. For stability condition, it is considered that the scheme is without the source term i.e. considering $rV=0$ . Therefore,

$\begin{array}{c}{C}_{i}^{n-1}={C}_{i}^{n}\left(1-{\sigma}^{2}\frac{\Delta t}{{\left(\Delta S\right)}^{2}}{S}_{i}^{2}\right)+{C}_{i+1}^{n}\left(r\frac{\Delta t}{2\Delta S}{S}_{i}+\frac{1}{2}{\sigma}^{2}\frac{\Delta t}{{\left(\Delta S\right)}^{2}}{S}_{i}^{2}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+{C}_{i-1}^{n}\left(\frac{1}{2}{\sigma}^{2}\frac{\Delta t}{{\left(\Delta S\right)}^{2}}{S}_{i}^{2}-r\frac{\Delta t}{2\Delta S}{S}_{i}\right)\end{array}$ (13)

Now Let ${S}_{\mathrm{max}}=\mathrm{max}Si$

${\alpha}_{1}={\sigma}^{2}\frac{\Delta t}{{\left(\Delta S\right)}^{2}}{\left({S}_{\mathrm{max}}\right)}^{2},$

${\alpha}_{2}=r\frac{\Delta t}{2\Delta S}{S}_{\mathrm{max}}$

Therefore

${C}_{i}^{n-1}={C}_{i}^{n}\left(1-{\alpha}_{1}\right)+{C}_{i+1}^{n}\left({\alpha}_{2}+\frac{{\alpha}_{1}}{2}\right)+{C}_{i-1}^{n}\left(\frac{{\alpha}_{1}}{2}-{\alpha}_{2}\right)$ (15)

For convex combination, the sum of coefficient is equal to 1. And also,

$\left(1-{\alpha}_{1}\right)\ge 0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left({\alpha}_{2}+\frac{{\alpha}_{1}}{2}\right)\ge 0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(\frac{{\alpha}_{1}}{2}-{\alpha}_{2}\right)\ge 0$

where, ${S}_{\mathrm{max}}=\mathrm{max}Si$ .

5. Numerical Results

Consider for purposes of graphical presentation the value of a call option with strike price K = 100. The risk-free interest rate r = 0.12, the time to expiration is T = 1 measured in years, and the volatility is = 0.10. The value of the call option is plotted over a range of stock prices $70\le S\le 130$ surrounding the strike price is illustrated in Figure 1.

Figure 1 is produced by MATLAB. The figure shows the price of the stock at different time.

Figure 1. Numerical solution for European call at different time steps with K = 100, r = 0.12, σ = 0.1, T = 1.

Figure 2 is the analytic solution of the Black-Scoles model for the same parameters as before is.

If the above two figures are combined into one, one can easily examine how the analytical and numerical solutions differ with each other. For identification, two different colors were used. Here the values of T, r, K, σ are same as previous and the temporal grid size ∆t = 0.01, and spatial grid size ∆S = 1.50.

From Figure 3, it is observed that the two solutions are almost superimposed. The red color represents the numerical solution and green represents the analytic solution. Two solutions cannot be distinguished by the naked eye. Hence the numerical scheme has a very small error.

Now to find the relative error for the scheme, the computation of the relative error in L_{1}-norm which is defined by

Figure 2. Analytic solution for European call at different time steps with K = 100, r = 0.12, σ = 0.1, T = 1.

Figure 3. Analytic solution and Numerical solution at initial time (K = 100, r = 0.12, σ = 0.1, T = 1).

$err=\frac{{\Vert {C}_{e}-{C}_{n}\Vert}_{1}}{{\Vert {C}_{e}\Vert}_{1}}$

where ${C}_{e}$ is the analytic solution and ${C}_{n}$ is the numerical solution of Black-Scholes model for European call option. We perform our numerical scheme for t = 0 to 1, r = 0.12, K = 100, σ = 0.10 with temporal grid size ∆t = 0.0500 and spatial grid size ∆S = 3 which satisfy the stability conditions.

The relative error is shown in Figure 4.

From Figure 4 it is seen that as time t goes from 0 to 1, the relative error becomes smaller and smaller. For explicit difference scheme, the relative error is less than 0.0035 which is quite acceptable. Figure 5 represents that the error is decreasing with respect to the smaller discretization parameters ∆t and ∆S, which shows the convergence of the explicit difference scheme.

Figure 4. Relative error for explicit difference scheme in the order of 10^{−3}.

Figure 5. Relative error of explicit difference scheme for different temporal and spatial grid size in the order of 10^{−3}.

Table 1. Comparison between the semi-implicit method and explicit method.

Table 2. Approximate value by explicit centered method and the exact value.

6. Comparison with Other Numerical Scheme

In this section, a comparison is made between our explicit difference scheme with another numerical scheme that was used in finding the approximate solution of Black-Scholes model which is Semi-implicit method [5] .

For a European call option with $0\le S\le 20$ , T = 0.25, K = 10, r = 0.1, σ = 0.4, with temporal grid size N = 2000 and spatial grid size M = 200, the Semi-implicit method and Explicit method were used to set the table below.

From Table 1 it is seen that explicit scheme gives better results than semi-implicit method. But the results obtained by the scheme that was developed here are not much close to the exact value. But if the temporal grid points are increased up to N = 41,000, spatial grid points up to M = 1000 and set the other parameters value as above then the following values (Table 2) were found for different stock prices.

7. Conclusion

This article concerned with the numerical solution of Black-Scholes model. The model is studied and derived an explicit finite difference scheme and established a stability condition of the scheme. The relative error of the explicit scheme was estimated by comparing the numerical solution with the analytical solution in L1-norms and presents the convergence features of the scheme graphically. The numerical simulation results are seen in good agreement with the well-known qualitative behavior of the Black-Scholes PDE. Also, a comparison is presented between the explicit method and the result obtained by another work using the semi-implicit method where our result is found more accurate than that of the semi-implicit method.

Cite this paper

Anwar, M.N. and Andallah, L.S. (2018) A Study on Numerical Solution of Black-Scholes Model. Journal of Mathematical Finance, 8, 372-381. https://doi.org/10.4236/jmf.2018.82024

References

- 1. Black, F. and Scholes, M. (1973) The Pricing of Options and Corporate Liabilities. Journal of Political Economy, 81, 637-654. https://doi.org/10.1086/260062
- 2. Merton, R.C. (1973) Theory of Rational Option Pricing. The Bell Journal of Economics and Management Science, 4, 141-183. https://doi.org/10.2307/3003143
- 3. Jodar, L., Sevilla-Peris, P., Cortes, J.C. and Sala, R. (2005) A New Direct Method for Solving the Black-Scholes Equation. Applied Mathematics Letters, 18, 29-32. https://doi.org/10.1016/j.aml.2002.12.016
- 4. Company, R., Gonzalez, A.L. and Jodar, L. (2006) Numerical Solution of Modified Black-Scholes Equation Pricing Stock Options with Discrete Dividend. Mathematical and Computer Modeling, 44, 1058-1068. https://doi.org/10.1016/j.mcm.2006.03.009
- 5. Dura, G. and Mosneagu, A.-M. (2010) Numerical Approximation of Black-Scholes Equation. An. Stiint. Univ.”Al. I. Cuza” Iasi. Mat, 56, 39-64.
- 6. Shinde, A.S. and Takale, K.C. (2012) Study of Black-Scholes Model and its Applications. Procedia Engineering, 38, 270-279. https://doi.org/10.1016/j.proeng.2012.06.035
- 7. Smith, G.D. (1985) Numerical Solution of Partial Differential Equations: Finite Difference Methods. Clarendon Press, Oxford.
- 8. Duffy, D.J. (2013) Finite Difference Methods in Financial Engineering: A Partial Differential Equation Approach. John Wiley & Sons, Hoboken.