@@ -41,7 +41,9 @@ While $\tau_t$ is always continuous, $\lambda$ can exhibit jumps. Since the time
  \Phi_{y_t, u} = {\mathbb E}\left[e^{i u x_{\tau_t}}\right] = {\mathbb E}\left[{\mathbb E}\left[\left.e^{i u x_s}\right|\tau_t=s\right]\right]
 \end{equation}
 
-where the inside expectation is taken on $x_{\tau_t}$ conditional on a fixed value of $\tau_t = s$ and the outside expectation is on all possible values of $\tau_t$. If the random time $\tau_t$ is independent of $x_t$, the randomness due to the Lévy process can be integrated out using the characteristic function of $x_t$:
+where the inside expectation is taken on $x_{\tau_t}$ conditional on a fixed value of $\tau_t = s$ and the outside expectation is on all possible values of $\tau_t$.
+
+If the random time $\tau_t$ is independent of $x_t$, the randomness due to the Lévy process can be integrated out using the characteristic function of $x_t$:
 
 \begin{equation}
 \Phi_{y_t, u} = {\mathbb E}\left[e^{-\tau_t \phi_{x,u}}\right] = {\mathbb L}_{\tau_t}\left(\phi_{x_1,u}\right)
 
@@ -204,7 +204,7 @@ def integrated_log_laplace(
         r"""The log-Laplace transform of the cumulative process:
 
         \begin{equation}
-            e^{\phi_{t, u}} = {\mathbb E} \left[e^{i u \int_0^t x_s ds}\right]
+            \phi_{t, u} = \log {\mathbb E} \left[e^{-u \int_0^t x_s ds}\right]
         \end{equation}
         """
 
 
@@ -101,18 +101,29 @@ def sample(self, n: int, time_horizon: float = 1, time_steps: int = 100) -> Path
         return self.sample_from_draws(Paths.normal_draws(n, time_horizon, time_steps))
 
     def sample_from_draws(self, path_dw: Paths, *args: Paths) -> Paths:
+        kappa = self.variance_process.kappa
+        time_steps = path_dw.time_steps
+        time_horizon = path_dw.t
+        dt = path_dw.dt
+
         if args:
             path_dz = args[0]
         else:
-            # generate the background driving process samples if not provided
+            # BDLP runs at kappa-rescaled time, so sample over [0, kappa*T]
             path_dz = self.variance_process.bdlp.sample(
-                path_dw.samples, path_dw.t, path_dw.time_steps
+                path_dw.samples, kappa * time_horizon, time_steps
             )
-        dt = path_dw.dt
-        # sample the activity rate process
-        v = self.variance_process.sample_from_draws(path_dz)
-        # create the time-changed Brownian motion
-        dw = path_dw.data * np.sqrt(v.data * dt)
-        paths = np.zeros(dw.shape)
-        paths[1:] = np.cumsum(dw[:-1], axis=0) + path_dz.data
-        return Paths(t=path_dw.t, data=paths)
+
+        # Variance via the OU recursion v_{i+1} = v_i * e^{-kappa dt} + Delta z_{i+1}
+        decay = np.exp(-kappa * dt)
+        dz = np.diff(path_dz.data, axis=0)
+        v = np.zeros_like(path_dz.data)
+        v[0] = self.variance_process.rate
+        for i in range(time_steps):
+            v[i + 1] = v[i] * decay + dz[i]
+
+        # Price: x_t = integral sqrt(v) dW + rho * z_{kappa t}
+        diffusion = np.sqrt(v[:-1] * dt) * path_dw.data[:-1]
+        paths = np.zeros_like(path_dw.data)
+        paths[1:] = np.cumsum(diffusion, axis=0) + self.rho * path_dz.data[1:]
+        return Paths(t=time_horizon, data=paths)
@@ -34,8 +34,10 @@ def support(self, mean: float, std: float, points: int) -> FloatArray:
         return np.linspace(0, points, points + 1)
 
     def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector:
-        phi = self.poisson.characteristic_exponent(t, u)
-        return -self.intensity.integrated_log_laplace(t, phi)
+        # phi_x is the per-unit-time Lévy exponent of the inner Poisson;
+        # the integrated log-Laplace already absorbs the time horizon
+        phi_x = self.poisson.characteristic_exponent(1, u)
+        return -self.intensity.integrated_log_laplace(t, phi_x)
 
     def arrivals(self, t: float = 1) -> FloatArray:
         paths = self.intensity.sample(1, t, math.ceil(100 * t)).integrate()
 
@@ -315,6 +315,13 @@ def characteristic_exponent(
             t, u
         ) + self.jumps.characteristic_exponent(t, u)
 
+    def sample_from_draws(self, path1: Paths, *args: Paths) -> Paths:
+        diffusion = super().sample_from_draws(path1, *args)
+        jump_path = self.jumps.sample(
+            diffusion.samples, diffusion.t, diffusion.time_steps
+        )
+        return Paths(t=diffusion.t, data=diffusion.data + jump_path.data)
+
 
 class DoubleHeston(StochasticProcess1D):
     r"""Double Heston stochastic volatility model.
 
@@ -5,24 +5,25 @@
 
 import numpy as np
 from pydantic import Field
-from scipy.optimize import Bounds
 from scipy.stats import gamma, norm
 
 from ..ta.paths import Paths
 from ..utils.distributions import Exponential
 from ..utils.types import Float, FloatArrayLike, Vector
-from .base import Im, IntensityProcess
+from .base import IntensityProcess, StochasticProcess1D
 from .poisson import CompoundPoissonProcess, D
 from .weiner import WeinerProcess
 
 
-class Vasicek(IntensityProcess):
-    r"""Gaussian OU process, also know as the
-    [Vasiceck model](https://en.wikipedia.org/wiki/Vasicek_model)
+class Vasicek(StochasticProcess1D):
+    r"""Gaussian OU process, also known as the
+    [Vasicek model](https://en.wikipedia.org/wiki/Vasicek_model).
 
     Historically, the Vasicek model was used to model the short rate, but it can be
     used to model any process that reverts to a mean level at a rate proportional to
-    the difference between the current level and the mean level.
+    the difference between the current level and the mean level. Unlike intensity
+    processes, the Vasicek process can take negative values, so it is not constrained
+    to a positive domain.
 
     \begin{equation}
         dx_t = \kappa (\theta - x_t) dt + \sigma dw_t
@@ -37,11 +38,15 @@ class Vasicek(IntensityProcess):
     \end{equation}
     """
 
+    rate: float = Field(default=1.0, description=r"Initial value $x_0$")
+    kappa: float = Field(
+        default=1.0, gt=0, description=r"Mean reversion speed $\kappa$"
+    )
+    theta: float = Field(default=1.0, description=r"Mean level $\theta$")
     bdlp: WeinerProcess = Field(
         default_factory=WeinerProcess,
-        description="Background driving Weiner process",
+        description="Background driving Wiener process",
     )
-    theta: float = Field(default=1.0, gt=0, description=r"Mean rate $\theta$")
 
     def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector:
         r"""The characteristic exponent of the Vasicek process is given by
@@ -54,9 +59,6 @@ def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector:
         var = self.analytical_variance(t)
         return u * (-1j * mu + 0.5 * var * u)
 
-    def integrated_log_laplace(self, t: FloatArrayLike, u: Vector) -> Vector:
-        raise NotImplementedError
-
     def sample(self, n: int, time_horizon: float = 1, time_steps: int = 100) -> Paths:
         paths = Paths.normal_draws(n, time_horizon, time_steps)
         return self.sample_from_draws(paths)
@@ -74,8 +76,8 @@ def sample_from_draws(self, draws: Paths, *args: Paths) -> Paths:
             paths[t + 1, :] = x + dx
         return Paths(t=draws.t, data=paths)
 
-    def domain_range(self) -> Bounds:
-        return Bounds(-np.inf, np.inf)
+    def ekt(self, t: FloatArrayLike) -> FloatArrayLike:
+        return np.exp(-self.kappa * t)
 
     def analytical_mean(self, t: FloatArrayLike) -> FloatArrayLike:
         r"""The analytical mean of the Vasicek process is given by
@@ -149,10 +151,12 @@ def intensity(self) -> float:
 
     @abstractmethod
     def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector:
-        r"""
+        r"""Characteristic exponent of a non-Gaussian OU process driven by a
+        Lévy process $Z$ with characteristic exponent $\phi_Z$:
+
         \begin{equation}
-            \phi_{x_t,u} = iu x_0 e^{-\kappa t} +
-                \int_{ue^{-\kappa t}}^{u} \frac{\psi_Z(v)}{v} , dv
+            \phi_{x_t, u} = - iu\, x_0\, e^{-\kappa t}
+                + \int_{u e^{-\kappa t}}^{u} \frac{\phi_Z(v)}{v}\,dv
         \end{equation}
         """
 
@@ -180,30 +184,82 @@ def create(cls, rate: float = 1, decay: float = 1, kappa: float = 1) -> Self:
         )
 
     def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector:
-        r"""
+        r"""Closed form of the [NGOU][quantflow.sp.ou.NGOU] characteristic
+        exponent when the BDLP is a compound Poisson process with intensity
+        $\lambda$ and exponential jumps of decay $\beta$:
+
+        \begin{equation}
+            \phi_{x_t, u} =
+                - iu\, x_0\, e^{-\kappa t}
+                + \lambda \log\!\left(
+                    \frac{\beta - iu}{\beta - iu\, e^{-\kappa t}}
+                \right)
+        \end{equation}
+
+        Derivation. The characteristic function of the
+        [Exponential][quantflow.utils.distributions.Exponential.characteristic]
+        jumps is
+
         \begin{equation}
-            \phi_{x_t,u} =
-                - iu x_0 e^{-\kappa t}
-                + \lambda \log\!\left(\frac{\beta - iu}{\beta - iu e^{-\kappa t}}\right)
+            \Phi_J(v) = \frac{\beta}{\beta - iv},
         \end{equation}
+
+        so the BDLP unit-time characteristic exponent
+        ([CompoundPoissonProcess][quantflow.sp.poisson.CompoundPoissonProcess]
+        at $t = 1$) is
+
+        \begin{equation}
+            \phi_Z(v) = \lambda\bigl(1 - \Phi_J(v)\bigr)
+                = -\frac{i v\, \lambda}{\beta - iv},
+        \end{equation}
+
+        and $\phi_Z(v)/v = -i\lambda / (\beta - iv)$. Substituting
+        $w = \beta - iv$ (so $dv = i\, dw$) in the NGOU integral gives
+
+        \begin{equation}
+            \int_{u e^{-\kappa t}}^{u} \frac{\phi_Z(v)}{v}\, dv
+                = \lambda \int_{\beta - iu\, e^{-\kappa t}}^{\beta - iu}
+                    \frac{dw}{w}
+                = \lambda \log\!\left(
+                    \frac{\beta - iu}{\beta - iu\, e^{-\kappa t}}
+                \right),
+        \end{equation}
+
+        which combined with the drift term $-iu\, x_0\, e^{-\kappa t}$ yields
+        the closed form above.
         """
         b = self.beta
-        iu = Im * u
+        iu = 1j * u
         c1 = iu * np.exp(-self.kappa * t)
         c0 = self.intensity * np.log((b - c1) / (b - iu))
         return -c0 - c1 * self.rate
 
     def integrated_log_laplace(self, t: FloatArrayLike, u: Vector) -> Vector:
+        r"""Log-Laplace transform of the time-integrated Gamma-OU process.
+
+        \begin{equation}
+        \begin{aligned}
+            \phi(t, u) &= \log E\!\left[e^{-u \int_0^t x_s\, ds}\right] \\
+            &= -u\, x_0\, \frac{1 - e^{-\kappa t}}{\kappa}
+              + \frac{\lambda}{u + \beta\kappa}
+              \!\left[\beta\kappa\,
+                \log\!\frac{\beta\kappa + u(1 - e^{-\kappa t})}{\beta\kappa}
+                - u\kappa t\right]
+        \end{aligned}
+        \end{equation}
+
+        where $x_0$ is the initial rate, $\lambda$ is the BDLP intensity and
+        $\beta$ the exponential jump decay.
+        """
         kappa = self.kappa
         b = self.beta
-        iu = Im * u
-        iuk = iu / kappa
+        muk = -u / kappa
         ekt = np.exp(-kappa * t)
-        c1 = iuk * (1 - ekt)
+        c1 = muk * (1 - ekt)
         c0 = self.intensity * (
-            b * np.log(b / (iuk + (b - iuk) / ekt)) / (iuk - b) - kappa * t
+            b * np.log(b / (muk + (b - muk) / ekt)) / (muk - b) - kappa * t
         )
-        return np.exp(c0 + c1 * self.rate)
+        return c0 + c1 * self.rate
 
     def sample(self, n: int, time_horizon: float = 1, time_steps: int = 100) -> Paths:
         dt = time_horizon / time_steps
@@ -235,28 +291,43 @@ def _advance(
     ) -> int:
         x = pp[i - 1]
         kappa = self.kappa
-        t0 = i * dt
-        t1 = t0 + dt
+        # step i fills pp[i] from pp[i-1] over the interval [(i-1) dt, i dt]
+        t0 = (i - 1) * dt
+        t1 = i * dt
         a = arrival or t1
         pp[i] = x - kappa * x * (a - t0) - kappa * (x + jump) * (t1 - a) + jump
         return i + 1
 
-    def cumulative_characteristic2(self, t: FloatArrayLike, u: Vector) -> Vector:
-        """Formula from a paper"""
-        kappa = self.kappa
-        b = self.beta
-        iu = Im * u
-        iuk = iu / kappa
-        ekt = np.exp(-kappa * t)
-        c1 = iuk * (1 - ekt)
-        c0 = self.intensity * (b * np.log(b / (b - c1)) - iu * t) / (iuk - b)
-        return np.exp(c0 + c1 * self.rate)
-
     def analytical_mean(self, t: FloatArrayLike) -> FloatArrayLike:
-        return self.intensity / self.beta
+        r"""Analytical mean of the Gamma-OU process at time $t$.
+
+        \begin{equation}
+            E[x_t] = x_0\, e^{-\kappa t}
+                   + \frac{\lambda}{\beta}\bigl(1 - e^{-\kappa t}\bigr)
+        \end{equation}
+
+        which converges to the stationary mean $\lambda / \beta$ as
+        $t \to \infty$.
+        """
+        ekt = self.ekt(t)
+        return self.rate * ekt + (self.intensity / self.beta) * (1 - ekt)
 
     def analytical_variance(self, t: FloatArrayLike) -> FloatArrayLike:
-        return self.intensity / self.beta / self.beta
+        r"""Analytical variance of the Gamma-OU process at time $t$.
+
+        \begin{equation}
+            \mathrm{Var}(x_t) = \frac{\lambda}{\beta^2}
+                \bigl(1 - e^{-2\kappa t}\bigr)
+        \end{equation}
+
+        which converges to the stationary variance $\lambda / \beta^2$ as
+        $t \to \infty$.
+        """
+        return self.intensity * (1 - self.ekt(2 * t)) / (self.beta * self.beta)
 
     def analytical_pdf(self, t: FloatArrayLike, x: FloatArrayLike) -> FloatArrayLike:
+        """Stationary marginal density of the Gamma-OU process. The transient
+        density is not available in closed form; use `pdf_from_characteristic`
+        for finite $t$.
+        """
         return gamma.pdf(x, self.intensity, scale=1 / self.beta)
@@ -32,10 +32,10 @@ def test_bns_horizon(bns: BNS) -> None:
 
 def test_bns_no_leverage_matches_integrated_laplace(bns: BNS) -> None:
     # at rho=0: E[exp(iu x_t)] = E[exp(-u^2/2 * int v_s ds)]
-    # which is the GammaOU integrated Laplace transform at i*u^2/2
-    u = np.array([0.1, 0.5, 1.0, 2.0, 5.0], dtype=complex)
+    # which is the GammaOU integrated Laplace transform evaluated at u^2/2
+    u = np.array([0.1, 0.5, 1.0, 2.0, 5.0])
     cf_bns = bns.characteristic(1, u)
-    cf_via_ou = bns.variance_process.integrated_log_laplace(1, 1j * u * u / 2)
+    cf_via_ou = np.exp(bns.variance_process.integrated_log_laplace(1, u * u / 2))
     np.testing.assert_allclose(cf_bns, cf_via_ou, atol=1e-12)
 
 
@@ -45,3 +45,16 @@ def test_bns_leverage(bns_leverage: BNS) -> None:
     assert bns_leverage.characteristic(1, 0) == 1
     # E[x_t] = rho * kappa * t * intensity / beta = -0.5 * 1 * 1 * 1.25 / 5
     assert m.mean_from_characteristic() == pytest.approx(-0.125, 1e-3)
+
+
+def test_bns_sample_moments(bns_leverage: BNS) -> None:
+    np.random.seed(42)
+    paths = bns_leverage.sample(5000, time_horizon=1.0, time_steps=200)
+    assert paths.data.shape == (201, 5000)
+    # Compare empirical moments at T=1 against the characteristic-function moments
+    m = bns_leverage.marginal(1)
+    terminal = paths.data[-1]
+    assert terminal.mean() == pytest.approx(
+        float(m.mean_from_characteristic()), abs=0.02
+    )
+    assert terminal.std() == pytest.approx(float(m.std_from_characteristic()), abs=0.02)
@@ -1,3 +1,4 @@
+import numpy as np
 import pytest
 
 from quantflow.sp.heston import DoubleHeston, DoubleHestonJ, Heston, HestonJ
@@ -83,3 +84,21 @@ def test_double_heston_jumps_characteristic(
     characteristic_tests(m)
     assert m.mean() == pytest.approx(0.0, abs=1e-6)
     assert m.std() == pytest.approx(0.5, rel=0.05)
+
+
+def test_heston_jumps_sample_includes_jumps(
+    heston_jumps: HestonJ[DoubleExponential],
+) -> None:
+    np.random.seed(42)
+    paths = heston_jumps.sample(5000, time_horizon=1.0, time_steps=200)
+    # If jumps were dropped, std would collapse to sqrt(1 - jump_fraction) * 0.5
+    # ~ 0.418, well below 0.5
+    assert paths.data[-1].std() == pytest.approx(0.5, abs=0.03)
+
+
+def test_double_heston_jumps_sample_includes_jumps(
+    double_heston_jumps: DoubleHestonJ[DoubleExponential],
+) -> None:
+    np.random.seed(42)
+    paths = double_heston_jumps.sample(5000, time_horizon=1.0, time_steps=200)
+    assert paths.data[-1].std() == pytest.approx(0.5, abs=0.03)