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;