ABSTRACT

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 inte- gral, 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.

Keywords:

Discrete Double Barrier Option, Merton Jump Diffusion Model, Discrete Convolution, Monte Carlo Method

1. 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 non- smooth 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.

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, ${\sigma }_J^2$ $=T$ 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 $\max \left(S_T-K,0\right),$ $L\le K\le U,L\le S_0\le U$ ; 5) the stock price can change continuously.

Merton jump diffusion model assumes that the stock price $S_t$ is in accordance with the following stochastic differential equation,

$\frac{\text{d}{S}_t}{S_{t-}}=\left(\mu -\lambda E\left[J_t\right]\right)\text{d}t+\sigma \text{d}{W}_t+J_t\text{d}{N}_t,$ (1)

where $\ln \left(1+J_t\right)~N\left({\mu }_J,{\sigma }_J^2\right)$ , W_t is standard Wiener process, $N_t$ is Poisson

process with intensity λ. We denote $Y_t=\ln \left(1+J_t\right)$ , $\gamma =\mu -\lambda E\left[J_t\right]-\frac{1}{2}{\sigma }^2$ . By

Itô Lemma, we obtain

$\text{d}\ln S_t=\gamma \text{d}t+\sigma \text{d}{W}_t+Y_t\text{d}{N}_t.$ (2)

We denote $A_i=\left\{S_{T_i}\in \left[L,U\right]\right\},i=1,2,\cdots ,m$ . In risk-neutral measure, if we assume risk-free interest rate r is constant. (Thus, the drift $\mu$ is replaced by the interest rate r.) By the definition of the option its price at the beginning time is equal to

${\text{e}}^{-rT}E\left[\max \left(S_T-K,0\right),1_{A_1}{1}_{A_2}\cdots 1_{A_m}\right].$ (3)

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

$V\left(S,T\right)={\displaystyle {\int }_{\left(x_1,x_2,\cdots ,x_m\right)\in D}\left(L{\text{e}}^{x_1+x_2+\cdots +x_m}-K\right)\tilde{f}\left(x_1,x_2,\cdots ,x_m\right)\text{d}{x}_m\cdots \text{d}{x}_2\text{d}{x}_1},$ (4)

where

$\begin{array}{l}D=\{x_1,x_2,\cdots ,x_m|x_1,x_1+x_2,\cdots ,x_1+x_2+\cdots +x_{m-1}\in \left[0,\ln \frac{U}{L}\right],\\ x_1+x_2+\cdots +x_m\in \left[\ln \frac{K}{L},\ln \frac{U}{L}\right]\},\end{array}$

$\tilde{f}\left(x_1,x_2,\cdots ,x_m\right)=h\left(x_1-\ln \frac{L}{S_0}\right){\displaystyle {\prod }_{i=2}^{m}h\left(x_i\right)},$

$h\left(x\right)={\displaystyle {\sum }_{i=0}^{\infty }\frac{{\left(\lambda T/m\right)}^i}{i!}{\text{e}}^{-\frac{\lambda T}{m}}\varphi \left(\frac{\lambda T}{m}+i{\mu }_J,\frac{{\sigma }^{2}T}{m}+i{\sigma }_J^2\right)}$ ,

$\varphi \left(\mu ,{\sigma }^2\right)$ is the p.d.f. of $N\left(\mu ,{\sigma }^2\right)$ .

Proof. Step 1. Assuming that the p.d.f. of random variable is known, we calculate the analytical result.

We denote $\widetilde{{\xi }_i}=\ln S_{T_i}-\ln S_{T_{i-1}},i=1,2,\cdots ,m$ , we have $S_{T_i}=S_0{\text{e}}^{{\displaystyle {\sum }_{j=1}^i\widetilde{{\xi }_j}}}$ . Let ${\xi }_1=$ Math_25#, We get $S_{T_i}=L{\text{e}}^{{\displaystyle {\sum }_{j=1}^i{\xi }_j}}$ . For $i=1,2,\cdots ,m$ ,

$A_i=\left\{L\le S_{T_i}\le U\right\}=\left\{\ln L\le \ln S_{T_i}\le \ln U\right\}=\left\{0\le {\displaystyle {\sum }_{j=1}^i{\xi }_j\le \ln \frac{U}{L}}\right\}$ .

From Equation (4), the value of this barrier option is,

$\begin{array}{l}V\left(S,T\right)\\ ={\text{e}}^{-rT}E\left[\max \left(S_T-K,0\right)1_{A_1}{1}_{A_2}\cdots 1_{A_m}\right]\\ ={\text{e}}^{-rT}E\left[\max \left(L{\text{e}}^{{\displaystyle {\sum }_{j=1}^m{\xi }_j}}-K,0\right)1_{\left\{0\le {\displaystyle {\sum }_{j=1}^1{\xi }_j}\le \ln \frac{U}{L}\right\}}{1}_{\left\{0\le {\displaystyle {\sum }_{j=1}^2{\xi }_j}\le \ln \frac{U}{L}\right\}}\cdots 1_{\left\{0\le {\displaystyle {\sum }_{j=1}^m{\xi }_j}\le \ln \frac{U}{L}\right\}}\right]\\ ={\text{e}}^{-rT}E\left[\left(L{\text{e}}^{{\displaystyle {\sum }_{j=1}^m{\xi }_j}}-K\right)1_{\left\{0\le {\displaystyle {\sum }_{j=1}^1{\xi }_j}\le \ln \frac{U}{L}\right\}}{1}_{\left\{0\le {\displaystyle {\sum }_{j=1}^2{\xi }_j}\le \ln \frac{U}{L}\right\}}\cdots 1_{\left\{\ln \frac{K}{L}\le {\displaystyle {\sum }_{j=1}^m{\xi }_j}\le \ln \frac{U}{L}\right\}}\right]\\ ={\text{e}}^{-rT}{\displaystyle {\int }_{\left(x_1,x_2,\cdots ,x_m\right)\in D}\left(L{\text{e}}^{x_1+x_2+\cdots +x_m}-K\right){\displaystyle {\prod }_{i=1}^m{h}_i\left(x_i\right)\text{d}{x}_m\cdots \text{d}{x}_2\text{d}{x}_1}},\end{array}$ (5)

where

$\begin{array}{l}D=\{x_1,x_2,\cdots ,x_m|x_1,x_1+x_2,\cdots ,x_1+x_2+\cdots +x_{m-1}\in \left[0,\ln \frac{U}{L}\right],\\ x_1+x_2+\cdots +x_m\in \left[\ln \frac{K}{L},\ln \frac{U}{L}\right]\},\end{array}$

$h\left(x_i\right)$ is p.d.f. of ${\xi }_i,i=2,3,\cdots ,m$ .

Step 2. Calculate $h\left(x_i\right)$ .

For $i=2,3,\cdots ,m$ from Equation (2), We know ${\xi }_i$ is sum of independent

random variables including a normal random variable $\eta ~N\left(\frac{\gamma T}{m},\frac{{\sigma }^{2}T}{m}\right)$ and a

series of normal random variables ${\eta }_k~N\left({\mu }_J,{\sigma }_J^2\right)$ . The length of the series is follow a Poisson distribution with intensity $\lambda T/m$ . As we know while the length

is j, sum of these variables ${\xi }_i~N\left(\frac{\gamma T}{m}+j{\mu }_J,\frac{{\sigma }^{2}T}{m}+j{\sigma }_J^2\right)$ , we get the p.d.f. of ${\xi }_i$ ,

$h_i={\displaystyle {\sum }_{j=0}^{\infty }\frac{{\left(\lambda T/m\right)}^j}{j!}{\text{e}}^{-\frac{\lambda T}{m}}\varphi \left(\frac{\gamma T}{m}+j{\mu }_J,\frac{{\sigma }^{2}T}{m}+j{\sigma }_J^2\right),}$ (6)

where $\varphi \left(\mu ,{\sigma }^2\right)$ is the p.d.f. of $N\left(\mu ,{\sigma }^2\right)$ . Because $h_i\left(x\right)$ has no relation with i, in $i=2,3,\cdots ,m$ , we denote $h_i\left(x\right)$ with $h\left(x\right)$ . In particular, as there is

a constant shift in ${\xi }_1$ , $h_1\left(x\right)=h\left(x-\ln \frac{L}{S_0}\right)$ .

Plug $h\left(x\right)$ into Equation (5), we finally obtain (4). □

Specially, if parameters is change to different constant in different monitor interval, $h_i\left(x\right)$ 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 ${\xi }_i,i=1,2,\cdots ,m$ decide the option value at t = T. As we have already known the distribution of ${\xi }_i$ , we can use Monte Carlo method to simulate ${\xi }_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.

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 $d=\ln \frac{U}{L}/n$ , where n is an in-

teger number. We used discrete random variable ${\xi }_{i,n}$ to approximate ${\xi }_i$ , ${\xi }_{i,n}$ satisfy

$P\left({\xi }_{i,n}=kd\right)=p_{k,n},$ (7)

where

$p_{k,n}=P\left(\left(k-0.5\right)d\le {\xi }_2<\left(k+0.5\right)d\right)$ , (8)

$k=\cdots ,-1,0,1,2,\cdots ,i=2,3,\cdots ,m$ . Specially, we discrete random variable ${\xi }_{1,n}$ to approximate ${\xi }_1$ , ${\xi }_{1,n}$ satisfy

$P\left({\xi }_{i,n}=kd\right)={\tilde{p}}_{k,n},$ (9)

where

${\tilde{p}}_{k,n}=P\left(\left(k-0.5\right)d\le {\xi }_1<\left(k+0.5\right)d\right).$ (10)

From Equation (5), we know

$\begin{array}{l}V\left(S,T\right)\\ ={\text{e}}^{-rT}E\left[\max \left(L{\text{e}}^{{\displaystyle {\sum }_{j=1}^m{\xi }_j}}-K,0\right)1_{\left\{0\le {\displaystyle {\sum }_{j=1}^1{\xi }_j}\le \ln \frac{U}{L}\right\}}{1}_{\left\{0\le {\displaystyle {\sum }_{j=1}^2{\xi }_j}\le \ln \frac{U}{L}\right\}}\cdots 1_{\left\{0\le {\displaystyle {\sum }_{j=1}^m{\xi }_j}\le \ln \frac{U}{L}\right\}}\right].\end{array}$ (11)

So we can use ${\tilde{V}}_n\left(S,T\right)$ to approximate $V\left(S,T\right)$ , where

(12)

We denote ${\zeta }_{k,n}={\displaystyle {\sum }_{j=1}^k{\xi }_{j,n}}{1}_{\left\{0<{\displaystyle {\sum }_{j=1}^1{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}\cdots 1_{\left\{0<{\displaystyle {\sum }_{j=1}^k{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}$ , $k=1,2,\cdots ,m$ . From

Equation (12), we know

${\tilde{V}}_n\left(S,T\right)={\text{e}}^{-rT}E\left(\max \left(L{\text{e}}^{{\zeta }_{m,n}}-K,0\right)\right).$ (13)

So we can calculate ${\tilde{V}}_n\left(S,T\right)$ by calculate ${\zeta }_{m,n}$ . From the definition of ${\zeta }_{1,n}$ , we know

${\zeta }_{1,n}={\xi }_{1,n}{1}_{\left\{0<{\xi }_{j,n}\le \ln \frac{U}{L}\right\}}.$ (14)

So the value of ${\zeta }_{1,n}$ may be $0,d,2d,\cdots ,nd$ , and satisfy

$P\left({\zeta }_{1,n}=kd\right)=\{\begin{array}{l}{\tilde{p}}_{k,n},k=1,2,\cdots ,n\\ {\widehat{p}}_{0,n},k=0\end{array}$ (15)

where ${\widehat{p}}_{0,n}=1-{\displaystyle {\sum }_{k=1}^n{\tilde{p}}_{k,n}}$ . Notice that ${\zeta }_{1,n}\ge 0$ . ${\zeta }_{1,n}>0$ is equivalent to

$1_{\left\{0<{\xi }_{1,n}\le \ln \frac{U}{L}\right\}}=1$ . ${\zeta }_{1,n}=0$ is equivalent to $1_{\left\{0<{\xi }_{1,n}\le \ln \frac{U}{L}\right\}}=0$ . In fact, we have the

following lemma:

Lemma 3.1. For $k=1,2,\cdots ,m,{\zeta }_{k,n}\ge 0$ , and

${\zeta }_{k,n}=0\iff 1_{\left\{0<{\displaystyle {\sum }_{j=1}^1{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}{1}_{\left\{0<{\displaystyle {\sum }_{j=1}^2{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}\cdots 1_{\left\{0<{\displaystyle {\sum }_{j=1}^k{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}=0,$ (16)

${\zeta }_{k,n}>0\iff 1_{\left\{0<{\displaystyle {\sum }_{j=1}^1{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}{1}_{\left\{0<{\displaystyle {\sum }_{j=1}^2{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}\cdots 1_{\left\{0<{\displaystyle {\sum }_{j=1}^k{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}=1.$ (17)

If ${\zeta }_{k,n}>0$ , then ${\zeta }_{k,n}={\displaystyle {\sum }_{j=1}^k{\xi }_{j,n}}$ .

Proof. From the definition of ${\zeta }_{k,n}$ , it can be easily proved.□

For $k=2,3,\cdots ,m$ , by lemma 3.1, we have

$\begin{array}{c}{\zeta }_{k,n}=\{\begin{array}{l}{\displaystyle {\sum }_{j=1}^k{\xi }_{j,n}}{1}_{\left\{0<{\displaystyle {\sum }_{j=1}^k{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}},1_{\left\{0<{\displaystyle {\sum }_{j=1}^1{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}\cdots 1_{\left\{0<{\displaystyle {\sum }_{j=1}^{k-1}{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}=1\\ 0,1_{\left\{0<{\displaystyle {\sum }_{j=1}^1{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}\cdots 1_{\left\{0<{\displaystyle {\sum }_{j=1}^{k-1}{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}=0\end{array}\\ =\{\begin{array}{l}\left({\zeta }_{k-1,n}+{\xi }_{k,n}\right)1_{\left\{0<{\displaystyle {\sum }_{j=1}^k{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}},1_{\left\{0<{\displaystyle {\sum }_{j=1}^1{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}\cdots 1_{\left\{0<{\displaystyle {\sum }_{j=1}^{k-1}{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}=1\\ 0,1_{\left\{0<{\displaystyle {\sum }_{j=1}^1{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}\cdots 1_{\left\{0<{\displaystyle {\sum }_{j=1}^{k-1}{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}=0\end{array}.\end{array}$ (18)

So we know $0\le {\zeta }_{k,n}\le \ln \frac{U}{L}=nd$ . As the possible value of ${\zeta }_{k,n}$ is an integer

multiple of d, the value of ${\zeta }_{k,n}$ may be $0,d,2d,\cdots ,nd$ . For $i=0,1,2,\cdots ,n$ , let

$q_{i,k,n}=P\left({\zeta }_{k,n}=id\right)$ . Notice that ${\zeta }_{k,n}>0\iff 1_{\left\{0<{\displaystyle {\sum }_{j=1}^1{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}\cdots 1_{\left\{0<{\displaystyle {\sum }_{j=1}^k{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}=1$ ,

from Equation (18),

$\begin{array}{c}{q}_{i,k,n}=P\left({\zeta }_{k,n}=id\right)\\ =P\left(\left({\zeta }_{k,n-1}+{\xi }_{k,n}\right)1_{\left\{0<{\displaystyle {\sum }_{j=1}^k{\xi }_{j,n}}\le \ln \frac{U}{L}\right\}}{1}_{\left\{{\zeta }_{k,n-1}\right\}}=id\right)\\ ={\displaystyle {\sum }_{j=1}^{n}P\left({\zeta }_{k,n-1}=jd\right)P\left({\xi }_{k,n}=\left(i-j\right)d\right)}\\ ={\displaystyle {\sum }_{j=1}^n{q}_{j,k-1,n}{p}_{i-j,n}},\end{array}$ (19)

$q_{0,k,n}=1-{\displaystyle {\sum }_{i=1}^n{q}_{i,k,n}}.$ (20)

Using Equation (19) and Equation (20) we can calculate all $q_{i,k,n}$ step by step. Notice that $-\left(n-1\right)\le i-j\le n-1$ , we just need to precompute $p_{k,n}$ for $k=-n+1,-n+2,\cdots ,n-1$ . Plug it into Equation (13), we have

$\begin{array}{c}{\tilde{V}}_n\left(S,T\right)={\text{e}}^{-rT}E\left(\max \left(L{\text{e}}^{{\zeta }_{m,n}}-K,0\right)\right)\\ {\text{=e}}^{-rT}{\displaystyle {\sum }_{i=1}^n{q}_{i,m,n}\max \left(L{\text{e}}^{id}-K,0\right)}.\end{array}$ (21)

It means we could calculate ${\tilde{V}}_n\left(S,T\right)$ by Equation (19), Equation (20) and Equation (21) after we got $p_{k,n}$ and ${\tilde{p}}_{k,n}$ . Actually, as we don’t need $q_{0,k,n}$ to calculate ${\tilde{V}}_n\left(S,T\right)$ 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 ${\xi }_i$ To get $p_{k,n}$ and ${\tilde{p}}_{k,n}$ . Denote

$g_j\left(x\right)=\frac{{\left(\lambda T/m\right)}^j}{j!}{\text{e}}^{-\frac{\lambda T}{m}}\varphi \left(\frac{\gamma T}{m}+j{\mu }_J,\frac{{\sigma }^{2}T}{m}+j{\sigma }_J^2\right)\left(x\right)$ .

Because ${\displaystyle {\sum }_{j=0}^{\infty }{g}_j\left(x\right)}$ is uniform convergence, for $i=2,3,\cdots ,m$ , by the definition of $p_{k,n}$ , we have

$\begin{array}{c}{p}_{k,n}=P\left(\left(k-0.5\right)d\le {\xi }_2<\left(k+0.5\right)d\right)={\displaystyle {\int }_{\left(k-0.5\right)d}^{\left(k+0.5\right)d}h\left(x\right)\text{d}x}\\ ={\displaystyle {\int }_{\left(k-0.5\right)d}^{\left(k+0.5\right)d}{\displaystyle {\sum }_{j=0}^{\infty }{g}_j\left(x\right)}\text{d}x}={\displaystyle {\sum }_{j=0}^{\infty }{\displaystyle {\int }_{\left(k-0.5\right)d}^{\left(k+0.5\right)d}{g}_j\left(x\right)\text{d}x}}\\ ={\displaystyle {\sum }_{j=0}^{\infty }\frac{{\left(\lambda T/m\right)}^j}{j!}{\text{e}}^{-\frac{\lambda T}{m}}(\Phi \left(\frac{\gamma T}{m}+j{\mu }_J,\frac{{\sigma }^{2}T}{m}+j{\sigma }_J^2\right)\left(\left(k+0.5\right)d\right)}\\ -\Phi \left(\frac{\gamma T}{m}+j{\mu }_J,\frac{{\sigma }^{2}T}{m}+j{\sigma }_J^2\right)\left(\left(k-0.5\right)d\right))\\ ={\displaystyle {\sum }_{j=0}^{\infty }\frac{{\left(\lambda T/m\right)}^j}{j!}{\text{e}}^{-\frac{\lambda T}{m}}(\Phi \left(0,1\right)\left(\frac{\left(k+0.5\right)d-\gamma T/m-j{\mu }_J}{\sqrt{{\sigma }^{2}T/m+j{\sigma }_J^2}}\right)}\\ -\Phi \left(0,1\right)\left(\frac{\left(k-0.5\right)d-\gamma T/m-j{\mu }_J}{\sqrt{{\sigma }^{2}T/m+j{\sigma }_J^2}}\right)),\end{array}$ (22)

where $\Phi \left(\mu ,{\sigma }^2\right)$ is c.d.f. of $N\left(\mu ,{\sigma }^2\right)$ , and $\Phi \left(0,1\right)$ is c.d.f. of standard normal distribution. Similarly, from the p.d.f. of ${\xi }_1$ , we can get

$\begin{array}{c}{\tilde{p}}_{k,n}={\displaystyle {\sum }_{j=0}^{\infty }\frac{{\left(\lambda T/m\right)}^j}{j!}{\text{e}}^{-\frac{\lambda T}{m}}(\Phi \left(0,1\right)\left(\frac{\left(k+0.5\right)d-\ln \left(L/S_0\right)-\gamma T/m-j{\mu }_J}{\sqrt{{\sigma }^{2}T/m+j{\sigma }_J^2}}\right)}\\ -\Phi \left(0,1\right)\left(\frac{\left(k-0.5\right)d-\ln \left(L/S_0\right)-\gamma T/m-j{\mu }_J}{\sqrt{{\sigma }^{2}T/m+j{\sigma }_J^2}}\right)).\end{array}$ (23)

Unfortunately, we can’t calculate $p_{k,n}$ and ${\tilde{p}}_{k,n}$ directly although we have got their analytic expression. Let

$\begin{array}{c}{p}_{k,n,n^{\prime }}={\displaystyle {\sum }_{j=0}^{n^{\prime }}\frac{{\left(\lambda T/m\right)}^j}{j!}{\text{e}}^{-\frac{\lambda T}{m}}(\Phi \left(0,1\right)\left(\frac{\left(k+0.5\right)d-\gamma T/m-j{\mu }_J}{\sqrt{{\sigma }^{2}T/m+j{\sigma }_J^2}}\right)}\\ -\Phi \left(0,1\right)\left(\frac{\left(k-0.5\right)d-\gamma T/m-j{\mu }_J}{\sqrt{{\sigma }^{2}T/m+j{\sigma }_J^2}}\right)),\end{array}$ (24)

$\begin{array}{c}{\tilde{p}}_{k,n,n^{\prime }}={\displaystyle {\sum }_{j=0}^{n^{\prime }}\frac{{\left(\lambda T/m\right)}^j}{j!}{\text{e}}^{-\frac{\lambda T}{m}}(\Phi \left(0,1\right)\left(\frac{\left(k+0.5\right)d-\ln \left(L/S_0\right)-\gamma T/m-j{\mu }_J}{\sqrt{{\sigma }^{2}T/m+j{\sigma }_J^2}}\right)}\\ -\Phi \left(0,1\right)\left(\frac{\left(k-0.5\right)d-\ln \left(L/S_0\right)-\gamma T/m-j{\mu }_J}{\sqrt{{\sigma }^{2}T/m+j{\sigma }_J^2}}\right)),\end{array}$ (25)

from above, we know $\underset{n^{\prime }\to \infty }{\lim }{p}_{k,n,n^{\prime }}=p_{k,n},\underset{n^{\prime }\to \infty }{\lim }{\tilde{p}}_{k,n,n^{\prime }}={\tilde{p}}_{k,n}$ , and $p_{k,n,n^{\prime }},{\tilde{p}}_{k,n,n^{\prime }}$ can

be directly calculate. By using $p_{k,n,n^{\prime }},{\tilde{p}}_{k,n,n^{\prime }}$ to approximate $p_{k,n},{\tilde{p}}_{k,n}$ , we finally get an approximation of ${\tilde{V}}_n\left(S,T\right)$ , and denote it with ${\tilde{V}}_{n,n^{\prime }}\left(S,T\right)$ , who is a computable approximation of $V\left(S,T\right)$ . If parameter is different in different monitor interval, we should just refresh $g_j\left(x\right)$ 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 S_t. Then we can use the mean price of option base on these S_t to estimate $V\left(S,T\right)$ according to the law of large numbers. While the price of option just bases on the value of $S_t$ on monitor date, we just need to simulation these values on monitor date. From Equation (2) we can conveniently simulate $\ln S_t$ then get $S_t$ indirectly. In other words, we simulate ${\xi }_i,i=1,2,\cdots ,m$ . As we know ${\xi }_i$ is sum of independent random variables including a normal random variable

$\eta ~N\left(\frac{\gamma T}{m},\frac{{\sigma }^{2}T}{m}\right)$ and a series of normal random variables ${\eta }_k~N\left({\mu }_J,{\sigma }_J^2\right)$ .

The length of the series is follow a Poisson distribution with intensity $\lambda T/m$ . We simulate $\eta$ , and use a while loop to simulate ${\eta }_k$ , and exit the loop with probability ${\text{e}}^{-\lambda T/m}$ each round. As usual, estimation error of Moto Carlo simulation is Proportional to $1/\sqrt{N}$ while estimated standard deviation is ${\sigma }_V/\sqrt{N}$ , ${\sigma }_V$ is standard deviation of $V\left(S,T\right)$ , $N$ is the number of simulation path. Because we can know the option value if $S_t$ 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 $\lambda =0$ in case 1, case 2 and case 3 to compare the result with other literatures. It shows when $\lambda$ 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 $S_0$ , 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 $S_0$ and parameters $K=100,$ $\sigma =0.25,$ $T=0.5,$ $r=0.05,$ $L=95,$ $U=110$ , see Table 2.

Case 3: Prices of discrete double knock-out call option monitored weekly (25 times) for different values of the underlying asset $S_0$ and parameters $K=100,$ $\sigma =0.25,$ $T=0.5,$ $r=0.05,$ $L=95,$ $U=110$ , see Table 3.

Secondly, we compare the efficiency of numerical method and simulation method in case 4 and case 5. Poisson intensity $\lambda >0$ , i.e., jump is considered. We set $\lambda =5$ , ${\mu }_J=0.05,{\sigma }_J^2=0.05$ , other parameters except U, L and $S_0$ 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 $\sigma$ and $r$ . $\sigma$ will increase from 0.1 to 0.25 in linearly step in monitor intervals. $r$ will decrease from 0.05 to 0.02 in linearly step in monitor intervals. In this case, as we should refresh $g_j\left(x\right)$ in each iteration, while n = 2000, CPU time mainly cost by generate $p_{k,n,n^{\prime }}$ whose time cost is $O\left(mn{n}^{\prime }\right)$ . 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.

5. 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.

Cite this paper

Li, M.J. (2017) Numerical Methods for Discrete Double Barrier Option Pricing Based on Merton Jump Diffusion Model. Open Journal of Statistics, 7, 446-458. https://doi.org/10.4236/ojs.2017.73032

References