Hi Euginio,
You code seems to error out as I can see some pitfalls which I have mentioned below:
DW = sqrt(dt)*randn(1,N)
This line forces identical brownian shocks on all paths.
XS(i+1,k)
This is outside the for loop hence it overwrites all but the last i, leaving mostly zeros.
Nt = poissrnd(lambda*dt, nsim, N);
... if Nt(i) >= lb
You create Nt as nsim × N, then treat it like a 1-D vector (Nt(i))
Below is the vectorised implementation of your code:
function S = simulateJDBernoulli(mu,sigma,lambda,a,b,p,T,N,nsim,S0)
dt = T/N;
nudt = (mu-0.5*sigma^2)*dt;
sqdt = sqrt(dt);
logS = zeros(N+1,nsim);
logS(1,:) = log(S0);
for k = 1:N
dW = sqdt*randn(1,nsim);
dN = poissrnd(lambda*dt,1,nsim);
% Jump magnitudes (vectorised)
J = (rand(1,nsim) < p).*a + (rand(1,nsim) >= p).*b;
J = J .* dN; % if dN>1 we add multiple jumps
% One-step update in log-space
logS(k+1,:) = logS(k,:) + nudt + sigma*dW + J;
end
S = exp(logS);
end
You can then call the function:
paths = simulateJDBernoulli(0.03,0.12,2.0,0.10,-0.05,0.6, ...
2.0,1000,50,100);
plot(0:1000,paths(:,1:5)), grid on
title('Five sample paths of Bernoulli jump–diffusion')
Output:
Hope this helps!
