I found this great post deriving the solution to the Merton Jump-Diffusion SDE
$$S_t = S_0\exp\left(\left(\mu - \frac{\sigma^2}{2}\right)t + \sigma W_t\right)\prod_{j=0}^{N_t}V_j$$
The first part of the equation (inside the exponent) just looks like Geometric Brownian Motion (GBM) to me, which I am ok with. However I am unsure about the product of $V$'s; it seems $N_t$ is a poisson process with intensity $\lambda$ and jump sizes $Y_j$ which are iid Normally distributed $Y_j ∼ N(\mu, \delta^2)$. But I don't know how to include these parameters in my simulation.
For GBM python code, I have followed the wiki page
import numpy as np
import math
dt = 1/252 # time increments
mu = 0.01 # stock price daily drift (1%)
sigma = 0.14 # stock price daily volatility (14%)
T = 400 # number of periods (dyas) to simulate
So = 100 # initial stock price
drift = mu - 0.5*sigma**2
diffusion = sigma*np.random.normal(loc=0, scale=math.sqrt(dt), size=(T,)) # normal distribution
S = np.exp(drift + diffusion)
S = So*S.cumprod() # GBM
I also found this post and this jupyter notebook; the latter has explicit code on the Merton jump model, however it is not very well commented and so I don't understand what is happening.
How can I interpret the product of $V$'s in the equation above, and how can I include these jumps in our simulation? (BONUS: preferably using numpy instead of explicit loops).
