Simulating Kou's Asymmetric Double Exponential Random Variable

In Python, I import NumPy module to generate Kou's a vector of Asymmetric Double Exponential Random Variables (ADERV). I attempt to apply Glasserman's method for simulating the aforementioned from his 2004 book Monte Carlo Methods in Financial Engineering on page 139, which is the following (with slight difference in notation):

• generate $N \sim \text{Poisson}(\lambda (t_{i+1}-t_{i}))$ where $\lambda$ is the number of jumps per annum
• generate $K \sim \text{Binomial(N,p)}$ where $p$ is the probability of upward jumps
• generate $R_{1} \sim \text{Exponential}(K / \eta_{1})$ and $R_{2} \sim \text{Exponential}((N-K) / \eta_{2})$ where $1/\eta_{1}$ and $1/\eta_{2}$ are the magnitude of up and down jumps, respectively.
• set $M=R_{1} - R_{2}$ where $M$ is the ADERV

Given $t_{i+1}-t_{i} = 10 / 252$, below is the code to generate ADERVs

import numpy as np import seaborn as sns import import matplotlib.pyplot as plt n = np.random.poisson(lam=20 * (10 / 252), size=100_000) k = np.random.binomial(n=n, p=0.50) r1 = np.random.exponential(scale=k / 25) r2 = np.random.exponential(scale=(n - k) / 25) m = r1 - r2 sns.displot(m, kind='kde') 

Question is is the above code a correct implementation?

-----Further discussion

I discretize $\frac{dS_{t}}{S_{t-}} = \left[r - \left(\frac{p}{\eta_{1}} - \frac{q}{\eta_{2}} \right) \lambda \right]dt + \sigma dW_{t}^{\mathbb{Q}} + d\left(\sum_{j=1}^{N_{t}} V_{j} - 1 \right)$ where $W_{t}^{\mathbb{Q}}$ is a Brownian motion under risk-neutral measure, $N_{t}$ is a Poisson process with intensity rate $\lambda > 0$, and $\{V_{j}\}$ is a series of IID non-negative r.v. s.t. $\Upsilon = \log(V_{j})$ has an asymmetric double exponential distribution to simulate MC paths to price a contingent claim $h(X)$.

As an easier pricing exercise and to sanity check, I define $h(X)$ to be an AtM Vanilla European put with $K=100$ and price it with $r=0$, $\lambda=20$, $p=q=0.5$, $\sigma=0.3$, $\eta_{1}=\eta_{2}=25$, and $T=10/252$. I use MC to simulate 100k paths. Lastly, I do the same thing but define $h(X)$ to be AtM Vanilla European call with same strike as the put. The random number generator is Mersenne Twister with seed number 20240101 in NumPy library.

The put price is 2.9566 and the call is 3.0876. Per put-call parity, I expect equivalent (or very close). However, if I assume $\lambda=0$, then the put price and call is roughly inline (call 2.3745 vs. put 2.3887). It could be I'm compensating for the jumps incorrectly (i.e., the drift).