diff --git a/docs/api/options/calibration.md b/docs/api/options/calibration.md
index 49fb5865..f385e66c 100644
--- a/docs/api/options/calibration.md
+++ b/docs/api/options/calibration.md
@@ -1,9 +1,15 @@
 # Vol Model Calibration
 
-::: quantflow.options.calibration.OptionEntry
+::: quantflow.options.calibration.base.OptionEntry
 
-::: quantflow.options.calibration.VolModelCalibration
+::: quantflow.options.calibration.base.VolModelCalibration
 
-::: quantflow.options.heston_calibration.HestonCalibration
+::: quantflow.options.calibration.heston.HestonCalibration
 
-::: quantflow.options.heston_calibration.HestonJCalibration
+::: quantflow.options.calibration.heston.HestonJCalibration
+
+::: quantflow.options.calibration.heston.DoubleHestonCalibration
+
+::: quantflow.options.calibration.heston.DoubleHestonJCalibration
+
+::: quantflow.options.calibration.bns.BNSCalibration
diff --git a/docs/examples/vol_surface_bns_calibration.py b/docs/examples/vol_surface_bns_calibration.py
new file mode 100644
index 00000000..6774bc48
--- /dev/null
+++ b/docs/examples/vol_surface_bns_calibration.py
@@ -0,0 +1,34 @@
+import json
+
+from docs.examples._utils import assets_path, print_model
+from quantflow.options.calibration import BNSCalibration
+from quantflow.options.pricer import OptionPricer
+from quantflow.options.surface import VolSurface, VolSurfaceInputs, surface_from_inputs
+from quantflow.sp.bns import BNS
+
+# Load a saved volatility surface snapshot and build the surface
+with open("docs/examples/volsurface.json") as fp:
+    surface: VolSurface = surface_from_inputs(VolSurfaceInputs(**json.load(fp)))
+
+surface.bs()
+surface.disable_outliers()
+
+# Create a BNS pricer with initial parameters
+pricer = OptionPricer(
+    model=BNS.create(vol=0.5, kappa=1.0, decay=10.0, rho=-0.2),
+)
+
+calibration: BNSCalibration[BNS] = BNSCalibration(
+    pricer=pricer,
+    vol_surface=surface,
+    moneyness_weight=0.5,
+)
+
+result = calibration.fit()
+print(result.message)
+print_model(calibration.model)
+
+# Plot the calibrated smile for all maturities and save as PNG
+fig = calibration.plot_maturities(max_moneyness=1.5, support=101)
+fig.update_layout(title="BNS Calibrated Smiles")
+fig.write_image(assets_path("bns_calibrated_smile.png"), width=1200)
diff --git a/docs/examples/vol_surface_heston_calibration.py b/docs/examples/vol_surface_heston_calibration.py
index 6a7c5cda..08b77d69 100644
--- a/docs/examples/vol_surface_heston_calibration.py
+++ b/docs/examples/vol_surface_heston_calibration.py
@@ -1,7 +1,7 @@
 import json
 
 from docs.examples._utils import assets_path, print_model
-from quantflow.options.heston_calibration import HestonCalibration
+from quantflow.options.calibration import HestonCalibration
 from quantflow.options.pricer import OptionPricer, OptionPricingMethod
 from quantflow.options.surface import VolSurfaceInputs, surface_from_inputs
 from quantflow.sp.heston import Heston
diff --git a/docs/examples/vol_surface_hestonj_calibration.py b/docs/examples/vol_surface_hestonj_calibration.py
index a6083832..94f0da2e 100644
--- a/docs/examples/vol_surface_hestonj_calibration.py
+++ b/docs/examples/vol_surface_hestonj_calibration.py
@@ -1,7 +1,7 @@
 import json
 
 from docs.examples._utils import assets_path, print_model
-from quantflow.options.heston_calibration import HestonJCalibration
+from quantflow.options.calibration import HestonJCalibration
 from quantflow.options.pricer import OptionPricer
 from quantflow.options.surface import VolSurface, VolSurfaceInputs, surface_from_inputs
 from quantflow.sp.heston import HestonJ
diff --git a/docs/glossary.md b/docs/glossary.md
index e53ff156..2a648e88 100644
--- a/docs/glossary.md
+++ b/docs/glossary.md
@@ -67,7 +67,7 @@ In the [Heston model][quantflow.sp.heston.Heston] the variance process $v_t$ is
 process, so the same condition applies with $\sigma$ being the vol of vol. The
 [CIR.is_positive][quantflow.sp.cir.CIR.is_positive] property checks whether the
 condition holds. The
-[HestonCalibration][quantflow.options.heston_calibration.HestonCalibration] class provides a
+[HestonCalibration][quantflow.options.calibration.heston.HestonCalibration] class provides a
 `feller_enforce` flag (default `True`) that imposes this as a hard inequality constraint
 during optimisation.
 
diff --git a/docs/javascripts/mathjax.js b/docs/javascripts/mathjax.js
index ace66202..1bad4ffd 100644
--- a/docs/javascripts/mathjax.js
+++ b/docs/javascripts/mathjax.js
@@ -9,11 +9,18 @@ window.MathJax = {
   options: {
     ignoreHtmlClass: ".*|",
     processHtmlClass: "arithmatex"
+  },
+  startup: {
+    typeset: false
   }
 };
 
 document$.subscribe(() => {
-  MathJax.startup.promise.then(() => {
-    MathJax.typesetPromise();
+  if (typeof MathJax === "undefined" || !MathJax.startup) return;
+  MathJax.startup.promise = MathJax.startup.promise.then(() => {
+    MathJax.startup.output.clearCache();
+    MathJax.typesetClear();
+    MathJax.texReset();
+    return MathJax.typesetPromise();
   });
 });
diff --git a/docs/tutorials/bns_calibration.md b/docs/tutorials/bns_calibration.md
new file mode 100644
index 00000000..8e4af104
--- /dev/null
+++ b/docs/tutorials/bns_calibration.md
@@ -0,0 +1,87 @@
+# BNS Volatility Model
+
+This tutorial calibrates the [BNS][quantflow.sp.bns.BNS] stochastic-volatility
+model (Barndorff-Nielsen and Shephard) to an implied volatility surface, using
+the same workflow as the Heston tutorial in
+[Volatility Surface](volatility_surface.md).
+
+BNS is structurally different from Heston. The variance process is a
+non-Gaussian Ornstein-Uhlenbeck process driven by a pure-jump Lévy process
+(Gamma-OU in this implementation), and the leverage effect is introduced by
+correlating the same jumps into the log-price.
+
+## Model Parameters
+
+[BNSCalibration][quantflow.options.calibration.bns.BNSCalibration] fits five
+parameters to the surface:
+
+| Parameter | Description |
+|---|---|
+| `v0` | Initial variance ($v_0$) |
+| `theta` | Long-run variance ($\theta = \lambda / \beta$) |
+| `kappa` | Mean reversion speed of the variance process |
+| `beta` | Exponential decay rate of the BDLP jump-size distribution |
+| `rho` | Leverage parameter (correlation between jumps in variance and log-price) |
+
+The BDLP intensity is set as $\lambda = \theta \beta$ so that the stationary
+mean of the Gamma-OU variance process equals $\theta$. This gives the same
+$(v_0, \theta)$ parameterisation as Heston.
+
+Because the variance is built from positive jumps and exponential mean
+reversion, it stays positive by construction. No Feller-style penalty is
+needed.
+
+## How BNS Fits the Surface
+
+The mechanism that produces a smile in BNS is structurally different from
+Heston. Heston relies on a diffusive volatility-of-variance $\sigma$ for the
+wings and a spot-variance correlation $\rho$ for the skew, both accumulating
+as $\sqrt{T}$.
+
+BNS instead injects discrete jumps directly into the variance process: each
+jump in $v_t$ is mirrored, scaled by $\rho$, into the log-price. The wing
+thickness is governed by the jump-size distribution (controlled by $\beta$)
+and the skew by $\rho$.
+
+A consequence of this structural difference is that the calibrator often
+settles at a small $\kappa$ together with a large $\theta$. The time scale of
+mean reversion is $1/\kappa$, so when $\kappa$ is small the variance process
+barely relaxes towards $\theta$ over the calibration horizon and stays close
+to $v_0$ throughout.
+
+In that regime $\theta$ is only weakly identified by the surface and the
+optimizer can move it freely as long as the jump-driven smile dynamics are
+preserved. The headline number to read in the output is $v_0$, which sets the
+at-the-money level.
+
+## Calibration
+
+The fit reuses [VolModelCalibration][quantflow.options.calibration.base.VolModelCalibration]
+two-stage optimiser from the Heston tutorial: L-BFGS-B for basin search,
+followed by trust-region reflective on the residual vector with parameter
+bounds.
+
+--8<-- "docs/examples/output/vol_surface_bns_calibration.out"
+
+## Calibrated Smile
+
+[![BNS calibrated smile](../assets/examples/bns_calibrated_smile.png)](../assets/examples/bns_calibrated_smile.png){target="_blank"}
+
+The fit is good for medium and long maturities and visibly off at the front
+expiries. This is the same short-maturity gap seen for Heston and
+Heston-jump-diffusion.
+
+The cause here is structural: BNS adds jumps, but they live in the variance
+process, not directly in the log-price. The jump-driven contribution to the
+log-price is bounded by the size of the variance jumps multiplied by $|\rho|$,
+which is small for short tenors.
+
+A model with explicit jumps in the log-price (such as
+[HestonJ][quantflow.sp.heston.HestonJ]) or a rough volatility model is better
+suited to the steep short-term skew observed in crypto markets.
+
+## Code
+
+```python
+--8<-- "docs/examples/vol_surface_bns_calibration.py"
+```
diff --git a/docs/tutorials/index.md b/docs/tutorials/index.md
index 8f68fdd8..4241fbf0 100644
--- a/docs/tutorials/index.md
+++ b/docs/tutorials/index.md
@@ -6,3 +6,4 @@ Step-by-step guides for common quantflow workflows.
 |---|---|
 | [Option Pricing](option_pricing.md) | Price a European option with the Black-Scholes and Heston-jump-diffusion models |
 | [Volatility Surface](volatility_surface.md) | Fetch live option data, build an implied volatility surface, and calibrate Heston and jump-diffusion models |
+| [BNS Volatility Model](bns_calibration.md) | Calibrate the Barndorff-Nielsen and Shephard stochastic-volatility model to an implied volatility surface |
diff --git a/docs/tutorials/volatility_surface.md b/docs/tutorials/volatility_surface.md
index c23e1683..9a5682e3 100644
--- a/docs/tutorials/volatility_surface.md
+++ b/docs/tutorials/volatility_surface.md
@@ -100,7 +100,7 @@ surface2 = surface_from_inputs(inputs)    # VolSurfaceInputs -> VolSurface
 
 ## Calibrating the Heston Model
 
-[HestonCalibration][quantflow.options.heston_calibration.HestonCalibration] fits the
+[HestonCalibration][quantflow.options.calibration.heston.HestonCalibration] fits the
 five Heston parameters ($v_0$, $\theta$, $\kappa$, $\sigma$, $\rho$) to the implied
 volatility surface using a two-stage optimisation:
 
@@ -131,7 +131,7 @@ down-weights near-expiry options.
 
 ### Plotting the Calibrated Smile
 
-Use [plot_maturities()][quantflow.options.calibration.VolModelCalibration.plot_maturities]
+Use [plot_maturities()][quantflow.options.calibration.base.VolModelCalibration.plot_maturities]
 to produce a Plotly figure overlaying market bid/ask implied vols against the model smile
 for all maturities at once:
 
@@ -178,7 +178,7 @@ the motivation for the Heston jump-diffusion model described in the next section
 
 ## Calibrating the Heston Jump-Diffusion Model
 
-[HestonJCalibration][quantflow.options.heston_calibration.HestonJCalibration] extends the
+[HestonJCalibration][quantflow.options.calibration.heston.HestonJCalibration] extends the
 Heston calibration with a compound Poisson jump component via the
 [HestonJ][quantflow.sp.heston.HestonJ] model. Jumps are drawn from a
 [DoubleExponential][quantflow.utils.distributions.DoubleExponential] distribution,
@@ -239,3 +239,4 @@ The calibrated parameter vector for the jump-diffusion model is:
 | `jump intensity` | Jump arrival rate (jumps per year) |
 | `jump variance` | Variance of a single jump |
 | `jump asymmetry` | Asymmetry of the jump distribution ([DoubleExponential][quantflow.utils.distributions.DoubleExponential]) |
+
diff --git a/mkdocs.yml b/mkdocs.yml
index b7bd1981..2260d10b 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -107,6 +107,7 @@ nav:
           - Types: api/utils/types.md
   - Tutorials:
       - tutorials/index.md
+      - BNS Volatility Model: tutorials/bns_calibration.md
       - CIR Process: tutorials/cir.md
       - Option Pricing: tutorials/option_pricing.md
       - Pricing Method Comparison: tutorials/pricing_method_comparison.md
@@ -145,5 +146,4 @@ markdown_extensions:
 
 extra_javascript:
   - javascripts/mathjax.js
-  - https://polyfill.io/v3/polyfill.min.js?features=es6
   - https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js
diff --git a/quantflow/options/calibration/__init__.py b/quantflow/options/calibration/__init__.py
new file mode 100644
index 00000000..2ba19d79
--- /dev/null
+++ b/quantflow/options/calibration/__init__.py
@@ -0,0 +1,23 @@
+from .base import (
+    ModelCalibrationEntryKey,
+    OptionEntry,
+    VolModelCalibration,
+)
+from .bns import BNSCalibration
+from .heston import (
+    DoubleHestonCalibration,
+    DoubleHestonJCalibration,
+    HestonCalibration,
+    HestonJCalibration,
+)
+
+__all__ = [
+    "BNSCalibration",
+    "DoubleHestonCalibration",
+    "DoubleHestonJCalibration",
+    "HestonCalibration",
+    "HestonJCalibration",
+    "ModelCalibrationEntryKey",
+    "OptionEntry",
+    "VolModelCalibration",
+]
diff --git a/quantflow/options/calibration.py b/quantflow/options/calibration/base.py
similarity index 99%
rename from quantflow/options/calibration.py
rename to quantflow/options/calibration/base.py
index c853612a..541296c5 100644
--- a/quantflow/options/calibration.py
+++ b/quantflow/options/calibration/base.py
@@ -13,8 +13,8 @@
 from quantflow.sp.base import StochasticProcess1D
 from quantflow.utils import plot
 
-from .pricer import OptionPricerBase
-from .surface import OptionPrice, VolSurface
+from ..pricer import OptionPricerBase
+from ..surface import OptionPrice, VolSurface
 
 M = TypeVar("M", bound=StochasticProcess1D)
 
diff --git a/quantflow/options/calibration/bns.py b/quantflow/options/calibration/bns.py
new file mode 100644
index 00000000..10c88732
--- /dev/null
+++ b/quantflow/options/calibration/bns.py
@@ -0,0 +1,55 @@
+from __future__ import annotations
+
+from typing import Generic, TypeVar
+
+import numpy as np
+from scipy.optimize import Bounds
+
+from quantflow.sp.bns import BNS
+
+from .base import VolModelCalibration
+
+B = TypeVar("B", bound=BNS)
+
+
+class BNSCalibration(VolModelCalibration[B], Generic[B]):
+    r"""Calibration of the [BNS][quantflow.sp.bns.BNS] stochastic volatility model.
+
+    The parameter vector is `[v0, theta, kappa, beta, rho]` where
+
+    - `v0` is the initial variance ($v_0 = \text{variance\_process.rate}$)
+    - `theta` is the long-run variance ($\theta = \lambda / \beta$)
+    - `kappa` is the mean-reversion speed of the variance process
+    - `beta` is the exponential decay of the BDLP jump-size distribution
+    - `rho` is the leverage parameter (correlation between jumps in variance and
+      jumps in log-price)
+
+    The BDLP intensity is set as $\lambda = \theta \beta$ so that the stationary
+    mean of the variance process equals $\theta$, mirroring the Heston
+    parameterisation. The Gamma-OU variance process is positive by construction,
+    so no Feller-style penalty is needed.
+    """
+
+    def get_bounds(self) -> Bounds:
+        vol_range = self.implied_vol_range()
+        vol_lb = 0.5 * vol_range.lb[0]
+        vol_ub = 1.5 * vol_range.ub[0]
+        v2 = vol_lb**2
+        v2u = vol_ub**2
+        return Bounds(
+            [v2, v2, 1e-3, 1.0, -0.9],
+            [v2u, v2u, np.inf, np.inf, 0.0],
+        )
+
+    def get_params(self) -> np.ndarray:
+        vp = self.model.variance_process
+        theta = vp.intensity / vp.beta
+        return np.asarray([vp.rate, theta, vp.kappa, vp.beta, self.model.rho])
+
+    def set_params(self, params: np.ndarray) -> None:
+        vp = self.model.variance_process
+        vp.rate = params[0]
+        vp.kappa = params[2]
+        vp.bdlp.jumps.decay = params[3]
+        vp.bdlp.intensity = params[1] * params[3]
+        self.model.rho = params[4]
diff --git a/quantflow/options/heston_calibration.py b/quantflow/options/calibration/heston.py
similarity index 97%
rename from quantflow/options/heston_calibration.py
rename to quantflow/options/calibration/heston.py
index 92a5a376..94952c20 100644
--- a/quantflow/options/heston_calibration.py
+++ b/quantflow/options/calibration/heston.py
@@ -9,8 +9,8 @@
 from quantflow.sp.heston import DoubleHeston, DoubleHestonJ, Heston, HestonJ
 from quantflow.sp.jump_diffusion import D
 
-from .calibration import VolModelCalibration
-from .pricer import OptionPricer
+from ..pricer import OptionPricer
+from .base import VolModelCalibration
 
 H = TypeVar("H", bound=Heston)
 
@@ -67,7 +67,7 @@ def penalize(self) -> float:
 class HestonJCalibration(HestonCalibration[HestonJ[D]], Generic[D]):
     """Calibration of the [HestonJ][quantflow.sp.heston.HestonJ] model with jumps.
 
-    Extends [HestonCalibration][quantflow.options.heston_calibration.HestonCalibration]
+    Extends [HestonCalibration][quantflow.options.calibration.heston.HestonCalibration]
     by appending jump parameters to the parameter vector and bounds.
     """
 
@@ -260,11 +260,11 @@ class DoubleHestonJCalibration(DoubleHestonCalibration[DoubleHestonJ[D]], Generi
     """Calibration of the [DoubleHestonJ][quantflow.sp.heston.DoubleHestonJ] model.
 
     Extends
-    [DoubleHestonCalibration][quantflow.options.heston_calibration.DoubleHestonCalibration]
+    [DoubleHestonCalibration][quantflow.options.calibration.heston.DoubleHestonCalibration]
     by appending the jump parameters of `heston1` to the parameter vector and bounds.
 
     Overrides `warm_start` to fit a full
-    [HestonJCalibration][quantflow.options.heston_calibration.HestonJCalibration]
+    [HestonJCalibration][quantflow.options.calibration.heston.HestonJCalibration]
     to the short-dated options, so that the jump parameters are also initialised
     before the joint optimisation.
     """
diff --git a/quantflow/sp/bns.py b/quantflow/sp/bns.py
index a167e2bc..a7fa1f35 100644
--- a/quantflow/sp/bns.py
+++ b/quantflow/sp/bns.py
@@ -4,11 +4,10 @@
 
 import numpy as np
 from pydantic import Field
-from scipy.special import xlogy
 
 from ..ta.paths import Paths
 from ..utils.types import FloatArrayLike, Vector
-from .base import Im, StochasticProcess1D
+from .base import StochasticProcess1D
 from .ou import GammaOU
 
 
@@ -55,7 +54,48 @@ def create(cls, vol: float, kappa: float, decay: float, rho: float) -> Self:
         )
 
     def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector:
-        return -self._zeta(t, 0.5 * Im * u * u, self.rho * u)
+        r"""Characteristic exponent of the BNS process with Gamma-OU variance.
+
+        \begin{equation}
+            \phi(t, u) = B v_0
+                - \lambda \left[
+                    \frac{\kappa t \, A}{\beta - A}
+                    + \frac{\beta}{\beta - A}
+                    \log\frac{\beta - i u \rho + B}{\beta - i u \rho}
+                \right]
+        \end{equation}
+
+        with
+
+        \begin{equation}
+            \begin{aligned}
+                A &= i u \rho - \frac{u^2}{2 \kappa}, \\
+                B &= \frac{u^2 (1 - e^{-\kappa t})}{2 \kappa}.
+            \end{aligned}
+        \end{equation}
+
+        where $v_0$ is the initial variance, $\kappa$ is the mean-reversion speed,
+        $\rho$ is the leverage parameter, and $(\lambda, \beta)$ are the intensity
+        and exponential-jump rate of the background driving Lévy process.
+        """
+        v = self.variance_process
+        k = v.kappa
+        beta = v.beta
+        intensity = v.intensity
+        v0 = v.rate
+        rho = self.rho
+
+        iur = 1j * u * rho
+        u2 = u * u
+        a = iur - 0.5 * u2 / k
+        b = 0.5 * u2 * (1 - np.exp(-k * t)) / k
+
+        diffusion = b * v0
+        bdlp = intensity * (
+            k * t * a / (beta - a)
+            + beta / (beta - a) * np.log((beta - iur + b) / (beta - iur))
+        )
+        return diffusion - bdlp
 
     def sample(self, n: int, time_horizon: float = 1, time_steps: int = 100) -> Paths:
         return self.sample_from_draws(Paths.normal_draws(n, time_horizon, time_steps))
@@ -76,20 +116,3 @@ def sample_from_draws(self, path_dw: Paths, *args: Paths) -> Paths:
         paths = np.zeros(dw.shape)
         paths[1:] = np.cumsum(dw[:-1], axis=0) + path_dz.data
         return Paths(t=path_dw.t, data=paths)
-
-    # Internal characteristics function methods (see docs)
-
-    def _zeta(self, t: Vector, a: Vector, b: Vector) -> Vector:
-        k = self.variance_process.kappa
-        c = a * (1 - np.exp(-k * t)) / k
-        g = (a + b) / self.variance_process.beta
-        return Im * c * self.variance_process.rate - self.variance_process.intensity * (
-            self._i(b + c, g) - self._i(b, g)
-        )
-
-    def _i(self, x: Vector, g: Vector) -> Vector:
-        k = self.variance_process.kappa
-        beta = self.variance_process.beta
-        l1 = xlogy(k - Im * g, x + Im * beta)
-        l2 = xlogy(g / (g + Im * k) / k, beta * g / k - x)
-        return l1 + l2
diff --git a/quantflow_tests/test_bns.py b/quantflow_tests/test_bns.py
new file mode 100644
index 00000000..92cf4f10
--- /dev/null
+++ b/quantflow_tests/test_bns.py
@@ -0,0 +1,47 @@
+import numpy as np
+import pytest
+
+from quantflow.sp.bns import BNS
+from quantflow_tests.utils import characteristic_tests
+
+
+@pytest.fixture
+def bns() -> BNS:
+    return BNS.create(vol=0.5, decay=5, kappa=1, rho=0)
+
+
+@pytest.fixture
+def bns_leverage() -> BNS:
+    return BNS.create(vol=0.5, decay=5, kappa=1, rho=-0.5)
+
+
+def test_bns(bns: BNS) -> None:
+    m = bns.marginal(1)
+    characteristic_tests(m)
+    assert bns.characteristic(1, 0) == 1
+    assert m.mean_from_characteristic() == pytest.approx(0.0, abs=1e-6)
+    # stationary variance = vol*vol = 0.25, so std on [0, 1] = 0.5
+    assert m.std_from_characteristic() == pytest.approx(0.5, 1e-3)
+
+
+def test_bns_horizon(bns: BNS) -> None:
+    m = bns.marginal(2)
+    # variance scales linearly with t when v0 equals stationary mean
+    assert m.std_from_characteristic() == pytest.approx(2**0.5 * 0.5, 1e-3)
+
+
+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)
+    cf_bns = bns.characteristic(1, u)
+    cf_via_ou = bns.variance_process.integrated_log_laplace(1, 1j * u * u / 2)
+    np.testing.assert_allclose(cf_bns, cf_via_ou, atol=1e-12)
+
+
+def test_bns_leverage(bns_leverage: BNS) -> None:
+    m = bns_leverage.marginal(1)
+    characteristic_tests(m)
+    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)
diff --git a/quantflow_tests/test_options.py b/quantflow_tests/test_options.py
index b83ad178..cbf6ad48 100644
--- a/quantflow_tests/test_options.py
+++ b/quantflow_tests/test_options.py
@@ -6,12 +6,13 @@
 import pytest
 
 from quantflow.options import bs
-from quantflow.options.calibration import ModelCalibrationEntryKey, OptionEntry
-from quantflow.options.heston_calibration import (
+from quantflow.options.calibration import (
     DoubleHestonCalibration,
     DoubleHestonJCalibration,
     HestonCalibration,
     HestonJCalibration,
+    ModelCalibrationEntryKey,
+    OptionEntry,
 )
 from quantflow.options.inputs import OptionInput
 from quantflow.options.pricer import OptionPricer
diff --git a/quantflow_tests/test_ou.py b/quantflow_tests/test_ou.py
index eba82a61..d3ef0759 100644
--- a/quantflow_tests/test_ou.py
+++ b/quantflow_tests/test_ou.py
@@ -1,6 +1,5 @@
 import pytest
 
-from quantflow.sp.bns import BNS
 from quantflow.sp.ou import GammaOU, Vasicek
 from quantflow_tests.utils import analytical_tests, characteristic_tests
 
@@ -15,11 +14,6 @@ def gamma_ou() -> GammaOU:
     return GammaOU.create(decay=10, kappa=5)
 
 
-@pytest.fixture
-def bns() -> BNS:
-    return BNS.create(vol=0.5, decay=5, kappa=1, rho=0)
-
-
 def test_marginal(gamma_ou: GammaOU) -> None:
     m = gamma_ou.marginal(1)
     assert m.mean() == 1
@@ -39,10 +33,3 @@ def test_vasicek(vasicek: Vasicek) -> None:
     assert m.variance() == pytest.approx(0.1)
     assert m.mean_from_characteristic() == pytest.approx(1.0, 1e-3)
     assert m.std_from_characteristic() == pytest.approx(m.std(), 1e-3)
-
-
-def test_bns(bns: BNS):
-    m = bns.marginal(1)
-    assert bns.characteristic(1, 0) == 1
-    assert m.mean() == 0.0
-    # assert pytest.approx(m.std(), 1e-3) == 0.5