# Understanding and simulating the jumps in Merton's Jump-Diffusion SDE?

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