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  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  on the down and out call option pricing. After this period, Reimer and Sandmann  use the binary tree model research barriers options. Then Gao, Huang and Subrahmanyam  pushed forward the pricing analysis of the American barrier option by using the decomposition technique, and Dai and Kwok  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  obtained pricing formula in the class of barrier options under jump diffusion model, and Xie  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  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 non- smooth solutions, Golbabai, Ballestraand Ahmadian used finite element method (FEM) to improve orders of convergence in  . In  , 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  , 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  . 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.
2. 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 maturity; 2) The time for observation has is set in advance, in the subsequent derivation, 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 ; 5) the stock price can change continuously.
Merton jump diffusion model assumes that the stock price is in accordance with the following stochastic differential equation,
where , Wt is standard Wiener process, is Poisson
process with intensity λ. We denote , . By
Itô Lemma, we obtain
We denote . In risk-neutral measure, if we assume 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
is the p.d.f. of .
Proof. Step 1. Assuming that the p.d.f. of random variable is known, we calculate the analytical result.
We denote , we have . Let Math_25#, We get . For ,
From Equation (4), the value of this barrier option is,
is p.d.f. of .
Step 2. Calculate .
For from Equation (2), We know 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 . As we know while the length
is j, sum of these variables , we get the p.d.f. of ,
where is the p.d.f. of . Because has no relation with i, in , we denote with . In particular, as there is
a constant shift in , .
Plug into Equation (5), we finally obtain (4). □
Specially, if parameters is change to different constant in different monitor interval, will be different while they can be calculate in the same method by different parameters.
3. Numerical Valuation Algorithm
From section 2 we know random variables decide the option value at t = T. As we have already known the distribution of , we can use Monte Carlo method to simulate 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.
We use discrete method to calculate the High dimensional multiple integral by calculate an approximate discrete convolution. Monte Carlo method is also used for compared.
Firstly, we introduce the discrete method. Let , where n is an in-
teger number. We used discrete random variable to approximate , satisfy
. Specially, we discrete random variable to approximate , satisfy
From Equation (5), we know
So we can use to approximate , where
We denote , . From
Equation (12), we know
So we can calculate by calculate . From the definition of , we know
So the value of may be , and satisfy
where . Notice that . is equivalent to
. is equivalent to . In fact, we have the
Lemma 3.1. For , and
If , then .
Proof. From the definition of , it can be easily proved.□
For , by lemma 3.1, we have
So we know . As the possible value of is an integer
multiple of d, the value of may be . For , let
. Notice that ,
from Equation (18),
Using Equation (19) and Equation (20) we can calculate all step by step. Notice that , we just need to precompute for . Plug it into Equation (13), we have
It means we could calculate by Equation (19), Equation (20) and Equation (21) after we got and . Actually, as we don’t need to calculate by Equation (21), there is no need to calculate Equation (20). Equation (19) is a discrete convolution we have to calculate m times.
We use the p.d.f. of To get and . Denote
Because is uniform convergence, for , by the definition of , we have
where is c.d.f. of , and is c.d.f. of standard normal distribution. Similarly, from the p.d.f. of , we can get
Unfortunately, we can’t calculate and directly although we have got their analytic expression. Let
from above, we know , and can
be directly calculate. By using to approximate , we finally get an approximation of , and denote it with , who is a computable approximation of . If parameter is different in different monitor interval, we should just refresh in each iteration to keep the algorithm run correctly.
Secondly, we introduce our Monte Carlo simulation method. The main idea Monte Carlo simulation method is generating thousands of simulations of St. Then we can use the mean price of option base on these St to estimate according to the law of large numbers. While the price of option just bases on the value of on monitor date, we just need to simulation these values on monitor date. From Equation (2) we can conveniently simulate then get indirectly. In other words, we simulate . As we know 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 . We simulate , and use a while loop to simulate , and exit the loop with probability each round. As usual, estimation error of Moto Carlo simulation is Proportional to while estimated standard deviation is , is standard deviation of , is the number of simulation path. Because we can know the option value if 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.
4. Numerical Results
We use some cases to compare the estimation accuracy of two methods introduced above.
Firstly, we use parameters that Poisson intensity 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 , the strike price is 100, the volatility 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.
Case 2: Prices of discrete double knock-out call option monitored daily (125 times) for different values of the underlying asset and parameters , see Table 2.
Table 1. 5 monitoring dates, S0, K = 100, σ = 0.25, T = 0.5, r = 0.05, L = 95, U = 110.
Table 2. 125 monitoring dates, S0, K = 100, σ = 0.25, T = 0.5, r = 0.05, L = 95, U = 110.
Case 3: Prices of discrete double knock-out call option monitored weekly (25 times) for different values of the underlying asset and parameters , see Table 3.
Secondly, we compare the efficiency of numerical method and simulation method in case 4 and case 5. Poisson intensity , i.e., jump is considered. We set , , other parameters except U, L and 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.
Table 3. 25 monitoring dates, S0, K = 100, σ = 0.25, T = 0.5, r = 0.05, L = 95, U = 110.
Table 4. λ = 5, μJ = 0.05, = 0.05, L = 95, U = 105.
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 . will increase from 0.1 to 0.25 in linearly step in monitor intervals. will decrease from 0.05 to 0.02 in linearly step in monitor intervals. In this case, as we should refresh in each iteration, while n = 2000, CPU time mainly cost by generate whose time cost is . 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.
Case 6: see Table 6.
Table 5. λ = 5, μJ = 0.05, = 0.05, L = 80, U = 140.
Table 6. σ, r are not constants.
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.
 Farnoosh, R., Sobhani, A., Rezazadeh, H. and Beheshti, M. (2016) Numerical Method for Discrete Double Barrier Option Pricing with Time-Dependent Parameters. Computational Economics, 70, 2006-2013.
 Golbabai, A., Ballestra, L.V. and Ahmadian, D. (2014) A Highly Accurate Finite Element Method to Price Discrete Double Barrier Options. Computational Economics, 44, 153-173.
 Ahmadian, D. and Ballestra, L.V. (2015) A Numerical Method to Price Discrete Double Barrier Options under a Constant Elasticity of Variance Model with Jump Diffusion. International Journal of Computer Mathematics, 92, 2310-2328.