Add calibration for hestonj (#26)

notebooks/api/sp/ou.rst

@@ -3,7 +3,7 @@ OU Processes
 ================
 
 These are the classes that implement gaussian and non-gaussian
-`Ornstein-Uhlenbeck <https://en.wikipedia.org/wiki/Ornstein%E2%80%93Uhlenbeck_process>` process.
+`Ornstein-Uhlenbeck <https://en.wikipedia.org/wiki/Ornstein%E2%80%93Uhlenbeck_process>`_ process.
 
 
 .. currentmodule:: quantflow.sp.ou

quantflow/options/calibration.py

@@ -1,19 +1,20 @@
 from __future__ import annotations
 
 from abc import ABC, abstractmethod
-from dataclasses import dataclass, field, replace
+from dataclasses import dataclass, field
 from datetime import datetime
 from decimal import Decimal
 from typing import Any, Generic, NamedTuple, Sequence, TypeVar
 
 import numpy as np
 import pandas as pd
+from pydantic import BaseModel, Field
 from scipy.optimize import Bounds, OptimizeResult, minimize
 
 from quantflow.sp.base import StochasticProcess1D
 from quantflow.sp.heston import Heston, HestonJ
+from quantflow.sp.jump_diffusion import D
 from quantflow.utils import plot
-from quantflow.utils.distributions import DoubleExponential
 
 from .pricer import OptionPricer
 from .surface import OptionPrice, VolSurface
@@ -55,26 +56,34 @@ def price_range(self) -> Bounds:
         return self._prince_range
 
 
-@dataclass
-class VolModelCalibration(ABC, Generic[M]):
+class VolModelCalibration(BaseModel, ABC, Generic[M], arbitrary_types_allowed=True):
     """Abstract class for calibration of a stochastic volatility model"""
 
     pricer: OptionPricer[M]
     """The :class:`.OptionPricer` for the model"""
-    vol_surface: VolSurface[Any]
+    vol_surface: VolSurface[Any] = Field(repr=False)
     """The :class:`.VolSurface` to calibrate the model with"""
     minimize_method: str | None = None
     """The optimization method to use - if None, the default is used"""
-    moneyness_weight: float = 0.0
-    """The weight for penalize moneyness as it moves away from 0
+    moneyness_weight: float = Field(default=0.0, ge=0.0)
+    """The weight for penalize options with moneyness as it moves away from 0
 
     The weight is applied as exp(-moneyness_weight * moneyness), therefore
     a value of 0 won't penalize moneyness at all
     """
-    options: dict[ModelCalibrationEntryKey, OptionEntry] = field(default_factory=dict)
+    ttm_weight: float = Field(default=0.0, ge=0.0, le=1.0)
+    """The weight for penalize options with ttm as it approaches 0
+
+    The weight is applied as `1 - ttm_weight*exp(-ttm)`, therefore
+    a value of 0 won't penalize ttm at all, a value of 1 will penalize
+    options with ttm->0 the most
+    """
+    options: dict[ModelCalibrationEntryKey, OptionEntry] = Field(
+        default_factory=dict, repr=False
+    )
     """The options to calibrate"""
 
-    def __post_init__(self) -> None:
+    def model_post_init(self, _ctx: Any) -> None:
         if not self.options:
             self.vol_surface.bs()
             for option in self.vol_surface.option_prices():
@@ -121,7 +130,7 @@ def remove_implied_above(self, quantile: float = 0.95) -> VolModelCalibration:
         for key, entry in self.options.items():
             if entry.implied_vol_range().ub <= exclude_above:
                 options[key] = entry
-        return replace(self, options=options)
+        return self.model_copy(update=dict(options=options))
 
     def get_bounds(self) -> Bounds | None:  # pragma: no cover
         """Get the parameter bounds for the calibration"""
@@ -193,7 +202,6 @@ def plot(
         )
 
 
-@dataclass
 class HestonCalibration(VolModelCalibration[Heston]):
     """A :class:`.VolModelCalibration` for the :class:`.Heston`
     stochastic volatility model"""
@@ -233,8 +241,7 @@ def penalize(self) -> float:
         return self.feller_penalize * neg * neg
 
 
-@dataclass
-class HestonJCalibration(VolModelCalibration[HestonJ[DoubleExponential]]):
+class HestonJCalibration(VolModelCalibration[HestonJ[D]], Generic[D]):
     """A :class:`.VolModelCalibration` for the :class:`.HestonJ`
     stochastic volatility model with :class:`.DoubleExponential`
     jumps
@@ -248,22 +255,24 @@ def get_bounds(self) -> Sequence[Bounds] | None:
         vol_ub = 1.5 * vol_range.ub[0]
         return Bounds(
             [vol_lb * vol_lb, vol_lb * vol_lb, 0.0, 0.0, -0.9],
-            [vol_ub * vol_ub, vol_ub * vol_ub, np.inf, np.inf, 0],
+            [vol_ub * vol_ub, vol_ub * vol_ub, np.inf, np.inf, 0.0],
         )
 
     def get_params(self) -> np.ndarray:
-        return np.asarray(
-            [
-                self.model.variance_process.rate,
-                self.model.variance_process.theta,
-                self.model.variance_process.kappa,
-                self.model.variance_process.sigma,
-                self.model.rho,
-                self.model.jumps.intensity,
-                self.model.jumps.jumps.decay,
-                self.model.jumps.jumps.log_kappa,
-            ]
-        )
+        params = [
+            self.model.variance_process.rate,
+            self.model.variance_process.theta,
+            self.model.variance_process.kappa,
+            self.model.variance_process.sigma,
+            self.model.rho,
+            self.model.jumps.intensity,
+            self.model.jumps.jumps.variance(),
+        ]
+        try:
+            params.append(self.model.jumps.jumps.asymmetry())
+        except NotImplementedError:
+            pass
+        return np.asarray(params)
 
     def set_params(self, params: np.ndarray) -> None:
         self.model.variance_process.rate = params[0]
@@ -272,8 +281,11 @@ def set_params(self, params: np.ndarray) -> None:
         self.model.variance_process.sigma = params[3]
         self.model.rho = params[4]
         self.model.jumps.intensity = params[5]
-        self.model.jumps.jumps.decay = params[6]
-        self.model.jumps.jumps.kappa = np.exp(params[7])
+        self.model.jumps.jumps.set_variance(params[6])
+        try:
+            self.model.jumps.jumps.set_asymmetry(params[7])
+        except IndexError:
+            pass
 
     def penalize(self) -> float:
         kt = 2 * self.model.variance_process.kappa * self.model.variance_process.theta

quantflow/options/pricer.py

@@ -1,10 +1,10 @@
 from __future__ import annotations
 
-from dataclasses import dataclass, field
 from typing import Any, Generic, NamedTuple, TypeVar, cast
 
 import numpy as np
 import pandas as pd
+from pydantic import BaseModel, Field
 
 from quantflow.sp.base import StochasticProcess1D
 from quantflow.utils import plot
@@ -106,13 +106,14 @@ def plot(self, series: str = "implied_vol", **kwargs: Any) -> Any:
         return plot.plot_vol_cross(self.df, series=series, **kwargs)
 
 
-@dataclass
-class OptionPricer(Generic[M]):
+class OptionPricer(BaseModel, Generic[M], arbitrary_types_allowed=True):
     """Pricer for options"""
 
     model: M
     """The stochastic process used for pricing"""
-    ttm: dict[int, MaturityPricer] = field(default_factory=dict, repr=False)
+    ttm: dict[int, MaturityPricer] = Field(
+        default_factory=dict, repr=False, exclude=True
+    )
     """Cache for :class:`.MaturityPricer` for different time to maturity"""
     n: int = 128
     """NUmber of discretization points for the marginal distribution"""

quantflow/sp/base.py

@@ -169,7 +169,9 @@ class IntensityProcess(StochasticProcess1D):
     """
 
     rate: float = Field(default=1.0, gt=0, description="Instantaneous initial rate")
+    r"""Instantaneous initial rate :math:`r_0`"""
     kappa: float = Field(default=1.0, gt=0, description="Mean reversion speed")
+    r"""Mean reversion speed :math:`\kappa`"""
 
     @abstractmethod
     def integrated_log_laplace(self, t: FloatArrayLike, u: Vector) -> Vector:

quantflow/sp/heston.py

@@ -160,6 +160,8 @@ def create(  # type: ignore [override]
                 v_0 &= {\tt rate}\cdot{\tt var}_d
             \end{align}
 
+        :param jump_distribution: The distribution of jump size (currently only
+            :class:`.Normal` and :class:`.DoubleExponential` are supported)
         :param rate: define the initial value of the variance process
         :param vol: The standard deviation of the price process, normalized by the
             square root of time, as time tends to infinity
@@ -174,8 +176,6 @@ def create(  # type: ignore [override]
         :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)
         """
         jd = JumpDiffusion.create(
             jump_distribution,

quantflow/sp/jump_diffusion.py

@@ -19,14 +19,16 @@ class JumpDiffusion(StochasticProcess1D, Generic[D]):
         dx_t = \sigma d w_t + d N_t
 
     where :math:`w_t` is a Weiner process with standard deviation :math:`\sigma`
-    and :math:`N_t` is a compound poisson process
-    with intensity :math:`\lambda` and jump distribution `D`
+    and :math:`N_t` is a :class:`.CompoundPoissonProcess`
+    with intensity :math:`\lambda` and generic jump distribution `D`
     """
 
     diffusion: WeinerProcess = Field(
         default_factory=WeinerProcess, description="diffusion"
     )
+    """The diffusion process is a standard :class:`.WeinerProcess`"""
     jumps: CompoundPoissonProcess[D] = Field(description="jump process")
+    """The jump process is a generic :class:`.CompoundPoissonProcess`"""
 
     def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector:
         return self.diffusion.characteristic_exponent(
@@ -63,6 +65,8 @@ def create(
         """Create a jump-diffusion model with a given jump distribution, volatility
         and jump fraction.
 
+        :param jump_distribution: The distribution of jump sizes (currently only
+            :class:`.Normal` and :class:`.DoubleExponential` are supported)
         :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)

quantflow/sp/ou.py

@@ -16,16 +16,28 @@
 
 
 class Vasicek(IntensityProcess):
-    """Gaussian OU process, also know as the Vasiceck model.
+    r"""Gaussian OU process, also know as the `Vasiceck model`_.
 
-    Strictly, it is not an intensity process since it is not positive
+    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.
+
+    .. math::
+        dx_t = \kappa (\theta - x_t) dt + \sigma dw_t
+
+    It derives from :class:`.IntensityProcess`, although, it is not strictly
+    an intensity process since it is not positive.
+
+    .. _`Vasiceck model`: https://en.wikipedia.org/wiki/Vasicek_model
     """
 
     bdlp: WeinerProcess = Field(
         default_factory=WeinerProcess,
         description="Background driving Weiner process",
     )
+    """Background driving Weiner process"""
     theta: float = Field(default=1.0, gt=0, description="Mean rate")
+    r"""Mean rate :math:`\theta`"""
 
     def characteristic_exponent(self, t: FloatArrayLike, u: Vector) -> Vector:
         mu = self.analytical_mean(t)

quantflow/sp/poisson.py

@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 from abc import abstractmethod
 from typing import Any, Generic, TypeVar
 
@@ -7,11 +9,12 @@
 from scipy.optimize import Bounds
 from scipy.stats import poisson
 
-from ..ta.paths import Paths
-from ..utils.distributions import Distribution1D
-from ..utils.functions import factorial
-from ..utils.transforms import TransformResult
-from ..utils.types import FloatArray, FloatArrayLike, Vector
+from quantflow.ta.paths import Paths
+from quantflow.utils.distributions import Distribution1D
+from quantflow.utils.functions import factorial
+from quantflow.utils.transforms import TransformResult
+from quantflow.utils.types import FloatArray, FloatArrayLike, Vector
+
 from .base import Im, StochasticProcess1D, StochasticProcess1DMarginal
 
 D = TypeVar("D", bound=Distribution1D)
@@ -144,7 +147,7 @@ def cdf_jacobian(self, t: FloatArrayLike, n: Vector) -> np.ndarray:
 class CompoundPoissonProcess(PoissonBase, Generic[D]):
     """A generic Compound Poisson process."""
 
-    intensity: float = Field(default=1.0, ge=0, description="intensity rate")
+    intensity: float = Field(default=1.0, gt=0, description="jump intensity rate")
     r"""Intensity rate :math:`\lambda` of the Poisson process
 
     It determines the number of jumps in the same way as the :class:`.PoissonProcess`
@@ -182,6 +185,32 @@ def analytical_variance(self, t: FloatArrayLike) -> FloatArrayLike:
         """Expected variance at a time horizon"""
         return self.intensity * t * (self.jumps.variance() + self.jumps.mean() ** 2)
 
+    @classmethod
+    def create(
+        cls,
+        jump_distribution: type[D],
+        *,
+        vol: float = 0.5,
+        jump_intensity: float = 100,
+        jump_asymmetry: float = 0.0,
+    ) -> CompoundPoissonProcess[D]:
+        """Create a Compound Poisson process with a given jump distribution, volatility,
+        jump intensity a nd jump asymmetry .
+
+        :param jump_distribution: The distribution of jump size (currently only
+            :class:`.Normal` and :class:`.DoubleExponential` are supported)
+        :param vol: Annualized standard deviation
+        :param jump_intensity: The average number of jumps per year
+        :param jump_asymmetry: The asymmetry of the jump distribution (0 for symmetric,
+            only used by distributions with asymmetry)
+        """
+        variance = vol * vol
+        jump_distribution_variance = variance / jump_intensity
+        jumps = jump_distribution.from_variance_and_asymmetry(
+            jump_distribution_variance, jump_asymmetry
+        )
+        return cls(intensity=jump_intensity, jumps=jumps)
+
 
 class MarginalDiscrete1D(StochasticProcess1DMarginal):
     def pdf_from_characteristic(

quantflow/utils/distributions.py

@@ -23,6 +23,18 @@ def from_variance_and_asymmetry(cls, variance: float, asymmetry: float) -> Self:
         """Create a distribution from variance and asymmetry"""
         raise NotImplementedError
 
+    def asymmetry(self) -> float:
+        """Asymmetry of the distribution"""
+        raise NotImplementedError
+
+    def set_variance(self, variance: float) -> None:
+        """Set the variance of the distribution"""
+        raise NotImplementedError
+
+    def set_asymmetry(self, asymmetry: float) -> None:
+        """Set the asymmetry of the distribution"""
+        raise NotImplementedError
+
 
 class Exponential(Distribution1D):
     r"""A :class:`.Marginal1D` for the `Exponential distribution`_
@@ -109,6 +121,10 @@ def support(self, points: int = 100, *, std_mult: float = 4) -> FloatArray:
             points,
         )
 
+    def set_variance(self, variance: float) -> None:
+        """Set the variance of the distribution"""
+        self.sigma = np.sqrt(variance)
+
 
 class DoubleExponential(Exponential):
     r"""A :class:`.Marginal1D` for the generalized double exponential distribution
@@ -204,3 +220,15 @@ def support(self, points: int = 100, *, std_mult: float = 4) -> FloatArray:
             self.mean() + std_mult * self.std(),
             points,
         )
+
+    def asymmetry(self) -> float:
+        """Asymmetry of the distribution"""
+        return np.log(self.kappa)
+
+    def set_variance(self, variance: float) -> None:
+        """Set the variance of the distribution"""
+        k2 = self.kappa * self.kappa
+        self.decay = np.sqrt((1 + k2 * k2) / (variance * k2))
+
+    def set_asymmetry(self, asymmetry: float) -> None:
+        self.kappa = np.exp(asymmetry)

quantflow_tests/test_options.py

@@ -22,7 +22,7 @@
 
 @pytest.fixture
 def heston() -> OptionPricer[Heston]:
-    return OptionPricer(Heston.create(vol=0.5, kappa=1, sigma=0.8, rho=0))
+    return OptionPricer(model=Heston.create(vol=0.5, kappa=1, sigma=0.8, rho=0))
 
 
 @pytest.fixture

quantflow_tests/test_options_pricer.py

@@ -9,7 +9,7 @@
 @pytest.fixture
 def pricer() -> OptionPricer[HestonJ[DoubleExponential]]:
     return OptionPricer(
-        HestonJ.create(DoubleExponential, vol=0.5, kappa=1, sigma=0.8, rho=0)
+        model=HestonJ.create(DoubleExponential, vol=0.5, kappa=1, sigma=0.8, rho=0)
     )