diff --git a/fypy/model/levy/BilateralGamma.py b/fypy/model/levy/BilateralGamma.py
index edf668a..9ca91f3 100644
--- a/fypy/model/levy/BilateralGamma.py
+++ b/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,21 +77,32 @@ 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]):
         """
@@ -90,12 +110,19 @@ def symbol(self, xi: Union[float, np.ndarray]):
         :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,7 +213,7 @@ 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]):
         """
@@ -184,7 +221,7 @@ def symbol(self, xi: Union[float, np.ndarray]):
         :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
diff --git a/fypy/model/levy/TemperedStable.py b/fypy/model/levy/TemperedStable.py
new file mode 100644
index 0000000..1cb9265
--- /dev/null
+++ b/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])
diff --git a/fypy/pricing/fourier/CarrMadanEuropeanPricer.py b/fypy/pricing/fourier/CarrMadanEuropeanPricer.py
index c79a9f0..0ae7ccf 100644
--- a/fypy/pricing/fourier/CarrMadanEuropeanPricer.py
+++ b/fypy/pricing/fourier/CarrMadanEuropeanPricer.py
@@ -6,32 +6,45 @@
 
 
 class CarrMadanEuropeanPricer(StrikesPricer):
-    def __init__(self,
-                 model: FourierModel,
-                 alpha: float = 0.75,
-                 eta: float = 0.1,
-                 N: int = 2 ** 9):
+    def __init__(
+        self, model: FourierModel, alpha: float = 0.75, eta: float = 0.1, N: int = 2**9
+    ):
+        """Carr-Madan method for Pricing European options under a Fourier model (i.e. using ChF)
+
+
+        Args:
+            model (FourierModel): FourierModel, model to price under
+            alpha (float, optional):  Defaults to 0.75.
+            eta (float, optional):  Defaults to 0.1.
+            N (int, optional):  Defaults to 2**9.
+        """
         self._model = model
-        self._alpha = alpha  # contour shift param (see Lee for recommendations)
+        self._alpha = alpha
         self._eta = eta
         self._N = N
         self._logS0 = np.log(self._model.spot())
 
-    def price_strikes_fill(self,
-                           T: float,
-                           K: np.ndarray,
-                           is_calls: np.ndarray,
-                           output: np.ndarray):
+    def price(self, T: float, K: float, is_call: bool) -> float:
+        """
+        Price a single strike of European option
+        :param T: float, time to maturity
+        :param K: float, strike of option
+        :param is_call: bool, indicator of if strike is call (true) or put (false)
+        :return: float, price of option
+        """
         lam = 2 * np.pi / (self._N * self._eta)
         b = self._N * lam / 2
 
         uv = np.arange(1, self._N + 1)  # TODO: check
         ku = -b + lam * (uv - 1)
         vj = (uv - 1) * self._eta
-        psij = self._chf(T=T, xi=vj - (self._alpha + 1) * 1j) \
-               / (self._alpha ** 2 + self._alpha - vj ** 2 + 1j * (2 * self._alpha + 1) * vj)
+
+        psij = self._chf(T=T, xi=vj - (self._alpha + 1) * 1j) / (
+            self._alpha**2 + self._alpha - vj**2 + 1j * (2 * self._alpha + 1) * vj
+        )
 
         disc = self._model.discountCurve(T)
+
         temp = (disc * self._eta / 3) * np.exp(1j * vj * b) * psij
         ind = np.zeros_like(uv)
         ind[0] = 1
@@ -40,10 +53,14 @@ def price_strikes_fill(self,
         Cku = np.real(np.exp(-self._alpha * ku) * fft(temp) / np.pi)
 
         logK = np.log(K)
-        # istrike = np.floor((logK + b) / lam + 1)
-        # xp = [ku(istrike) ku(istrike + 1)]
-        # yp = [Cku(istrike) Cku(istrike + 1)];
-        price = real(interp1(xp, yp, logK));
+        istrike = int(np.floor((logK + b) / lam + 1)) - 1
+
+        xp = [ku[istrike], ku[istrike + 1]]
+        yp = [Cku[istrike], Cku[istrike + 1]]
+        price = float(np.interp(logK, xp, yp))
+        if not is_call:
+            price = price - float((self._model.forwardCurve(T) * disc - K * disc))
+        return price
 
     def _chf(self, T: float, xi: np.ndarray):
         return self._model.chf(T, xi) * np.exp(1j * self._logS0 * xi)
diff --git a/fypy/pricing/fourier/HilbertEuropeanPricer.py b/fypy/pricing/fourier/HilbertEuropeanPricer.py
new file mode 100644
index 0000000..5dbeb26
--- /dev/null
+++ b/fypy/pricing/fourier/HilbertEuropeanPricer.py
@@ -0,0 +1,64 @@
+from fypy.termstructures.EquityForward import EquityForward
+from fypy.model.FourierModel import FourierModel
+import numpy as np
+from scipy.fft import fft
+from fypy.pricing.StrikesPricer import StrikesPricer
+
+
+class HilbertEuropeanPricer(StrikesPricer):
+    def __init__(
+        self,
+        model: FourierModel,
+        alpha: float = 0.75,
+        eta: float = 0.1,
+        N: int = 2**9,
+        Nh: int = 2**5,
+    ):
+        """Hilbert method for Pricing European options under a Fourier model (i.e. using ChF)
+
+
+        Args:
+            model (FourierModel): FourierModel, model to price under
+            alpha (float, optional):  Defaults to 0.75.
+            eta (float, optional):  Defaults to 0.1.
+            N (int, optional):  Defaults to 2**9.
+            Nh (int, optional):  Defaults to 2**9.
+        """
+        self._model = model
+        self._alpha = alpha
+        self._eta = eta
+        self._N = N
+        self._Nh = Nh
+        self._h = 2 * np.pi / Nh
+        self._logS0 = np.log(self._model.spot())
+
+    def price(self, T: float, K: float, is_call: bool) -> float:
+        """
+        Price a single strike of European option
+        :param T: float, time to maturity
+        :param K: float, strike of option
+        :param is_call: bool, indicator of if strike is call (true) or put (false)
+        :return: float, price of option
+        """
+        gridL = np.arange(-int(self._N / 2), 0)
+        gridR = -gridL[::-1]
+        H = (
+            np.sum(
+                self._g(self._h * gridL, T, K) * (np.cos(np.pi * gridL) - 1) / gridL
+                + self._g(self._h * gridR, T, K) * (np.cos(np.pi * gridR) - 1) / gridR
+            )
+            / np.pi
+        )
+        disc = self._model.discountCurve(T)
+        price = 0.5 * np.real(
+            self._model.forwardCurve(T) * disc - K * disc + 1j * disc * H
+        )
+        if not is_call:
+            price = price - (self._model.forwardCurve(T) * disc - K * disc)
+        return price
+
+    def _g(self, xi: np.ndarray, T: float, K: float):
+        return np.exp(-1j * xi * np.log(K / self._model.spot())) * (
+            self._model.spot() * self._model.chf(T, xi - 1j)
+            - K * self._model.chf(T, xi)
+        )
diff --git a/tests/pricing/fourier/Test_Carr_Madan_European.py b/tests/pricing/fourier/Test_Carr_Madan_European.py
new file mode 100644
index 0000000..5732c3a
--- /dev/null
+++ b/tests/pricing/fourier/Test_Carr_Madan_European.py
@@ -0,0 +1,77 @@
+import unittest
+
+from fypy.model.levy import VarianceGamma, NIG, CMGY, KouJD
+from fypy.model.levy.MertonJD import MertonJD
+from fypy.pricing.fourier.CarrMadanEuropeanPricer import CarrMadanEuropeanPricer
+from fypy.model.levy.BlackScholes import *
+from fypy.termstructures.DiscountCurve import DiscountCurve_ConstRate
+from fypy.termstructures.EquityForward import EquityForward
+
+
+class Test_Carr_Madan_European(unittest.TestCase):
+
+    def test_levy_models(self):
+        S0 = 100
+        r = 0.05
+        q = 0.01
+        T = 1
+
+        disc_curve = DiscountCurve_ConstRate(rate=r)
+        div_disc = DiscountCurve_ConstRate(rate=q)
+        fwd = EquityForward(S0=S0, discount=disc_curve, divDiscount=div_disc)
+
+        # 1) Black Scholes
+        model = BlackScholes(sigma=0.15, forwardCurve=fwd, discountCurve=disc_curve)
+        pricer = CarrMadanEuropeanPricer(model=model, N=2**20)
+        price = pricer.price(T=T, K=S0, is_call=True)
+
+        self.assertAlmostEqual(price, 7.94871378854164, 5)
+
+        # 2) Variance Gamma
+        model = VarianceGamma(
+            sigma=0.2, theta=0.1, nu=0.85, forwardCurve=fwd, discountCurve=disc_curve
+        )
+        pricer = CarrMadanEuropeanPricer(model=model, N=2**20)
+        price = pricer.price(T=T, K=S0, is_call=True)
+
+        self.assertAlmostEqual(price, 10.13935062748614, 5)
+
+        # 3) NIG
+        model = NIG(
+            alpha=15, beta=-5, delta=0.5, forwardCurve=fwd, discountCurve=disc_curve
+        )
+        pricer = CarrMadanEuropeanPricer(model=model, N=2**20)
+        price = pricer.price(T=T, K=S0, is_call=True)
+
+        self.assertAlmostEqual(price, 9.63000693130414, 5)
+
+        # 4) MJD
+        model = MertonJD(
+            sigma=0.12,
+            lam=0.4,
+            muj=-0.12,
+            sigj=0.18,
+            forwardCurve=fwd,
+            discountCurve=disc_curve,
+        )
+        pricer = CarrMadanEuropeanPricer(model=model, N=2**20)
+        price = pricer.price(T=T, K=S0, is_call=True)
+
+        self.assertAlmostEqual(price, 8.675684635426279, 5)
+
+        # 5) CGMY
+        model = CMGY(forwardCurve=fwd, discountCurve=disc_curve)
+        pricer = CarrMadanEuropeanPricer(model=model, N=2**20)
+        price = pricer.price(T=T, K=S0, is_call=True)
+        self.assertAlmostEqual(price, 5.80222163947386, 5)
+
+        # 6) Kou's Jump Diffusion
+        model = KouJD(forwardCurve=fwd, discountCurve=disc_curve)
+        pricer = CarrMadanEuropeanPricer(model=model, N=2**20)
+        price = pricer.price(T=T, K=S0, is_call=True)
+
+        self.assertAlmostEqual(price, 11.92430307601936, 5)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/pricing/fourier/Test_Hilbert_European.py b/tests/pricing/fourier/Test_Hilbert_European.py
new file mode 100644
index 0000000..157e784
--- /dev/null
+++ b/tests/pricing/fourier/Test_Hilbert_European.py
@@ -0,0 +1,77 @@
+import unittest
+
+from fypy.model.levy import VarianceGamma, NIG, CMGY, KouJD
+from fypy.model.levy.MertonJD import MertonJD
+from fypy.pricing.fourier.HilbertEuropeanPricer import HilbertEuropeanPricer
+from fypy.model.levy.BlackScholes import *
+from fypy.termstructures.DiscountCurve import DiscountCurve_ConstRate
+from fypy.termstructures.EquityForward import EquityForward
+
+
+class Test_Hilbert_European(unittest.TestCase):
+
+    def test_levy_models(self):
+        S0 = 100
+        r = 0.05
+        q = 0.01
+        T = 1
+
+        disc_curve = DiscountCurve_ConstRate(rate=r)
+        div_disc = DiscountCurve_ConstRate(rate=q)
+        fwd = EquityForward(S0=S0, discount=disc_curve, divDiscount=div_disc)
+
+        # 1) Black Scholes
+        model = BlackScholes(sigma=0.15, forwardCurve=fwd, discountCurve=disc_curve)
+        pricer = HilbertEuropeanPricer(model=model, N=2**20)
+        price = pricer.price(T=T, K=S0, is_call=True)
+
+        self.assertAlmostEqual(price, 7.94871378854164, 5)
+
+        # 2) Variance Gamma
+        model = VarianceGamma(
+            sigma=0.2, theta=0.1, nu=0.85, forwardCurve=fwd, discountCurve=disc_curve
+        )
+        pricer = HilbertEuropeanPricer(model=model, N=2**20)
+        price = pricer.price(T=T, K=S0, is_call=True)
+
+        self.assertAlmostEqual(price, 10.13935062748614, 5)
+
+        # 3) NIG
+        model = NIG(
+            alpha=15, beta=-5, delta=0.5, forwardCurve=fwd, discountCurve=disc_curve
+        )
+        pricer = HilbertEuropeanPricer(model=model, N=2**20)
+        price = pricer.price(T=T, K=S0, is_call=True)
+
+        self.assertAlmostEqual(price, 9.63000693130414, 5)
+
+        # 4) MJD
+        model = MertonJD(
+            sigma=0.12,
+            lam=0.4,
+            muj=-0.12,
+            sigj=0.18,
+            forwardCurve=fwd,
+            discountCurve=disc_curve,
+        )
+        pricer = HilbertEuropeanPricer(model=model, N=2**20)
+        price = pricer.price(T=T, K=S0, is_call=True)
+
+        self.assertAlmostEqual(price, 8.675684635426279, 5)
+
+        # 5) CGMY
+        model = CMGY(forwardCurve=fwd, discountCurve=disc_curve)
+        pricer = HilbertEuropeanPricer(model=model, N=2**20)
+        price = pricer.price(T=T, K=S0, is_call=True)
+        self.assertAlmostEqual(price, 5.80222163947386, 5)
+
+        # 6) Kou's Jump Diffusion
+        model = KouJD(forwardCurve=fwd, discountCurve=disc_curve)
+        pricer = HilbertEuropeanPricer(model=model, N=2**20)
+        price = pricer.price(T=T, K=S0, is_call=True)
+
+        self.assertAlmostEqual(price, 11.92430307601936, 5)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/pricing/fourier/suite.py b/tests/pricing/fourier/suite.py
index a625d22..f260e6b 100644
--- a/tests/pricing/fourier/suite.py
+++ b/tests/pricing/fourier/suite.py
@@ -2,15 +2,23 @@
 from tests.pricing.fourier.Test_Proj_European import Test_Proj_European
 from tests.pricing.fourier.Test_Lewis_European import Test_Lewis_European
 from tests.pricing.fourier.Test_GilPeleaz_European import Test_GilPeleaz_European
+from tests.pricing.fourier.Test_Carr_Madan_European import Test_Carr_Madan_European
+from tests.pricing.fourier.Test_Hilbert_European import Test_Hilbert_European
 
 
 def test_suite():
     suite = unittest.TestSuite()
-    for test in (Test_Proj_European, Test_Lewis_European, Test_GilPeleaz_European):
+    for test in (
+        Test_Proj_European,
+        Test_Lewis_European,
+        Test_GilPeleaz_European,
+        Test_Carr_Madan_European,
+        Test_Hilbert_European,
+    ):
         suite.addTest(unittest.TestLoader().loadTestsFromTestCase(test))
 
     return suite
 
 
-if __name__ == '__main__':
-    unittest.TextTestRunner(verbosity=2).run(test_suite())
\ No newline at end of file
+if __name__ == "__main__":
+    unittest.TextTestRunner(verbosity=2).run(test_suite())