Pricing of average strike Asian call option using numerical PDE methods
aa r X i v : . [ q -f i n . C P ] J un Pricing of average strike Asian call option using numerical PDE methodsAbhishek Kumar , Ashwin Waikos and Siddhartha P. Chakrabarty Department of Mathematics, Indian Institute of Technology Guwahati, Guwahati 781039, Assam, India
Abstract
In this paper, a standard PDE for the pricing of arithmetic average strike Asian call option is pre-sented. A Crank-Nicolson Implicit Method and a Higher Order Compact finite difference scheme forthis pricing problem is derived. Both these schemes were implemented for various values of risk freerate and volatility. The option prices for the same set of values of risk free rate and volatility was alsocomputed using Monte Carlo simulation. The comparative results of the two numerical PDE methodsshows close match with the Monte Carlo results, with the Higher Order Compact scheme exhibiting abetter match. To the best of our knowledge, this is the first work to use the numerical PDE approach forpricing Asian call options with average strike.
Key Words : Asian Option; Crank Nicolson Implicit Method; Higher Order Compact; Monte Carlo Simulation
Options are one of the most common financial derivatives that are traded both in exchanges as well as overthe counter [1, 2, 3, 4]. The two most common options are the European and the American options. Thepurchase or sale of the underlying asset for the option takes place at the discretion of the holder of the optionfor a price called the exercise price or the strike price at a fixed time called expiration time (denoted by T ) incase of European option and any time in [0 , T ] in case of American option. Options can be further classifiedas call or put depending whether the holder has the right to buy or the right to sell the underlying asset. Inaddition there are numerous other options that can broadly be classified as exotic options [1, 2, 5, 6]. Thespecialty of these kinds of options is that the final payoff is more sophisticated and sometimes depends onsome function of the path of prices of underlying asset. One of the common path dependent exotic optionsis the Asian option [1, 2, 5, 6, 7], whose payoff depends on the average historical stock prices. In this paper,we will focus on pricing of an Asian call option with arithmetic average strike.We begin with the assumption that the stock prices S ( t ) follow a geometric Brownian motion given bythe stochastic differential equation [1, 2, 5, 6, 7], dS ( t ) = rS ( t ) dt + σS ( t ) dW ( t ) , where r is the average rate of growth of asset prices or drift, σ is the volatility and W ( t )(0 ≤ t ≤ T ) is aBrownian motion or Wiener process under the risk neutral measure P . The payoff for an Asian call option [1, 2, 5, 6, 7], with arithmetic average strike is given by, V ( T ) = max S ( T ) − T T Z S ( u ) du, . he price of the option at time t ∈ [0 , T ] (with filtration F t ) is given by the risk-neutral pricing formula [7], V ( t ) = E (cid:16) e − r ( T − t ) V ( T ) |F t (cid:17) . In order to deal with the challenge of the payoff V ( T ) being path-dependent we introduce a stochasticprocess I ( t ) [2, 7] given by, I ( t ) = t Z S ( u ) du. The stochastic differential equation for evolution of I ( t ) is thus given by, dI ( t ) = S ( t ) dt. Thus the Asian call option price V ( S, I, t ) for continuous arithmetic average strike satisfies the backwardPDE [2, 6, 7], ∂V∂t + rS ∂V∂S + 12 σ S ∂ V∂S + S ∂V∂I − rV = 0 . Note that the above problem is three dimensional which leads to greater computational costs. This moti-vates the reduction of the problem into lower dimension. [2]. For this purpose, a new variable R ( t ) = S ( t ) R t S ( u ) du is defined [2, 6]. This in turn motivates the ansatz V ( S, I, t ) = S · H ( R, t ) for some func-tion H ( R, t ) . It can be shown that the SDE satisfied by R ( t ) is [6], dR ( t ) = (cid:0) σ − µ ) R ( t ) (cid:1) dt − σR ( t ) dW ( t ) . Also the function H ( R, t ) satisfies the PDE [6], ∂H∂t + 12 σ R ∂ H∂R + (1 − rR ) ∂H∂R = 0 . The solution of this backward PDE requires a final condition and two boundary conditions which are out-lined below [2, 6].i.
Final Condition :
The final payoff for the option gives the final condition, H ( R ( T ) , T ) = max (cid:18) − T R ( T ) , (cid:19) . ii. Right Hand Boundary Condition :
The right hand boundary condition for R → ∞ can be obtained byobserving that since the integral R ( t ) is bounded, so S → for R → ∞ . For S → the option is notexercised rendering it’s value to be . Hence, H ( R, t ) = 0 for R → ∞ . iii. Left Hand Boundary Condition :
The left hand boundary condition for R → can be obtained fromthe similarity reduction equation. The term R∂H/∂R → as R → . Assuming that H is bounded itfollows that the term R ∂ H/∂R → as R → . This leads to the boundary condition, ∂H∂t + ∂H∂R = 0 for R → . he problem to be solved now reduces to, ∂H∂t + 12 σ R ∂ H∂R + (1 − rR ) ∂H∂R = 0 . subject to, H ( R ( T ) , T ) = max (cid:18) − R ( T ) T , (cid:19) ∂H∂t + ∂H∂R = 0 for R → (1) H = 0 , for R → ∞ . Once the solution H ( R, t ) is obtained the price of the Asian option is determined by, V ( S (0) , R (0) ,
0) = S (0) H ( R (0) , , where S (0) is the initial stock price. For pricing of options (especially exotic options) the standard methodused is Monte Carlo simulation [8, 9, 10]. This involves simulating the paths of the underlying asset andcalculating the option price based on this path. A large number of such simulations are run and the averageof the option prices from each simulation is taken to be the option price. Several methodologies have beenadopted for pricing of options using the numerical PDE approach. Some of the most commonly used onesare finite differences of lower order [11, 12, 13, 14] and higher order compact schemes for American options[15, 16] and option pricing in stochastic volatility model [17]. The problem of pricing the average strike Asian call option essentially entails solving for equation (1). Whilegeometric mean Asian option admits closed form solutions [8], the same is not true in case of arithmeticaverage Asian options. As such one has to seek a solution through numerical methods for PDEs [12]. Thereare several articles in literature which dwell upon numerical PDE approach to Asian option pricing. One ofthe first papers to deal with numerical PDE pricing of options is by Rogers and Shi [11]. In this paper, theauthors first reduce the problem of solving a parabolic PDE in two variables and present a highly accuratelower bound. Zvan et al. [13] in their technical report, do an extensive analysis of numerical PDE methodsof Asian options. They discuss the shortcomings of applying the usual numerical PDE techniques used forstandard options in case of Asian options. In particular they adapt flux limiting techniques from computa-tional fluid dynamics (CFD) to tackle the problem of spurious oscillations that arise in Asian options. Vecer[12] provided a numerical implementation of the Asian option pricing problem using the θ method. Duboisand Lelievre [14] extend the approach by Rogers and Shi [11] and propose a scheme which produced fastand accurate results. While all these papers [11, 12, 13, 14] do examine the pricing problem from the nu-merical PDE point of view, the focus is mostly on fixed strike options. Rogers and Shi [11] and Zvan et al.[13] present some results on average strike put options.Our main objective in this paper will be to use Higher Order Compact (HOC) scheme for this purpose.We will postpone the discussion on this until the next section. In this section we will present the Crank-icolson Implicit Method (CNIM) for solving equation (1) and compare the results with those obtained byMonte Carlo simulation [18].CNIM is obtained by taking the average between Forward-Time Centered-Space method (FTCS) andBackward-Time Centered-Space method (BTCS). For this purpose, let us define a finite difference dis-cretization of the PDE (equation(1)) with the uniform grid t n = t + ( n − · ∆ t , n = 1 : N + 1 and R i = R + ( i − · ∆ R , i = 1 : M + 1 , where ∆ t and ∆ R are the temporal and spatial mesh size respec-tively. The values used in this paper are t = 0 and t N +1 = T = 1 (i.e year option) with R = 0 and R M +1 = 5 . Let us also define the variables c ( R ) = σ R and d ( R ) = (1 − rR ) . H ( R, t ) at the point ( R i , t n ) is denoted by H ni . The CNIM discretization of equation (1) is then given by, H n +1 i − H ni ∆ t + c i ∆ R " H n +1 i +1 + H ni +1 − (cid:18) H n +1 i + H ni (cid:19) + H n +1 i − + H ni − + d i R " H n +1 i +1 + H ni +1 − H n +1 i − + H ni − = 0 (2)The above can now be rewritten as, G i H ni +1 + K i H ni + J i H ni − = D i H n +1 i +1 + E i H n +1 i + F i H n +1 i − (3)where, G i = − c i R − d i R K i = c i ∆ R + 1∆ t J i = − c i R + d i R D i = c i R + d i R E i = − c i ∆ R + 1∆ t F i = c i R − d i R . Let us define the vector, H ( n ) = ( H n , H n , H n , H n , . . . , H nM ) ⊤ for n = 1 : N + 1 . The CNIM can now be written in the matrix form as, BH ( n ) = AH ( n +1) + b ( n ) , (4)where, A = E D . . . F E D . . . ... . . . . . . . . . ... . . . . . . . . . D M − . . . F M E M . nd B = K G . . . J K G . . . ... . . . . . . . . . ... . . . . . . . . . G M − . . . J M K M .b ( n ) = F H n +11 − J H n . . .. . .D M H n +1 M +1 − G M H nM +1 . From second order finite differences, the one-sided difference is given by, ∂H∂R (cid:12)(cid:12)(cid:12)(cid:12) i = − H i +2 + 4 H i +1 − H i R + O (∆ R ) . Therefore, applying the same on the left boundary along with backward time approximation we get, H n = (1 − k ) H n +11 + 4 kH n +12 − kH n +13 , where k = ∆ t R . The right boundary H nM +1 = 0 follows from equation (1). The final condition is given by, H N +1 i = (cid:18) − R i T (cid:19) + . The CNIM formulation above was solved using MatLab
SSTM . The solution was obtained using an iterativeprocess which involved a tolerance criterion of | H new (1 , − H old (1 , | < ǫ . For the purpose of thisimplementation the number of time and space grid points were taken to be and respectively. Thetolerance level was taken to be ǫ = 10 − . Finite difference methods have been used for solving ODE and PDE problems for quite a long time. Theyare relatively easier to set up and solve, but require structured mesh [19, 20]. Finite element methods onthe contrary are more sophisticated, works well with irregular domains and are amenable to unstructuredmeshes. Finite element methods are relatively more challenging to implement. Standard finite differenceschemes like the CNIM (used in the previous section) are second order accurate. However, in financialapplications like option pricing a higher level of accuracy is desirable. A direct extension of the centraldifference schemes to achieve higher order accuracy would involve more node points on the stencil. Aninnovative way of achieving higher accuracy with lesser number of nodes in the stencil came by the way ofhigher order compact (HOC) schemes. Spotz and Carey [19] and Spotz [20] provide an excellent discourseon application of this scheme in case of viscous flow and computational mechanics. → .
06 0 .
06 0 . . . . σ ↓ CNIM MC CNIM MC CNIM MC S (0) = 100 .Despite it’s enormous potential of application to finance problems, HOC schemes have not been usedmuch in this area. Zhao et al. [21] presented a compact scheme for American option pricing with secondorder accuracy in space. Tangman et al. [15, 16] applied a HOC scheme for the pricing of American putoption. They do a comparative [15] analysis with a non-compact fourth order scheme. In their subsequentpaper they describe an improvement of a method suggested by Han and Wu [22]. In this chapter we shallapply HOC scheme to the setup (1) which is written in the following form [18], c ( R ) ∂ H∂R + d ( R ) ∂H∂R = g (5)where c ( R ) and d ( R ) are as defined in the previous section and g = − ∂H∂t . We now define notation δ [19, 20, 23] as follows, δ R f := ∂f∂R , δ R f := ∂ f∂R and so on.Thus rewriting Equation (5) in terms of finite difference discretization we get, c i δ R H i + d i δ R H i = g i (6)In the HOC scheme we derive the leading truncation error terms in terms of finite difference equation makinguse of the original equation. Denoting these terms by τ , the HOC scheme is obtained by subtracting τ backto the original finite difference discretization. Thus we have, c i δ R H i + d i δ R H i − τ i = g i (7)where, τ i = ∆ R (cid:18) c ∂ H∂R + 2 d ∂ H∂R (cid:19) + O (∆ R ) (8)rom the initial PDE (1) we get, ∂c∂R ∂ H∂R + c ∂ H∂R + ∂d∂R ∂H∂R + d ∂ H∂R = ∂g∂R (9)Therefore this can be rewritten in δ notation as, ∂ H∂R (cid:12)(cid:12)(cid:12)(cid:12) i = − c i (cid:20) δ R c i δ R H i + δ R d i δ R H i + d i δ R H i (cid:21) + 1 c i δ R g i (10)Differentiating the PDE (9) w.r.t. R once more and simplifying we get, ∂ c∂R ∂ H∂R + 2 ∂c∂R ∂ H∂R + c ∂ H∂R + ∂ d∂R ∂H∂R + 2 ∂d∂R ∂ H∂R + d ∂ H∂R = ∂ g∂R (11)which can be written in δ -notation as, − c i ∂ H∂R (cid:12)(cid:12)(cid:12)(cid:12) i = − δ R g i + δ R c i δ R H i + δ R d i δ R H i + 2 δ R d i δ R H i + (cid:0) δ R g i − δ R c i δ R H i − δ R d i δ R H i − d i δ R H i (cid:1) δ R c i c i (12) + (cid:0) δ R g i − δ R c i δ R H i − δ R d i δ R H i − d i δ R H i (cid:1) d i c i Making use of approximations in equation (10) and equation (12) in the truncation error ( τ i ) (equation (8))and hence substituting τ i in equation (7) we obtain, δ R H i (cid:20) c i + ∆ R d i c i + σ R i ∆ R d i c i + ∆ R σ − r ∆ R − ∆ R c i σ R i − ∆ R c i d i σ R i (cid:21) + δ R H i (cid:20) d i − r ∆ R d i c i + r ∆ R σ R i c i (cid:21) (13) = (cid:20) R δ R + ∆ R d i c i δ R − ∆ R c i ( σ R i ) δ R (cid:21) g i Note that, c i = σ R i ⇒ δ R c i = σ R i , δ R c i = σ and d i = 1 − rR i ⇒ δ R d i = − r, δ R d i = 0 .We define, A i = (cid:20) c i + ∆ R d i c i + σ R i ∆ R d i c i + ∆ R σ − r ∆ R − ∆ R c i σ R i − ∆ R c i d i σ R i (cid:21) B i = (cid:20) d i − r ∆ R d i c i + r ∆ R σ R i c i (cid:21) F i = ∆ R d i c i − ∆ R c i ( σ R i ) , k t R , k t R . We now apply the HOC scheme to the above equation (recalling that g = − ∂H∂t ) and obtain, k (cid:0) H n +1 i +1 − H n +1 i + H n +1 i − + H ni +1 − H ni + H ni − (cid:1) A i + k (cid:0) H n +1 i +1 − H n +1 i − + H ni +1 − H ni − (cid:1) B i = (cid:0) − H n +1 i + H ni (cid:1) − ∆ R " H n +1 i +1 − H n +1 i + H n +1 i − ∆ R + ∆ R (cid:20) H ni +1 − H ni + H ni − ∆ R (cid:21) − F i (cid:2) H n +1 i +1 − H n +1 i − (cid:3) + F i (cid:2) H ni +1 − H ni − (cid:3) (14)he above can now be rewritten as, G i H ni +1 + K i H ni + J i H ni − = D i H n +1 i +1 + E i H n +1 i + F i H n +1 i − (15)where, G i = − k A i − k B i + 112 + F i K i = 2 k A i + 56 J i = − k A i + k B i + 112 − F i D i = k A i + k B i + 112 + F i E i = − k A i + 56 F i = k A i − k B i + 112 − F i . The implicit method can be written in the matrix form, BH ( n ) = AH ( n +1) + b ( n ) , where H ( n ) , A, B and b ( n ) has already been defined in the previous section. As with the case of CNIM thenumber of time and space grid points were taken to be and respectively along with the tolerancelevel of ǫ = 10 − . The scheme was implemented using MatLab SSTM . r → .
06 0 .
06 0 . . . . σ ↓ HOC MC HOC MC HOC MC . . . . . S (0) = 100 . In this section we discuss the results obtained by using the CNIM and the HOC schemes as outlined in theprevious two sections. As already noted we could not find any results for average strike Asian call optionusing numerical PDE methods. For the purpose of comparison we used the Monte Carlo (MC) simulationas the benchmark value. We generated the path of a stock prices using the geometric Brownian motionrocess. We generated such paths and determined the option price from each of the paths generated.The average of all these option prices was taken to be the option price, for the purpose of comparisonwith the PDE methods. We generated the option prices using all the three methods for three values of r = 0 . , . , . and five values of σ = 0 . , . , . , . , . .A comparative study of results from the CNIM and the MC methods showed a close match. The com-parative results are presented in Table 1 along with the CPU time in seconds. The option prices for the threevalues of r against the five values of σ have been presented in the graphical form in Figures (1), (2) and (3).For r = 0 . (Figure (1)), the match was very close except the case where σ = 0 . . This slight differencein the option price is reduced when r = 0 . (Figure (2)). The other values for r = 0 . showed a close match.The results were similar for the case r = 0 . except for a very minimal difference in the case of σ = 0 . (Figure (3)). The CPU time in case of CNIM was however considerably lower ( < . seconds) as comparedwith the Monte Carlo simulation ( > seconds).The results and comparison of the HOC scheme and the MC method indicates excellent consonance.A comparison of the results from these two methods in terms of values and CPU time in seconds havebeen presented in a tabular form in Table 2. The option prices from both the methods are very close toeach other. In fact the results obtained from the HOC scheme show a better match with the MC simulationresults as compared with the CNIM method. This holds for all the values of r and σ and is evident fromthe comparative figures (Figures (4), (5), (6)) of HOC and MC. As was the case with CNIM, the CPU timetaken in case of the HOC scheme is significantly less ( < . seconds) in contrast to Monte Carlo simulationwhich requires at least seconds. In this paper we examined several ways of computing the price of an average strike Asian call option,namely Monte Carlo simulation and the numerical PDE approach. In case of option pricing, the benchmarkgenerally used is Monte Carlo simulation which suffers from some severe drawbacks like computationalcosts and a certain amount of uncertainty of pricing. In contrast, the usage of numerical PDE approachesthat we have taken results in lesser computational costs and also provides an unique answer. The numericalPDE approach in pricing the average strike Asian call option is by and large an unexplored area, since thisapproach applied to Asian option is mostly concentrated on the case of fixed strike. In this paper, we takethe PDE approach to the pricing problem and present two schemes to accomplish this numerically. Firstly,we use the Crank-Nicolson Implicit Method (which is second order) to solve the PDE and hence price theoption. Then, we present a Higher Order Compact scheme (fourth order) to solve the problem. Finallywe make a comparison of results obtained from the PDE approach with that of Monte Carlo. The resultsobtained using the two PDE techniques were in excellent agreement with the Monte Carlo results. Theresults obtained using Higher Order Scheme are closer to the Monte Carlo results as opposed to Crank-Nicolson Implicit method vis-a-vis Monte Carlo. This is more so in case of lower values of σ . To the best ofour knowledge, this is the first work to use the numerical PDE approach for pricing Asian call options withaverage strike. We believe this work would find more applications in the area of option pricing through thePDE approach. Acknowledgments
The authors express their deep gratitude to Dr. Jiten C. Kalita (Department of Mathematics, IIT Guwa-hati) and Prof. Anoop K. Dass (Department of Mechanical Engineering, IIT Guwahati) for their valuableguidance and suggestions during the preparation of the manuscript.
References [1] J. C. Hull,
Options, Futures and Other Derivatives , Prentice Hall of India, 2006.[2] P. Wilmott, S. Howison & J. Dewynne,
The Mathematics of Financial Derivatives , Cambridge Univer-sity Press, 1995.[3] S. Roman,
Introduction to the Mathematics of Finance: From Risk Management to Options Pricing ,Springer, 2004.[4] M. Capinski & T. Zastawniak,
Mathematics for Finance: An Introduction to Financial Engineering ,Springer, 2003.[5] D. J. Higham,
An Introduction to Financial Option Valuation: Mathematics, Stochastics and Compu-tation
Cambridge University Press, 2004.[6] R. Seydel,
Tools for Computational Finance , Springer, 2006.[7] S. Shreve,
Stochastic Calculus for Finance, Volume II: Continuous-Time Models , Springer-Verlag,2004.[8] P. Glasserman,
Monte Carlo Methods in Financial Engineering
Springer, 2003.[9] P. P. Boyle, Options: A Monte Carlo approach,
Journal of Financial Economics , vol. 4(3), pp. 323–338, 1977.[10] M. Broadie and P. Glasserman, Estimating security price derivatives using simulation,
ManagementScience , vol. 42(2), pp. 269–285, 1996.[11] L. C. G. Rogers & Z. Shi, The Value of an Asian option,
Journal of Applied Probability , vol. 32(4),pp. 1077–1088, 1995.[12] J. Vecer, A new PDE approach for pricing arithmetic average Asian options,
Journal of ComputationalFinance , vol. 4(4), pp. 105–113, 2001.[13] R. Zvan, P. A. Forsyth & K. Vetzal, Robust Numerical Methods for PDE Models of Asian Options,Technical report, Cheriton School of Computer Science, University of Waterloo, 1996.[14] F. Dubois & T. Lelievr, Efficient Pricing of Asian Options by the PDE Approach,
Journal of Compu-tational Finance , Vol. 8(2), pp. 55–64, Winter 2004/05.15] D. Y. Tangman, A. Gopaul & M. Bhuruth, Numerical pricing of options using high-order compactfinite difference schemes,
Journal of Computational and Applied Mathematics , vol. 218(2), pp. 270–280, 2008.[16] D. Y. Tangman, A. Gopaul & M. Bhuruth, A fast high-order finite difference algorithm for pricingAmerican options,
Journal of Computational and Applied Mathematics , vol. 222(1), pp. 17–29, 2008.[17] B. During & M. Fournie, High-order compact finite difference scheme for option pricing in stochasticvolatility models,
SSRN-1646885 .[18] A. Kumar & A. Waikos,
Pricing of Asian options using higher order compact scheme , B.Tech project,Department of Mathematics, Indian Institute of Technology Guwahati, 2011.[19] W. F. Spotz & G. F. Carey, High-order compact finite difference methods with applications to viscousflows, Technical report, Texas Institute for Computational and Applied Mathematics, 1994.[20] W. F. Spotz,
High order compact finite difference schemes for computational mechanics , Ph.D thesis,University of Texas at Austin, 1995.[21] J. Zhao, M. Davison & R. M. Corless, Compact finite difference method for American option pricing,
Journal of Computational and Applied Mathematics , vol. 206(1), pp.306–321, 2007.[22] H. Han & X. Wu, A fast numerical method for the BlackScholes equation of American options,
SIAMJournal of Numerical Analysis , vol. 41(6), pp. 2081-2095, 2003.[23] J. C. Kalita, D. C. Dalal & A. K. Dass, A class of higher order compact schemes for the unsteady two-dimensional convectiondiffusion equation with variable convection coefficients,
International Journalfor Numerical Methods in Fluids , vol. 38(12), pp. 1111–1131, 2002. .05 0.1 0.15 0.2 0.25 0.3 0.35 0.434567891011 σ O p t i on P r i ce Comparison between Crank Nicolson and Monte Carlo for r=0.06
Crank NicolsonMonte Carlo
Figure 1: Comparison between Monte Carlo Simulation (MC) and Crank Nicolson Implicit Method (CNIM)for r=0.06 σ O p t i on P r i ce Comparison between Crank Nicolson and Monte Carlo for r=0.1
Crank NicolsonMonte Carlo
Figure 2: Comparison between Monte Carlo Simulation (MC) and Crank Nicolson Implicit Method (CNIM)for r=0.1 .05 0.1 0.15 0.2 0.25 0.3 0.35 0.49101112131415 σ O p t i on P r i ce Comparison between Crank Nicolson and Monte Carlo for r=0.2
Crank NicolsonMonte Carlo
Figure 3: Comparison between Monte Carlo Simulation (MC) and Crank Nicolson Implicit Method (CNIM)for r=0.2 σ O p t i on P r i ce Comparison between Higher Order Compact and Monte Carlo for r=0.06
Higher Order CompactMonte Carlo
Figure 4: Comparison between Monte Carlo Simulation (MC) and Higher Order Compact Scheme (HOC)for r=0.06 .05 0.1 0.15 0.2 0.25 0.3 0.35 0.4456789101112 σ O p t i on P r i ce Comparison between Higher Order Compact and Monte Carlo for r=0.1
Higher Order CompactMonte Carlo
Figure 5: Comparison between Monte Carlo Simulation (MC) and Higher Order Compact Scheme (HOC)for r=0.1 σ O p t i on P r i ce Comparison between Higher Order Compact and Monte Carlo for r=0.2