An Accurate Numerical Integrator for the Solution of Black Scholes Financial Model Equation

In this paper the Black Scholes differential equation is transformed into a parabolic heat equation by appropriate change in variables. The transformed equation is semi-discretized by the Method of Lines (MOL). The evolving system of ordinary differential equations (ODEs) is integrated numerically by an L-stable trapezoidal-like integrator. Results show accuracy of relative maximum error of order 10−10.


Introduction
Financial derivative, in particular options, became very popular financial contracts in the last few decades.Options can be used to hedge assets and portfolios in order to control the risk due to movements in the share price.A European call (put) option provides the right share price to buy or sell a fixed number of assets at the fixed exercised price E, at the expiry time t 0 [1].In the early 1970's Fisher Black and Myron Scholes made a major breakthrough by deriving a partial differential equation that must be satisfied by the price of any derivative security dependent on a non-dividend-paying stock [2].According to [3] their work had a huge impact on how options were viewed in the financial world.In an idealized financial market, the price of a European option can be obtained as the solution of the Black-Scholes equation [4] [5].However, the Black Scholes equation has been derived under quite restrictive assumptions such as frictionless, liquid and complete market.In recent years nonlinear Black Scholes equations have been derived in order to model transaction costs arising in the hedging of portfolios [1] [6] [7] and feedback effects due to large traders [8]- [13].
In seeking the solution of the Black Scholes equation, emphasis is always laid on derivation of formula or equation for the price of the option of interest and computation of the price of the option.This calls for the usage of numerical methods because explicit theoretical solutions for the price of the option do not exist.From a binomial model, [6] derived an option price that takes into account transaction costs which approximates a Black Scholes price but with modified volatility.[14] [15] computed the option price of the Black Scholes equation as the solution of a nonlinear quasi-variational inequality.This approach has the disadvantage that the option price depends on the choice of the utility function.Seeking the analytical solution of the Black-Scholes equation, [16] used the Adomain decomposition method.Adomain decomposition can provide analytical approximations to a wide class of linear and nonlinear equations without perturbation, closure approximations or discretization.Solution for the Black-Scholes equation as a semigroup on spaces of continuous functions on (0, ∞) is presented in [17].
In the mathematical literature, only a few results can be found on the numerical discretization of Black Scholes equations.The numerical approaches vary from binomial approximations for American options in stochastic framework [18], Monte-Carlo methods [19], finite element discretization [20], and finite difference approximations [1].The numerical discretization of the Black Scholes equations with nonlinear volatilities has been performed using explicit finite difference schemes [21] [22].Explicit numerical schemes have the disadvantage that restrictive conditions on the discretization parameters, time and space steps, are needed to obtain stable, convergent schemes [23].Moreover the convergence order is the only one in time.
The Method of Lines (MOL) is a general procedure for the solution of time-dependent partial differential equations (PDEs) [24].The basic idea of the MOL is to replace the spatial (boundary value) derivatives in the PDEs with algebraic approximations [25].Ones this is done, the spatial derivatives are no longer stated explicitly in terms of the spatial dependent variables.Thus, only the initial value variable, typically time in a physical problem, remains, which results in a system of ODEs that approximates the original PDE.An accurate integration algorithm for initial value ODEs to compute an approximate numerical solution to the PDE can then be used for the numerical integration.One of the salient features of the MOL is the use of existing, and generally well established, numerical methods for ODEs [26].
This paper is organized thus: In Section 2 we transform the black Scholes equation into a heat equation by change in variables.In Section 3 we introduce an L-stable trapezoidal like integrator for the numerical integration of the transformed Black Scholes equation.Section 4 is devoted to the numerical test of the method on the transformed Black Scholes equation.Section 5 explains the computation of the errors and relative errors of the method while results are discussed in Section 6.

Transforming the Black Scholes Equation into a Parabolic Heat Equation
Given the Black Scholes equation: Subject to: ( ) where: Following [27] to transform the diffusion-advection-reaction Equation (1) into a parabolic heat PDE, we make the following change of variables: By taking appropriate partial derivatives of ( ) , V S t in Equation (4c) and substituting in Equation (1) yields Substituting e x K for S in Equation ( 5) and dividing by Defining And substituting for k in Equation ( 6) we obtain ( ) ( ) Making , x xx u u and u τ the subjects of Equations ( 10), ( 11) and ( 12) respectively we obtain: e ax b u w bu Substituting Equations ( 13), ( 14) and (15) into Equation ( 8) we obtain By letting the coefficients of w and x w in Equation ( 16) vanish identically Under the condition that and ( ) the Black Scholes equation given in Equation ( 1) is transformed into the parabolic heat equation PDE 2 , , 0, ,0 max e e ,0 , According to [28] European Call option as the solution to the Black Scholes equation on 0 ,0 , ~e as , Equating both hand sides of Equation ( 21) yields Substituting for S from Equation (4) in Equation ( 22) gives Substituting 2 2τ σ for ( ) , e e By equating the RHS of Equations ( 2b) and ( 24) , e e Substituting k for 2 2r σ from Equation (7) in Equation ( 25) we obtain ( ) Substituting Equation ( 26) for ( ) , u x τ in Equation ( 9) yields

Remark
By appropriate change in variables, Equation ( 1) is transformed into Equation ( 19) which is a parabolic heat equation to be discretized by the MOL.Equation ( 27) is the derived approximate theoretical solution to the transformed Black Scholes equation.

L-Stable Implicit Trapezoidal-Like Integrator
The trapezoidal-like integrator ) 1 e e e 1 e e e e e l l l l and λ λ are the first and second eigenvalues of the discretization matrix respectively; l is the time step; and ' denotes differentiation with respect to time.The derivation of the method (28) and the analysis of the order of accuracy are as discussed in [29], while the stability properties of the method are discussed in [30].

Computation of Absolute and Relative Errors
In this section we explain how the absolute errors and relative errors of the methods shown in Table 1 and Table 2 were obtained.

Figure 1 .
Figure 1.The graphs of the numerical solution of the new scheme and the theoretical solution.

Table 1 .
Solution approximations and errors of the new scheme.

Table 2 .
Solution approximations, errors and relative errors of the new scheme.