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

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?


1 Answer 1


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);
        Z = normrnd(0,1);
        S(t,s) = S(t-1,s)*exp(drift*dt + sigma*sqrt(dt)*Z + J);

for more details please check volopta

  • $\begingroup$ 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}$. $\endgroup$
    – user39039
    Commented Aug 13, 2016 at 17:11
  • $\begingroup$ Ok. You can use Cumsum . $\endgroup$
    – user16651
    Commented Aug 13, 2016 at 17:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.