# 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;

-

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.

-
I understand your point that randn() produces a random number each time. I have made changes (see code) so that instead of calculating $xT$, I use $xbefore$ for the payoff calculation but still the value I get is less than 9.013. – Moneyness Sep 3 '14 at 17:09
@Moneyness I ran your new script couple times and got values between 9.5-13.5, can you recheck just to make sure we have same code here? – emcor Sep 4 '14 at 10:19
I have just done it a few times I get values like 9.4552 8.7864 9.2066 11.4286 9.6279 8.7172 1.4979e5 16.6024 58.9300 1.6e7. There is too much variability. Also jumpTimes = 0:exprnd(lam):T seems to produce too many jump times for a process with $\lambda = 2$. Also check the values of $p$ in the code. Are also getting 1 sometimes during loop. – Moneyness Sep 4 '14 at 15:48

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.

-
".........jumps repeating at the same interval", Not true, it does not repeat, exprnd(lam) is the increment which produces different output each time. So it's a random increase. – Moneyness Sep 5 '14 at 18:06
Humor me. Try printing the jumpTimes vector and see what it contains. – q.t.f. Sep 5 '14 at 19:01
Alright, I really don't know why its does not work. Anywhere I have changed the calculations of jumpTimes. The results are now quite stable around 7.7 but still not 9.013. – Moneyness Sep 5 '14 at 19:19