Multigrid Method for a Two Dimensional Fully Nonlinear Black-Scholes Equation with a Nonlinear Volatility Function

Aicha Driouchand Hassan Al Moatassime

1LAMFA UMR 7352 CNRS,University of Picardie Jules Verne,Amiens,France;

2Department of Mathematics,Faculty of science and technology,University Cadi Ayyad,Marrakesh,Morocco.

Abstract.This paper deals with the task of pricing European basket options in the presence of transaction costs.We develop a model that incorporates the illiquidity of the market into the classical two-assets Black-Scholes framework.We perform a numerical simulation using finite difference method.We consider a nonlinear multigrid method in order to reduce computational costs.The objective of this paper is to investigate a deterministic extension for the Barles’and Soner’s model and to demonstrate the effectiveness of multigrid approach to solving a fully nonlinear two dimensional Black-Scholes problem.

Key words:Fully nonlinear equation,multigrid method,black-scholes equation,finite difference method,FAS algorithm.

1 Introduction

Basket option is a financial derivative where the underlying asset is in the form of a sum or an average of different assets gathered together in a basket.Pricing basket options for nonlinear Black-Scholes equations is an interesting problem in financial mathematics,since an analytical solution of the basket option does not exist even in the linear case.Therefore,finding numerical solutions becomes an essential approach for solving the governing pricing equations(see[2]).In this work we are concerned with pricing European options in incomplete markets,more precisely,markets that are subjected to transaction costs.In this case the Black-Scholes theory given in[7]fails.To overcome this problem,several researchers have developped some extensions to the Black-Scholes equation by incorporating more realistic assumptions on financial markets(see[10,11,13,20]for more details).We focus our attention on a nonlinear deterministic Black-Scholes extension given by Barles and Soner in[6].They proposed an enlarged volatility to capture the influence of transaction costs,

The unknownVis the European call option,Sis the underlying asset,Tis the expiration date of the optionV,σis the constant volatility given in the linear Black-Scholes model,r≥0 is the risk-free interest rate constant,ais a nonnegative parameter that summarizes risk aversion and transaction costs and Ψ denotes the solution to the nonlinear ordinary differential equation,

Figure 1 is a representation of the solution of ODE(1.3)using a fourth order Runge-Kutta method.The problem(1.1)has been studied theoretically in[6]and the existence of a unique continuous viscosity solution has been shown using the theory of stochastic optimal control[14].Recently,a constructive theoretical approach for the problem(1.1)was given in[1].For numerical simulations of option pricing problems,Monte-Carlo methods are widely used.They are easy to implement and can handle multiple dimensions but they suffer from slow convergence(see[15,16]).Over the last two decades,numerical methods based on PDEs have also gained popularity since they allow an easy computation of Greeks which represent a basic instrument for risk management.However,PDE methods cannot compete with Monte-Carlo methods in two or more dimensions in terms of computing time,consequently,they are limited to one or two dimensions.The aim of our work is to push forward the limit of what can be done using PDE methods.There is a considerable literature that reviews numerical simulations of the option price when markets are subject to transaction costs using PDEs,we quote[4,19,22],where several nonlinear Black-Scholes models in the presence of transaction costs have been studied for European and American options,[30]presented a positivity preserving scheme for the Barles’and Soner’s model and they studied its convergence.[21]used a grid stretching technique on non-uniform meshes for the computation and[24,25]reviewed nonlinear parabolic equations arising from markets under transaction costs.However,there are only few works for the computation of nonlinear Black-Scholes equations depending on more than one underlying asset,among them we cite[29],they gave a numerical analysis for two dimensional nonlinear Black-Scholes equation for Spread option in incomplete markets and[3]where the authors gave a theoretical and numerical study of multi-asset nonlinear Black-Scholes equation arising from models with general transaction costs.Another problem is the boundary conditions for Black-Scholes equations,they have a boundary condition at infinity,which proves problematic for numerical simulations.One can overcome this difficulty by choosing a far field boundary condition similar to conditions provided in[18].The main scope of our work is to give an extended version of the Barles’and Soner’s model for solving numerically the two dimensional deterministic pricing problem,the difficulty resides in the nonlinear term,since it depends on different parameters such as the volatility of the two assets,the correlation coefficient and second derivatives of the option price.

Figure 1:Solution Ψ to(1.3).

This paper is organized as follows.In Section 2 we present the two dimensional Barles’and Soner’s model for a basket option,we give in Section 3 the discrete scheme and boundary conditions for our model.We review the multigrid full approximation storage technique for nonlinear problems in Section 4 then numerical results are given in the section thereafter.Finally a conclusion is in order.

2 2D nonlinear deterministic Black-Scholes model

The purpose of our study is to analyze the impact of transaction costs on the option price when dealing with two underlying assets.We assume that the stock valuesS1(t)andS2(t)are following a geometric Brownian motion governed by

whereσiandμi(i=1,2)are nonnegative constants that represent the volatility and the drift for theithasset respectively andBi(t)are correlated Brownian motions.Our aim is to price the basket option,the payoff in this case is given by

whereKis the strike price forS1andS2anda1+a2=1.The literature on hedging portfolios with correlated risky assets is very scanty,in[5]the value of the option is obtained by solving optimization problems and using asymptotic analysis.In our work we use a completely deterministic approach,we give a generalization of the Barles’and Soner’s model for pricing a European basket call,

For(S1,S2)∈R2+,t∈(0,T),whereρi,j=ρifi/=jandρij=1 otherwise,ρis the correlation coefficient between the two assetsS1andS2and Ψ is the solution of(1.3).Barles and Soner stated that Ψ has an asymptotic behaviour similar to the identity for largex.This propriety allows us to simulate Eq.(2.3)with this simplification.

Figure 2:Graph of the payoff of basket option(K=100,a1=a2=0.5).

3 Finite difference discretization

3.1 Boundary condition

It is known that the Black-Scholes equation is defined on an infinite domain[0,+∞)×[0,+∞).Therefore,in order to solve it numerically it is necessary to introduce artificial boundary conditions to make the computational domain bounded.Through rigorous error estimates,[18]shows that the optimal domain size for the far field boundary condition is given by

wherei={1,2}.We denote our computational domain by Ω =[0,S1max]×[0,S2max].A very common type of boundary condition in finance is the linearity condition[12];such a condition on a boundary?Ω reads

Eq.(3.2)means that the option price is nearly linear with respect to the asset price at a neighborhood of the infinity.However,a direct numerical simulation of Black-Scholes equation supplemented with this type of boundary conditions may cause some inaccuracies.Nevertheless,in[27,28]the authors show that Black-Scholes problems with linearity conditions can be treated with sufficient accuracy using a theoretical stability result.In our work we use both linearity conditions and Dirichlet boundary conditions,i.e.ifS1andS2are both equal to 0 then the option priceV(0,0,t)=0 for allt∈[0,T].If only one asset vanishes,i.eSi=0(i=1,2)then Eq.(2.3)withSi=0 leads us to the one dimensional Barles and Soner equation for the assetSj,j={1,2}i,for example,ifS1=0 then we solve

with the Dirichlet boundary condition atS2max,

IfSi=Simaxfor(i=1,2)then we use the linearity condition,this leads us to the following equation

forj∈{1,2}i.Finally,for the far field boundariesS1=S1maxandS2=S2maxwe use the linearity condition,

3.2 The numerical scheme

In the following,we adopt the two dimensional Crank-Nicolson finite difference scheme for Eq.(2.3)with nine point stencil.We use a variable transformation in time in order to have a forward parabolic equation,then Eq.(2.3)becomes

For the sake of simplicity we setx=S1andy=S2.Equation(3.6)becomes

Remark 3.1.One can notice that since we deal with nonsmooth initial data some inaccuracies may occur when approximating the derivatives using finite difference method.The accuracy can be improved using the grid stretching technique presented in[23].For the sake of simplicity of exposure,we do not use this technique and we restrict ourselves to uniform grids.

4 The multigrid method

In this section,we introduce the concepts of the multigrid method established in[8].The idea of this method is based on the fact that,low frequencies of components of the error on a fine grid become high frequencies on a coarser grid,therefore,we can accelerate the convergence of the solution of a fine grid problem by computing corrections on a coarser grid and then interpolating them back to the fine grid problem.In the following,we consider uniform two dimensional grids.

There are different multigrid approaches.The most known are the V-cycle,W-cycle,and full multigrid cycle[26].In our numerical simulations,we use the V-cycle method,which consists on recursively solving a system of equations on a sequence of grids of different levels from finest to coarsest.Each iteration transfers information through a series of restrictions(from a fine to a coarser grid)and prolongations(from a coarse to afiner grid).

Figure 3:Multigrid sketch on uniform 2D meshes.

4.1 Prolongation and restriction in two dimensions

4.1.1 Prolongation

We use a bilinear interpolation in order to have the solution values on the fine grid.We denote byi,jthe position indices in the coarse grid.The bilinear interpolation is given vertically by

whereucandufare the solutions on the coarse grid and the fine grid respectively.

4.1.2 Restriction

We use a combination of restriction using only the identity in some points of the grid and the full weighting method[9].For example,we can choose our restriction operator vertically as the following

Without loss of generality we use the identity operator when restricting points with dirichlet condition

4.2 The full approximation storage algorithm

An effective approach of the multigrid method to nonlinear problems is the FAS algorithm which is well established as a fast solver for nonlinear PDEs.In fact FAS algorithm is a generalization of the coarse-grid correction principle to nonlinear problems[17].Let(Ωk)1 ≤k≤nbe a sequence ofngrids Ωkof the domain Ω where Ωnis the finest grid and Ω1is the coarsest grid.Suppose that on every grid we have an iterative methodGkfor solving the problemAk(Uk)=fk,whereAkis the nonlinear operator.We denotethe final result afterνkiterations.The FAS algorithm is described as follows

It is worth noting that FAS reduces to the multigrid algorithm whenAkis a linear operator[8,17].However,FAS requires storing the values of the solutions over the entire grids.

4.2.1 Application to the Barles’and Soner’s model

5 Numerical results

5.1 The multigrid cycle

We begin our simulations with a demonstration of the effectiveness of the multigrid method.All the numerical experiments are performed on uniform grids.For the computations,we used the following parametersS1max=300,S2max=300,T=1(one year),a=0.02,K=100,σ1=0.2,σ2=0.2,ρ=0.5,r=0.1.We used 3 levels of grids withν=2 for presmoothing andν=3 for coarse grid corrections.

Figure 4 represents a rate of convergence comparison for both the multigrid and monogrid(standard Gauss-Seidel method on fine grids)in the linear case,Figure 5 and 6 represent the rate of convergence comparison for multigrid and monogrid when treating Ψ as the identity and in the case of Barles’and Soner’s model respectively.We have compared the total CPU time elapsed using the multigrid approach on uniform meshes and monogrid.All the numerical experiments were implemented in Mac OS X,2.7GHz Intel core i5 processor.

In Table 1 we compute the total CPU time elapsed for the two methods(monogrid and multigrid).The multigrid approach becomes significantly faster than monogrid(s-tandard Gauss-Seidel method)for higher mesh resolutions.In the case of meshes having a discretisation sizes 5122and 10242,we found that our implementation of the multigrid solver(Barles’and Soner’s equation)becomes 40 and 59 faster than a standard Gauss-Seidel method respectively.

Table 1:Comparison of the elapsed CPU time for the multigrid method for different discretization steps in the linear and nonlinear case.

Figure 4:A comparison of rate of convergence for multigrid and monogrid methods in the linear case on different mesh sizes(1282,2562and 5122)for 3 levels of grids.

Figure 5:A comparison of rate of convergence for multigrid and monogrid when Ψ(x)=x on different mesh sizes(1282,2562and 5122)for 3 levels of grids.

5.2 Back to the Barles’and Soner’s 2D model:the impact of transaction costs

Figure 6:A comparison of rate of convergence for the multigrid and monogrid in the case of Barles’and Soner’s model on different mesh sizes(1282,2562and 5122)for 3 levels of grids.

Figure 7:Comparison of the difference between the linear and the nonlinear case for differents transaction costs parameters a.

We end our numerical experiments by demonstrating the effect of the transaction cost parameteraon the option price. In Figure 7 we illustrate the influence of transaction costs on the numerical solution,by plotting the difference between the solution obtained with the Barles’and Soner’s model and the solution of the standard Black-Scholes equation.We see that,the higherais,the greater the option price becomes.Note that increasing the parameterahas a considerable impact on the stability of the Crank-Nicolson scheme;for the space step ΔSi=2.3i={1,2}anda=0.02,the time step Δt=0.001 was sufficient,contrary to the casea=0.4 where it has been necessary to consider a smaller time step(Δt=0.0001)to overcome the instability of the scheme and to get a stable numerical solution.

In what follows,we analyze different aspects of the proposed transaction cost model(two dimensional model of Barles and Soner).For numerical tests we used the following parameters.

Table 2:The chosen parameters for the numerical tests.

We denote byVTC,the option price using the Barles’and Soner’s model,VLinearthe option price using the linear Black-Scholes model and Γ=ΔVTC,the Gamma of the option that is a financial measure used by traders for hedging tasks.The Gamma represents the convexity of the option price with respect to the underlying asset,the higher it is,the greater the convexity and the price of the option varies rapidly in value.

Figures 8,9 and 10 illustrate the impact of transaction costs and the option price att=0.It should be noted that transaction costs are influenced by several factors,including the volatility of each underlying and the second derivatives of the assets.In Figure 8(Test 1),we observe that the maximum of these transaction costs is reached at the neighborhood of 2K,this is justified by the behaviour of Gamma(The ’Gamma’is maximum for values close to 2K).Moreover,the volatility effect on the transaction costs is observed in Figures 8 and 9.The ’Gamma’has an opposite behaviour with transaction costs,as it tends to decrease when the volatility increases.Finally,in Figure 10(Test3),we notice a symmetry in(10(a)-10(c))due to the fact that the two assets have the same volatility,again,the maximum of the transaction costs is reached when the two assets are close to 2Kand these costs are canceled whenS1andS2are far from 2K.

6 Conclusion

The major contribution of this work focuses on two aspects,a deterministic extension of the Barles’and Soner’s model for two assets and a numerical investigation of the effect of the Multigrid method.The obtained results confirmed the relevance of the extension proposed.In fact,we have shown numerically that these costs depend primarily on asset prices,the magnitude of the second derivatives(i.e.the ‘Gamma’)and the volatility of each asset.We also examined the effectiveness of the multigrid method when dealing with a two dimensional problem,we showed that it gives better results in terms of the execution time in contrast to the standard iterative method on a fine grid.

Figure 8:Test 1:Asset S1is more volatile than S2.

Figure 9:Test 2:Assets S1and S2are low volatiles.

Figure 10:Test 3:Assets S1and S2having the same volatility.

Acknowledgments

This work was supported by the ARC-Mathematics research federation FR3399.We would like to thank Olivier Goubet,Professor at university of Picardie Jules Verne and Mostafa Abounouh,Professor at Cadi Ayyad university for their insights and expertise that substantially improved the manuscript.We also would like to thank the anonymous reviewers for their precious advices and comments.