Numerical Methods for Discrete Double Barrier Option Pricing Based on Merton Jump Diffusion Model

As a kind of weak-path dependent options, barrier options are an important kind of exotic options. Because the pricing formula for pricing barrier options with discrete observations cannot avoid computing a high dimensional integral, numerical calculation is time-consuming. In the current studies, some scholars just obtained theoretical derivation, or gave some simulation calculations. Others impose underlying assets on some strong assumptions, for example, a lot of calculations are based on the Black-Scholes model. This thesis considers Merton jump diffusion model as the basic model to derive the pricing formula of discrete double barrier option; numerical calculation method is used to approximate the continuous convolution by calculating discrete convolution. Then we compare the results of theoretical calculation with simulation results by Monte Carlo method, to verify their efficiency and accuracy. By comparing the results of degeneration constant parameter model with the results of previous models we verified the calculation method is correct indirectly. Compared with the Monte Carlo simulation method, the numerical results are stable. Even if we assume the simulation results are accurate, the time consumed by the numerical method to achieve the same accuracy is much less than the Monte Carlo simulation method.


Introduction
Options as a kind of important financial derivatives have been developed rapidly in recent decades.With the increase of investment demands and the progress of theoretical level, all kinds of exotic options were born.As an important path dependent exotic option, barrier option was accounted for nearly half of the share.
Their research has a very important practical significance.So, new exotic option pricing problem has become an important research topic.
Black and Scholes [1] set up the B-S model and get the pricing formula of the call option.Common barrier options pricing has a long history.It was first studied by Merton [2] on the down and out call option pricing.After this period, Reimer and Sandmann [3] use the binary tree model research barriers options.
Then Gao, Huang and Subrahmanyam [4] pushed forward the pricing analysis of the American barrier option by using the decomposition technique, and Dai and Kwok [5] has got the pricing formula of American down and in call option.
The model of the underlying asset is also developed from the original B-S model to a variety of more complex models.Such as Wang, Du [6] obtained pricing formula in the class of barrier options under jump diffusion model, and Xie [7] put forward the constant elasticity of variance (CEV) pricing barrier option process.However, these developments are based on the constant parameter model.Farnoosh, Sobhani, Rezazadeh [8] proposed a method based on the time dependent parametric B-S model of the barrier options to extend the research of this topic to another direction.
Looking at the development of barrier option pricing, the barrier option pricing theory is from a study based on the continuous observation idealized model to the effective discrete observation model, and the stochastic model of the target stock is also developed from B-S model with constant parameters to model with time-dependent parameter, or constant parameter jump diffusion model, CEV model, SVJD model and so on.In order to make the scope of application more extensive, the goal of this paper is the further development for the numerical algorithm of the option pricing based on Merton jump diffusion model.
In the numerical calculation section, having a wide range of applications, Monte Carlo simulation method is certainly one of the alternative calculation methods.But it is often accompanied by a wide range of accuracy and low computational efficiency.As traditional numerical algorithm performs badly in nonsmooth solutions, Golbabai, Ballestraand Ahmadian used finite element method (FEM) to improve orders of convergence in [9].In [10], Ahmadian and Ballestra found it also performs well under a constant elasticity of variance model with jump diffusion.For some models, the other way to get numerical method is directly starting from the analytic solution.In [8], Farnoosh, Sobhani, Rezazadeh is proposed that for time dependent numerical parameters of barrier options based on B-S model an analytical solution can be obtained.They started from solving partial differential equations.By constructing a heat conduction equation, they finally got the analytical solution.However, it still can't avoid calculate a multiple integral of a high dimension.Using the numerical solution of the Romberg extrapolation algorithm to calculate is taking time and precision into account.But consider from the efficiency of the integral calculation, it is not as good as quadrature method which been used by Milev and Tagliani in general B-S model [11].The difference is probably cause by the influence of boundary.
The quadrature method used on general B-S model is able to extend to more complex model.Furthermore, for time dependent parameters, it also can achieve good results.
In the second section, we will consider knock-out call option as an example to introduce risk neutral pricing based on Merton jump diffusion model.Finally, we obtain an analytical formula of the discrete double barrier option.
In the third section, we will introduce the application of quadrature method to calculate the high dimensional integrals generation generated in the second section.We also introduce Monte Carlo simulation method to compute the result for comparison.
In the fourth section, the computational efficiency of the numerical method is compared with samples in literatures and results of simulation method.
In the conclusion, we note that this calculation method is significantly better in accuracy and efficiency, and there is a certain promotion potential.

Discrete Double Barrier Options Model
We assume the European discrete knock-out double barrier option has the following conditions: 1) European options, i.e., option can only be exercised at We denote In risk-neutral measure, if we as- sume risk-free interest rate r is constant.(Thus, the drift µ is replaced by the interest rate r.)By the definition of the option its price at the beginning time is equal to ( ) We can calculate this formula to get the analytical formula of the discrete double barrier option.
Theorem 2.1.The value of option with above assume is given by the following integral where φ µ σ is the p.d.f. of ( ) From Equation ( 4), the value of this barrier option is, ( ( ) where Step 2. Calculate ( ) .As we know while the length is j, sum of these variables where ( ) In particular, as there is a constant shift in 1 ξ , ( ) Plug ( ) h x into Equation ( 5), we finally obtain (4).□ Specially, if parameters is change to different constant in different monitor interval,

( )
i h x will be different while they can be calculate in the same method by different parameters.

Numerical Valuation Algorithm
From section 2 we know random variables , 1, 2, , i i m ξ =  decide the option value at t = T.As we have already known the distribution of i ξ , we can use Monte Carlo method to simulate i ξ and get some possible option value at t = T.
By weak law of large numbers, the means of the Discounting of these price is converging to the option price at t = 0.However, we can't expect it as an efficient selection.
Since we have got the analytical formula of the discrete double barrier option.
We can also calculate the numerical result by calculate the analytical formula.
While there is a High dimensional multiple integral, we should anticipate some tricks.As the integration has a break at the boundary of the integral area, we compare and find that Milev's method to discrete the integration, solve directly perform a better precision than the extrapolation method, such as the Romberg extrapolation algorithm.
Using Equation ( 19) and Equation (20) we can calculate all , , i k n q step by step.
Notice that ( ) . Plug it into Equation (13), we have =e max e , 0 . where Unfortunately, we can't calculate from above, we know  Secondly, we introduce our Monte Carlo simulation method.The main idea Monte Carlo simulation method is generating thousands of simulations of S t .
Then we can use the mean price of option base on these S t to estimate ( ) according to the law of large numbers.While the price of option just bases on the value of t S on monitor date, we just need to simulation these values on monitor date.From Equation (2) we can conveniently simulate ln t S then get t S indirectly.In other words, we simulate , 1, 2, , i i m ξ =  .As we know i ξ is sum of independent random variables including a normal random variable  and a series of normal random variables The length of the series is follow a Poisson distribution with intensity T m λ .We simulate η , and use a while loop to simulate k η , and exit the loop with probability e T m λ − each round.As usual, estimation error of Moto Carlo simulation is Proportional to 1 N while estimated standard deviation is V N σ , V σ is standard deviation of ( ) , V S T , N is the number of simulation path.Because we can know the option value if t S touch the barrier in early monitor date, so in simulation, we can also stop a simulation round in advance.However, the probability to stop is different for different parameters.So we use no-stop method to simulate and record the time in worst situation.

Numerical Results
We use some cases to compare the estimation accuracy of two methods introduced above.
Firstly, we use parameters that Poisson intensity 0 λ = in case 1, case 2 and case 3 to compare the result with other literatures.It shows when λ tends to 0, results tends to the results base on B-S model.Using the estimate value and estimate standard error calculate by Monte Carlo Simulation as the real value, the t statistic of estimate error of numerical algorithm with n = 1000 is around 1. It means its estimate error is at the same level with the estimate standard error of Monte Carlo Simulation, i.e., the accuracy of numerical algorithm with n = 1000 is like the accuracy of Monte Carlo Simulation with simulation times = 100000.However, CPU time cost by the latter is thousands times of the former.To estimate the numerical error on the option price, we will compare the numerical result with the accurate result with n = 10000 and calculate the error.
Case 1: Prices of discrete double knock-out call option in 5 monitoring dates.The current price of the underlying asset is 0 S , the strike price is 100, the vola- tility is 20% per annum, the call option has six months remaining to maturity, the risk-free rate is 10% per annum (compounded continuously), the lower barrier is placed at 95, and the upper barrier is imposed at 110.Numerical method estimate numerical error (written in brackets).Monte Carlo Simulation also estimate standard error (written in brackets), Milev's result is copy from his paper.Because they are generated by different computers, so CPU time in present methods and CPU time in Milev's paper cannot be compared.Table 1 shows the result., see Table 2.    3.
Secondly, we compare the efficiency of numerical method and simulation method in case 4 and case 5. Poisson intensity 0 λ > , i.e., jump is considered.
We set 5 λ = , 2 0.05, 0.05 , other parameters except U, L and 0 S are the same with case 2. We could learn that while U and L are quite close, option is Worthless.So standard deviation of option value at maturity is so large relative to the option value at now.Both of the two methods have a large relative error.
While U and L are not quite close, although calculate result's accuracy is a bit worse than Simulation result's, it is still much better in efficiency as it cost much less CPU time.
Case 4: U and L are close so that the option has a high probability to knock-out, see Table 4.  Case 5: U and L are not quite close so that the option has a lower probability to knock-out, see Table 5.
Finally, if we predict parameters are not constants, we could set different value of parameters in different monitor interval (interval between two adjacent monitor date).In case 6, parameters are the same with case 5, except σ and r .σ will increase from 0.1 to 0.25 in linearly step in monitor intervals.r will de- crease from 0.05 to 0.02 in linearly step in monitor intervals.In this case, as we should refresh ( ) Case 6: see Table 6.

Conclusion
While mainstream methods to price European discrete knock-out double barrier option are not better than Monte Carlo simulation method in estimated error order, the advantage of the presented numerical algorithm is that it costs less compute time than simulation method.It is benefit from its direct calculate thought.Because of the main idea of this method that is not complicated, it can be extended to the case stock price obeys Merton jump diffusion model and work well on it.However, as this numerical method is based on the analytical solution, it is applicable only to those stochastic models for which the transition probability can be computed in closed-form.

maturity; 2 ); 5 )
The time for observation has is set in advance, we assume observation time interval is consistent.For different observation intervals, the subsequent theory still holds; 3) All the observation time have the same U and L as barrier.When the stock price at the observation time higher than U or lower than L, the option knock-out; 4) Ignore the transaction costs and taxes.It means if the option hasn't knocked out at observation time, at the time of maturity, the option value is the stock price can change continuously.

For 2 ,
3, , i m =  from Equation (2), We know i ξ is sum of independent random variables including a normal random variable The length of the series is follow a Poisson distribution with intensity T m λ

Firstly, we
introduce the discrete method.Let ln U d n L =, where n is an integer number.We used discrete random variable , convergence, for 2, 3, , i m =  , by the definition of , k n p , we have( . If parameter is different in different monitor interval, we should just refresh ( ) j g x in each iteration to keep the algorithm run correctly.

Case 3 :
Prices of discrete double knock-out call option monitored weekly (25 times) for different values of the underlying asset 0

j g x
in each iteration, while n = 2000, CPU time mainly cost by generate , , k n n p ′ whose time cost is ( ) O mnn′ .From the result we could learn that while n = 2000, the accuracy of numerical algorithm is a bit worse than Monte Carlo simulation.But while it just cost 1.4% CPU time of simulation method, as CPU time cost is proportional to n which means calculate result's convergence rate is faster than simulation while n = 2000.It still has a much better efficiency.