Simulating Merton Jump-Diffusion Paths for Option Pricing

# How to simulate a jump-diffusion process?


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

## Answer by user16651 (score 3)

https://quant.stackexchange.com/a/29610

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