Debugging Monte Carlo Barrier Option Pricing with Jump Diffusion

# Valuation of barrier options in Jump diffusion model


# Valuation of barrier options in Jump diffusion model












I am trying to evaluate the value of a Barrier option using Monte carlo method. The stock follows a jump diffusion model. I am using the method described in Metwally and Atiya. The authors describe the steps so writing the algorithm in matlab say, should be easy. I have implemented the the first algorithm in matlab, described in this paper but my results are not the same as those of the authors. For example, my code below gives 5.1 but according to the authors results it should be 9.013.

The other problem I have is that the probability $P_i$ is negative or more than 1 sometimes during simulation. Could the formula in the paper be wrong?. How can it be coded to avoid this. I have used it as it is in the paper.

```
clc 
clear all
t = cputime;

%%%%%%%%%%%%%%%%%%% Parameters %%%%%%%%%%%%%%%%
S0 = 100.0;
X = 110.0;
H = 85.0;
R = 1.0;
r = 0.05;
sigma = 0.25;
T = 1.0;

%%%%%%%%%%%%%%%% Jump Parameters %%%%%%%%%%%%%%
lam = 2.0;
muA = 0.0;
 sigmaA = 0.1;

%%%%%%%%%%%%%%% calculated parameters %%%%%%%%%%
k = exp(muA+0.5*sigmaA*sigmaA)-1;
c = r-0.5*sigma^2-lam*k;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 1e5; % Monte carlo runs
DP = zeros(N,1);
for n = 1:N
I = 1;
jumpTimes = 0:exprnd(lam):T; %interjump times Exp(lam)
K = size(jumpTimes,2)-1;
jumpTimes(end+1) = T;
x = log(S0);
for i = 1:K+1
    tau = jumpTimes(i+1)-jumpTimes(i);
    xbefore = x + c*tau + sigma*sqrt(tau)*randn();

    p = 1.0-exp(-2.0*(log(H)-x)*(log(H)-xbefore)/(tau*sigma^2));
    p = p*(xbefore > log(H));
    b = (jumpTimes(i+1)-jumpTimes(i))/(1.0-p);
    s = jumpTimes(i)+b*rand();

    if s <= jumpTimes(i+1) && s >= jumpTimes(i)
    gamma = exp(-(x-xbefore+c*tau)^2/(2*sigma^2*tau))/(sigma*sqrt(2*pi*tau));
    g = (x-log(H))/(2*gamma*pi*sigma^2)*(s-jumpTimes(i))^(-1.5)*(jumpTimes(i+1)-s)^(-0.5)*...
        exp(-((xbefore-log(H)-c*(jumpTimes(i+1)-s))^2/(2*(jumpTimes(i+1)-s)*sigma^2)+...
        (x-log(H)+c*(s-jumpTimes(i)))^2/(2*(s-jumpTimes(i))*sigma^2)));
    DP(n)= R*b*g*exp(-r*s);
    I = 0;
    break
    end
    A = muA + sigmaA*randn();
    xafter = xbefore + A;
    if xafter <= log(H)
    DP(n) = R*exp(-r*jumpTimes(i+1));
    I = 0;
    break
    end
    x = xafter;
end
if I==1 % no crossing happened
    DP(n) = exp(-r*T)*max(exp(xbefore) - X, 0.0);
end

end

DOC = mean(DP)
e = (cputime-t)/60;
```

## Answer by emcor (score 3)

https://quant.stackexchange.com/a/14608

The error is, you are not storing the random numbers for the same path at the end:

```
xbefore = x + c*tau + sigma*sqrt(tau)*randn()

A = muA + sigmaA*randn();

xafter = xbefore + A;
```

But then at end you set a different path here by creating a new random number:

```
xT = log(S0)+(c+muA*lambda)*T+sqrt((sigma^2+(muA^2+sigmaA^2)*lambda)*T)*randn();
```

randn() generates a new random variable each time, so you need to store the specific path to make sure that it takes the same path at xT, and not take a new path by using another *randn() for xT.

e.g. use another dummy variables:

```
rand = randn();
```

to use the same number at each loop, and

```
path = path+randn();
```

to store the total path for xT at the end.

## Answer by q.t.f. (score 0)

https://quant.stackexchange.com/a/14641

Your code has a bug that it picks one random jump time by exprnd(lam), then makes a vector of jump times with jumps repeating at the same interval. It should pick a new random time for each jump.