# How to simulate a jump-diffusion process?

I would like to price Asian and Digital options under Merton's jump-diffusion model. To that end, I will have to simulate from a jump diffusion process.

In general, the stock price process is given by $$S(t) = S(0)e^{(r-q-\omega)t+X(t)},$$ where $\omega$ is the drift term that makes the discounted stock price process a martingale and where $X(t) = \sigma W(t)+\sum_{i=1}^{N(t)} Y_i$ is a jump-diffusion process. The process $N(t)$ is a Poisson process with intensity $\lambda$ and the jump sizes $Y_i$ are iid Normally distributed $Y_i \sim N(\mu, \delta^2)$.

So, the issue here is the simulation of $X(t)$. I would like to simulate $X(t)$ as follows: first, I construct a grid $[0,T]$ with grid step size $\Delta t = 1/252$ (number of trading days in one year). I now want to compute $X$ on this grid stepwise. Therefore: $$X_{n\Delta t} = X_{(n-1) \Delta t}+\sigma \sqrt{\Delta t} \epsilon_n + p_n,$$ where $\epsilon_n$ are standard normal random variables (easy to simulate in Matlab) and where $p_n$ are Compound Poisson random variables independent of $\epsilon_n$. The question is, how to simulate from these $p_n$?

I was thinking to write a compound Poisson generator myself:

function [p] = CPoissonGenerator(mu,delta,lambda,dt,rows,columns)
N = poissrnd((1/lambda)*dt,rows,columns);
Y = normrnd(mu,delta,1,N);
cP = cumsum([0,Y]);
p = cP(end);
end


However, I do not know how to generate a matrix of Poisson random variables without using a loop in the above code? This can be very time consuming. Furthermore, I'm not sure if the program is correct. Is there a better way to do this?

you don't need generate a matrix of Poisson random variables:

% S0 = spot price
% K  = strike price
% r = risk free rate
% q  = dividend yield
% sigma = volatility
% T = maturity
% lambda_J= jump frequency
% mu_J    = jump mean parameter
% sigma_J = jump volatility parameter
% N_T = number of time steps
% N_S = number of stock price paths
%----------------------------------------------------------------------------
dt = T/N_T;
kappa = exp(mu_J) - 1;
drift = r - q - lambda*kappa - 0.5*sigma^2;
S = zeros(N_T,N_S);
for s=1:N_S
S(1,s) = S0;
for t=2:N_T
J = 0;
if lambda_J ~= 0
Nt = poissrnd(lambda_J*dt);
if Nt > 0
for i=1:Nt
J = J + normrnd(mu_J - sigma_J^2/2,sigma_J);
end
end
end
Z = normrnd(0,1);
S(t,s) = S(t-1,s)*exp(drift*dt + sigma*sqrt(dt)*Z + J);
end


for more details please check volopta

• Code can be made easier by using 'Cumsum' command in Matlab. Still, I do not particularly like the use of the loops. It is quite time consuming, especially if the number of paths is $10^{5}$. – user39039 Aug 13 '16 at 17:11
• Ok. You can use Cumsum . – user16651 Aug 13 '16 at 17:15