From d748fd7df5ddfe2012c98811629906665e9b3047 Mon Sep 17 00:00:00 2001
From: Luca <luca@quantmind.com>
Date: Sun, 12 Jan 2025 12:30:09 +0000
Subject: [PATCH] Jump diffusion and hestonj use the same parameters

---
 .vscode/launch.json                      |  2 +-
 notebooks/api/sp/jump_diffusion.rst      |  8 ----
 notebooks/examples/heston_vol_surface.md | 15 ++++---
 notebooks/models/bns.md                  |  2 +-
 notebooks/models/heston.md               | 28 ++++++------
 notebooks/models/heston_jumps.md         | 20 ++++++---
 notebooks/models/jump_diffusion.md       | 57 ++++++++++++++++--------
 notebooks/models/ou.md                   | 32 ++++++-------
 notebooks/reference/glossary.md          |  2 +-
 notebooks/theory/levy.md                 |  2 +-
 quantflow/sp/heston.py                   | 55 ++++++++++++-----------
 quantflow/sp/jump_diffusion.py           | 57 +++++++++++++-----------
 quantflow/utils/distributions.py         | 14 ++++++
 quantflow_tests/test_heston.py           |  3 +-
 quantflow_tests/test_jump_diffusion.py   | 13 +++---
 quantflow_tests/test_options_pricer.py   |  4 +-
 16 files changed, 182 insertions(+), 132 deletions(-)

diff --git a/.vscode/launch.json b/.vscode/launch.json
index a185bdd..53da042 100644
--- a/.vscode/launch.json
+++ b/.vscode/launch.json
@@ -14,7 +14,7 @@
             "args": [
                 "-x",
                 "-vvv",
-                "quantflow_tests/test_utils.py",
+                "quantflow_tests/test_jump_diffusion.py",
             ]
         },
     ]
diff --git a/notebooks/api/sp/jump_diffusion.rst b/notebooks/api/sp/jump_diffusion.rst
index 4bff41b..e1293eb 100644
--- a/notebooks/api/sp/jump_diffusion.rst
+++ b/notebooks/api/sp/jump_diffusion.rst
@@ -14,11 +14,3 @@ The most famous jump-diffusion model is the Merton model, which was introduced b
    :member-order: groupwise
    :autosummary:
    :autosummary-nosignatures:
-
-
-
-.. autoclass:: Merton
-   :members:
-   :member-order: groupwise
-   :autosummary:
-   :autosummary-nosignatures:
diff --git a/notebooks/examples/heston_vol_surface.md b/notebooks/examples/heston_vol_surface.md
index 9f24d40..3001152 100644
--- a/notebooks/examples/heston_vol_surface.md
+++ b/notebooks/examples/heston_vol_surface.md
@@ -11,21 +11,24 @@ kernelspec:
   name: python3
 ---
 
-# Heston Volatility Surface
+# HestonJ Volatility Surface
 
-Here we study the Implied volatility surface of the Heston model. The Heston model is a stochastic volatility model that is widely used in the finance industry to price options.
+Here we study the Implied volatility surface of the Heston model with jumps.
+The Heston model is a stochastic volatility model that is widely used in the finance industry to price options.
 
 ```{code-cell} ipython3
 from quantflow.sp.heston import HestonJ
+from quantflow.utils.distributions import DoubleExponential
 from quantflow.options.pricer import OptionPricer
 
-pricer = OptionPricer(model=HestonJ.exponential(
+pricer = OptionPricer(model=HestonJ.create(
+    DoubleExponential,
     vol=0.5,
     kappa=2,
-    rho=0.0,
-    sigma=0.6,
+    rho=-0.2,
+    sigma=0.8,
     jump_fraction=0.5,
-    jump_asymmetry=1.0
+    jump_asymmetry=0.2
 ))
 pricer
 ```
diff --git a/notebooks/models/bns.md b/notebooks/models/bns.md
index cc30bf6..1b57bb9 100644
--- a/notebooks/models/bns.md
+++ b/notebooks/models/bns.md
@@ -28,7 +28,7 @@ This means that the characteristic function of $y_t$ can be represented as
     &= {\mathbb E}\left[\exp{\left(-\tau_t \phi_{w, u} + i u \rho z_{\kappa t}\right)}\right]
 \end{align}
 
-$\phi_{w, u}$ is the characteristic exponent of $w_1$. The second equivalence is a consequence of $w$ and $\tau$ being independent, as discussed in [the time-changed Lévy](./levy.md) process section.
+$\phi_{w, u}$ is the characteristic exponent of $w_1$. The second equivalence is a consequence of $w$ and $\tau$ being independent, as discussed in [the time-changed Lévy](../theory/levy.md) process section.
 
 ```{code-cell}
 from quantflow.sp.bns import BNS
diff --git a/notebooks/models/heston.md b/notebooks/models/heston.md
index f012aa9..252923e 100644
--- a/notebooks/models/heston.md
+++ b/notebooks/models/heston.md
@@ -7,7 +7,7 @@ jupytext:
     format_version: 0.13
     jupytext_version: 1.16.6
 kernelspec:
-  display_name: Python 3 (ipykernel)
+  display_name: .venv
   language: python
   name: python3
 ---
@@ -31,20 +31,20 @@ This means that the characteristic function of $y_t=x_{\tau_t}$ can be represent
      &= e^{-a_{t,u} - b_{t,u} \nu_0}
 \end{align}
 
-```{code-cell}
+```{code-cell} ipython3
 from quantflow.sp.heston import Heston
 pr = Heston.create(vol=0.6, kappa=2, sigma=1.5, rho=-0.1)
 pr
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 # check that the variance CIR process is positive
 pr.variance_process.is_positive, pr.variance_process.marginal(1).std()
 ```
 
 ## Characteristic Function
 
-```{code-cell}
+```{code-cell} ipython3
 from quantflow.utils import plot
 m = pr.marginal(0.1)
 plot.plot_characteristic(m)
@@ -58,26 +58,26 @@ The immaginary part of the characteristic function is given by the correlation c
 
 Here we compare the marginal distribution at a time in the future $t=1$ with a normal distribution with the same standard deviation.
 
-```{code-cell}
+```{code-cell} ipython3
 plot.plot_marginal_pdf(m, 128, normal=True, analytical=False)
 ```
 
 Using log scale on the y axis highlighs the probability on the tails much better
 
-```{code-cell}
+```{code-cell} ipython3
 plot.plot_marginal_pdf(m, 128, normal=True, analytical=False, log_y=True)
 ```
 
 ## Option pricing
 
-```{code-cell}
+```{code-cell} ipython3
 from quantflow.options.pricer import OptionPricer
 from quantflow.sp.heston import Heston
 pricer = OptionPricer(Heston.create(vol=0.6, kappa=2, sigma=0.8, rho=-0.2))
 pricer
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 import plotly.express as px
 import plotly.graph_objects as go
 from quantflow.options.bs import black_call
@@ -89,7 +89,7 @@ fig.add_trace(go.Scatter(x=r.moneyness_ttm, y=b.time_value, name=b.name, line=di
 fig.show()
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 fig = None
 for ttm in (0.05, 0.1, 0.2, 0.4, 0.6, 1):
     fig = pricer.maturity(ttm).plot(fig=fig, name=f"t={ttm}")
@@ -102,17 +102,17 @@ The simulation of the Heston model is heavily dependent on the simulation of the
 
 The code implements algorithms from {cite:p}heston-simulation
 
-```{code-cell}
+```{code-cell} ipython3
 from quantflow.sp.heston import Heston
 pr = Heston.create(vol=0.6, kappa=2, sigma=0.8, rho=-0.4)
 pr
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 pr.sample(20, time_horizon=1, time_steps=1000).plot().update_traces(line_width=0.5)
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 import pandas as pd
 from quantflow.utils import plot
 
@@ -122,12 +122,12 @@ df = pd.DataFrame(mean, index=paths.time)
 plot.plot_lines(df)
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 std = dict(std=pr.marginal(paths.time).std(), simulated=paths.std())
 df = pd.DataFrame(std, index=paths.time)
 plot.plot_lines(df)
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 
 ```
diff --git a/notebooks/models/heston_jumps.md b/notebooks/models/heston_jumps.md
index 8679a85..51e872b 100644
--- a/notebooks/models/heston_jumps.md
+++ b/notebooks/models/heston_jumps.md
@@ -30,37 +30,41 @@ where $j_t$ is a double exponential Compound Poisson process which adds three ad
 * the jump percentage (fraction) contribution to the total variance
 * the jump asymmetry is defined as a parameter greater than 0; 1 means jump are symmetric
 
-The jump process is independent of the other Brownian motions.
+The jump process is independent from the Brownian motions. See [HestonJ](../api/sp/heston.rst#quantflow.sp.heston.HestonJ) for python
+API documentation.
 
 ```{code-cell} ipython3
 from quantflow.sp.heston import HestonJ
-pr = HestonJ.exponential(
+from quantflow.utils.distributions import DoubleExponential
+pr = HestonJ.create(
+    DoubleExponential,
     vol=0.6,
     kappa=2,
     sigma=0.8,
-    rho=-0.2,
+    rho=-0.0,
     jump_intensity=50,
     jump_fraction=0.2,
-    jump_asymmetry=1.2
+    jump_asymmetry=0.0
 )
 pr
 ```
 
 ```{code-cell} ipython3
 from quantflow.utils import plot
-plot.plot_marginal_pdf(pr.marginal(0.1), 128, normal=True, analytical=False)
+plot.plot_marginal_pdf(pr.marginal(0.5), 128, normal=True, analytical=False)
 ```
 
 ```{code-cell} ipython3
 from quantflow.options.pricer import OptionPricer
-from quantflow.sp.heston import Heston
 pricer = OptionPricer(pr)
 pricer
 ```
 
 ```{code-cell} ipython3
+
 fig = None
 for ttm in (0.05, 0.1, 0.2, 0.4, 0.6, 1):
+
     fig = pricer.maturity(ttm).plot(fig=fig, name=f"t={ttm}")
 fig.update_layout(title="Implied black vols", height=500)
 ```
@@ -68,3 +72,7 @@ fig.update_layout(title="Implied black vols", height=500)
 ```{code-cell} ipython3
 
 ```
+
+```{code-cell} ipython3
+
+```
diff --git a/notebooks/models/jump_diffusion.md b/notebooks/models/jump_diffusion.md
index 23dd3fd..7ac0ba1 100644
--- a/notebooks/models/jump_diffusion.md
+++ b/notebooks/models/jump_diffusion.md
@@ -6,7 +6,7 @@ jupytext:
     format_version: 0.13
     jupytext_version: 1.16.6
 kernelspec:
-  display_name: Python 3 (ipykernel)
+  display_name: .venv
   language: python
   name: python3
 ---
@@ -17,21 +17,26 @@ The library allows to create a vast array of jump-diffusion models. The most fam
 
 ## Merton Model
 
-```{code-cell}
-from quantflow.sp.jump_diffusion import Merton
+```{code-cell} ipython3
+from quantflow.sp.jump_diffusion import JumpDiffusion
+from quantflow.utils.distributions import Normal
 
-pr = Merton.create(diffusion_percentage=0.2, jump_intensity=50, jump_skew=-0.5)
-pr
+merton = JumpDiffusion.create(Normal, jump_fraction=0.8, jump_intensity=50)
 ```
 
 ### Marginal Distribution
 
-```{code-cell}
-m = pr.marginal(0.02)
+```{code-cell} ipython3
+m = merton.marginal(0.02)
 m.std(), m.std_from_characteristic()
 ```
 
-```{code-cell}
+```{code-cell} ipython3
+m2 = jd.marginal(0.02)
+m2.std(), m2.std_from_characteristic()
+```
+
+```{code-cell} ipython3
 from quantflow.utils import plot
 
 plot.plot_marginal_pdf(m, 128, normal=True, analytical=False, log_y=True)
@@ -39,7 +44,7 @@ plot.plot_marginal_pdf(m, 128, normal=True, analytical=False, log_y=True)
 
 ### Characteristic Function
 
-```{code-cell}
+```{code-cell} ipython3
 plot.plot_characteristic(m)
 ```
 
@@ -47,17 +52,17 @@ plot.plot_characteristic(m)
 
 We can price options using the `OptionPricer` tooling.
 
-```{code-cell}
+```{code-cell} ipython3
 from quantflow.options.pricer import OptionPricer
-pricer = OptionPricer(pr)
+pricer = OptionPricer(merton)
 pricer
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 fig = None
 for ttm in (0.05, 0.1, 0.2, 0.4, 0.6, 1):
     fig = pricer.maturity(ttm).plot(fig=fig, name=f"t={ttm}")
-fig.update_layout(title="Implied black vols", height=500)
+fig.update_layout(title="Implied black vols - Merton", height=500)
 ```
 
 This term structure of volatility demostrates one of the principal weakness of the Merton's model, and indeed of all jump diffusion models based on Lévy processes, namely the rapid flattening of the volatility surface as time-to-maturity increases.
@@ -67,14 +72,30 @@ For very short time-to-maturities, however, the model has no problem in producin
 
 ### MC paths
 
-```{code-cell}
-pr.sample(20, time_horizon=1, time_steps=1000).plot().update_traces(line_width=0.5)
+```{code-cell} ipython3
+merton.sample(20, time_horizon=1, time_steps=1000).plot().update_traces(line_width=0.5)
 ```
 
 ## Exponential Jump Diffusion
 
-This is a variation of the Mertoin model, where the jump distribution is a double exponential, one for the negative jumps and one for the positive jumps.
+This is a variation of the Mertoin model, where the jump distribution is a double exponential.
+The advantage of this model is that it allows for an asymmetric jump distribution, which can be useful in some cases, for example options prices with a skew.
+
+```{code-cell} ipython3
+from quantflow.utils.distributions import DoubleExponential
+
+jd = JumpDiffusion.create(DoubleExponential, jump_fraction=0.8, jump_intensity=50, jump_asymmetry=0.2)
+pricer = OptionPricer(jd)
+pricer
+```
+
+```{code-cell} ipython3
+fig = None
+for ttm in (0.05, 0.1, 0.2, 0.4, 0.6, 1):
+    fig = pricer.maturity(ttm).plot(fig=fig, name=f"t={ttm}")
+fig.update_layout(title="Implied black vols - Double-exponential Jump Diffusion ", height=500)
+```
+
+```{code-cell} ipython3
 
-```{code-cell}
-from 
 ```
diff --git a/notebooks/models/ou.md b/notebooks/models/ou.md
index 3a2ce25..e0aafd1 100644
--- a/notebooks/models/ou.md
+++ b/notebooks/models/ou.md
@@ -6,7 +6,7 @@ jupytext:
     format_version: 0.13
     jupytext_version: 1.16.6
 kernelspec:
-  display_name: Python 3 (ipykernel)
+  display_name: .venv
   language: python
   name: python3
 ---
@@ -19,7 +19,7 @@ The general definition of an Ornstein-Uhlebeck (OU) process is as the solution o
     d x_t = -\kappa x_t dt + d z_t
 \end{equation}
 
-where $z$, with $z_0 = 0$, is a [Lévy](./levy.md) process. As $z_t$ drives the OU process, it is usually referred to as a background driving Lévy process (**BDLP**).
+where $z$, with $z_0 = 0$, is a [Lévy](../theory/levy.md) process. As $z_t$ drives the OU process, it is usually referred to as a background driving Lévy process (**BDLP**).
 
 The OU process can be integrated into the formula (see Appendix below).
 
@@ -57,13 +57,13 @@ which means the process admits a stationary probability distribution equal to
     x_t \sim N\left(\theta, \frac{\sigma^2}{2\kappa}\right)\ \ t\rightarrow\infty
 \end{equation}
 
-```{code-cell}
+```{code-cell} ipython3
 from quantflow.sp.ou import Vasicek
-pr = Vasicek()
+pr = Vasicek(kappa=2)
 pr
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 pr.sample(20, time_horizon=1, time_steps=1000).plot().update_traces(
     line_width=0.5
 ).update_layout(
@@ -71,12 +71,12 @@ pr.sample(20, time_horizon=1, time_steps=1000).plot().update_traces(
 )
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 m = pr.marginal(1)
 m.mean(), m.std()
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 m.mean_from_characteristic(), m.std_from_characteristic()
 ```
 
@@ -105,7 +105,7 @@ where $m_n$ are the jump times of the Poisson process $N_{\kappa t} and $j_n$ ar
 
 ### Integrated Intensity
 
-One of the advantages of these OU processes is that they offer a great deal of analytical tractability. For example, the integrated value of the process, which can be used as a time change for [Lévy processes](./levy.md), is given by
+One of the advantages of these OU processes is that they offer a great deal of analytical tractability. For example, the integrated value of the process, which can be used as a time change for [Lévy processes](../theory/levy.md), is given by
 
 \begin{align}
    \int_0^t x_s ds &= \epsilon_t x_0 + \int_0^t \epsilon_{t-s} d z_{\kappa s} = \frac{z_{\kappa t} - x_t + x_0}{\kappa}\\
@@ -125,7 +125,7 @@ The library provides an implementation of the non-gaussian OU process in the for
 
 In this case, the BDLP is an exponential compound Poisson process with Lévy density $\lambda\beta e^{-\beta x}$, in other words, the [exponential compound Poisson](./poisson.md) process with intensity $\lambda$ and decay $\beta$.
 
-```{code-cell}
+```{code-cell} ipython3
 from quantflow.sp.ou import GammaOU
 
 pr = GammaOU.create(decay=10, kappa=5)
@@ -140,11 +140,11 @@ The charatecristic exponent of the $\Gamma$-OU process is given by, see {cite:p}
     \phi_{u, t} = -x_{0} i u e^{-\kappa t} - \lambda\ln\left(\frac{\beta-iue^{-\kappa t}}{\beta -iu}\right)
 \end{equation}
 
-```{code-cell}
+```{code-cell} ipython3
 pr.marginal(1).mean(), pr.marginal(1).std()
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 import numpy as np
 from quantflow.utils import plot
 
@@ -152,14 +152,14 @@ m = pr.marginal(5)
 plot.plot_marginal_pdf(m, 128)
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 from quantflow.utils.plot import plot_characteristic
 plot_characteristic(m)
 ```
 
 ### Sampling Gamma OU
 
-```{code-cell}
+```{code-cell} ipython3
 from quantflow.sp.ou import GammaOU
 pr = GammaOU.create(decay=10, kappa=5)
 
@@ -170,7 +170,7 @@ pr.sample(50, time_horizon=1, time_steps=1000).plot().update_traces(line_width=0
 
 Test the simulated meand and stadard deviation against the values from the invariant gamma distribution.
 
-```{code-cell}
+```{code-cell} ipython3
 import pandas as pd
 from quantflow.utils import plot
 
@@ -180,7 +180,7 @@ df = pd.DataFrame(mean, index=paths.time)
 plot.plot_lines(df)
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 std = dict(std=pr.marginal(paths.time).std(), simulated=paths.std())
 df = pd.DataFrame(std, index=paths.time)
 plot.plot_lines(df)
@@ -198,6 +198,6 @@ The integration of the OU process can be achieved by multiplying both sides of t
     x_t &= x_0 e^{-\kappa t} + \int_0^t e^{-\kappa\left(t - s\right)} d z_s
 \end{align}
 
-```{code-cell}
+```{code-cell} ipython3
 
 ```
diff --git a/notebooks/reference/glossary.md b/notebooks/reference/glossary.md
index b1d4839..f7ad50b 100644
--- a/notebooks/reference/glossary.md
+++ b/notebooks/reference/glossary.md
@@ -33,7 +33,7 @@ of a real-valued random variable $X$ is the function given by
 
 The Hurst exponent is a measure of the long-term memory of time series. The Hurst exponent is a measure of the relative tendency of a time series either to regress strongly to the mean or to cluster in a direction.
 
-Check this study on the [Hurst exponent with OHLC data](./applications/hurst).
+Check this study on the [Hurst exponent with OHLC data](../applications/hurst).
 
 ## Moneyness
 
diff --git a/notebooks/theory/levy.md b/notebooks/theory/levy.md
index a45664a..7ad9e2d 100644
--- a/notebooks/theory/levy.md
+++ b/notebooks/theory/levy.md
@@ -36,7 +36,7 @@ The independence and stationarity of the increments of the Lévy process imply t
 
 where the **characteristic exponent** $\phi_{x_1,u}$ is given by the [Lévy–Khintchine formula](https://en.wikipedia.org/wiki/L%C3%A9vy_process).
 
-There are several Lévy processes in the literature, including, importantly, the [Poisson process](./poisson.md), the compound Poisson process, and the Brownian motion.
+There are several Lévy processes in the literature, including, importantly, the [Poisson process](../models/poisson.md), the compound Poisson process, and the Brownian motion.
 
 +++
 
diff --git a/quantflow/sp/heston.py b/quantflow/sp/heston.py
index b8a944b..bcbd6e3 100644
--- a/quantflow/sp/heston.py
+++ b/quantflow/sp/heston.py
@@ -6,11 +6,11 @@
 from pydantic import Field
 
 from quantflow.ta.paths import Paths
-from quantflow.utils.distributions import DoubleExponential
 from quantflow.utils.types import FloatArrayLike, Vector
 
 from .base import StochasticProcess1D
 from .cir import CIR
+from .jump_diffusion import JumpDiffusion
 from .poisson import CompoundPoissonProcess, D
 
 
@@ -112,24 +112,31 @@ def sample_from_draws(self, path1: Paths, *args: Paths) -> Paths:
 
 
 class HestonJ(Heston, Generic[D]):
-    r"""The Heston stochastic volatility model with jumps
+    r"""A generic Heston stochastic volatility model with jumps
 
     The Heston model with jumps is an extension of the classical square-root
     stochastic volatility model of Heston (1993) with the addition of jump
-    processes. The jumps are modeled as compound Poisson processes
+    processes. The jumps are modeled via a compound Poisson process
 
     .. math::
         d x_t &= d w^1_t + d N_t\\
         d v_t &= \kappa (\theta - v_t) dt + \nu \sqrt{v_t} dw^2_t \\
         \rho dt &= {\tt E}[dw^1 dw^2]
+
+    This model is generic and therefore allows for different types of jumps
+    distributions **D**.
+
+    The Bates model is obtained by using the
+    :class:`.Normal` distribution for the jump sizes.
     """
 
     jumps: CompoundPoissonProcess[D] = Field(description="jumps process")
-    """Jump process driven by a compound Poisson process"""
+    """Jump process driven by a :class:`.CompoundPoissonProcess`"""
 
     @classmethod
-    def exponential(
+    def create(  # type: ignore [override]
         cls,
+        jump_distribution: type[D],
         *,
         rate: float = 1.0,
         vol: float = 0.5,
@@ -139,8 +146,8 @@ def exponential(
         theta: float | None = None,
         jump_intensity: float = 100,  # number of jumps per year
         jump_fraction: float = 0.1,  # percentage of variance due to jumps
-        jump_asymmetry: float = 1,
-    ) -> HestonJ[DoubleExponential]:
+        jump_asymmetry: float = 0.0,
+    ) -> HestonJ[D]:
         r"""Create an Heston model with :class:`.DoubleExponential` jumps.
 
         To understand the parameters lets introduce the following notation:
@@ -163,22 +170,23 @@ def exponential(
             variance and price processes
         :param theta: The long-term mean of the variance process, if `None`, it
             defaults to the diffusion variance given by :math:`{\tt var}_d`
-        :param jump_intensity: The number of jumps per year
-        :param jump_fraction: The percentage of variance due to jumps
-        :param jump_asymmetry: The asymmetry of the jump distribution (1 for symmetric)
-        :param jump_distribution: The distribution of the jumps, either 'normal'
-            or 'double_exponential' (default)
+        :param jump_intensity: The average number of jumps per year
+        :param jump_fraction: The fraction of variance due to jumps (between 0 and 1)
+        :param jump_asymmetry: The asymmetry of the jump distribution
+            (0 for symmetric jumps)
+        :param jump_distribution: The distribution of the jumps
+            (normal distribution for the Merton model)
         """
-        if jump_fraction <= 0 or jump_fraction >= 1:
-            raise ValueError("jump_percentage must be between 0 and 1")
-        total_variance = vol * vol
-        jump_variance = total_variance * jump_fraction
-        diffusion_variance = total_variance - jump_variance
-        jump_distribution_variance = jump_variance / jump_intensity
-        jumps = DoubleExponential.from_moments(
-            variance=jump_distribution_variance, kappa=jump_asymmetry
+        jd = JumpDiffusion.create(
+            jump_distribution,
+            vol=vol,
+            jump_intensity=jump_intensity,
+            jump_fraction=jump_fraction,
+            jump_asymmetry=jump_asymmetry,
         )
-        return HestonJ[DoubleExponential](
+        total_variance = vol * vol
+        diffusion_variance = total_variance * (1 - jump_fraction)
+        return cls(
             variance_process=CIR(
                 rate=rate * diffusion_variance,
                 kappa=kappa,
@@ -186,10 +194,7 @@ def exponential(
                 theta=theta if theta is not None else diffusion_variance,
             ),
             rho=rho,
-            jumps=CompoundPoissonProcess(
-                intensity=jump_intensity,
-                jumps=jumps,
-            ),
+            jumps=jd.jumps,
         )
 
     def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector:
diff --git a/quantflow/sp/jump_diffusion.py b/quantflow/sp/jump_diffusion.py
index d9b84ba..3a18f10 100644
--- a/quantflow/sp/jump_diffusion.py
+++ b/quantflow/sp/jump_diffusion.py
@@ -1,12 +1,11 @@
 from __future__ import annotations
 
-from typing import Generic, Self
+from typing import Generic
 
 import numpy as np
 from pydantic import Field
 
 from ..ta.paths import Paths
-from ..utils.distributions import Normal
 from ..utils.types import FloatArrayLike, Vector
 from .base import StochasticProcess1D
 from .poisson import CompoundPoissonProcess, D
@@ -52,35 +51,39 @@ def analytical_mean(self, t: FloatArrayLike) -> FloatArrayLike:
     def analytical_variance(self, t: FloatArrayLike) -> FloatArrayLike:
         return self.diffusion.analytical_variance(t) + self.jumps.analytical_variance(t)
 
-
-class Merton(JumpDiffusion[Normal]):
-    """Merton jump-diffusion model"""
-
     @classmethod
     def create(
         cls,
+        jump_distribution: type[D],
         vol: float = 0.5,
-        diffusion_percentage: float = 0.5,
         jump_intensity: float = 100,
-        jump_skew: float = 0.0,
-    ) -> Self:
-        """Create a Merton jump-diffusion model with a given volatility,
-        diffusion percentage, jump skewness, and jump intensity"""
+        jump_fraction: float = 0.5,
+        jump_asymmetry: float = 0.0,
+    ) -> JumpDiffusion[D]:
+        """Create a jump-diffusion model with a given jump distribution, volatility
+        and jump fraction.
+
+        :param vol: total annualized standard deviation
+        :param jump_intensity: The average number of jumps per year
+        :param jump_fraction: The fraction of variance due to jumps (between 0 and 1)
+        :param jump_asymmetry: The asymmetry of the jump distribution (0 for symmetric,
+            only used by distributions with asymmetry)
+
+        If the jump distribution is set to the :class:`.Normal` distribution, the
+        model reduces to a Merton jump-diffusion model.
+        """
         variance = vol * vol
-        jump_std = 1.0
-        jump_mean = 0.0
-        if diffusion_percentage > 1:
-            raise ValueError("diffusion_percentage must be less than 1")
-        elif diffusion_percentage < 0:
-            raise ValueError("diffusion_percentage must be greater than 0")
-        elif diffusion_percentage == 1:
-            jump_intensity = 0
+        if jump_fraction >= 1:
+            raise ValueError("jump_fraction must be less than 1")
+        elif jump_fraction <= 0:
+            raise ValueError("jump_fraction must be greater than 0")
         else:
-            jump_std = np.sqrt(variance * (1 - diffusion_percentage) / jump_intensity)
-            jump_mean = jump_skew / jump_intensity
-        return cls(
-            diffusion=WeinerProcess(sigma=np.sqrt(variance * diffusion_percentage)),
-            jumps=CompoundPoissonProcess[Normal](
-                intensity=jump_intensity, jumps=Normal(mu=jump_mean, sigma=jump_std)
-            ),
-        )
+            jump_variance = variance * jump_fraction
+            jump_distribution_variance = jump_variance / jump_intensity
+            jumps = jump_distribution.from_variance_and_asymmetry(
+                jump_distribution_variance, jump_asymmetry
+            )
+            return cls(
+                diffusion=WeinerProcess(sigma=np.sqrt(variance * (1 - jump_fraction))),
+                jumps=CompoundPoissonProcess(intensity=jump_intensity, jumps=jumps),
+            )
diff --git a/quantflow/utils/distributions.py b/quantflow/utils/distributions.py
index 94ab94d..5b9dd2a 100644
--- a/quantflow/utils/distributions.py
+++ b/quantflow/utils/distributions.py
@@ -18,6 +18,11 @@ class Distribution1D(Marginal1D):
     def sample(self, n: int) -> np.ndarray:
         """Sample from the distribution"""
 
+    @classmethod
+    def from_variance_and_asymmetry(cls, variance: float, asymmetry: float) -> Self:
+        """Create a distribution from variance and asymmetry"""
+        raise NotImplementedError
+
 
 class Exponential(Distribution1D):
     r"""A :class:`.Marginal1D` for the `Exponential distribution`_
@@ -76,6 +81,11 @@ class Normal(Distribution1D):
     mu: float = Field(default=0, description="mean")
     sigma: float = Field(default=1, gt=0, description="standard deviation")
 
+    @classmethod
+    def from_variance_and_asymmetry(cls, variance: float, asymmetry: float) -> Self:
+        """The normal distribution is symmetric, so the asymmetry is ignored"""
+        return cls(mu=0, sigma=np.sqrt(variance))
+
     @property
     def sigma2(self) -> float:
         return self.sigma**2
@@ -135,6 +145,10 @@ def log_kappa(self) -> float:
         """The log of the :attr:`.kappa` parameter"""
         return np.log(self.kappa)
 
+    @classmethod
+    def from_variance_and_asymmetry(cls, variance: float, asymmetry: float) -> Self:
+        return cls.from_moments(variance=variance, kappa=np.exp(asymmetry))
+
     @classmethod
     def from_moments(
         cls,
diff --git a/quantflow_tests/test_heston.py b/quantflow_tests/test_heston.py
index a477407..74c7997 100644
--- a/quantflow_tests/test_heston.py
+++ b/quantflow_tests/test_heston.py
@@ -12,7 +12,8 @@ def heston() -> Heston:
 
 @pytest.fixture
 def heston_jumps() -> HestonJ[DoubleExponential]:
-    return HestonJ.exponential(
+    return HestonJ.create(
+        DoubleExponential,
         vol=0.5,
         kappa=1,
         sigma=0.5,
diff --git a/quantflow_tests/test_jump_diffusion.py b/quantflow_tests/test_jump_diffusion.py
index b39589a..c6fb00c 100644
--- a/quantflow_tests/test_jump_diffusion.py
+++ b/quantflow_tests/test_jump_diffusion.py
@@ -1,16 +1,17 @@
 import pytest
 
-from quantflow.sp.jump_diffusion import Merton
+from quantflow.sp.jump_diffusion import JumpDiffusion
+from quantflow.utils.distributions import Normal
 
 
 @pytest.fixture
-def merton() -> Merton:
-    return Merton.create(diffusion_percentage=0.2, jump_skew=-0.1)
+def merton() -> JumpDiffusion[Normal]:
+    return JumpDiffusion.create(Normal, jump_fraction=0.8)
 
 
-def test_characteristic(merton: Merton) -> None:
+def test_characteristic(merton: JumpDiffusion[Normal]) -> None:
     m = merton.marginal(1)
-    assert m.mean() < 0
+    assert m.mean() == 0
     assert pytest.approx(m.std()) == pytest.approx(0.5, 1.0e-3)
     pdf = m.pdf_from_characteristic(128)
     assert pdf.x[0] < 0
@@ -18,7 +19,7 @@ def test_characteristic(merton: Merton) -> None:
     assert -pdf.x[0] != pdf.x[-1]
 
 
-def test_sampling(merton: Merton) -> None:
+def test_sampling(merton: JumpDiffusion[Normal]) -> None:
     paths = merton.sample(1000, time_horizon=1, time_steps=1000)
     mean = paths.mean()
     assert mean[0] == 0
diff --git a/quantflow_tests/test_options_pricer.py b/quantflow_tests/test_options_pricer.py
index 0da7067..761044f 100644
--- a/quantflow_tests/test_options_pricer.py
+++ b/quantflow_tests/test_options_pricer.py
@@ -8,7 +8,9 @@
 
 @pytest.fixture
 def pricer() -> OptionPricer[HestonJ[DoubleExponential]]:
-    return OptionPricer(HestonJ.exponential(vol=0.5, kappa=1, sigma=0.8, rho=0))
+    return OptionPricer(
+        HestonJ.create(DoubleExponential, vol=0.5, kappa=1, sigma=0.8, rho=0)
+    )
 
 
 @pytest.mark.skipif(not has_plotly, reason="Plotly not installed")