Carr-Madan and Hilbert method

fypy/model/levy/BilateralGamma.py

@@ -7,15 +7,19 @@
 
 
 class _BilateralGammaBase(LevyModel):
-    def __init__(self,
-                 forwardCurve: ForwardCurve,
-                 discountCurve: DiscountCurve,
-                 params: np.ndarray):
+    def __init__(
+        self,
+        forwardCurve: ForwardCurve,
+        discountCurve: DiscountCurve,
+        params: np.ndarray,
+    ):
         """
         Bilateral Gamma (BG) model
         :param forwardCurve: ForwardCurve term structure
         """
-        super().__init__(forwardCurve=forwardCurve, discountCurve=discountCurve, params=params)
+        super().__init__(
+            forwardCurve=forwardCurve, discountCurve=discountCurve, params=params
+        )
 
     # ================
     # Model Parameters
@@ -52,7 +56,12 @@ def cumulants(self, T: float) -> Cumulants:
         :param T: float, time to maturity (time at which cumulants are evaluated)
         :return: Cumulants object
         """
-        alpha_p, lam_p, alpha_m, lam_m = self.alpha_p, self.lambda_p, self.alpha_m, self.lambda_m
+        alpha_p, lam_p, alpha_m, lam_m = (
+            self.alpha_p,
+            self.lambda_p,
+            self.alpha_m,
+            self.lambda_m,
+        )
 
         rn_drift = self.risk_neutral_log_drift()
 
@@ -68,34 +77,52 @@ def factorial(n):
         # Using equation (2.8)
 
         def cumulant(n: int):
-            return factorial(n - 1) * (alpha_p / lam_p ** n + (-1) ** n * alpha_m / lam_m ** n)
-
-        return Cumulants(T=T,
-                         rn_drift=rn_drift,
-                         c1=T * rn_drift,  # + cumulant(1)
-                         c2=T * cumulant(2),
-                         c4=T * cumulant(4))
+            return factorial(n - 1) * (
+                alpha_p / lam_p**n + (-1) ** n * alpha_m / lam_m**n
+            )
+
+        return Cumulants(
+            T=T,
+            rn_drift=rn_drift,
+            c1=T * rn_drift,  # + cumulant(1)
+            c2=T * cumulant(2),
+            c4=T * cumulant(4),
+        )
 
     def convexity_correction(self) -> float:
         """
         Computes the convexity correction for the Levy model, added to log process drift to ensure
         risk neutrality
         """
-        alpha_p, lam_p, alpha_m, lam_m = self.alpha_p, self.lambda_p, self.alpha_m, self.lambda_m
-        return -np.log((lam_p / (lam_p - 1)) ** alpha_p * (lam_m / (lam_m + 1)) ** alpha_m)
+        alpha_p, lam_p, alpha_m, lam_m = (
+            self.alpha_p,
+            self.lambda_p,
+            self.alpha_m,
+            self.lambda_m,
+        )
+        return -np.log(
+            (lam_p / (lam_p - 1)) ** alpha_p * (lam_m / (lam_m + 1)) ** alpha_m
+        )
 
     def symbol(self, xi: Union[float, np.ndarray]):
         """
         Levy symbol, uniquely defines Characteristic Function via: chf(T,xi) = exp(T*symbol(xi)),  for all T>=0
         :param xi: np.ndarray or float, points in frequency domain
         :return: np.ndarray or float, symbol evaluated at input points in frequency domain
         """
-        alpha_p, lam_p, alpha_m, lam_m = self.alpha_p, self.lambda_p, self.alpha_m, self.lambda_m
+        alpha_p, lam_p, alpha_m, lam_m = (
+            self.alpha_p,
+            self.lambda_p,
+            self.alpha_m,
+            self.lambda_m,
+        )
 
         rn_drift = self.risk_neutral_log_drift()
 
-        return 1j * xi * rn_drift \
-               + np.log((lam_p / (lam_p - 1j * xi)) ** alpha_p * (lam_m / (lam_m + 1j * xi)) ** alpha_m)
+        return 1j * xi * rn_drift + np.log(
+            (lam_p / (lam_p - 1j * xi)) ** alpha_p
+            * (lam_m / (lam_m + 1j * xi)) ** alpha_m
+        )
 
     # =============================
     # Calibration Interface Implementation
@@ -112,30 +139,37 @@ def default_params(self) -> Optional[np.ndarray]:
 
 
 class BilateralGamma(_BilateralGammaBase):
-    def __init__(self,
-                 forwardCurve: ForwardCurve,
-                 discountCurve: DiscountCurve,
-                 alpha_p: float,
-                 lambda_p: float,
-                 alhpa_m: float,
-                 lambda_m: float):
+    def __init__(
+        self,
+        forwardCurve: ForwardCurve,
+        discountCurve: DiscountCurve,
+        alpha_p: float,
+        lambda_p: float,
+        alhpa_m: float,
+        lambda_m: float,
+    ):
         """
         Bilateral Gamma (BG) model
         :param forwardCurve: ForwardCurve term structure
         """
-        super().__init__(forwardCurve=forwardCurve, discountCurve=discountCurve,
-                         params=np.asarray([alpha_p, lambda_p, alhpa_m, lambda_m]))
+        super().__init__(
+            forwardCurve=forwardCurve,
+            discountCurve=discountCurve,
+            params=np.asarray([alpha_p, lambda_p, alhpa_m, lambda_m]),
+        )
 
 
 class BilateralGammaMotion(_BilateralGammaBase):
-    def __init__(self,
-                 forwardCurve: ForwardCurve,
-                 discountCurve: DiscountCurve,
-                 alpha_p: float,
-                 lambda_p: float,
-                 alhpa_m: float,
-                 lambda_m: float,
-                 sigma: float):
+    def __init__(
+        self,
+        forwardCurve: ForwardCurve,
+        discountCurve: DiscountCurve,
+        alpha_p: float,
+        lambda_p: float,
+        alhpa_m: float,
+        lambda_m: float,
+        sigma: float,
+    ):
         """
         Bilateral Gamma Motion (BGM) model, an extension of Bilateral Gamma model to include Brownian motion
         component, resulting in significantly better calibration and elimimation of smile kinks produced by
@@ -145,8 +179,11 @@ def __init__(self,
 
         :param forwardCurve: ForwardCurve term structure
         """
-        super().__init__(forwardCurve=forwardCurve, discountCurve=discountCurve,
-                         params=np.asarray([alpha_p, lambda_p, alhpa_m, lambda_m, sigma]))
+        super().__init__(
+            forwardCurve=forwardCurve,
+            discountCurve=discountCurve,
+            params=np.asarray([alpha_p, lambda_p, alhpa_m, lambda_m, sigma]),
+        )
 
     # ================
     # Model Parameters
@@ -167,7 +204,7 @@ def cumulants(self, T: float) -> Cumulants:
         :return: Cumulants object
         """
         cumulants = super().cumulants(T)
-        cumulants.c2 += self.sigma ** 2
+        cumulants.c2 += T * self.sigma**2
 
         return cumulants
 
@@ -176,15 +213,15 @@ def convexity_correction(self) -> float:
         Computes the convexity correction for the Levy model, added to log process drift to ensure
         risk neutrality
         """
-        return super().convexity_correction() - 0.5 * self.sigma ** 2
+        return super().convexity_correction() - 0.5 * self.sigma**2
 
     def symbol(self, xi: Union[float, np.ndarray]):
         """
         Levy symbol, uniquely defines Characteristic Function via: chf(T,xi) = exp(T*symbol(xi)),  for all T>=0
         :param xi: np.ndarray or float, points in frequency domain
         :return: np.ndarray or float, symbol evaluated at input points in frequency domain
         """
-        return super().symbol(xi=xi) - 0.5 * self.sigma ** 2 * xi * xi
+        return super().symbol(xi=xi) - 0.5 * self.sigma**2 * xi * xi
 
     # =============================
     # Calibration Interface Implementation

fypy/model/levy/TemperedStable.py

@@ -0,0 +1,144 @@
+from fypy.model.levy.LevyModel import LevyModel
+from fypy.model.FourierModel import Cumulants
+from fypy.termstructures.ForwardCurve import ForwardCurve
+from fypy.termstructures.DiscountCurve import DiscountCurve
+import numpy as np
+from typing import List, Tuple, Optional, Union
+from scipy.special import gamma
+
+
+class TemperedStable(LevyModel):
+    def __init__(
+        self,
+        forwardCurve: ForwardCurve,
+        discountCurve: DiscountCurve,
+        alpha_p: float = 0.2,
+        beta_p: float = 0.5,
+        lambda_p: float = 1,
+        alpha_m: float = 0.3,
+        beta_m: float = 0.3,
+        lambda_m: float = 2,
+    ):
+        """
+        Carr-Geman-Madan-Yor (CGMY) model.  When Y=0, this model reduces to VG
+        :param forwardCurve: ForwardCurve term structure
+        :param C: float, viewed as a measure of the overall level of activity, and influences kurtosis
+        :param G: float, rate of exponential decay on the right tail
+        :param M: float, rate of exponential decay on the left tail. Typically for equities G < M, ie the left
+            tail is then heavier than the right (more down risk)
+        :param Y: float, controls the "fine structure" of the process
+        """
+        super().__init__(
+            forwardCurve=forwardCurve,
+            discountCurve=discountCurve,
+            params=np.asarray([alpha_p, beta_p, lambda_p, alpha_m, beta_m, lambda_m]),
+        )
+
+    # ================
+    # Model Parameters
+    # ================
+
+    @property
+    def alpha_p(self) -> float:
+        """Model Parameter"""
+        return self._params[0]
+
+    @property
+    def beta_p(self) -> float:
+        """Model Parameter"""
+        return self._params[1]
+
+    @property
+    def lambda_p(self) -> float:
+        """Model Parameter"""
+        return self._params[2]
+
+    @property
+    def alpha_m(self) -> float:
+        """Model Parameter"""
+        return self._params[3]
+
+    @property
+    def beta_m(self) -> float:
+        """Model Parameter"""
+        return self._params[4]
+
+    @property
+    def lambda_m(self) -> float:
+        """Model Parameter"""
+        return self._params[5]
+
+    # =============================
+    # Fourier Interface Implementation
+    # =============================
+
+    def cumulants(self, T: float) -> Cumulants:
+        """
+        Evaluate the cumulants of the model at a given time. This is useful e.g. to figure out integration bounds etc
+        during pricing
+        :param T: float, time to maturity (time at which cumulants are evaluated)
+        :return: Cumulants object
+        """
+        alpha_p, beta_p, lambda_p = self.alpha_p, self.beta_p, self.lambda_p
+        alpha_m, beta_m, lambda_m = self.alpha_m, self.beta_m, self.lambda_m
+
+        rn_drift = self.risk_neutral_log_drift()
+
+        def cumulants_gen(n: int):
+
+            return gamma(n - beta_p) * alpha_p / (lambda_p ** (n - beta_p)) + (
+                -1
+            ) ** n * gamma(n - beta_m) * alpha_m / (lambda_m ** (n - beta_m))
+
+        return Cumulants(
+            T=T,
+            rn_drift=rn_drift,
+            c1=T * (rn_drift + cumulants_gen(1)),
+            c2=T * cumulants_gen(2),
+            c4=T * cumulants_gen(3),
+        )
+
+    def symbol(self, xi: Union[float, np.ndarray]):
+        """
+        Levy symbol, uniquely defines Characteristic Function via: chf(T,xi) = exp(T*symbol(xi)),  for all T>=0
+        :param xi: np.ndarray or float, points in frequency domain
+        :return: np.ndarray or float, symbol evaluated at input points in frequency domain
+        """
+        alpha_p, beta_p, lambda_p = self.alpha_p, self.beta_p, self.lambda_p
+        alpha_m, beta_m, lambda_m = self.alpha_m, self.beta_m, self.lambda_m
+        rn_drift = self.risk_neutral_log_drift()
+
+        return (
+            1j * xi * rn_drift
+            + alpha_p
+            * gamma(-beta_p)
+            * ((lambda_p - 1j * xi) ** beta_p - lambda_p**beta_p)
+            + alpha_m
+            * gamma(-beta_m)
+            * ((lambda_m + 1j * xi) ** beta_m - lambda_m**beta_m)
+        )
+
+    def convexity_correction(self) -> float:
+        """
+        Computes the convexity correction for the Levy model, added to log process drift to ensure
+        risk neutrality
+        """
+        alpha_p, beta_p, lambda_p = self.alpha_p, self.beta_p, self.lambda_p
+        alpha_m, beta_m, lambda_m = self.alpha_m, self.beta_m, self.lambda_m
+        return -(
+            alpha_p * gamma(-beta_p) * ((lambda_p - 1) ** beta_p - lambda_p**beta_p)
+            + alpha_m * gamma(-beta_m) * ((lambda_m + 1) ** beta_m - lambda_m**beta_m)
+        )
+
+    # =============================
+    # Calibration Interface Implementation
+    # =============================
+
+    def num_params(self) -> int:
+        return 5
+
+    def param_bounds(self) -> Optional[List[Tuple]]:
+        return [(0, np.inf), (0, 1), (0, np.inf), (0, np.inf), (0, 1), (0, np.inf)]
+
+    def default_params(self) -> Optional[np.ndarray]:
+        return np.asarray([1, 0.5, 1, 1, 0.3, 1])