Add BNS2

dev/quantflow.dockerfile

@@ -20,6 +20,7 @@ COPY mkdocs.yml ./
 COPY dev/ ./dev/
 COPY docs/ ./docs/
 COPY quantflow/ ./quantflow/
+COPY app/ ./app/
 RUN uv run ./dev/build-examples
 RUN uv run mkdocs build
 
@@ -31,14 +32,11 @@ WORKDIR /app
 # Copy virtualenv from builder
 COPY --from=builder /build/.venv /app/.venv
 
-# Copy application code
+# Copy application code (app/ from builder includes built docs)
 COPY quantflow/ ./quantflow/
-COPY app/ ./app/
+COPY --from=builder /build/app ./app
 COPY pyproject.toml ./
 
-# Copy built documentation
-COPY --from=builder /build/app/docs ./app/docs
-
 ENV PYTHONPATH=/app
 ENV PYTHONUNBUFFERED=1
 ENV PATH="/app/.venv/bin:$PATH"

docs/api/options/calibration.md

@@ -13,3 +13,5 @@
 ::: quantflow.options.calibration.heston.DoubleHestonJCalibration
 
 ::: quantflow.options.calibration.bns.BNSCalibration
+
+::: quantflow.options.calibration.bns.BNS2Calibration

docs/api/sp/bns.md

@@ -1,3 +1,6 @@
 # Barndorff-Nielson & Shephard process
 
 ::: quantflow.sp.bns.BNS
+
+
+::: quantflow.sp.bns.BNS2

docs/api/sp/index.md

@@ -16,7 +16,8 @@ This page gives an overview of all stochastic processes available in the library
 |---|---|
 | [CIR][quantflow.sp.cir.CIR] | Cox-Ingersoll-Ross square-root diffusion |
 | [Vasicek][quantflow.sp.ou.Vasicek] | Gaussian Ornstein-Uhlenbeck process |
-| [GammaOU][quantflow.sp.ou.GammaOU] | Non-Gaussian OU process driven by a Gamma subordinator |
+| [NGOU][quantflow.sp.ou.NGOU] | Generic non-Gaussian OU process driven by a pure-jump Lévy process |
+| [GammaOU][quantflow.sp.ou.GammaOU] | Non-Gaussian OU process with Gamma stationary marginal |
 
 ### Jump processes
 
@@ -33,14 +34,24 @@ This page gives an overview of all stochastic processes available in the library
 | [Heston][quantflow.sp.heston.Heston] | Classical square-root stochastic volatility model |
 | [HestonJ][quantflow.sp.heston.HestonJ] | Heston model with compound Poisson jumps |
 | [DoubleHeston][quantflow.sp.heston.DoubleHeston] | Two independent Heston variance processes |
-| [DoubleHestonJ][quantflow.sp.heston.DoubleHestonJ] | Double Heston with compound Poisson jumps on the second component |
+| [DoubleHestonJ][quantflow.sp.heston.DoubleHestonJ] | Double Heston with compound Poisson jumps on the first component |
+| [BNS][quantflow.sp.bns.BNS] | Barndorff-Nielsen and Shephard model with Gamma-OU variance |
+| [BNS2][quantflow.sp.bns.BNS2] | Two-factor BNS with convex-combination variance |
 
 ### Jump diffusion
 
 | Process | Description |
 |---|---|
 | [JumpDiffusion][quantflow.sp.jump_diffusion.JumpDiffusion] | Diffusion with compound Poisson jumps |
 
+### Copulas
+
+| Process | Description |
+|---|---|
+| [Copula][quantflow.sp.copula.Copula] | Abstract base class for bivariate copulas |
+| [IndependentCopula][quantflow.sp.copula.IndependentCopula] | Independence copula $C(u, v) = u v$ |
+| [FrankCopula][quantflow.sp.copula.FrankCopula] | Archimedean Frank copula |
+
 ## Base classes
 
 ::: quantflow.sp.base.StochasticProcess

docs/examples/vol_surface_bns2_calibration.py

@@ -0,0 +1,39 @@
+import json
+
+from docs.examples._utils import assets_path, print_model
+from quantflow.options.calibration import BNS2Calibration
+from quantflow.options.pricer import OptionPricer
+from quantflow.options.surface import VolSurface, VolSurfaceInputs, surface_from_inputs
+from quantflow.sp.bns import BNS, BNS2
+
+# 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()
+
+# Two-factor BNS: a fast factor for short maturities and a slow one for long.
+# Opposite-sign leverages lets one factor lift the OTM-call wing (rho>0) while
+# the other carries the equity-style downside skew (rho<0).
+pricer = OptionPricer(
+    model=BNS2(
+        bns1=BNS.create(vol=0.4, kappa=20.0, decay=20.0, rho=-0.6),
+        bns2=BNS.create(vol=0.5, kappa=0.3, decay=5.0, rho=0.3),
+        weight=0.3,
+    )
+)
+
+calibration: BNS2Calibration[BNS2] = BNS2Calibration(
+    pricer=pricer,
+    vol_surface=surface,
+    moneyness_weight=0.5,
+)
+
+result = calibration.fit()
+print(result.message)
+print_model(calibration.model)
+
+fig = calibration.plot_maturities(max_moneyness=1.5, support=101)
+fig.update_layout(title="BNS2 Calibrated Smiles")
+fig.write_image(assets_path("bns2_calibrated_smile.png"), width=1200)

docs/index.md

@@ -1,4 +1,4 @@
-# <a href="https://quantmind.github.io/quantflow"><img src="https://raw.githubusercontent.com/quantmind/quantflow/main/docs/assets/quantflow-light.svg" width=300 /></a>
+# <a href="https://quantmind.github.io/quantflow"><img src="https://raw.githubusercontent.com/quantmind/quantflow/main/docs/assets/logos/quantflow-lockup.svg" width=300 /></a>
 
 [![PyPI version](https://badge.fury.io/py/quantflow.svg)](https://badge.fury.io/py/quantflow)
 [![Python versions](https://img.shields.io/pypi/pyversions/quantflow.svg)](https://pypi.org/project/quantflow)

docs/tutorials/bns_calibration.md

@@ -1,16 +1,16 @@
 # 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).
+model (Barndorff-Nielsen and Shephard) and its two-factor extension
+[BNS2][quantflow.sp.bns.BNS2] 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
+## Single-factor BNS
 
 [BNSCalibration][quantflow.options.calibration.bns.BNSCalibration] fits five
 parameters to the surface:
@@ -25,46 +25,38 @@ parameters to the surface:
 
 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.
+$(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.
 
-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
+### 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$.
+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.
+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
+### Calibrated parameters
 
-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.
+The fit uses the
+[VolModelCalibration][quantflow.options.calibration.base.VolModelCalibration]
+two-stage optimiser: 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
@@ -74,14 +66,62 @@ 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.
+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
+### Code
 
 ```python
 --8<-- "docs/examples/vol_surface_bns_calibration.py"
 ```
+
+## Two-factor BNS
+
+The original multi-factor BNS extends the single-factor model by replacing the
+variance with a convex combination of independent Gamma-OU processes. With
+weight $w \in [0, 1]$ and a single Brownian motion driving the diffusion,
+
+\begin{equation}
+    \begin{aligned}
+        \sigma^2_t &= w\, v^1_t + (1 - w)\, v^2_t \\
+        dx_t &= \sigma_t\, dw_t
+            + \rho_1\, dz^1_{\kappa_1 t}
+            + \rho_2\, dz^2_{\kappa_2 t}
+    \end{aligned}
+\end{equation}
+
+Pairing a fast-mean-reverting factor with a slow one decouples the
+short-maturity skew from the long-maturity level, in the same spirit as the
+[DoubleHeston][quantflow.sp.heston.DoubleHeston] extension of Heston.
+
+[BNS2Calibration][quantflow.options.calibration.bns.BNS2Calibration] fits
+eleven parameters:
+
+`[v01, theta1, kappa_delta, beta1, rho1, v02, theta2, kappa2, beta2, rho2, w]`
+
+with `kappa1 = kappa2 + kappa_delta` enforcing that the first factor
+mean-reverts faster than the second. Both leverage parameters are free in
+$[-0.9, 0.9]$: a positive $\rho_i$ produces up-jumps in the log-price that
+lift the OTM call wing, while a negative one produces equity-style downside
+skew. There is no warm start, so the optimiser starts from the user-supplied
+initial parameters; pick distinct timescales for `bns1` and `bns2` (and
+consider opposite-sign leverages) to seed a meaningful two-factor fit.
+
+### Calibrated parameters
+
+--8<-- "docs/examples/output/vol_surface_bns2_calibration.out"
+
+[![BNS2 calibrated smile](../assets/examples/bns2_calibrated_smile.png)](../assets/examples/bns2_calibrated_smile.png){target="_blank"}
+
+The two-factor variant adds flexibility on the term structure: the fast factor
+absorbs short-dated skew while the slow factor anchors the long end. The
+remaining short-maturity gap is structural in the same way as the single-factor
+case: BNS2 still injects jumps only through the variance process, so the
+log-price wings are bounded by the jump sizes scaled by $|\rho_i|$.
+
+### Code
+
+```python
+--8<-- "docs/examples/vol_surface_bns2_calibration.py"
+```

quantflow/options/calibration/__init__.py

@@ -3,7 +3,7 @@
     OptionEntry,
     VolModelCalibration,
 )
-from .bns import BNSCalibration
+from .bns import BNS2Calibration, BNSCalibration
 from .heston import (
     DoubleHestonCalibration,
     DoubleHestonJCalibration,
@@ -12,6 +12,7 @@
 )
 
 __all__ = [
+    "BNS2Calibration",
     "BNSCalibration",
     "DoubleHestonCalibration",
     "DoubleHestonJCalibration",

quantflow/options/calibration/bns.py

@@ -5,7 +5,7 @@
 import numpy as np
 from scipy.optimize import Bounds
 
-from quantflow.sp.bns import BNS
+from quantflow.sp.bns import BNS, BNS2
 
 from .base import VolModelCalibration
 
@@ -38,7 +38,7 @@ def get_bounds(self) -> Bounds:
         v2u = vol_ub**2
         return Bounds(
             [v2, v2, 1e-3, 1.0, -0.9],
-            [v2u, v2u, np.inf, np.inf, 0.0],
+            [v2u, v2u, np.inf, np.inf, 0.9],
         )
 
     def get_params(self) -> np.ndarray:
@@ -53,3 +53,85 @@ def set_params(self, params: np.ndarray) -> None:
         vp.bdlp.jumps.decay = params[3]
         vp.bdlp.intensity = params[1] * params[3]
         self.model.rho = params[4]
+
+
+B2 = TypeVar("B2", bound=BNS2)
+
+
+class BNS2Calibration(VolModelCalibration[B2], Generic[B2]):
+    r"""Calibration of the [BNS2][quantflow.sp.bns.BNS2] two-factor BNS model.
+
+    The parameter vector is
+
+    `[v01, theta1, kappa_delta, beta1, rho1, v02, theta2, kappa2, beta2, rho2, w]`
+
+    where `kappa1 = kappa2 + kappa_delta` with `kappa_delta > 0`, enforcing that
+    the first (short-maturity) factor mean-reverts faster than the second, and
+    `w` is the convex-combination weight of the first variance factor. The same
+    $(v_0, \theta)$ parameterisation as
+    [BNSCalibration][quantflow.options.calibration.bns.BNSCalibration] is used
+    for each factor: the BDLP intensity is set as $\lambda_i = \theta_i \beta_i$
+    so the stationary mean of $v^i$ equals $\theta_i$.
+
+    Both leverage parameters are free in $[-0.9, 0.9]$: a positive $\rho_i$
+    produces up-jumps in the log-price that lift the OTM call wing, while a
+    negative one produces equity-style downside skew. The joint fit relies on
+    the user-supplied initial parameters: pick distinct timescales for `bns1`
+    and `bns2` (and consider opposite-sign leverages) to give the optimiser a
+    meaningful two-factor starting point.
+
+    TODO: improve this calibration. The 11-parameter fit is slow (finite-diff
+    Jacobian dominates) and tends to collapse the two timescales into a near
+    single-factor solution unless the initial conditions force them apart.
+    Candidate improvements: analytic Jacobian of the characteristic exponent,
+    a smarter warm start that does not bias the kappas to merge, and tighter
+    bounds on `kappa1` and `beta_i`.
+    """
+
+    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-4, 1.0, -0.9, v2, v2, 1e-3, 1.0, -0.9, 0.0],
+            [v2u, v2u, np.inf, np.inf, 0.9, v2u, v2u, 5.0, np.inf, 0.9, 1.0],
+        )
+
+    def get_params(self) -> np.ndarray:
+        vp1 = self.model.bns1.variance_process
+        vp2 = self.model.bns2.variance_process
+        kappa_delta = max(vp1.kappa - vp2.kappa, 1e-4)
+        theta1 = vp1.intensity / vp1.beta
+        theta2 = vp2.intensity / vp2.beta
+        return np.asarray(
+            [
+                vp1.rate,
+                theta1,
+                kappa_delta,
+                vp1.beta,
+                self.model.bns1.rho,
+                vp2.rate,
+                theta2,
+                vp2.kappa,
+                vp2.beta,
+                self.model.bns2.rho,
+                self.model.weight,
+            ]
+        )
+
+    def set_params(self, params: np.ndarray) -> None:
+        vp1 = self.model.bns1.variance_process
+        vp1.rate = params[0]
+        vp1.bdlp.jumps.decay = params[3]
+        vp1.bdlp.intensity = params[1] * params[3]
+        self.model.bns1.rho = params[4]
+        vp2 = self.model.bns2.variance_process
+        vp2.rate = params[5]
+        vp2.kappa = params[7]
+        vp2.bdlp.jumps.decay = params[8]
+        vp2.bdlp.intensity = params[6] * params[8]
+        self.model.bns2.rho = params[9]
+        vp1.kappa = vp2.kappa + params[2]  # kappa2 + kappa_delta
+        self.model.weight = params[10]