Haiming Song,Ran Zhang,?and WenYi Tian2Department of Mathematics,Jilin University,Changchun,3002,P.R.China.

2Department of Mathematics,Hong Kong Baptist University,Hong Kong.

Received 15 January 2014;Accepted 23 February 2014

Spectral Method for the Black-Scholes Model of American Options Valuation

Haiming Song1,Ran Zhang1,?and WenYi Tian21Department of Mathematics,Jilin University,Changchun,130012,P.R.China.

2Department of Mathematics,Hong Kong Baptist University,Hong Kong.

Received 15 January 2014;Accepted 23 February 2014

.In this paper,we devote ourselves to the research of numerical methods for American option pricing problems under the Black-Scholes model.The optimal exercise boundary which satisfies a nonlinear Volterra integral equation is resolved by a high-order collocation method based on graded meshes.For the other spatialdomain boundary,an artificial boundary condition is applied to the pricing problem for the effective truncation of the semi-infinite domain.Then,the front-fixing and stretching transformations are employed to change the truncated problem in an irregular domain into a one-dimensional parabolic problem in[-1,1].The Chebyshev spectral method coupled with fourth-order Runge-Kutta method is proposed for the resulting parabolic problem related to the options.The stability of the semi-discrete numerical method is established for the parabolic problem transformed from the original model.Numerical experiments are conducted to verify the performance of the proposed methods and compare them with some existing methods.

AMS subject classifications:35A35,90A09,65K10,65M12,65M60

American option pricing,Black-Scholes model,optimal exercise boundary,frontfixing,Chebyshev spectral method,Runge-Kutta method

1 Introduction

Options as a kind of important financial derivatives have a wide range of applications.The option pricing problem,especially the American option pricing problem,attracts the attention of more and more financial practitioners.The distinctive feature of an American option is its early exercise privilege,that is,the holder is endowed with the additional right to exercise the option at any time prior to the date of expiration,which makes an American option worth more than its European counterpart.Because of its early exercise privilege,American option pricing problem exists an optimal exercise boundary[3,26],which makes American pricing problem a nonlinear problem,and no analytical solution exists comparing to that of European options.

The researches of American option pricing problems have been extensively developed in recent decades.Traditionally,there are two ways-the analytical method and the numerical method.For analytical results,we refer to[1,9,10,13,25,33].They managed to present the solutions in a closed form depending on the optimal exercise boundary,then the option value was determined as long as the optimal exercise boundary was given.However it was not known actually in practice.For numerical aspects,it can be divided into two categories in general,one is based on Monte Carlo approach[11,28,30],the other is based on partial differential equation(PDE)approach[2,3,19,20,22].We would like to further adopt the PDE approach in this paper,since Monte Carlo method requires demanding computational resource due to its slow convergence.

The Black-Scholes equation is one of the most effective PDE models[5,21].And numerical methods are of popular use and frequently resorted for the Black-Scholes model among financial practitioners.For example,lattice tree methods,finite difference methods and finite element methods have been developed and extensively studied in recent decades.Cox,Ross and Rubinstein firstly introduced the binomial model to price American options in their seminalpaper[15].Later,Amin and Khanna showed the convergence of the binomial method in[4].Finite difference methods have been developed and discussed for a long time for American option pricing problem[8,19].One may refer to[22]for the convergence analysis.Recently,finite element methods[3,20]have attracted more interest in this field for its solid theoretical framework,in particularly its efficiency and variability.Interested readers may refer to[26]and the references therein for a complete survey.In this paper,we discuss the Chebyshev spectral method to the same problem,which turns to be efficient and comparable to finite element method.

There exist four main challenges for the numerical treatment of the American option pricing problem:

·The optimalexercise boundary is unknown,which satisfies a highly nonlinearequation,so it is not easy to get a fine resolution of optimal exercise boundary.

·For the other boundary,since we cannot adopt numerical methods directly to unbounded domain,how to truncate the unbounded domain and control the truncated errors are the key issues.

·The solving domain is an irregular domain,and the partition has to be reconstructed for each time order in order to resolve the option values exactly.

·The last one is how to choose an efficient numerical method to solve the problem fast and accurately.

About the first challenge,Cox[15,23]has proved that the optimal exercise boundary satisfied a nonlinear Volterra integral equation,and Ma et al.[29]have solved the same nonlinear Volterra integral equation by a high-order collocation method and presented numerical results.We shall follow the idea of[29]to deal with the first difficulty.For the second challenge,Holmes and Yang have introduced an artificial boundary condition to truncate the unbounded problem related to American call options in[20],we will use this method to deal with the put options.For the third challenge,front-fixing transformation technique[19,20,26,27,34],a very effective method for free boundary problems,is applied to convert the irregular domain into a regular one.Then the truncated problem can be transformed into a one-dimensional parabolic problem in[-1,1]by stretching transformation.

The numerical method for solving the one-dimensional parabolic problem in[-1,1]is our main task.To our knowledge,there are no pervious efforts about the spectral method applied to the American option pricing problems.Here,the Chebyshev spectral method[16,32]with fourth-order Runge-Kutta method[31]is adopted to solve the parabolic problem related to option values.Numerical simulations verify the effectiveness of our algorithms and appear to be more accurate and efficient than the finite element method.Then the problem can be solved quickly by the spectral method achieving the same accuracy.

The rest of the paper is organized as follows.In Section 2,we shall describe the Black-Scholes model for American put option,and present a high-order collocation method to solve the optimal exercise boundary.Artificial boundary condition,front-fixing and stretching transformations are employed to formulate a one-dimensional parabolic problem in[-1,1]to replace the original problem in the unbounded domain In section 3.In section 4,the Chebyshev spectral method is applied to the one-dimensional parabolic problem in[-1,1],and the stability of the method is also presented in this section.In Section 5,numerical simulations are implemented to test the performance of the proposed method and compare with some existing methods.Some concluded remark is collected in Section 6.

2 The optimal exercise boundary

In this section,we mainly introduce the optimal exercise boundary of American options.A derivation of this part can be found in[35].For completeness,we outline the main process as follows.Assume thatSandtare the underlying asset price and time respectively,then the put option valueP=P(S,t)satisfies the following free boundary problem

wherer,σ,d,KandTare the interest rate,volatility,dividend rate,exercise price and maturity date.B(t)is the optimal exercise boundary and F+=max(F,0).In[15],Cox,Ross,and Rubinstein proved that the optimal exercise boundary satisfied the following nonlinear Volterra integral equation[23,29]:

In order to simplify the equation,lett=T-τ,y(τ)=ln eB(τ)=lnB(T-τ),and combine with 1-N(x)=N(-x),we obtain

Here,λ is a constant that satisfies the following transcendental equation

Setting Ym,j=Yτ(τm+cj(τm+1-τm))(j=1,...,J+1),we can express Ymτ(the restriction of Yτon intervalkm)by interpolation

with Lagrange interpolation polynomials

Therefore,the global collocation solution Yτon[0,T]is given by

where Φm(τ)is the characteristic function onkm.

Substituting Yτinto(2.3),we can derive the following block nonlinear systems

for eachm=0,1,...,M-1,where Fm,jare given by

Here,we adopt simplified Newton’s method to solve the nonlinear systems(2.7),where the Jacobian is calculated as follows:

where

3 Pricing models via transformation and truncation

With the optimal exercise boundary computed in Section 2,American put option pricing problem(2.1)becomes a parabolic problem on a known unbounded domain.In this section,artificial boundary condition and front-fixing transformation will be discussed for the Black-Scholes equation to reduce the unbounded domain to the desired rectangular bounded computational domain.

3.1 Artificial boundary condition

In this subsection,we shall present an artificial boundary condition[14,24],which can truncate the infinite domain problem into a compact one and ensure the accuracy of calculation.

Observing equation(2.1),we can find that it is a backward equation with variational coefficients.Using the standard variable transforms[26]

whereαandβare constants to be determined,the Black-Scholes equation(2.1)with Dirichlet boundary conditions can be rewritten as a forward diffusion problem with constant coefficients in some variable domain

where the constants,the right hand function and the free boundary after transformation read as follows

In the following,we collect some important facts about the optimal exercise boundary before and after the change of variables.

Lemma 3.1(Property of the optimal exercise boundary,see[23,26]).The free boundary B(t)is a nondecreasing and bounded function with

where BP∞=KX is the price of Permanent American put options,and

Let B=ln(BP∞/K)=ln(X),then we get B≤b(τ)≤0,τ∈[0,T].

Next,we introduce the far field estimate,which plays a pivotal role to show the effectiveness of numerically solving the option pricing Eq.(2.1)on some truncated domains.The first result is due to Holmes and Yang[20],who have shown the estimate for American call options.For other numerical methods to truncate the unbounded domain we refer to[3]and[19].

Lemma 3.2(see[20]).For a given positive number ε∈(0,1),we have

and due to Lemma 3.1 and 3.2,we can easily see that

Given some sufficiently small thresholdε,we can derive corresponding upper boundYforS(or equivalentlya(τ)fory)from(3.5).In such a way,the bounded variable domain approximating problem(3.2)can be reformulated by

Note that two parallel free boundaries are introduced in(3.6).

3.2 Front-fixing transformation

In this subsection,we will firstly introduce the front-fixing method,which can transform the variable domain problem(3.6)into a rectangular bounded domain problem.And then,by some standard transformations,we can get a heat equation with homogeneous boundary conditions in[-1,1].

By using the front-fixing transforms[34]

we can reformulate variable domain problem(3.6)into the following problem in a rectangular domain:

Next,using the standard variable transforms

problem(3.8)can be rewritten as the following equations with the spatial variable in[-1,1]:

We can see that problem(3.10)is a heat equation with homogeneous boundary conditions in[-1,1],which can be solved with high accuracy by numerical methods.

4 Spectral methods

In this section,numerical methods will be presented for pricing American options.In spatial direction,the Chebyshev spectral method[32]will be used to approximate equation(3.10),and a fourth-order Runge-Kutta method[31]will be presented in temporal direction.In additional,the stability of the semi-discrete spectral numerical scheme is established.

4.1 Notations

Before designing the spectral method for equation(3.10),we introduce some definitions,and recall some basic results to be used in the sequel.

Firstofall,we review the Chebyshev polynomials and their specialproperties.Tn(z)=cos(ncos-1(z)),n=0,1,...,N,are the Chebyshev polynomials,which satisfy the following three-term recurrence relation

and the special property related to their derivatives

Settingω(z)=(1-z2)-1/2as the Chebyshev weightfunction,we introduce the weighted inner product are weighted norm,respectively

And the weighted Sobolev spaces are defined by

4.2 Chebyshev spectral method

In this subsection,we mainly describe the Chebyshev spectral method for solving the equation(3.10)and prove the stability of this method.

Firstly,we approximate the unknown functionu(z,τ)by Chebyshev polynomials

where the bilinear form is defined as

Theorem 4.1.Assume b∈H1([0,T]),then the Chebyshev spectral method(4.4)is stable,i.e.

where C is a constant depending on T.

Proof.LetV=Uin equation(4.4),we have

In[12],the following estimates

holds and there exists a constantδ＞0 such that

Then we have the upper bound for the third term on the left hand side of(4.5)by(4.7)and Cauchy inequality withε(Page 706 in[17]),as follows

Substituting(4.6)and(4.8)into(4.5),it obtains

From the assumption on functionb,it indicates thatC(τ)is integrable over(0,T].Then using Gronwall inequality,we see that the above inequality is simplified as

which obtains the conclusion.

4.3 Implementation procedure

In this subsection,we mainly describe the implementation procedure of Chebyshev spectral method[31,32]for solving equation(3.10).

Letznbe the Chebyshev-Gauss-Lobatto pointszn=cos(πn/N),0≤n≤N,andτm=mT/M,0≤m≤Mbe the temporal partition grids.Notice that{amn=an(τm)}can be expressed explicitly by{Um(zn)=U(xn,τm)}from(4.3),that is

wherec0,cN=2 andcn=1,for 0＜n＜N.

The first order derivative of the unknown functionu(z,τm)can be approximated by

Using the property(4.2)of Chebyshev polynomials,we obtain

By comparing the coefficients of(4.11)and(4.12),we obtain

In a similar way,we get the second-order derivatives in the following form

then we have

Therefore,replacing the unknown function and its derivatives by their corresponding approximations(4.3),(4.11)and(4.13)in the equation(3.10)produces a ODE systems aboutan(τ).And the fourth-order Runge-Kutta(RK4)method[31]is utilized with equation(4.9)to solve the resulting ODEs problem as follows

whereUm=(Um(z1),...,Um(zN))and

5 Numerical experiments

In this section,we present numerical simulations to verify the theoretic analysis in Section 2 and Section 4,and check the efficiency of the collocation method and spectral method for American option pricing problems.

Here,we consider the model of one-year(T=1)American put options,setσ=0.2,K=10 in equation(2.1),and vary parametersranddfor following three cases

·Case I:r＜dwithr=0.005 andd=0.01;

·Case II:r=dwithr=d=0.01;

·Case III:r＞dwithr=0.05 andd=0.01.

Before verifying the efficiency of our constructive method,we firstly declare thatε=10-6is chosen in Theorem 3.2,then take the solution obtained by binomialmethod[23]with 40000 points in temporal direction as the numerical approximation of the exact solution.

In the following examples,we show the efficiency of the collocation method for solving the optimal exercise boundary.

Example 5.1.The Eq.(2.2)which satisfied by the optimal exercise boundary with above three cases is considered.We takeJ=3 and Tol=10-6in Eq.(2.8),and letM=25 in temporal direction.Then the approximation of optimal exercise boundary can be obtained by

Figure 1:The optimal exercise boundary computed by collocation method with M=25 and the binomial method with 40000 points for three cases.Case I:r＜d(left);Case II:r=d(middle);Case III:r＞d(right).

Table 1:L2 error and time cost about optimal exercise boundary computed by collocation method(CM)with M=25 and binomial method(BM)with M=10000 for three cases.

The results of the collocation method for these three cases are shown in Fig.1 with notation ”+”,and the solid lines represent the optimal exercise boundary computed by the binomial method with 40000 points in temporal direction.From Fig.1,we can see that the optimal exercise boundary computed by the collocation method approximates the exact solution well.We also give a comparison between the collocation method and binomial method in Table 5.1,from which,we can find that the collocation method is much faster than binomial method in the same accuracy.The numerical results confirm the general fact that binomial method needs more computational cost to guarantee the accuracy.

With the optimal exercise boundary obtained by the collocation method,the free boundary problem(2.1)is converted into a regular parabolic problem(3.10),which will be solved by Chebyshev spectral method.

Example 5.2.We consider to solve Eq.(3.10)by using spectral method for the three casesr＜d,r=dandr＞d.In spatial direction approximations,N=43,N=33 andN=15 are taken as the highest degrees of Chebyshev polynomials,respectively.And in temporal direction,we chooseM=10000,and use RK4 to solve(3.10).Then we can obtain the approximations of the option values by

Figure 2:The option values computed by the Chebyshev spectral method.Case I:r＜d(left);Case II:r=d(middle);Case III:r＞d(right).

whereα=α0,β=γα2+r.

Fig.2 lists the results of the spectral method for these three cases.From the numerical simulations,we can see that the option values computed by the spectral method is a good approximation of the exact solution.The comparison between the spectral method and finite element method is also made in Table 5.2,from which,we can notice that the spectral method is faster than finite element method with the same accuracy,especially the case ofr＞d.The results about option values in Table 5.2 also show that the bigger ratior/dis,the less time consuming and the more accurate option values are simulated by spectral method.This observation confirms with the truth that the spectral method works more effectively for smoother functions,in fact,withr/dbeing bigger,the initial value becomes smoother.In general,we can conclude that spectral method is an effective method for American option pricing problems.

6 Conclusions

In this paper,we introduce a collocation method based on the graded meshes for solving the nonlinear Volterra integral equation,which is derived for the optimal exercise boundary arising in pricing American option.With this boundary,we reformulated the Black-Scholes equation into a one-dimensional parabolic problem in[-1,1]by the artif icial boundary condition truncation,front-fixing and the stretching transformations.The resulting problem for the valuation of American option is solved by spectral method coupled with Runge-Kutta method.We also give the stability analysis of spectral method,and numerical simulations are provided to verify the efficiency of these numerical algorithms.

Table 2:L2 error,L∞error and time cost about option values computed by spectral method(SM)and finite element method(FEM)at t=0 for three cases.

Acknowledgments

The authors would like to thank the anonymous referees for their careful reading of the manuscript and their valuable comments.The authors also wish to thank the High Performance Computing Center of Jilin University and Computing Center of Jilin Province for essential support.This work was supported by the National Natural Science Foundation of China(Grant Nos.11271157,11371171),the Open Project Program of the State Key Lab of CAD&CG(A1302)of Zhejiang University and the Scientific Research Foundation for Returned Scholars,Ministry of Education of China.

[1]G.Barone-Adesi and R.E.Whaley,Efficient analytic approximation of American option values,J.Finance,42(1987),301-320.

[2]W.Allegretto,G.Barone-Adesi,and R.J.Elliott,Numerical evaluation of the critical price and American options,European J.Finance,1(1995),69-78.

[3]W.Allegretto,Y.Lin and H.Yang,Finite element error estimates for a nonlocal problem in American option valuation,SIAM J.Numer.Anal.,39(2001),834-857.

[4]K.Amin and A.Khanna,Convergence of American option values from discrete-to continuoustime financial models,Math.Finance,4(1994),289-304.

[5]F.Black and M.Scholes,The pricing of options and corporate liabilities,J.Pol.Econ.,81(1973),637-659.

[6]H.Brunner,Collocation Methods for Volterra Integral and Related Functional Equations,Cambridge University Press,2004.

[7]H.Brunner,The numerical solution of weakly singular Volterra integral equations by collocation on graded meshes,Math.Comp.,45(1985),417-437.

[8]M.Brennan and E.Schwartz,Finite difference methods and jump processes arising in the pricing of contingent claims:A synthesis,J.Financ.Quant.Anal.,13(1978),461-474.

[9]L.Badea and J.Wang,A new formulation for the valuation of American options,I:Solution uniqueness,Anal.Sci.Comput.,2000,3-16.

[10]L.Badea and J.Wang,A new formulation for the valuation of American options,I:Solution existence,Anal.Sci.Comput.,2000,17-33.

[11]P.Boyle,M.Broadie and P.Glasserman,Monte Carlo methods for security pricing,J.Econom.Dynam.Control,21(1997),1267-1321.

[12]N.Bressan and A.Quarteroni,Analysis of chebyshev collocation methods for parabolic equations,SIAM J.Numer.Anal.,23(1986),1138-1154.

[13]P.Carr,R.Jarrow and R.Myneni,Alternative characterizations of American put options,Math.Finance,2(1992),87-106.

[14]J.R.Cannon,The One-Dimensional Heat Equation,Cambridge University Press&Reading,MA,Addison Wesley,2004.

[15]J.C.Cox,S.A.Ross and M.Rubinstein,Option pricing:A simplified approach,J.Financ.Econ.,7(1979),229-263.

[16]C.Canuto,M.Y.Hussaini,A.Quarteroni and T.A.Zang,Spectral Methods,Fundamentals in Single Domains,Springer Press,Berlin,2006.

[17]L.C.Evans,Partial differential equations,American Mathematical Society,Providence,RI,second edition,2010.

[18]J.D.Evans,R.Kuske and J.B.Keller,American options on assets with dividends near expiry,Math.Finance,12(2002),219-237.

[19]H.Han and X.Wu,A fast numerical method for the Black-Scholes equation of American options,SIAM J.Numer.Anal.,41(2003),2081-2095.

[20]A.D.Holmes and H.Yang,A front-fixing finite element method for the valuation of American options,SIAM J.Sci.Comput.,30(2008),2158-2180.

[21]J.Hull,Fundamentals of Futures and Options Markets,Prentice Hall,2007.

[22]P.Jaillet,D.Lamberton and B.Lapeyre,Variational inequalities and the pricing of American options,Acta Appl.Math.,21(1990),263-289.

[23]L.Jiang,Mathematical Modeling and Methods of Option Pricing,World Scientific Press,Singapore,2005.

[24]P.Kangro and R.Nicolaides,Far field boundary conditions for Black-Scholes equations,SIAM J.Numer.Anal.,38(2000),1357-1368.

[25]I.J.Kim,The analytic valuation of American puts,Rev.Financial Studies,3(1990),547-572.

[26]Y.K.Kwok,Mathematical Models of Financial Derivatives,2nd edition,Springer Finance,Berlin Heidelberg,2008.

[27]H.G.Landau,Heat conduction in a melting solid,Quart.J.Appl.Math.,8(1950),81-94.

[28]F.A.Longstaff and E.S.Schwartz,Valuing American options by simulation:A simple leastsquares approach,Rev.Financial Studies,14(2001),113-147.

[29]J.Ma,K.Xiang and Y.Jiang,An Integral Equation Method with High-Order Collocation Implementations for Pricing American Put Options,Int.J.Econ.Finance,2(2010),102-112.

[30]L.C.G.Rogers,Monte Carlo valuation of American options,Math.Finance,12(2002),271-286.

[31]J.Shen and T.Tang,Spectral and High-Order Methods with Applications,Science Press,Beijing,2006.

[32]J.Shen,T.Tang and L.L.Wang,Spectral Methods,Algorithms,Analysis and Applications,Springer Press,Berlin,2011.

[33]W.K.Wong and K.Xu,Refining the quadratic approximation formula for an American option,Int.J.Theor.Appl.Finance,4(2001),773-781.

[34]L.X.Wu and Y.K.Kwok,A front-fixing finite difference method for the valuation of American options,J.Financ.Eng.,6(1997),83-97.

[35]R.Zhang,H.M.Song and N.N.Luan,Weak Galerkin finite element method for valuation of American options,Front.Math.China,9(2014),455-476.

?Corresponding author.Email addresses:songhm11@mails.jlu.edu.cn(H.Song),zhangran@jlu.edu.cn(R.Zhang),twymath@gmail.com(W.Tian)