Repository: PyFE/PyFENG (language: python)

Issue:
build a package



Discussion / hints:
We don't need to upload the build/dist folders into git.

Fix patch:
diff --git a/docs/source/conf.py b/docs/source/conf.py
index cd35a81..71e582a 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -12,17 +12,18 @@
 #
 import os
 import sys
-sys.path.insert(0, os.path.abspath('../..'))
+
+sys.path.insert(0, os.path.abspath("../.."))
 
 add_module_names = False
 # -- Project information -----------------------------------------------------
 
-project = 'pyfeng'
-copyright = '2021, Jaehyuk Choi'
-author = 'Jaehyuk Choi'
+project = "pyfeng"
+copyright = "2021, Jaehyuk Choi"
+author = "Jaehyuk Choi"
 
 # The full version, including alpha/beta/rc tags
-release = '0.1.1'
+release = "0.1.1"
 
 
 # -- General configuration ---------------------------------------------------
@@ -38,15 +39,15 @@
     "sphinx.ext.mathjax",
     "sphinx.ext.viewcode",
     "sphinx.ext.napoleon",
-    'm2r2',
+    "m2r2",
 ]
 
 # source_suffix = '.rst'
-source_suffix = ['.rst', '.md']
+source_suffix = [".rst", ".md"]
 
 
 # Add any paths that contain templates here, relative to this directory.
-templates_path = ['_templates']
+templates_path = ["_templates"]
 
 # List of patterns, relative to source directory, that match files and
 # directories to ignore when looking for source files.
@@ -64,4 +65,4 @@
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
-html_static_path = ['_static']
\ No newline at end of file
+html_static_path = ["_static"]
diff --git a/pyfeng/__init__.py b/pyfeng/__init__.py
index b7e04cf..6dfa640 100644
--- a/pyfeng/__init__.py
+++ b/pyfeng/__init__.py
@@ -1,15 +1,33 @@
-from .norm import Norm  # the order is sensitive because of `price_barrier` method. Put it before .bsm
+from .norm import (
+    Norm,
+)  # the order is sensitive because of `price_barrier` method. Put it before .bsm
 from .bsm import Bsm, BsmDisp
 from .cev import Cev
 from .gamma import InvGam, InvGauss
-from .sabr import SabrHagan2002, SabrNorm, SabrLorig2017, SabrChoiWu2021H, SabrChoiWu2021P
+from .sabr import (
+    SabrHagan2002,
+    SabrNorm,
+    SabrLorig2017,
+    SabrChoiWu2021H,
+    SabrChoiWu2021P,
+)
 from .sabr_int import SabrUncorrChoiWu2021
 from .sabr_mc import SabrCondMc
 from .nsvh import Nsvh1, NsvhMc
-from .multiasset import BsmSpreadKirk, BsmSpreadBjerksund2014, NormBasket, NormSpread, BsmBasketLevy1992, BsmMax2, \
-    BsmBasketMilevsky1998, BsmBasket1Bm, BsmBasketLowerBound, BsmBasketJsu
+from .multiasset import (
+    BsmSpreadKirk,
+    BsmSpreadBjerksund2014,
+    NormBasket,
+    NormSpread,
+    BsmBasketLevy1992,
+    BsmMax2,
+    BsmBasketMilevsky1998,
+    BsmBasket1Bm,
+    BsmBasketLowerBound,
+    BsmBasketJsu,
+)
 from .multiasset_mc import BsmNdMc, NormNdMc
+
 # Basket-Asian from ASP 2021
 from .multiasset_Ju2002 import BsmBasketAsianJu2002, BsmContinuousAsianJu2002
 from .asian import BsmAsianLinetsky2004, BsmAsianJsu
-from . import _ex as ex
diff --git a/pyfeng/asian.py b/pyfeng/asian.py
index 289a8a4..b8f756b 100644
--- a/pyfeng/asian.py
+++ b/pyfeng/asian.py
@@ -48,20 +48,20 @@ def moments(self, spot, texp, n=1):
 
         """
         lam = self.sigma[0]
-        v = (self.intr - self.divr - lam**2 / 2) / lam
+        v = (self.intr - self.divr - lam ** 2 / 2) / lam
         beta = v / lam
 
         ds = []
-        for j in range(n+1):
+        for j in range(n + 1):
             item0 = 2 ** n
-            for i in range(n+1):
+            for i in range(n + 1):
                 if i != j:
-                    item0 *= ((beta + j)**2 - (beta + i)**2)**(-1)
+                    item0 *= ((beta + j) ** 2 - (beta + i) ** 2) ** (-1)
             ds.append(item0)
         item1 = 0
-        for i in range(n+1):
-            item1 += (ds[i] * np.exp((lam**2 * i**2 / 2 + lam * i * v) * texp))
-        moment = (spot / texp)**n * math.factorial(n) / lam ** (2 * n) * item1
+        for i in range(n + 1):
+            item1 += ds[i] * np.exp((lam ** 2 * i ** 2 / 2 + lam * i * v) * texp)
+        moment = (spot / texp) ** n * math.factorial(n) / lam ** (2 * n) * item1
 
         return moment
 
@@ -84,8 +84,8 @@ def moment_mvsk(self, spot, texp):
         m1, m2, m3, m4 = moments[0], moments[1], moments[2], moments[3]
         mu = m1
         var = m2 - m1 ** 2
-        skew = (m3 - m1 ** 3 - 3 * m2 * m1 + 3 * m1 ** 3) / var**(3 / 2)
-        kurt = (m4 - 3 * m1 ** 4 - 4 * m3 * m1 + 6 * m2 * m1 ** 2) / var**2
+        skew = (m3 - m1 ** 3 - 3 * m2 * m1 + 3 * m1 ** 3) / var ** (3 / 2)
+        kurt = (m4 - 3 * m1 ** 4 - 4 * m3 * m1 + 6 * m2 * m1 ** 2) / var ** 2
 
         return mu, var, skew, kurt
 
@@ -106,16 +106,13 @@ def price(self, strike, spot, texp, cp=1):
         mu, var, skew, kurt = self.moment_mvsk(spot, texp)
 
         m = nsvh.Nsvh1(sigma=self.sigma)
-        m.calibrate_vsk(var, skew, kurt-3, texp, setval=True)
+        m.calibrate_vsk(var, skew, kurt - 3, texp, setval=True)
         price = m.price(strike, mu, texp, cp)
 
         return df * price
 
 
-
-
 class BsmAsianLinetsky2004:
-    
     def __init__(self, intr, divr, vol, texp, strike, spot, b, call):
         self.intr = intr
         self.divr = divr
@@ -125,21 +122,21 @@ def __init__(self, intr, divr, vol, texp, strike, spot, b, call):
         self.spot = spot
         self.b = b
         self.call = call
-        self.nu = 2*(self.intr-self.divr)/(self.vol**2)-1
-        self.tau = self.vol**2*self.texp/4
-        self.k = self.tau*self.strike/self.spot
-        
-    def find_zeros_real(self,nu):
+        self.nu = 2 * (self.intr - self.divr) / (self.vol ** 2) - 1
+        self.tau = self.vol ** 2 * self.texp / 4
+        self.k = self.tau * self.strike / self.spot
+
+    def find_zeros_real(self, nu):
         eigenval = []
         if nu <= -2:
-            for n in range(2,math.floor(np.abs(nu)/2)+1+1):
-                eigenvalue = np.abs(nu)-2*n+2
+            for n in range(2, math.floor(np.abs(nu) / 2) + 1 + 1):
+                eigenvalue = np.abs(nu) - 2 * n + 2
                 eigenval.append(eigenvalue)
         else:
             pass
         return np.array(eigenval)
 
-    def positive_eigenvalue(self,eigval):
+    def positive_eigenvalue(self, eigval):
         eigval_p = []
         for i in range(len(eigval)):
             real_i = np.array(eigval[i], dtype=complex).real[0]
@@ -149,41 +146,49 @@ def positive_eigenvalue(self,eigval):
                 pass
         return eigval_p
 
-    def find_zeros_imag(self,n,nu):
-        p = sympy.symbols('p')
+    def find_zeros_imag(self, n, nu):
+        p = sympy.symbols("p")
         eigenval = []
-        for j in range(1,n):
-            p_tilde = sympy.solve(p*(sympy.log(4*p)-1)-2*np.pi*(j+nu/4-1/2),p)
+        for j in range(1, n):
+            p_tilde = sympy.solve(
+                p * (sympy.log(4 * p) - 1) - 2 * np.pi * (j + nu / 4 - 1 / 2), p
+            )
             eigenval.append(p_tilde)
         eigenval = self.positive_eigenvalue(np.real(eigenval))
         eigenval = np.array(eigenval).transpose()
         return eigenval
-    
-    def eta_q(self,nu,eigenval,b):
-        f = lambda x: m.whitw((1-nu)/2, x/2, 1/(2*b))
+
+    def eta_q(self, nu, eigenval, b):
+        f = lambda x: m.whitw((1 - nu) / 2, x / 2, 1 / (2 * b))
         return complex(-derivative(f, eigenval, dx=1e-12))
 
-    def xi_p(self,nu,eigenval,b):
-        f = lambda x: m.whitw((1-nu)/2,complex(0,x/2),1/(2*b))
-        return complex(derivative(f,eigenval,dx=1e-12))
-        
-    def price_element_imag(self,nu,tau,p_value):
-        p1 = m.exp(-(nu**2+p_value**2)*tau/2)
-        p2 = (p_value*m.gamma(complex(nu/2,p_value/2)))/(4*self.xi_p(nu=nu,eigenval=p_value,b=self.b)*m.gamma(complex(1,p_value)))
-        const = (2*self.k)**((nu+3)/2)*m.exp(-1/(4*self.k))
-        w1 = m.whitw(-(nu+3)/2,complex(0,p_value/2),1/(2*self.k))
-        m1 = m.whitm((1-nu)/2,complex(0,p_value/2),1/(2*self.b))
+    def xi_p(self, nu, eigenval, b):
+        f = lambda x: m.whitw((1 - nu) / 2, complex(0, x / 2), 1 / (2 * b))
+        return complex(derivative(f, eigenval, dx=1e-12))
+
+    def price_element_imag(self, nu, tau, p_value):
+        p1 = m.exp(-(nu ** 2 + p_value ** 2) * tau / 2)
+        p2 = (p_value * m.gamma(complex(nu / 2, p_value / 2))) / (
+            4
+            * self.xi_p(nu=nu, eigenval=p_value, b=self.b)
+            * m.gamma(complex(1, p_value))
+        )
+        const = (2 * self.k) ** ((nu + 3) / 2) * m.exp(-1 / (4 * self.k))
+        w1 = m.whitw(-(nu + 3) / 2, complex(0, p_value / 2), 1 / (2 * self.k))
+        m1 = m.whitm((1 - nu) / 2, complex(0, p_value / 2), 1 / (2 * self.b))
         return p1 * p2 * const * w1 * m1
-        
-    def price_element_real(self,nu,tau,q_value):
-        p1 = m.exp(-(nu**2-q_value**2)*tau/2)
-        p2 = (q_value*m.gamma((nu+q_value)/2))/(4*self.eta_q(nu=nu,eigenval=q_value,b=self.b)*m.gamma(1+q_value))
-        const = (2*self.k)**((nu+3)/2)*m.exp(-1/(4*self.k))
-        w1 = m.whitw(-(nu+3)/2,q_value/2,1/(2*self.k))
-        m1 = m.whitm((1-nu)/2,q_value/2,1/(2*self.b))
+
+    def price_element_real(self, nu, tau, q_value):
+        p1 = m.exp(-(nu ** 2 - q_value ** 2) * tau / 2)
+        p2 = (q_value * m.gamma((nu + q_value) / 2)) / (
+            4 * self.eta_q(nu=nu, eigenval=q_value, b=self.b) * m.gamma(1 + q_value)
+        )
+        const = (2 * self.k) ** ((nu + 3) / 2) * m.exp(-1 / (4 * self.k))
+        w1 = m.whitw(-(nu + 3) / 2, q_value / 2, 1 / (2 * self.k))
+        m1 = m.whitm((1 - nu) / 2, q_value / 2, 1 / (2 * self.b))
         return p1 * p2 * const * w1 * m1
-        
-    def exact_asian_price(self,n_eig):
+
+    def exact_asian_price(self, n_eig):
         intr = self.intr
         divr = self.divr
         vol = self.vol
@@ -195,26 +200,31 @@ def exact_asian_price(self,n_eig):
         nu = self.nu
         tau = self.tau
         k = self.k
-        p = self.find_zeros_imag(n=n_eig+1,nu=nu)
+        p = self.find_zeros_imag(n=n_eig + 1, nu=nu)
         q = self.find_zeros_real(nu=nu)
-        
+
         imaginary_terms = []
         for i in range(len(p)):
-            imaginary_term = self.price_element_imag(nu=nu,tau=tau,p_value=p[i])
+            imaginary_term = self.price_element_imag(nu=nu, tau=tau, p_value=p[i])
             imaginary_terms.append(complex(imaginary_term))
-        
+
         real_terms = []
         for i in range(len(q)):
-            real_term = self.price_element_real(nu=nu,tau=tau,q_value=q[i])
+            real_term = self.price_element_real(nu=nu, tau=tau, q_value=q[i])
             real_terms.append(complex(real_term))
-            
-        P = np.sum(imaginary_terms)+np.sum(real_terms)
-        put_price = np.exp(-intr*texp)*(4*spot/(texp*vol**2))* P
-        call_price = put_price + (spot*(np.exp(-divr*texp)-np.exp(-intr*texp))/((intr-divr)*texp) - strike*np.exp(-intr*texp)) 
-        
+
+        P = np.sum(imaginary_terms) + np.sum(real_terms)
+        put_price = np.exp(-intr * texp) * (4 * spot / (texp * vol ** 2)) * P
+        call_price = put_price + (
+            spot
+            * (np.exp(-divr * texp) - np.exp(-intr * texp))
+            / ((intr - divr) * texp)
+            - strike * np.exp(-intr * texp)
+        )
+
         if call == False:
-            return(put_price.real)
+            return put_price.real
         elif call == True:
-            return(call_price.real)
+            return call_price.real
         else:
-            return("Wrong")
+            return "Wrong"
diff --git a/pyfeng/bsm.py b/pyfeng/bsm.py
index 4af7940..b8a8314 100644
--- a/pyfeng/bsm.py
+++ b/pyfeng/bsm.py
@@ -28,7 +28,9 @@ class Bsm(opt.OptAnalyticABC):
     """
 
     @staticmethod
-    def price_formula(strike, spot, sigma, texp, cp=1, intr=0.0, divr=0.0, is_fwd=False):
+    def price_formula(
+        strike, spot, sigma, texp, cp=1, intr=0.0, divr=0.0, is_fwd=False
+    ):
         """
         Black-Scholes-Merton model call/put option pricing formula (static method)
 
@@ -53,8 +55,8 @@ def price_formula(strike, spot, sigma, texp, cp=1, intr=0.0, divr=0.0, is_fwd=Fa
 
         # don't directly compute d1 just in case sigma_std is infty
         d1 = np.log(fwd / strike) / sigma_std
-        d2 = d1 - 0.5*sigma_std
-        d1 += 0.5*sigma_std
+        d2 = d1 - 0.5 * sigma_std
+        d1 += 0.5 * sigma_std
 
         cp = np.array(cp)
         price = fwd * spst.norm.cdf(cp * d1) - strike * spst.norm.cdf(cp * d2)
@@ -68,7 +70,9 @@ def vega(self, strike, spot, texp, cp=1):
         sigma_std = np.maximum(self.sigma * np.sqrt(texp), np.finfo(float).eps)
         d1 = np.log(fwd / strike) / sigma_std + 0.5 * sigma_std
 
-        vega = df * fwd * spst.norm.pdf(d1) * np.sqrt(texp)  # formula according to wikipedia
+        vega = (
+            df * fwd * spst.norm.pdf(d1) * np.sqrt(texp)
+        )  # formula according to wikipedia
         return vega
 
     def delta(self, strike, spot, texp, cp=1):
@@ -95,25 +99,29 @@ def gamma(self, strike, spot, texp, cp=1):
 
         fwd, df, divf = self._fwd_factor(spot, texp)
 
-        sigma_std = np.maximum(self.sigma * np.sqrt(texp), 100*np.finfo(float).eps)
+        sigma_std = np.maximum(self.sigma * np.sqrt(texp), 100 * np.finfo(float).eps)
         d1 = np.log(fwd / strike) / sigma_std + 0.5 * sigma_std
 
-        gamma = df * spst.norm.pdf(d1)/fwd/sigma_std  # formula according to wikipedia
+        gamma = (
+            df * spst.norm.pdf(d1) / fwd / sigma_std
+        )  # formula according to wikipedia
         if not self.is_fwd:
-            gamma *= (divf/df)**2
+            gamma *= (divf / df) ** 2
         return gamma
 
     def theta(self, strike, spot, texp, cp=1):
 
         fwd, df, divf = self._fwd_factor(spot, texp)
 
-        sigma_std = np.maximum(self.sigma * np.sqrt(texp), 100*np.finfo(float).eps)
+        sigma_std = np.maximum(self.sigma * np.sqrt(texp), 100 * np.finfo(float).eps)
         d1 = np.log(fwd / strike) / sigma_std
-        d2 = d1 - 0.5*sigma_std
-        d1 += 0.5*sigma_std
+        d2 = d1 - 0.5 * sigma_std
+        d1 += 0.5 * sigma_std
 
         # still not perfect; need to consider the derivative w.r.t. divr and is_fwd = True
-        theta = -0.5*spst.norm.pdf(d1)*fwd*self.sigma/np.sqrt(texp) - cp*self.intr*strike*spst.norm.cdf(cp*d2)
+        theta = -0.5 * spst.norm.pdf(d1) * fwd * self.sigma / np.sqrt(
+            texp
+        ) - cp * self.intr * strike * spst.norm.cdf(cp * d2)
         theta *= df
         return theta
 
@@ -145,31 +153,35 @@ def _impvol_newton(self, price, strike, spot, texp, cp=1, setval=False):
 
         # Exclude optoin price below intrinsic value or above max value (1 for call or k for put)
         # ind_solve can be scalar or array. scalar can be fine in np.abs(p_err[ind_solve])
-        ind_solve = (price_std - p_min > Bsm.IMPVOL_TOL) & (p_max - price_std > Bsm.IMPVOL_TOL)
+        ind_solve = (price_std - p_min > Bsm.IMPVOL_TOL) & (
+            p_max - price_std > Bsm.IMPVOL_TOL
+        )
 
         # initial guess = inflection point in sigma (volga=0)
-        _sigma = np.ones_like(ind_solve) * np.sqrt(2*np.abs(np.log(strike_std))/texp)
+        _sigma = np.ones_like(ind_solve) * np.sqrt(
+            2 * np.abs(np.log(strike_std)) / texp
+        )
 
         bsm_model.sigma = _sigma
         p_err = bsm_model.price(strike_std, 1.0, texp, cp) - price_std
-        #print(np.sign(p_err), _sigma)
-        
+        # print(np.sign(p_err), _sigma)
+
         if np.any(ind_solve):
             for k in range(32):  # usually iteration ends less than 10
                 vega = bsm_model.vega(strike_std, 1.0, texp, cp)
-                _sigma -= p_err/vega
+                _sigma -= p_err / vega
                 bsm_model.sigma = _sigma
                 p_err = bsm_model.price(strike_std, 1.0, texp, cp) - price_std
                 p_err_max = np.amax(np.abs(p_err[ind_solve]))
-                #print(k, p_err_max, _sigma)
+                # print(k, p_err_max, _sigma)
 
                 # ignore the error of the elements with ind_solve = False
                 if p_err_max < Bsm.IMPVOL_TOL:
                     break
-            #print(k)
+            # print(k)
 
             if p_err_max >= Bsm.IMPVOL_TOL:
-                warn_msg = f'impvol_newton did not converged within {k} iterations: max error = {p_err_max}'
+                warn_msg = f"impvol_newton did not converged within {k} iterations: max error = {p_err_max}"
                 warnings.warn(warn_msg, Warning)
 
         # Put Nan for the out-of-bound option prices
@@ -177,11 +189,17 @@ def _impvol_newton(self, price, strike, spot, texp, cp=1, setval=False):
 
         # Though not error is above tolerance, if the price is close to min or max, set 0 or inf
         _sigma = np.where(
-            (np.abs(p_err) >= Bsm.IMPVOL_TOL) & (np.abs(price_std - p_min) <= Bsm.IMPVOL_TOL),
-            0, _sigma)
+            (np.abs(p_err) >= Bsm.IMPVOL_TOL)
+            & (np.abs(price_std - p_min) <= Bsm.IMPVOL_TOL),
+            0,
+            _sigma,
+        )
         _sigma = np.where(
-            (np.abs(p_err) >= Bsm.IMPVOL_TOL) & (np.abs(price_std - p_max) <= Bsm.IMPVOL_TOL),
-            np.inf, _sigma)
+            (np.abs(p_err) >= Bsm.IMPVOL_TOL)
+            & (np.abs(price_std - p_max) <= Bsm.IMPVOL_TOL),
+            np.inf,
+            _sigma,
+        )
 
         if scalar_output:
             _sigma = np.asscalar(_sigma)
@@ -194,7 +212,7 @@ def _impvol_newton(self, price, strike, spot, texp, cp=1, setval=False):
     #### Class method
     impvol = _impvol_newton
 
-    def vol_smile(self, strike, spot, texp, cp=1, model='norm'):
+    def vol_smile(self, strike, spot, texp, cp=1, model="norm"):
         """
         Equivalent volatility smile for a given model
 
@@ -208,36 +226,46 @@ def vol_smile(self, strike, spot, texp, cp=1, model='norm'):
         Returns:
             volatility smile under the specified model
         """
-        if model.lower() == 'bsm':
+        if model.lower() == "bsm":
             return self.sigma * np.ones_like(strike + spot + texp + cp)
-        if model.lower() == 'norm':
+        if model.lower() == "norm":
             price = self.price(strike, spot, texp, cp=cp)
             return norm.Norm(None).impvol(price, strike, spot, texp, cp=cp)
-        elif model.lower() == 'norm-approx' or model.lower() == 'norm-grunspan':
+        elif model.lower() == "norm-approx" or model.lower() == "norm-grunspan":
             fwd, _, _ = self._fwd_factor(spot, texp)
             kk = strike / fwd
             lnk = np.log(kk)
-            if model.lower() == 'norm-approx':
-                return self.sigma * fwd * np.sqrt(kk) * (1 + lnk**2/24) / (1 + self.sigma**2 * texp/24)
+            if model.lower() == "norm-approx":
+                return (
+                    self.sigma
+                    * fwd
+                    * np.sqrt(kk)
+                    * (1 + lnk ** 2 / 24)
+                    / (1 + self.sigma ** 2 * texp / 24)
+                )
             else:
-                with np.errstate(divide='ignore', invalid='ignore'):
-                    term1 = np.where(np.fabs(lnk) > 1e-8, (kk - 1)/lnk, 2/(3 - kk))
-                    term2 = np.where(np.fabs(lnk) > 1e-8, (np.log(term1) - lnk/2) / lnk**2, 1/24)
-                return self.sigma * fwd * term1 * (1 - term2 * self.sigma**2 * texp)
+                with np.errstate(divide="ignore", invalid="ignore"):
+                    term1 = np.where(np.fabs(lnk) > 1e-8, (kk - 1) / lnk, 2 / (3 - kk))
+                    term2 = np.where(
+                        np.fabs(lnk) > 1e-8,
+                        (np.log(term1) - lnk / 2) / lnk ** 2,
+                        1 / 24,
+                    )
+                return self.sigma * fwd * term1 * (1 - term2 * self.sigma ** 2 * texp)
         else:
-            raise ValueError(f'Unknown model: {model}')
+            raise ValueError(f"Unknown model: {model}")
 
     def _price_suboptimal(self, strike, spot, texp, cp=1, strike2=None):
         strike2 = strike if strike2 is None else strike2
         fwd, df, _ = self._fwd_factor(spot, texp)
 
         sigma_std = np.maximum(self.sigma * np.sqrt(texp), np.finfo(float).eps)
-        d1 = np.log(fwd/strike2)/sigma_std
-        d2 = d1 - 0.5*sigma_std
-        d1 += 0.5*sigma_std
+        d1 = np.log(fwd / strike2) / sigma_std
+        d2 = d1 - 0.5 * sigma_std
+        d1 += 0.5 * sigma_std
 
         price = fwd * spst.norm.cdf(cp * d1) - strike * spst.norm.cdf(cp * d2)
-        price *= cp*df
+        price *= cp * df
         return price
 
     def _barrier_params(self, barrier, spot):
@@ -251,7 +279,9 @@ def _barrier_params(self, barrier, spot):
         Returns:
             barrier option pricing parameters (psi, spot_mirror)
         """
-        psi = np.power(barrier / spot, 2 * (self.intr - self.divr) / self.sigma ** 2 - 1)
+        psi = np.power(
+            barrier / spot, 2 * (self.intr - self.divr) / self.sigma ** 2 - 1
+        )
         spot_reflected = barrier ** 2 / spot
         return psi, spot_reflected
 
@@ -273,37 +303,45 @@ def price_barrier(self, strike, barrier, spot, texp, cp=1, io=-1):
 
         psi, spot_reflected = self._barrier_params(barrier, spot)
 
-        '''
+        """
         `mirror_sign` is +1/-1 if call/put remains same/flipped in the reflection principle
         +1 if (barrier < spot AND call) or (barrier > spot AND put), -1 otherwise
-        '''
+        """
         mirror_sign = np.where(barrier < spot, 1, -1) * cp
-        '''
+        """
         This is a trick to handle a trivial case:
         Knock-out call with spot < strike < barrier is worth zero.
         Knock-out put with barrier < strike < spot is worth zero.
         Without explicit adjustment, Knock-out price is negative, Knock-in price is higher than vanilla.
         In both scenario (mirror_sign = -1), we set strike = barrier, which will do the adjustment.
-        '''
+        """
         barrier = np.where(
-            mirror_sign > 0, barrier,
-            cp * np.maximum(cp*strike, cp*barrier)
+            mirror_sign > 0, barrier, cp * np.maximum(cp * strike, cp * barrier)
         )
 
         p_euro1 = np.where(
-            mirror_sign > 0, 0,
-            self._price_suboptimal(strike, spot, texp, cp=cp, strike2=barrier))
+            mirror_sign > 0,
+            0,
+            self._price_suboptimal(strike, spot, texp, cp=cp, strike2=barrier),
+        )
 
-        p_euro2 = self._price_suboptimal(strike, spot_reflected, texp, cp=mirror_sign*cp)
+        p_euro2 = self._price_suboptimal(
+            strike, spot_reflected, texp, cp=mirror_sign * cp
+        )
         p_euro2 -= np.where(
-            mirror_sign > 0, 0,
-            self._price_suboptimal(strike, spot_reflected, texp,
-                                   cp=mirror_sign*cp, strike2=barrier))
+            mirror_sign > 0,
+            0,
+            self._price_suboptimal(
+                strike, spot_reflected, texp, cp=mirror_sign * cp, strike2=barrier
+            ),
+        )
 
         p = p_euro1 + psi * p_euro2  # knock-in price
         p = np.where(
-            io > 0, p,  # knock-in type
-            self._price_suboptimal(strike, spot, texp, cp=cp) - p)
+            io > 0,
+            p,  # knock-in type
+            self._price_suboptimal(strike, spot, texp, cp=cp) - p,
+        )
 
         return p
 
@@ -320,10 +358,10 @@ def moments_vsk(self, texp=1):
         References:
             https://en.wikipedia.org/wiki/Log-normal_distribution
         """
-        ww = np.exp(texp*self.sigma**2)
+        ww = np.exp(texp * self.sigma ** 2)
         var = ww - 1
-        skew = (ww + 2)*np.sqrt(ww-1)
-        exkurt = ww**2*(ww*(ww + 2) + 3) - 6  # ww**4 + 2*ww**3 + 3*ww - 6
+        skew = (ww + 2) * np.sqrt(ww - 1)
+        exkurt = ww ** 2 * (ww * (ww + 2) + 3) - 6  # ww**4 + 2*ww**3 + 3*ww - 6
         return var, skew, exkurt
 
 
@@ -352,7 +390,7 @@ class BsmDisp(smile.OptSmileABC, Bsm):
     beta = 1  # equivalent to Black-Scholes
     pivot = 0
     sigma_disp = None
-    #_m_bsm = None
+    # _m_bsm = None
 
     def __init__(self, sigma, beta=1, pivot=0, *args, **kwargs):
         """
@@ -366,7 +404,7 @@ def __init__(self, sigma, beta=1, pivot=0, *args, **kwargs):
         """
         self.pivot = pivot
         self.beta = beta
-        #self._m_bsm = Bsm(sigma=beta*sigma, *args, **kwargs)
+        # self._m_bsm = Bsm(sigma=beta*sigma, *args, **kwargs)
         super().__init__(sigma, *args, **kwargs)
 
     @property
@@ -387,7 +425,7 @@ def disp_spot(self, spot):
         Returns:
             Displaces spot
         """
-        return self.beta*spot + (1-self.beta)*self.pivot
+        return self.beta * spot + (1 - self.beta) * self.pivot
 
     def disp_strike(self, strike, texp):
         """
@@ -400,12 +438,12 @@ def disp_strike(self, strike, texp):
         Returns:
             Displace strike
         """
-        return self.beta*strike + (1-self.beta)*self.forward(self.pivot, texp)
+        return self.beta * strike + (1 - self.beta) * self.forward(self.pivot, texp)
 
     def price(self, strike, spot, texp, cp=1):
         spot = self.disp_spot(spot)
         strike = self.disp_strike(strike, texp)
-        return (1/self.beta) * super().price(strike, spot, texp, cp=cp)
+        return (1 / self.beta) * super().price(strike, spot, texp, cp=cp)
 
     def delta(self, strike, spot, texp, cp=1):
         spot = self.disp_spot(spot)
@@ -431,18 +469,19 @@ def gamma(self, strike, spot, texp, cp=1):
     def theta(self, strike, spot, texp, cp=1):
         spot = self.disp_spot(spot)
         strike = self.disp_strike(strike, texp)
-        return (1/self.beta)*super().theta(
-            strike, spot, texp, cp=cp)
+        return (1 / self.beta) * super().theta(strike, spot, texp, cp=cp)
 
     def impvol(self, price_in, strike, spot, texp, cp=1, setval=False):
         spot = self.disp_spot(spot)
         strike = self.disp_strike(strike, texp)
-        sigma = (1/self.beta)*super().impvol(self.beta*price_in, strike, spot, texp, cp=cp, setval=False)
+        sigma = (1 / self.beta) * super().impvol(
+            self.beta * price_in, strike, spot, texp, cp=cp, setval=False
+        )
         if setval:
             self.sigma = sigma
         return sigma
 
-    def vol_smile(self, strike, spot, texp, cp=1, model='bsm'):
+    def vol_smile(self, strike, spot, texp, cp=1, model="bsm"):
         """
         Equivalent volatility smile for a given model
 
@@ -456,13 +495,19 @@ def vol_smile(self, strike, spot, texp, cp=1, model='bsm'):
         Returns:
             volatility smile under the specified model
         """
-        if model.lower() == 'norm-approx':
+        if model.lower() == "norm-approx":
             fwdd = self.forward(self.disp_spot(spot), texp)
             kkd = self.disp_strike(strike, texp) / fwdd
             lnkd = np.log(kkd)
             # self.sigma actually means self.beta * self._sigma
-            vol = self.sigma_disp * fwdd * np.sqrt(kkd) * (1 + lnkd**2/24) / (1 + self.sigma**2 * texp/24)
-        elif model.lower() == 'bsm-approx':
+            vol = (
+                self.sigma_disp
+                * fwdd
+                * np.sqrt(kkd)
+                * (1 + lnkd ** 2 / 24)
+                / (1 + self.sigma ** 2 * texp / 24)
+            )
+        elif model.lower() == "bsm-approx":
             fwd = self.forward(spot, texp)
             kk = strike / fwd
             lnk = np.log(kk)
@@ -472,19 +517,25 @@ def vol_smile(self, strike, spot, texp, cp=1, model='bsm'):
             lnkd = np.log(kkd)
 
             # self.sigma actually means self.beta * self.sigma_disp
-            vol = self.sigma_disp * (fwdd/fwd) * np.sqrt(kkd/kk)
-            vol *= (1 + lnkd**2/24) / (1 + lnk**2/24) * (1 + vol**2 * texp/24) / (1 + self.sigma**2 * texp/24)
+            vol = self.sigma_disp * (fwdd / fwd) * np.sqrt(kkd / kk)
+            vol *= (
+                (1 + lnkd ** 2 / 24)
+                / (1 + lnk ** 2 / 24)
+                * (1 + vol ** 2 * texp / 24)
+                / (1 + self.sigma ** 2 * texp / 24)
+            )
         else:
             vol = super().vol_smile(strike, spot, texp, model=model, cp=cp)
 
         return vol
 
     def price_barrier(self, strike, barrier, spot, *args, **kwargs):
-        return (1/self.beta)*self.bsm_model.price_barrier(
-            self.disp(strike), self.disp(barrier), self.disp(spot), *args, **kwargs)
+        return (1 / self.beta) * self.bsm_model.price_barrier(
+            self.disp(strike), self.disp(barrier), self.disp(spot), *args, **kwargs
+        )
 
     def moments_vsk(self, texp=1):
         rv = super().moments_vsk(self, texp)
         rv[0] /= self.beta
-        rv[1] /= self.beta**3
-        rv[2] /= self.beta**4
+        rv[1] /= self.beta ** 3
+        rv[2] /= self.beta ** 4
diff --git a/pyfeng/cev.py b/pyfeng/cev.py
index e2845e6..5925141 100644
--- a/pyfeng/cev.py
+++ b/pyfeng/cev.py
@@ -19,6 +19,7 @@ class Cev(opt.OptAnalyticABC, smile.OptSmileABC, smile.MassZeroABC):
         >>> m.price(np.arange(80, 121, 10), 100, 1.2)
         array([16.11757214, 10.00786871,  5.64880408,  2.89028476,  1.34128656])
     """
+
     sigma = None
     beta = 0.5
     is_bsm_sigma = False
@@ -37,7 +38,7 @@ def __init__(self, sigma, beta=0.5, intr=0.0, divr=0.0, is_fwd=False):
 
     def params_kw(self):
         params = super().params_kw()
-        extra = {'beta': self.beta}
+        extra = {"beta": self.beta}
         return {**params, **extra}  # Py 3.9, params | extra
 
     def mass_zero(self, spot, texp, log=False):
@@ -45,13 +46,20 @@ def mass_zero(self, spot, texp, log=False):
 
         betac = 1.0 - self.beta
         a = 0.5 / betac
-        sigma_std = np.maximum(self.sigma / np.power(fwd, betac) * np.sqrt(texp), np.finfo(float).eps)
+        sigma_std = np.maximum(
+            self.sigma / np.power(fwd, betac) * np.sqrt(texp), np.finfo(float).eps
+        )
         x = 0.5 / np.square(betac * sigma_std)
 
         if log:
-            log_mass = (a - 1)*np.log(x) - x - np.log(spsp.gamma(a))
-            log_mass += np.log(1 + (a-1)/x*(1 + (a-2)/x*(1 + (a-3)/x*(1 + (a-4)/x))))
-            with np.errstate(divide='ignore'):
+            log_mass = (a - 1) * np.log(x) - x - np.log(spsp.gamma(a))
+            log_mass += np.log(
+                1
+                + (a - 1)
+                / x
+                * (1 + (a - 2) / x * (1 + (a - 3) / x * (1 + (a - 4) / x)))
+            )
+            with np.errstate(divide="ignore"):
                 log_mass = np.where(x > 100, log_mass, np.log(spst.gamma.sf(x=x, a=a)))
             return log_mass
         else:
@@ -69,12 +77,14 @@ def mass_zero_t0(self, spot, texp):
         """
         fwd = self.forward(spot, texp)
         betac = 1.0 - self.beta
-        alpha = self.sigma/np.power(fwd, betac)
-        t0 = 0.5/(betac*alpha)**2
+        alpha = self.sigma / np.power(fwd, betac)
+        t0 = 0.5 / (betac * alpha) ** 2
         return t0
 
     @staticmethod
-    def price_formula(strike, spot, texp, sigma=None, cp=1, beta=0.5, intr=0.0, divr=0.0, is_fwd=False):
+    def price_formula(
+        strike, spot, texp, sigma=None, cp=1, beta=0.5, intr=0.0, divr=0.0, is_fwd=False
+    ):
         """
 
         Args:
@@ -97,15 +107,15 @@ def price_formula(strike, spot, texp, sigma=None, cp=1, beta=0.5, intr=0.0, divr
 
         betac = 1.0 - beta
         betac_inv = 1.0 / betac
-        alpha = sigma/np.power(fwd, betac)
-        sigma_std = np.maximum(alpha*np.sqrt(texp), np.finfo(float).eps)
+        alpha = sigma / np.power(fwd, betac)
+        sigma_std = np.maximum(alpha * np.sqrt(texp), np.finfo(float).eps)
         kk = strike / fwd
         x = 1.0 / np.square(betac * sigma_std)
         y = np.power(kk, 2 * betac) * x
 
         # Need to clean up the case beta > 0
         if beta > 1.0:
-            raise ValueError('Cannot handle beta value higher than 1.0')
+            raise ValueError("Cannot handle beta value higher than 1.0")
 
         ncx2_sf = spst.ncx2.sf
         ncx2_cdf = spst.ncx2.cdf
@@ -114,7 +124,7 @@ def price_formula(strike, spot, texp, sigma=None, cp=1, beta=0.5, intr=0.0, divr
         price = np.where(
             cp > 0,
             fwd * ncx2_sf(y, 2 + betac_inv, x) - strike * ncx2_cdf(x, betac_inv, y),
-            strike * ncx2_sf(x, betac_inv, y) - fwd * ncx2_cdf(y, 2 + betac_inv, x)
+            strike * ncx2_sf(x, betac_inv, y) - fwd * ncx2_cdf(y, 2 + betac_inv, x),
         )
         return disc_fac * price
 
@@ -127,11 +137,29 @@ def delta(self, strike, spot, texp, cp=1):
         y = k_star * np.power(strike, 2 / betac_inv)
 
         if self.beta < 1.0:
-            delta = 0.5 * (cp - 1) + spst.ncx2.sf(y, 2 + betac_inv, x) + 2 * x / betac_inv * \
-                    (spst.ncx2.pdf(y, 4 + betac_inv, x) - strike / fwd * spst.ncx2.pdf(x, betac_inv, y))
+            delta = (
+                0.5 * (cp - 1)
+                + spst.ncx2.sf(y, 2 + betac_inv, x)
+                + 2
+                * x
+                / betac_inv
+                * (
+                    spst.ncx2.pdf(y, 4 + betac_inv, x)
+                    - strike / fwd * spst.ncx2.pdf(x, betac_inv, y)
+                )
+            )
         else:
-            delta = 0.5 * (cp - 1) + spst.ncx2.sf(x, -betac_inv, y) - 2 * x / betac_inv * \
-                    (spst.ncx2.pdf(x, -betac_inv, y) - strike / fwd * spst.ncx2.pdf(y, 4 - betac_inv, x))
+            delta = (
+                0.5 * (cp - 1)
+                + spst.ncx2.sf(x, -betac_inv, y)
+                - 2
+                * x
+                / betac_inv
+                * (
+                    spst.ncx2.pdf(x, -betac_inv, y)
+                    - strike / fwd * spst.ncx2.pdf(y, 4 - betac_inv, x)
+                )
+            )
 
         delta *= df if self.is_fwd else divf
         return delta
@@ -141,13 +169,15 @@ def cdf(self, strike, spot, texp, cp=1):
 
         betac = 1.0 - self.beta
         betac_inv = 1.0 / betac
-        alpha = self.sigma/np.power(fwd, betac)
-        sigma_std = np.maximum(alpha*np.sqrt(texp), np.finfo(float).eps)
+        alpha = self.sigma / np.power(fwd, betac)
+        sigma_std = np.maximum(alpha * np.sqrt(texp), np.finfo(float).eps)
         kk = strike / fwd
         x = 1.0 / np.square(betac * sigma_std)
         y = np.power(kk, 2 * betac) * x
 
-        cdf = np.where(cp > 0, spst.ncx2.cdf(x, betac_inv, y), spst.ncx2.sf(x, betac_inv, y))
+        cdf = np.where(
+            cp > 0, spst.ncx2.cdf(x, betac_inv, y), spst.ncx2.sf(x, betac_inv, y)
+        )
         return cdf
 
     def gamma(self, strike, spot, texp, cp=1):
@@ -159,17 +189,29 @@ def gamma(self, strike, spot, texp, cp=1):
         y = k_star * np.power(strike, 2 / betac_inv)
 
         if self.beta < 1.0:
-            gamma = (2 + betac_inv - x) * spst.ncx2.pdf(y, 4 + betac_inv, x) + x * spst.ncx2.pdf(y, 6 + betac_inv, x) + \
-                    strike / fwd * (x * spst.ncx2.pdf(x, betac_inv, y) - y * spst.ncx2.pdf(x, 2 + betac_inv, y))
+            gamma = (
+                (2 + betac_inv - x) * spst.ncx2.pdf(y, 4 + betac_inv, x)
+                + x * spst.ncx2.pdf(y, 6 + betac_inv, x)
+                + strike
+                / fwd
+                * (
+                    x * spst.ncx2.pdf(x, betac_inv, y)
+                    - y * spst.ncx2.pdf(x, 2 + betac_inv, y)
+                )
+            )
         else:
-            gamma = (x * spst.ncx2.pdf(x, -betac_inv, y) - y * spst.ncx2.pdf(x, 2 - betac_inv, y)) + \
-                    strike / fwd * ((2 - betac_inv - x) * spst.ncx2.pdf(y, 4 - betac_inv, x) + \
-                                    x * spst.ncx2.pdf(y, 6 - betac_inv, x))
+            gamma = (
+                x * spst.ncx2.pdf(x, -betac_inv, y)
+                - y * spst.ncx2.pdf(x, 2 - betac_inv, y)
+            ) + strike / fwd * (
+                (2 - betac_inv - x) * spst.ncx2.pdf(y, 4 - betac_inv, x)
+                + x * spst.ncx2.pdf(y, 6 - betac_inv, x)
+            )
 
-        gamma *= 2 * (divf/betac_inv)**2 / df * x / fwd
+        gamma *= 2 * (divf / betac_inv) ** 2 / df * x / fwd
 
         if self.is_fwd:
-            gamma *= (df/divf)**2
+            gamma *= (df / divf) ** 2
 
         return gamma
 
@@ -184,9 +226,13 @@ def vega(self, strike, spot, texp, cp=1):
         y = k_star * np.power(strike, 2 / betac_inv)
 
         if self.beta < 1.0:
-            vega = - fwd * spst.ncx2.pdf(y, 4 + betac_inv, x) + strike * spst.ncx2.pdf(x, betac_inv, y)
+            vega = -fwd * spst.ncx2.pdf(y, 4 + betac_inv, x) + strike * spst.ncx2.pdf(
+                x, betac_inv, y
+            )
         else:
-            vega = fwd * spst.ncx2.pdf(x, -betac_inv, y) - strike * spst.ncx2.pdf(y, 4 - betac_inv, x)
+            vega = fwd * spst.ncx2.pdf(x, -betac_inv, y) - strike * spst.ncx2.pdf(
+                y, 4 - betac_inv, x
+            )
 
         sigma = self.sigma * spot ** (self.beta - 1)
 
diff --git a/pyfeng/_ex.py b/pyfeng/ex.py
similarity index 100%
rename from pyfeng/_ex.py
rename to pyfeng/ex.py
diff --git a/pyfeng/gamma.py b/pyfeng/gamma.py
index 0e397ab..197d346 100644
--- a/pyfeng/gamma.py
+++ b/pyfeng/gamma.py
@@ -19,20 +19,25 @@ class InvGam(smile.OptSmileABC):
         >>> m.price(np.arange(80, 121, 10), 100, 1.2)
         array([15.49803779,  9.53595458,  5.49889751,  3.02086661,  1.60505654])
     """
+
     sigma = None
 
     @staticmethod
-    def price_formula(strike, spot, texp, alpha, beta, cp=1, intr=0.0, divr=0.0, is_fwd=False):
+    def price_formula(
+        strike, spot, texp, alpha, beta, cp=1, intr=0.0, divr=0.0, is_fwd=False
+    ):
         disc_fac = np.exp(-texp * intr)
         fwd_scale = spot * (1.0 if is_fwd else np.exp(-texp * divr) / disc_fac) / beta
-        kk = strike/beta
+        kk = strike / beta
 
         price = np.where(
             cp > 0,
-            fwd_scale*spst.gamma.cdf(x=1/kk, a=alpha-1) - kk*spst.gamma.cdf(x=1/kk, a=alpha),
-            kk*spst.gamma.sf(x=1/kk, a=alpha) - fwd_scale*spst.gamma.sf(x=1/kk, a=alpha-1)
+            fwd_scale * spst.gamma.cdf(x=1 / kk, a=alpha - 1)
+            - kk * spst.gamma.cdf(x=1 / kk, a=alpha),
+            kk * spst.gamma.sf(x=1 / kk, a=alpha)
+            - fwd_scale * spst.gamma.sf(x=1 / kk, a=alpha - 1),
         )
-        return disc_fac*beta * price
+        return disc_fac * beta * price
 
     def alpha_beta(self, spot, texp):
         """
@@ -51,18 +56,22 @@ def alpha_beta(self, spot, texp):
         """
 
         fwd = self.forward(spot, texp)
-        alpha = 1/(np.exp(self.sigma**2*texp)-1) + 2
-        beta = (alpha-1)*fwd
+        alpha = 1 / (np.exp(self.sigma ** 2 * texp) - 1) + 2
+        beta = (alpha - 1) * fwd
         return alpha, beta
 
     def price(self, strike, spot, texp, cp=1):
         alpha, beta = self.alpha_beta(spot, texp)
-        return self.price_formula(strike, spot, texp, alpha, beta, cp=cp, intr=self.intr, divr=self.divr)
+        return self.price_formula(
+            strike, spot, texp, alpha, beta, cp=cp, intr=self.intr, divr=self.divr
+        )
 
     def cdf(self, strike, spot, texp, cp=1):
         alpha, beta = self.alpha_beta(spot, texp)
-        x = strike/beta
-        cdf = np.where(cp > 0, spst.gamma.cdf(1/x, a=alpha), spst.gamma.sf(1/x, a=alpha))
+        x = strike / beta
+        cdf = np.where(
+            cp > 0, spst.gamma.cdf(1 / x, a=alpha), spst.gamma.sf(1 / x, a=alpha)
+        )
         return cdf
 
 
@@ -82,24 +91,25 @@ class InvGauss(smile.OptSmileABC):
         >>> m.price(np.arange(80, 121, 10), 100, 1.2)
         array([15.71924064,  9.70753358,  5.54459412,  2.95300168,  1.48019682])
     """
+
     sigma = None
 
     def price(self, strike, spot, texp, cp=1):
         fwd, df, _ = self._fwd_factor(spot, texp)
-        sig2_inv = np.exp(self.sigma**2 * texp) - 1
-        ig = spst.invgauss(mu=sig2_inv, scale=1/sig2_inv)
+        sig2_inv = np.exp(self.sigma ** 2 * texp) - 1
+        ig = spst.invgauss(mu=sig2_inv, scale=1 / sig2_inv)
         kk = strike / fwd
         price = np.where(
             cp > 0,
-            ig.cdf(1/kk) - kk * ig.sf(kk),
-            kk * ig.cdf(kk) - ig.sf(1/kk),
+            ig.cdf(1 / kk) - kk * ig.sf(kk),
+            kk * ig.cdf(kk) - ig.sf(1 / kk),
         )
-        return df*fwd*price
+        return df * fwd * price
 
     def cdf(self, strike, spot, texp, cp=1):
         fwd, df, _ = self._fwd_factor(spot, texp)
-        sig2_inv = np.exp(self.sigma**2 * texp) - 1
-        ig = spst.invgauss(mu=sig2_inv, scale=1/sig2_inv)
+        sig2_inv = np.exp(self.sigma ** 2 * texp) - 1
+        ig = spst.invgauss(mu=sig2_inv, scale=1 / sig2_inv)
         x = strike / fwd
         cdf = np.where(cp > 0, ig.sf(x), ig.cdf(x))
         return cdf
diff --git a/pyfeng/garch.py b/pyfeng/garch.py
index b6c9147..0c9b597 100644
--- a/pyfeng/garch.py
+++ b/pyfeng/garch.py
@@ -41,8 +41,16 @@ def vol_paths(self, tobs, n_path=None):
         w_t[0, :] = 2 * np.log(self.sigma)
 
         for i in range(1, n_dt + 1):
-            w_t[i, :] = w_t[i - 1, :] + (self.mr * self.theta * np.exp(-w_t[i - 1, :]) - self.mr - self.vov ** 2 / 2) * \
-                        self.dt + self.vov * np.sqrt(self.dt) * rn_norm[i - 1, :]
+            w_t[i, :] = (
+                w_t[i - 1, :]
+                + (
+                    self.mr * self.theta * np.exp(-w_t[i - 1, :])
+                    - self.mr
+                    - self.vov ** 2 / 2
+                )
+                * self.dt
+                + self.vov * np.sqrt(self.dt) * rn_norm[i - 1, :]
+            )
 
         v_t = np.exp(w_t)
 
@@ -50,16 +58,16 @@ def vol_paths(self, tobs, n_path=None):
 
     def cond_fwd_vol(self, texp):
         """
-            Returns new forward and volatility conditional on volatility path (e.g., sigma_T, integrated variance)
-            The forward and volatility are standardized in the sense that F_0 = 1 and sigma_0 = 1
-            Therefore, they should be scaled by the original F_0 and sigma_0 values
+        Returns new forward and volatility conditional on volatility path (e.g., sigma_T, integrated variance)
+        The forward and volatility are standardized in the sense that F_0 = 1 and sigma_0 = 1
+        Therefore, they should be scaled by the original F_0 and sigma_0 values
 
-            Args:
-                theta: the long term average
-                mr: coefficient of dt
-                texp: time-to-expiry
+        Args:
+            theta: the long term average
+            mr: coefficient of dt
+            texp: time-to-expiry
 
-            Returns: (forward, volatility)
+        Returns: (forward, volatility)
         """
         rhoc = np.sqrt(1.0 - self.rho ** 2)
         tobs = self.tobs(texp)
@@ -67,12 +75,19 @@ def cond_fwd_vol(self, texp):
         v_paths = self.vol_paths(tobs)
         sigma_paths = np.sqrt(v_paths)
         sigma_final = sigma_paths[-1, :]
-        int_sigma = scint.simps(sigma_paths, dx=texp/n_dt, axis=0)
-        int_var = scint.simps(sigma_paths ** 2, dx=texp/n_dt, axis=0)
-        int_sigma_inv = scint.simps(1 / sigma_paths, dx=texp/n_dt, axis=0)
-
-        fwd_cond = np.exp(self.rho * (2 * (sigma_final - self.sigma) / self.vov - self.mr * self.theta * int_sigma_inv /
-                                      self.vov + (self.mr / self.vov + self.vov / 4) * int_sigma - self.rho * int_var / 2))  # scaled by initial value
+        int_sigma = scint.simps(sigma_paths, dx=texp / n_dt, axis=0)
+        int_var = scint.simps(sigma_paths ** 2, dx=texp / n_dt, axis=0)
+        int_sigma_inv = scint.simps(1 / sigma_paths, dx=texp / n_dt, axis=0)
+
+        fwd_cond = np.exp(
+            self.rho
+            * (
+                2 * (sigma_final - self.sigma) / self.vov
+                - self.mr * self.theta * int_sigma_inv / self.vov
+                + (self.mr / self.vov + self.vov / 4) * int_sigma
+                - self.rho * int_var / 2
+            )
+        )  # scaled by initial value
 
         vol_cond = rhoc * np.sqrt(int_var / texp)
 
diff --git a/pyfeng/heston_mc.py b/pyfeng/heston_mc.py
index a88d756..ef616bb 100644
--- a/pyfeng/heston_mc.py
+++ b/pyfeng/heston_mc.py
@@ -18,37 +18,48 @@ def vol_paths(self, tobs):
         n_dt = len(dt)
 
         dB_t = self._bm_incr(tobs, cum=False)  # B_t (0 <= s <= 1)
-        vv_path = np.empty([n_dt+1, self.n_path])  # variance series: V0, V1,...,VT
-        vv = self.sigma**2
+        vv_path = np.empty([n_dt + 1, self.n_path])  # variance series: V0, V1,...,VT
+        vv = self.sigma ** 2
         vv_path[0, :] = vv
         for i in range(n_dt):
-            vv = vv + self.mr * (self.theta ** 2 - vv) * dt[i] + np.sqrt(vv) * self.vov * dB_t[i, :]  # Euler method
-            vv = vv + 0.25 * self.vov**2 * (dB_t[i, :]**2 - dt[i])  # Milstein method
+            vv = (
+                vv
+                + self.mr * (self.theta ** 2 - vv) * dt[i]
+                + np.sqrt(vv) * self.vov * dB_t[i, :]
+            )  # Euler method
+            vv = vv + 0.25 * self.vov ** 2 * (
+                dB_t[i, :] ** 2 - dt[i]
+            )  # Milstein method
             vv[vv < 0] = 0  # variance should be larger than zero
-            vv_path[i+1, :] = vv
+            vv_path[i + 1, :] = vv
 
         # return normalized sigma, e.g., sigma(0) = 1
-        return np.sqrt(vv_path)/self.sigma
+        return np.sqrt(vv_path) / self.sigma
 
     def cond_fwd_vol(self, texp):
 
         tobs = self.tobs(texp)
         n_steps = len(tobs)
         sigma_paths = self.vol_paths(tobs)
-        vv0 = self.sigma**2
+        vv0 = self.sigma ** 2
         vv_ratio = sigma_paths[-1, :]
-        int_var_std = scint.simps(sigma_paths**2, dx=1, axis=0) / n_steps
+        int_var_std = scint.simps(sigma_paths ** 2, dx=1, axis=0) / n_steps
 
-        int_sig_dw = ((vv_ratio - 1) * vv0 - self.mr * texp * (self.theta ** 2 - int_var_std * vv0)) / self.vov
-        fwd_cond = np.exp(self.rho * int_sig_dw - 0.5*self.rho**2 * int_var_std * vv0 * texp)
-        vol_cond = np.sqrt((1 - self.rho**2) * int_var_std)
+        int_sig_dw = (
+            (vv_ratio - 1) * vv0
+            - self.mr * texp * (self.theta ** 2 - int_var_std * vv0)
+        ) / self.vov
+        fwd_cond = np.exp(
+            self.rho * int_sig_dw - 0.5 * self.rho ** 2 * int_var_std * vv0 * texp
+        )
+        vol_cond = np.sqrt((1 - self.rho ** 2) * int_var_std)
 
         # return normalized forward and volatility
         return fwd_cond, vol_cond
 
 
 class HestonCondMcQE:
-    '''
+    """
     Conditional MC for Heston model based on QE discretization scheme by Andersen(2008)
 
     Underlying price is assumed to follow a geometric Brownian motion.
@@ -72,10 +83,10 @@ class HestonCondMcQE:
                [13.12057087,  0.29501411, 35.90207168],
                [13.1042753 ,  0.29476289, 35.89245755],
                [13.09047939,  0.29547721, 35.86410028]])
-    '''
+    """
 
     def __init__(self, vov=1, kappa=0.5, rho=-0.9, theta=0.04):
-        '''
+        """
         Initiate a Heston model
 
         Args:
@@ -83,7 +94,7 @@ def __init__(self, vov=1, kappa=0.5, rho=-0.9, theta=0.04):
             kappa: speed of variance's mean-reversion, strictly positive
             rho: correlation between BMs of price and vol
             theta: long-term mean (equilibirum level) of the variance, strictly positive
-        '''
+        """
         self.vov = vov
         self.kappa = kappa
         self.rho = rho
@@ -93,8 +104,21 @@ def __init__(self, vov=1, kappa=0.5, rho=-0.9, theta=0.04):
         self.rx_results = None
         self.dis = 1e-3
 
-    def price(self, strike, spot, texp, sigma, delta, intr=0, divr=0, psi_c=1.5, path=10000, scheme='QE', seed=None):
-        '''
+    def price(
+        self,
+        strike,
+        spot,
+        texp,
+        sigma,
+        delta,
+        intr=0,
+        divr=0,
+        psi_c=1.5,
+        path=10000,
+        scheme="QE",
+        seed=None,
+    ):
+        """
         Conditional MC routine for Heston model
         Generate paths for vol only using QE discretization scheme.
         Compute integrated variance and get BSM prices vector for all strikes.
@@ -114,7 +138,7 @@ def price(self, strike, spot, texp, sigma, delta, intr=0, divr=0, psi_c=1.5, pat
 
         Return:
             BSM price vector for all strikes
-        '''
+        """
         self.sigma = sigma
         self.bsm_model = bsm.Bsm(self.sigma, intr=intr, divr=divr)
         self.delta = delta
@@ -124,15 +148,18 @@ def price(self, strike, spot, texp, sigma, delta, intr=0, divr=0, psi_c=1.5, pat
         vt = self.sigma ** 2 * np.ones([self.path, self.step + 1])
         np.random.seed(seed)
 
-        if scheme == 'QE':
+        if scheme == "QE":
             u = np.random.uniform(size=(self.path, self.step))
 
             expo = np.exp(-self.kappa * self.delta)
             for i in range(self.step):
                 # compute m, s_square, psi given vt(i)
                 m = self.theta + (vt[:, i] - self.theta) * expo
-                s2 = vt[:, i] * (self.vov ** 2) * expo * (1 - expo) / self.kappa + self.theta * (self.vov ** 2) * \
-                     ((1 - expo) ** 2) / (2 * self.kappa)
+                s2 = vt[:, i] * (self.vov ** 2) * expo * (
+                    1 - expo
+                ) / self.kappa + self.theta * (self.vov ** 2) * ((1 - expo) ** 2) / (
+                    2 * self.kappa
+                )
                 psi = s2 / m ** 2
 
                 # compute vt(i+1) given psi
@@ -142,18 +169,20 @@ def price(self, strike, spot, texp, sigma, delta, intr=0, divr=0, psi_c=1.5, pat
                 b = np.sqrt(b2)
                 a = m[below] / (1 + b2)
                 z = spst.norm.ppf(u[below, i])
-                vt[below, i+1] = a * (b + z) ** 2
+                vt[below, i + 1] = a * (b + z) ** 2
 
                 above = np.where(psi > psi_c)[0]
                 p = (psi[above] - 1) / (psi[above] + 1)
                 beta = (1 - p) / m[above]
                 for k in range(len(above)):
                     if u[above[k], i] > p[k]:
-                        vt[above[k], i+1] = beta[k] ** -1 * np.log((1 - p[k]) / (1 - u[above[k], i]))
+                        vt[above[k], i + 1] = beta[k] ** -1 * np.log(
+                            (1 - p[k]) / (1 - u[above[k], i])
+                        )
                     else:
-                        vt[above[k], i+1] = 0
+                        vt[above[k], i + 1] = 0
 
-        elif scheme == 'TG':
+        elif scheme == "TG":
             if np.all(self.rx_results) == None:
                 self.psi_points, self.rx_results = self.prepare_rx()
 
@@ -161,8 +190,11 @@ def price(self, strike, spot, texp, sigma, delta, intr=0, divr=0, psi_c=1.5, pat
             for i in range(self.step):
                 # compute m, s_square, psi given vt(i)
                 m = self.theta + (vt[:, i] - self.theta) * expo
-                s2 = vt[:, i] * (self.vov ** 2) * expo * (1 - expo) / self.kappa + self.theta * (self.vov ** 2) * \
-                     ((1 - expo) ** 2) / (2 * self.kappa)
+                s2 = vt[:, i] * (self.vov ** 2) * expo * (
+                    1 - expo
+                ) / self.kappa + self.theta * (self.vov ** 2) * ((1 - expo) ** 2) / (
+                    2 * self.kappa
+                )
                 psi = s2 / m ** 2
 
                 rx = np.array([self.find_rx(j) for j in psi])
@@ -171,65 +203,91 @@ def price(self, strike, spot, texp, sigma, delta, intr=0, divr=0, psi_c=1.5, pat
                 mu_v = np.zeros_like(z)
                 sigma_v = np.zeros_like(z)
                 mu_v[:, i] = rx * m / (spst.norm.pdf(rx) + rx * spst.norm.cdf(rx))
-                sigma_v[:, i] = np.sqrt(s2) * psi ** (-0.5) / (spst.norm.pdf(rx) + rx * spst.norm.cdf(rx))
+                sigma_v[:, i] = (
+                    np.sqrt(s2)
+                    * psi ** (-0.5)
+                    / (spst.norm.pdf(rx) + rx * spst.norm.cdf(rx))
+                )
 
-                vt[:, i+1] = np.fmax(mu_v[:, i] + sigma_v[:, i] * z[:, i], 0)
+                vt[:, i + 1] = np.fmax(mu_v[:, i] + sigma_v[:, i] * z[:, i], 0)
 
-        elif scheme == 'Euler':
+        elif scheme == "Euler":
             z = np.random.normal(size=(self.path, self.step))
             for i in range(self.step):
-                vt[:, i+1] = vt[:, i] + self.kappa * (self.theta - np.max(vt[:, i], 0)) * self.delta + \
-                             self.vov * np.sqrt(np.max(vt[:, i], 0) * self.delta) * z[:, i]
+                vt[:, i + 1] = (
+                    vt[:, i]
+                    + self.kappa * (self.theta - np.max(vt[:, i], 0)) * self.delta
+                    + self.vov * np.sqrt(np.max(vt[:, i], 0) * self.delta) * z[:, i]
+                )
             below_0 = np.where(vt < 0)
             vt[below_0] = 0
 
-        elif scheme == 'Milstein':
+        elif scheme == "Milstein":
             z = np.random.normal(size=(self.path, self.step))
             for i in range(self.step):
-                vt[:, i+1] = vt[:, i] + self.kappa * (self.theta - np.max(vt[:, i], 0)) * self.delta + self.vov * \
-                             np.sqrt(np.max(vt[:, i], 0) * self.delta) * z[:, i] + \
-                             self.vov**2 * 0.25 * (z[:, i]**2 - 1) * self.delta
+                vt[:, i + 1] = (
+                    vt[:, i]
+                    + self.kappa * (self.theta - np.max(vt[:, i], 0)) * self.delta
+                    + self.vov * np.sqrt(np.max(vt[:, i], 0) * self.delta) * z[:, i]
+                    + self.vov ** 2 * 0.25 * (z[:, i] ** 2 - 1) * self.delta
+                )
             below_0 = np.where(vt < 0)
             vt[below_0] = 0
 
-        elif scheme == 'KJ':
+        elif scheme == "KJ":
             z = np.random.normal(size=(self.path, self.step))
             for i in range(self.step):
-                vt[:, i+1] = (vt[:, i] + self.kappa * self.theta * self.delta + self.vov * \
-                             np.sqrt(np.max(vt[:, i], 0) * self.delta) * z[:, i] + \
-                             self.vov**2 * 0.25 * (z[:, i]**2 - 1) * self.delta) / (1 + self.kappa * self.delta)
+                vt[:, i + 1] = (
+                    vt[:, i]
+                    + self.kappa * self.theta * self.delta
+                    + self.vov * np.sqrt(np.max(vt[:, i], 0) * self.delta) * z[:, i]
+                    + self.vov ** 2 * 0.25 * (z[:, i] ** 2 - 1) * self.delta
+                ) / (1 + self.kappa * self.delta)
             below_0 = np.where(vt < 0)
             vt[below_0] = 0
 
         # compute integral of vt, equivalent spot and vol
         vt_int = spint.simps(vt, dx=self.delta)
-        spot_cmc = spot * np.exp(self.rho * (vt[:, -1] - vt[:, 0] - self.kappa * (self.theta * texp - vt_int))
-                                 / self.vov - self.rho ** 2 * vt_int / 2)
+        spot_cmc = spot * np.exp(
+            self.rho
+            * (vt[:, -1] - vt[:, 0] - self.kappa * (self.theta * texp - vt_int))
+            / self.vov
+            - self.rho ** 2 * vt_int / 2
+        )
         vol_cmc = np.sqrt((1 - self.rho ** 2) * vt_int / texp)
 
         # compute bsm price vector for the given strike vector
         price_cmc = np.zeros_like(strike)
         for j in range(len(strike)):
-            price_cmc[j] = np.mean(self.bsm_model.price_formula(strike[j], spot_cmc, vol_cmc, texp, intr=intr, divr=divr))
+            price_cmc[j] = np.mean(
+                self.bsm_model.price_formula(
+                    strike[j], spot_cmc, vol_cmc, texp, intr=intr, divr=divr
+                )
+            )
 
         return price_cmc
 
     def prepare_rx(self):
-        '''
+        """
         Pre-calculate r(x) and store the result
         for TG scheme only
-        '''
-        fx = lambda rx: rx * spst.norm.pdf(rx) + spst.norm.cdf(rx) * (1 + rx ** 2) / \
-                        ((spst.norm.pdf(rx) + rx * spst.norm.cdf(rx)) ** 2) - 1
+        """
+        fx = (
+            lambda rx: rx * spst.norm.pdf(rx)
+            + spst.norm.cdf(rx)
+            * (1 + rx ** 2)
+            / ((spst.norm.pdf(rx) + rx * spst.norm.cdf(rx)) ** 2)
+            - 1
+        )
         rx_results = np.linspace(-2, 100, 10 ** 5)
         psi_points = fx(rx_results)
 
         return psi_points, rx_results
 
     def find_rx(self, psi):
-        '''
+        """
         Return r(psi) according to the pre_calculated results
-        '''
+        """
 
         if self.rx_results[self.psi_points >= psi].size == 0:
             print("Caution: input psi too large")
@@ -238,7 +296,10 @@ def find_rx(self, psi):
             print("Caution: input psi too small")
             return self.rx_results[0]
         else:
-            return (self.rx_results[self.psi_points >= psi][0] + self.rx_results[self.psi_points <= psi][-1])/2
+            return (
+                self.rx_results[self.psi_points >= psi][0]
+                + self.rx_results[self.psi_points <= psi][-1]
+            ) / 2
 
 
 class HestonMcAe:
@@ -256,6 +317,7 @@ class HestonMcAe:
         >>> price_ae
         8.946951375550809
     """
+
     def __init__(self, vov=1, kappa=0.5, rho=-0.9, theta=0.04, r=0):
         """
         Initiate a Heston model
@@ -286,22 +348,43 @@ def ch_f(self, texp, sigma_0, sigma_t, chi_dim):
             ch_f: characteristic function of the distribution of integral sigma_t
 
         """
-        gamma_f = lambda a: np.sqrt(self.kappa ** 2 - 2 * self.vov**2 * a * 1j)
+        gamma_f = lambda a: np.sqrt(self.kappa ** 2 - 2 * self.vov ** 2 * a * 1j)
 
         temp_f = lambda a: gamma_f(a) * texp
 
-        ch_f_part_1 = lambda a: gamma_f(a) * np.exp(-0.5 * (temp_f(a) - self.kappa * texp)) \
-                                * (1 - np.exp(-self.kappa * texp)) / (self.kappa * (1 - np.exp(-temp_f(a))))
-
-        ch_f_part_2 = lambda a: np.exp((sigma_0**2 + sigma_t**2) / self.vov ** 2 * \
-                                       (self.kappa * (1 + np.exp(-self.kappa * texp)) / (1 - np.exp(-self.kappa * texp))
-                                        - gamma_f(a) * (1 + np.exp(-temp_f(a))) / (1 - np.exp(-temp_f(a)))))
-
-        ch_f_part_3 = lambda a: spsp.iv(0.5 * chi_dim - 1, np.sqrt(sigma_0**2 * sigma_t**2) * 4 * gamma_f(a) *
-                                   np.exp(-0.5 * temp_f(a)) / (self.vov ** 2 * (1 - np.exp(-temp_f(a))))) / \
-                                spsp.iv(0.5 * chi_dim - 1, np.sqrt(sigma_0**2 * sigma_t**2) * 4 * self.kappa *
-                                   np.exp(-0.5 * self.kappa * texp) / (
-                                           self.vov ** 2 * (1 - np.exp(- self.kappa * texp))))
+        ch_f_part_1 = (
+            lambda a: gamma_f(a)
+            * np.exp(-0.5 * (temp_f(a) - self.kappa * texp))
+            * (1 - np.exp(-self.kappa * texp))
+            / (self.kappa * (1 - np.exp(-temp_f(a))))
+        )
+
+        ch_f_part_2 = lambda a: np.exp(
+            (sigma_0 ** 2 + sigma_t ** 2)
+            / self.vov ** 2
+            * (
+                self.kappa
+                * (1 + np.exp(-self.kappa * texp))
+                / (1 - np.exp(-self.kappa * texp))
+                - gamma_f(a) * (1 + np.exp(-temp_f(a))) / (1 - np.exp(-temp_f(a)))
+            )
+        )
+
+        ch_f_part_3 = lambda a: spsp.iv(
+            0.5 * chi_dim - 1,
+            np.sqrt(sigma_0 ** 2 * sigma_t ** 2)
+            * 4
+            * gamma_f(a)
+            * np.exp(-0.5 * temp_f(a))
+            / (self.vov ** 2 * (1 - np.exp(-temp_f(a)))),
+        ) / spsp.iv(
+            0.5 * chi_dim - 1,
+            np.sqrt(sigma_0 ** 2 * sigma_t ** 2)
+            * 4
+            * self.kappa
+            * np.exp(-0.5 * self.kappa * texp)
+            / (self.vov ** 2 * (1 - np.exp(-self.kappa * texp))),
+        )
 
         ch_f = lambda a: ch_f_part_1(a) * ch_f_part_2(a) * ch_f_part_3(a)
         return ch_f
@@ -320,7 +403,9 @@ def gen_vov_t(self, chi_dim, chi_lambda, texp, n_paths):
 
         """
         cof = self.vov ** 2 * (1 - np.exp(-self.kappa * texp)) / (4 * self.kappa)
-        sigma_t = np.sqrt(cof * np.random.noncentral_chisquare(chi_dim, chi_lambda, n_paths))
+        sigma_t = np.sqrt(
+            cof * np.random.noncentral_chisquare(chi_dim, chi_lambda, n_paths)
+        )
         return sigma_t
 
     def gen_s_t(self, spot, sigma_t, sigma_0, texp, integral_sigma_t, n_paths):
@@ -338,15 +423,32 @@ def gen_s_t(self, spot, sigma_t, sigma_0, texp, integral_sigma_t, n_paths):
             s_t: stock price at time T
         """
 
-        integral_sqrt_sigma_t = (sigma_t**2 - sigma_0**2 - self.kappa * self.theta * texp + self.kappa * integral_sigma_t)\
-                                / self.vov
-        mean = np.log(spot) + (self.r * texp - 0.5 * integral_sigma_t + self.rho * integral_sqrt_sigma_t)
+        integral_sqrt_sigma_t = (
+            sigma_t ** 2
+            - sigma_0 ** 2
+            - self.kappa * self.theta * texp
+            + self.kappa * integral_sigma_t
+        ) / self.vov
+        mean = np.log(spot) + (
+            self.r * texp - 0.5 * integral_sigma_t + self.rho * integral_sqrt_sigma_t
+        )
         sigma_2 = np.sqrt((1 - self.rho ** 2) * integral_sigma_t)
         s_t = np.exp(mean + sigma_2 * np.random.normal(size=n_paths))
         return s_t
 
-    def price(self, strike, spot, texp, sigma_0, intr=0, divr=0, n_paths=10000, seed=None,
-              dis_can="Inverse-Gaussian", call=1):
+    def price(
+        self,
+        strike,
+        spot,
+        texp,
+        sigma_0,
+        intr=0,
+        divr=0,
+        n_paths=10000,
+        seed=None,
+        dis_can="Inverse-Gaussian",
+        call=1,
+    ):
         """
         Args:
             strike: strike price
@@ -364,8 +466,11 @@ def price(self, strike, spot, texp, sigma_0, intr=0, divr=0, n_paths=10000, seed
             np.random.seed(seed)
 
         chi_dim = (4 * self.theta * self.kappa) / (self.vov ** 2)
-        chi_lambda = (4 * self.kappa * np.exp(-self.kappa * texp)) / \
-                     ((self.vov ** 2) * (1 - np.exp(-self.kappa * texp))) * sigma_0**2
+        chi_lambda = (
+            (4 * self.kappa * np.exp(-self.kappa * texp))
+            / ((self.vov ** 2) * (1 - np.exp(-self.kappa * texp)))
+            * sigma_0 ** 2
+        )
 
         sigma_t = self.gen_vov_t(chi_dim, chi_lambda, texp, n_paths)
 
@@ -375,7 +480,7 @@ def price(self, strike, spot, texp, sigma_0, intr=0, divr=0, n_paths=10000, seed
         moment_2st = (derivative(ch_f, 0, n=2, dx=1e-5) / (1j ** 2)).real
 
         if dis_can == "Inverse-Gaussian":
-            scale_ig = moment_1st**3 / (moment_2st - moment_1st**2)
+            scale_ig = moment_1st ** 3 / (moment_2st - moment_1st ** 2)
             miu_ig = moment_1st / scale_ig
             integral_sigma_t = spst.invgauss.rvs(miu_ig, scale=scale_ig)
             s_t = self.gen_s_t(spot, sigma_t, sigma_0, texp, integral_sigma_t, n_paths)
@@ -394,4 +499,4 @@ def price(self, strike, spot, texp, sigma_0, intr=0, divr=0, n_paths=10000, seed
         else:
             price_ae = np.fmax(strike - s_t, 0).mean()
 
-        return np.exp(- self.r * texp) * price_ae
+        return np.exp(-self.r * texp) * price_ae
diff --git a/pyfeng/multiasset.py b/pyfeng/multiasset.py
index 5407e99..b80117f 100644
--- a/pyfeng/multiasset.py
+++ b/pyfeng/multiasset.py
@@ -29,7 +29,9 @@ class BsmSpreadKirk(opt.OptMaABC):
 
     def price(self, strike, spot, texp, cp=1):
         df = np.exp(-texp * self.intr)
-        fwd = np.array(spot) * (1.0 if self.is_fwd else np.exp(-texp * np.array(self.divr)) / df)
+        fwd = np.array(spot) * (
+            1.0 if self.is_fwd else np.exp(-texp * np.array(self.divr)) / df
+        )
         assert fwd.shape[-1] == self.n_asset
 
         fwd1 = fwd[..., 0] - np.minimum(strike, 0)
@@ -37,7 +39,7 @@ def price(self, strike, spot, texp, cp=1):
 
         sig1 = self.sigma[0] * fwd[..., 0] / fwd1
         sig2 = self.sigma[1] * fwd[..., 1] / fwd2
-        sig_spd = np.sqrt(sig1*(sig1 - 2.0*self.rho*sig2) + sig2**2)
+        sig_spd = np.sqrt(sig1 * (sig1 - 2.0 * self.rho * sig2) + sig2 ** 2)
         price = bsm.Bsm.price_formula(fwd2, fwd1, sig_spd, texp, cp=cp, is_fwd=True)
         return df * price
 
@@ -62,27 +64,32 @@ class BsmSpreadBjerksund2014(opt.OptMaABC):
 
     def price(self, strike, spot, texp, cp=1):
         df = np.exp(-texp * self.intr)
-        fwd = np.array(spot) * (1.0 if self.is_fwd else np.exp(-texp * np.array(self.divr)) / df)
+        fwd = np.array(spot) * (
+            1.0 if self.is_fwd else np.exp(-texp * np.array(self.divr)) / df
+        )
         assert fwd.shape[-1] == self.n_asset
 
         fwd1 = fwd[..., 0]
         fwd2 = fwd[..., 1]
 
-        std11 = self.sigma[0]**2 * texp
-        std12 = self.sigma[0]*self.sigma[1] * texp
-        std22 = self.sigma[1]**2 * texp
+        std11 = self.sigma[0] ** 2 * texp
+        std12 = self.sigma[0] * self.sigma[1] * texp
+        std22 = self.sigma[1] ** 2 * texp
 
         aa = fwd2 + strike
-        bb = fwd2/aa
-        std = np.sqrt(std11 - 2*bb*self.rho*std12 + bb**2*std22)
+        bb = fwd2 / aa
+        std = np.sqrt(std11 - 2 * bb * self.rho * std12 + bb ** 2 * std22)
 
-        d3 = np.log(fwd1/aa)
-        d1 = (d3 + 0.5*std11 - bb*(self.rho*std12 - 0.5*bb*std22)) / std
-        d2 = (d3 - 0.5*std11 + self.rho*std12 + bb*(0.5*bb - 1)*std22) / std
-        d3 = (d3 - 0.5*std11 + 0.5*bb**2*std22) / std
+        d3 = np.log(fwd1 / aa)
+        d1 = (d3 + 0.5 * std11 - bb * (self.rho * std12 - 0.5 * bb * std22)) / std
+        d2 = (d3 - 0.5 * std11 + self.rho * std12 + bb * (0.5 * bb - 1) * std22) / std
+        d3 = (d3 - 0.5 * std11 + 0.5 * bb ** 2 * std22) / std
 
-        price = cp*(fwd1*spst.norm.cdf(cp*d1) - fwd2*spst.norm.cdf(cp*d2)
-                    - strike*spst.norm.cdf(cp*d3))
+        price = cp * (
+            fwd1 * spst.norm.cdf(cp * d1)
+            - fwd2 * spst.norm.cdf(cp * d2)
+            - strike * spst.norm.cdf(cp * d3)
+        )
 
         return df * price
 
@@ -119,14 +126,17 @@ def __init__(self, sigma, cor=None, weight=None, intr=0.0, divr=0.0, is_fwd=Fals
 
     def price(self, strike, spot, texp, cp=1):
         df = np.exp(-texp * self.intr)
-        fwd = np.array(spot) * (1.0 if self.is_fwd else np.exp(-texp * np.array(self.divr)) / df)
+        fwd = np.array(spot) * (
+            1.0 if self.is_fwd else np.exp(-texp * np.array(self.divr)) / df
+        )
         assert fwd.shape[-1] == self.n_asset
 
         fwd_basket = fwd @ self.weight
         vol_basket = np.sqrt(self.weight @ self.cov_m @ self.weight)
 
         price = norm.Norm.price_formula(
-            strike, fwd_basket, vol_basket, texp, cp=cp, is_fwd=True)
+            strike, fwd_basket, vol_basket, texp, cp=cp, is_fwd=True
+        )
         return df * price
 
 
@@ -142,6 +152,7 @@ class NormSpread(opt.OptMaABC):
         >>> m.price(np.arange(-2, 3) * 10, [100, 120], 1.3)
         array([17.95676186, 13.74646821, 10.26669936,  7.47098719,  5.29057157])
     """
+
     weight = np.array([1, -1])
 
     price = NormBasket.price
@@ -166,18 +177,20 @@ class BsmBasketLevy1992(NormBasket):
                28.05196621, 24.70229571, 21.800801  , 19.28360474, 17.09570196,
                15.19005654])
     """
+
     def price(self, strike, spot, texp, cp=1):
         df = np.exp(-texp * self.intr)
-        fwd = np.array(spot) * (1.0 if self.is_fwd else np.exp(-texp * np.array(self.divr)) / df)
+        fwd = np.array(spot) * (
+            1.0 if self.is_fwd else np.exp(-texp * np.array(self.divr)) / df
+        )
         assert fwd.shape[-1] == self.n_asset
 
         fwd_basket = fwd * self.weight
         m1 = np.sum(fwd_basket, axis=-1)
         m2 = np.sum(fwd_basket @ np.exp(self.cov_m * texp) * fwd_basket, axis=-1)
 
-        sig = np.sqrt(np.log(m2/(m1**2))/texp)
-        price = bsm.Bsm.price_formula(
-            strike, m1, sig, texp, cp=cp, is_fwd=True)
+        sig = np.sqrt(np.log(m2 / (m1 ** 2)) / texp)
+        price = bsm.Bsm.price_formula(strike, m1, sig, texp, cp=cp, is_fwd=True)
         return df * price
 
 
@@ -202,19 +215,24 @@ class BsmBasketMilevsky1998(NormBasket):
                24.49577509, 21.38543199, 18.77356434, 16.56909804, 14.69831445,
                13.10186928])
     """
+
     def price(self, strike, spot, texp, cp=1):
         df = np.exp(-texp * self.intr)
-        fwd = np.array(spot) * (1.0 if self.is_fwd else np.exp(-texp * np.array(self.divr)) / df)
+        fwd = np.array(spot) * (
+            1.0 if self.is_fwd else np.exp(-texp * np.array(self.divr)) / df
+        )
         assert fwd.shape[-1] == self.n_asset
 
         fwd_basket = fwd * self.weight
         m1 = np.sum(fwd_basket, axis=-1)
         m2 = np.sum(fwd_basket @ np.exp(self.cov_m * texp) * fwd_basket, axis=-1)
 
-        alpha = 1/(m2/m1**2-1) + 2
-        beta = (alpha-1)*m1
+        alpha = 1 / (m2 / m1 ** 2 - 1) + 2
+        beta = (alpha - 1) * m1
 
-        price = gamma.InvGam.price_formula(strike, m1, texp, alpha, beta, cp=cp, is_fwd=True)
+        price = gamma.InvGam.price_formula(
+            strike, m1, texp, alpha, beta, cp=cp, is_fwd=True
+        )
         return df * price
 
 
@@ -243,25 +261,29 @@ def __init__(self, sigma, cor=None, weight=None, intr=0.0, divr=0.0, is_fwd=Fals
     def price(self, strike, spot, texp, cp=1):
         sig = self.sigma
         df = np.exp(-texp * self.intr)
-        fwd = np.array(spot) * (1.0 if self.is_fwd else np.exp(-texp * np.array(self.divr)) / df)
+        fwd = np.array(spot) * (
+            1.0 if self.is_fwd else np.exp(-texp * np.array(self.divr)) / df
+        )
         assert fwd.shape[-1] == self.n_asset
 
         sig_std = sig * np.sqrt(texp)
-        spd_rho = np.sqrt(np.dot(sig, sig) - 2*self.rho*sig[0]*sig[1])
+        spd_rho = np.sqrt(np.dot(sig, sig) - 2 * self.rho * sig[0] * sig[1])
         spd_std = spd_rho * np.sqrt(texp)
 
         # -x and y as rows
         # supposed to be -log(fwd/strike) but strike is added later
-        xx = -np.log(fwd)/sig_std - 0.5*sig_std
+        xx = -np.log(fwd) / sig_std - 0.5 * sig_std
 
-        fwd_ratio = fwd[0]/fwd[1]
-        yy = np.log([fwd_ratio, 1/fwd_ratio])/spd_std + 0.5*spd_std
+        fwd_ratio = fwd[0] / fwd[1]
+        yy = np.log([fwd_ratio, 1 / fwd_ratio]) / spd_std + 0.5 * spd_std
 
-        rho12 = np.array([self.rho*sig[1]-sig[0], self.rho*sig[0]-sig[1]])/spd_rho
+        rho12 = (
+            np.array([self.rho * sig[1] - sig[0], self.rho * sig[0] - sig[1]]) / spd_rho
+        )
 
         mu0 = np.zeros(2)
-        cor_m1 = rho12[0] + (1-rho12[0])*np.eye(2)
-        cor_m2 = rho12[1] + (1-rho12[1])*np.eye(2)
+        cor_m1 = rho12[0] + (1 - rho12[0]) * np.eye(2)
+        cor_m2 = rho12[1] + (1 - rho12[1]) * np.eye(2)
 
         strike_isscalar = np.isscalar(strike)
         strike = np.atleast_1d(strike)
@@ -273,13 +295,21 @@ def price(self, strike, spot, texp, cp=1):
         parity = 0 if np.all(cp > 0) else self.m_switch.price(0, fwd, texp) + fwd[1]
         price = np.zeros_like(strike, float)
         for k in range(n_strike):
-            xx_ = xx + np.log(strike[k])/sig_std
-            term1 = fwd[0] * (spst.norm.cdf(yy[0]) - spst.multivariate_normal.cdf(np.array([xx_[0], yy[0]]), mu0, cor_m1))
-            term2 = fwd[1] * (spst.norm.cdf(yy[1]) - spst.multivariate_normal.cdf(np.array([xx_[1], yy[1]]), mu0, cor_m2))
-            term3 = strike[k] * np.array(spst.multivariate_normal.cdf(xx_ + sig_std, mu0, self.cor_m))
-
-            assert(term1 + term2 + term3 >= strike[k])
-            price[k] = (term1 + term2 + term3 - strike[k])
+            xx_ = xx + np.log(strike[k]) / sig_std
+            term1 = fwd[0] * (
+                spst.norm.cdf(yy[0])
+                - spst.multivariate_normal.cdf(np.array([xx_[0], yy[0]]), mu0, cor_m1)
+            )
+            term2 = fwd[1] * (
+                spst.norm.cdf(yy[1])
+                - spst.multivariate_normal.cdf(np.array([xx_[1], yy[1]]), mu0, cor_m2)
+            )
+            term3 = strike[k] * np.array(
+                spst.multivariate_normal.cdf(xx_ + sig_std, mu0, self.cor_m)
+            )
+
+            assert term1 + term2 + term3 >= strike[k]
+            price[k] = term1 + term2 + term3 - strike[k]
 
             if cp[k] < 0:
                 price[k] += strike[k] - parity
@@ -334,7 +364,7 @@ def root(fac, std, strike):
         scalar_output = np.isscalar(np.sum(fac * std, axis=-1) - strike)
         strike = np.atleast_1d(strike)
 
-        with np.errstate(divide='ignore', invalid='ignore'):
+        with np.errstate(divide="ignore", invalid="ignore"):
             log_k = np.where(strike > 0, np.log(strike), 1)
 
             # Initial guess with linearlized assmption
@@ -358,8 +388,16 @@ def root(fac, std, strike):
 
         for k in range(32):
             y_vec = fac * np.exp(std * x[ind, None])
-            y = np.log(np.sum(y_vec, axis=-1)) - log_k[ind] if log else np.sum(y_vec, axis=-1) - strike[ind]
-            dy = np.sum(std * y_vec, axis=-1) / np.sum(y_vec, axis=-1) if log else np.sum(std * y_vec, axis=-1)
+            y = (
+                np.log(np.sum(y_vec, axis=-1)) - log_k[ind]
+                if log
+                else np.sum(y_vec, axis=-1) - strike[ind]
+            )
+            dy = (
+                np.sum(std * y_vec, axis=-1) / np.sum(y_vec, axis=-1)
+                if log
+                else np.sum(std * y_vec, axis=-1)
+            )
             x[ind] -= y / dy
             if len(y) == 0:
                 print(ind, y_vec, y)
@@ -368,7 +406,9 @@ def root(fac, std, strike):
                 break
 
         if y_err_max > BsmBasket1Bm.IMPVOL_TOL:
-            warn_msg = f'root did not converge within {k} iterations: max error = {y_err_max}'
+            warn_msg = (
+                f"root did not converge within {k} iterations: max error = {y_err_max}"
+            )
             warnings.warn(warn_msg, Warning)
 
         return x[0] if scalar_output else x
@@ -380,21 +420,22 @@ def price(self, strike, spot, texp, cp=1):
         fwd_basket = fwd * self.weight
         sigma_std = self.sigma * np.sqrt(texp)
         cp = np.array(cp)
-        d2 = -cp * self.root(fwd_basket * np.exp(-sigma_std**2/2), sigma_std, strike)
+        d2 = -cp * self.root(
+            fwd_basket * np.exp(-(sigma_std ** 2) / 2), sigma_std, strike
+        )
 
         if np.isscalar(d2):
-            d1 = d2 + cp*sigma_std
+            d1 = d2 + cp * sigma_std
         else:
-            d1 = d2[:, None] + np.atleast_1d(cp)[:, None]*sigma_std
+            d1 = d2[:, None] + np.atleast_1d(cp)[:, None] * sigma_std
 
-        price = np.sum(fwd_basket*spst.norm.cdf(d1), axis=-1)
+        price = np.sum(fwd_basket * spst.norm.cdf(d1), axis=-1)
         price -= strike * spst.norm.cdf(d2)
-        price *= cp*df
+        price *= cp * df
         return price
 
 
 class BsmBasketLowerBound(NormBasket):
-
     def V_func(self, fwd):
         """
         Generate factor matrix V
@@ -407,7 +448,9 @@ def V_func(self, fwd):
         gg /= np.linalg.norm(gg)
 
         # equation 22，generate Q_1 and V_1
-        Q1 = self.chol_m.T @ gg / np.sqrt(gg @ self.cov_m @ gg)  # in py, this is a row vector
+        Q1 = (
+            self.chol_m.T @ gg / np.sqrt(gg @ self.cov_m @ gg)
+        )  # in py, this is a row vector
         V1 = self.chol_m @ Q1
         V1 = V1[:, None]
 
@@ -466,18 +509,52 @@ def moment_vsk(self, fwd, texp):
         m1 = sum(self.weight[i] * fwd[i] for i in range(n))
 
         m2_index = [i for i in product(np.arange(n), repeat=2)]
-        m2 = sum(self.weight[i] * self.weight[j] * fwd[i] * fwd[j] *
-                 np.exp(self.sigma[i] * self.sigma[j] * self.cor_m[i][j] * texp) for i, j in m2_index)
+        m2 = sum(
+            self.weight[i]
+            * self.weight[j]
+            * fwd[i]
+            * fwd[j]
+            * np.exp(self.sigma[i] * self.sigma[j] * self.cor_m[i][j] * texp)
+            for i, j in m2_index
+        )
 
         m3_index = [i for i in product(np.arange(n), repeat=3)]
-        m3 = sum(self.weight[i] * self.weight[j] * self.weight[l] * fwd[i] * fwd[j] * fwd[l] *
-                 np.exp(sum(self.sigma[ii] * self.sigma[jj] * self.cor_m[ii][jj] for ii, jj in
-                            combinations(np.array([i, j, l]), 2)) * texp) for i, j, l in m3_index)
+        m3 = sum(
+            self.weight[i]
+            * self.weight[j]
+            * self.weight[l]
+            * fwd[i]
+            * fwd[j]
+            * fwd[l]
+            * np.exp(
+                sum(
+                    self.sigma[ii] * self.sigma[jj] * self.cor_m[ii][jj]
+                    for ii, jj in combinations(np.array([i, j, l]), 2)
+                )
+                * texp
+            )
+            for i, j, l in m3_index
+        )
 
         m4_index = [i for i in product(np.arange(n), repeat=4)]
-        m4 = sum(self.weight[i] * self.weight[j] * self.weight[l] * self.weight[k] * fwd[i] * fwd[j] * fwd[l] *
-                 fwd[k] * np.exp(sum(self.sigma[ii] * self.sigma[jj] * self.cor_m[ii][jj] for ii, jj in
-                                      combinations(np.array([i, j, l, k]), 2)) * texp) for i, j, l, k in m4_index)
+        m4 = sum(
+            self.weight[i]
+            * self.weight[j]
+            * self.weight[l]
+            * self.weight[k]
+            * fwd[i]
+            * fwd[j]
+            * fwd[l]
+            * fwd[k]
+            * np.exp(
+                sum(
+                    self.sigma[ii] * self.sigma[jj] * self.cor_m[ii][jj]
+                    for ii, jj in combinations(np.array([i, j, l, k]), 2)
+                )
+                * texp
+            )
+            for i, j, l, k in m4_index
+        )
 
         var = m2 - m1 ** 2
         skew = (m3 - m1 ** 3 - 3 * m2 * m1 + 3 * m1 ** 3) / var ** (3 / 2)
@@ -505,7 +582,7 @@ def price(self, strike, spot, texp, cp=1):
         var, skew, kurt = self.moment_vsk(fwd, texp)
 
         m = nsvh.Nsvh1(sigma=self.sigma)
-        m.calibrate_vsk(var, skew, kurt-3, texp, setval=True)
+        m.calibrate_vsk(var, skew, kurt - 3, texp, setval=True)
         price = m.price(strike, fwd_basket, texp, cp)
 
         return df * price
diff --git a/pyfeng/multiasset_Ju2002.py b/pyfeng/multiasset_Ju2002.py
index e867c0c..f8d0bec 100644
--- a/pyfeng/multiasset_Ju2002.py
+++ b/pyfeng/multiasset_Ju2002.py
@@ -24,221 +24,392 @@ def __init__(self, sigma, cor=None, weight=None, intr=0.0, divr=0.0, is_fwd=Fals
             divr: vector of dividend/convenience yield (foreign interest rate) 0-D or (n_asset, ) array
             is_fwd: if True, treat `spot` as forward price. False by default.
         """
-        super().__init__(sigma, cor=cor, weight = weight, intr=intr, divr=divr, is_fwd=is_fwd)
+        super().__init__(
+            sigma, cor=cor, weight=weight, intr=intr, divr=divr, is_fwd=is_fwd
+        )
         global num_asset
         num_asset = len(self.weight)
-    
-    def average_s(self, spot, texp, basket = True):
-        #cal the forward price of asset num in the basket
-        if(basket):
+
+    def average_s(self, spot, texp, basket=True):
+        # cal the forward price of asset num in the basket
+        if basket:
             if np.isscalar(spot):
-                 spot = np.full(num_asset,spot)
+                spot = np.full(num_asset, spot)
             if np.isscalar(self.divr):
-                 self.divr = np.full(num_asset,self.divr)
+                self.divr = np.full(num_asset, self.divr)
             av_s = np.zeros(num_asset)
             for num in range(num_asset):
-                av_s[num] = (self.weight[num]*spot[num]*np.exp((self.intr-self.divr[num])*texp))
+                av_s[num] = (
+                    self.weight[num]
+                    * spot[num]
+                    * np.exp((self.intr - self.divr[num]) * texp)
+                )
         else:
             if np.isscalar(spot):
-                 spot = np.full(num_asset,spot)
+                spot = np.full(num_asset, spot)
             if np.isscalar(self.divr):
-                 self.divr = np.full(num_asset,self.divr)
+                self.divr = np.full(num_asset, self.divr)
             av_s = np.zeros(num_asset)
             for num in range(num_asset):
-                av_s[num] = (self.weight[num]*spot[num]*np.exp((self.intr-self.divr[num])*texp/(num_asset-1)*num))
+                av_s[num] = (
+                    self.weight[num]
+                    * spot[num]
+                    * np.exp(
+                        (self.intr - self.divr[num]) * texp / (num_asset - 1) * num
+                    )
+                )
         self.av_s = av_s
-    
-    def average_rho(self, texp, basket = True):
-        #cal the rho between asset i and j
-        if(basket):
-            av_rho = np.zeros((num_asset,num_asset))
+
+    def average_rho(self, texp, basket=True):
+        # cal the rho between asset i and j
+        if basket:
+            av_rho = np.zeros((num_asset, num_asset))
             for i in range(num_asset):
                 for j in range(num_asset):
-                    av_rho[i,j] = self.cor_m[i,j]*self.sigma[i]*self.sigma[j]*texp
+                    av_rho[i, j] = (
+                        self.cor_m[i, j] * self.sigma[i] * self.sigma[j] * texp
+                    )
         else:
-            av_rho = np.zeros((num_asset,num_asset))
+            av_rho = np.zeros((num_asset, num_asset))
             for i in range(num_asset):
-                for j in range(i,num_asset):
-                    av_rho[i,j] = self.sigma[0]**2*texp/(num_asset-1)*i
-                    av_rho[j,i] = av_rho[i,j]
+                for j in range(i, num_asset):
+                    av_rho[i, j] = self.sigma[0] ** 2 * texp / (num_asset - 1) * i
+                    av_rho[j, i] = av_rho[i, j]
         self.av_rho = av_rho
-    
+
     def u1(self, spot, texp):
-        #the first momentum of log normal distribution#
+        # the first momentum of log normal distribution#
         u1_value = self.av_s.sum()
         # u1_value = self.weight @ (spot * np.exp((self.intr-self.divr)*texp))
         return u1_value
-    
+
     def u2(self, z):
-        #the second momentum of log normal distribution#
+        # the second momentum of log normal distribution#
         u2_value = 0
         for i in range(num_asset):
             for j in range(num_asset):
-              u2_value += self.av_s[i] * self.av_s[j] *np.exp(z*z*self.av_rho[i,j])
+                u2_value += (
+                    self.av_s[i] * self.av_s[j] * np.exp(z * z * self.av_rho[i, j])
+                )
         return u2_value
-    
+
     def u2_1st_der(self):
         u2_1st_value = 0
         for i in range(num_asset):
             for j in range(num_asset):
-              u2_1st_value += self.av_s[i] * self.av_s[j] *self.av_rho[i,j]
+                u2_1st_value += self.av_s[i] * self.av_s[j] * self.av_rho[i, j]
         return u2_1st_value
-    
+
     def u2_2nd_der(self):
         u2_2nd_value = 0
         for i in range(num_asset):
             for j in range(num_asset):
-              u2_2nd_value += self.av_s[i] * self.av_s[j] *pow(self.av_rho[i,j],2)
+                u2_2nd_value += self.av_s[i] * self.av_s[j] * pow(self.av_rho[i, j], 2)
         return u2_2nd_value
-    
+
     def u2_3rd_der(self):
         u2_3rd_value = 0
         for i in range(num_asset):
             for j in range(num_asset):
-              u2_3rd_value += self.av_s[i] * self.av_s[j] *pow(self.av_rho[i,j],3)
+                u2_3rd_value += self.av_s[i] * self.av_s[j] * pow(self.av_rho[i, j], 3)
         return u2_3rd_value
-    
+
     def ak_bar(self):
-        #calculate the average a, save to self#
+        # calculate the average a, save to self#
         av_a = self.av_rho @ self.av_s
         self.av_a = av_a
-    
+
     def e_a12_a2(self):
-        return 2 * self.av_s @ pow(self.av_a,2)
-    
+        return 2 * self.av_s @ pow(self.av_a, 2)
+
     def e_a12_a22(self):
         value = 0
         for i in range(num_asset):
             for j in range(num_asset):
-                value += self.av_a[i]*self.av_s[i]*self.av_rho[i,j]*self.av_a[j]*self.av_s[j]
-        value *=8
-        value += 2*self.u2_1st_der()*self.u2_2nd_der()
+                value += (
+                    self.av_a[i]
+                    * self.av_s[i]
+                    * self.av_rho[i, j]
+                    * self.av_a[j]
+                    * self.av_s[j]
+                )
+        value *= 8
+        value += 2 * self.u2_1st_der() * self.u2_2nd_der()
         return value
-    
+
     def e_a13_a3(self):
-        return 6* self.av_s @ pow(self.av_a,3)
-    
+        return 6 * self.av_s @ pow(self.av_a, 3)
+
     def e_a1_a2_a3(self):
         value = 0
         for i in range(num_asset):
             for j in range(num_asset):
-                value += self.av_s[i]*pow(self.av_rho[i,j],2)*self.av_a[j]*self.av_s[j]
+                value += (
+                    self.av_s[i]
+                    * pow(self.av_rho[i, j], 2)
+                    * self.av_a[j]
+                    * self.av_s[j]
+                )
         value *= 6
         return value
-    
+
     def e_a23(self):
         value = 0
-        temp = np.zeros((num_asset,num_asset))
+        temp = np.zeros((num_asset, num_asset))
         for i in range(num_asset):
             for j in range(num_asset):
-                temp[i,j] = pow(self.av_s[i],0.5)*self.av_rho[i,j]*pow(self.av_s[j],0.5)
+                temp[i, j] = (
+                    pow(self.av_s[i], 0.5) * self.av_rho[i, j] * pow(self.av_s[j], 0.5)
+                )
         for i in range(num_asset):
             for j in range(num_asset):
                 for k in range(num_asset):
-                    value += temp[i,j]*temp[j,k]*temp[k,i]
+                    value += temp[i, j] * temp[j, k] * temp[k, i]
         value *= 8
         return value
-    
+
     def func_a1(self, z):
-        return -pow(z,2)*self.u2_1st_der()/2/self.u2(0)
-    
+        return -pow(z, 2) * self.u2_1st_der() / 2 / self.u2(0)
+
     def func_a2(self, z):
-        return 2*pow(self.func_a1(z),2)-pow(z,4)*self.u2_2nd_der()/2/self.u2(0)
-    
+        return 2 * pow(self.func_a1(z), 2) - pow(
+            z, 4
+        ) * self.u2_2nd_der() / 2 / self.u2(0)
+
     def func_a3(self, z):
-        return 6*self.func_a1(z)*self.func_a2(z)-4*pow(self.func_a1(z),3)-pow(z,6)*self.u2_3rd_der()/2/self.u2(0)
-    
+        return (
+            6 * self.func_a1(z) * self.func_a2(z)
+            - 4 * pow(self.func_a1(z), 3)
+            - pow(z, 6) * self.u2_3rd_der() / 2 / self.u2(0)
+        )
+
     def func_b1(self, spot, texp, z):
-        return pow(z,4)*self.e_a12_a2()/4/pow(self.u1(spot,texp),3)
-        
+        return pow(z, 4) * self.e_a12_a2() / 4 / pow(self.u1(spot, texp), 3)
+
     def func_b2(self, z):
-        return pow(self.func_a1(z),2)-self.func_a2(z)/2
-    
+        return pow(self.func_a1(z), 2) - self.func_a2(z) / 2
+
     def func_c1(self, spot, texp, z):
-        return -self.func_a1(z)*self.func_b1(spot, texp, z)
-    
+        return -self.func_a1(z) * self.func_b1(spot, texp, z)
+
     def func_c2(self, spot, texp, z):
-        return pow(z,6)*(9*self.e_a12_a22()+4*self.e_a13_a3())/144/pow(self.u1(spot,texp),4)
-    
+        return (
+            pow(z, 6)
+            * (9 * self.e_a12_a22() + 4 * self.e_a13_a3())
+            / 144
+            / pow(self.u1(spot, texp), 4)
+        )
+
     def func_c3(self, spot, texp, z):
-        return pow(z,6)*(4*self.e_a1_a2_a3()+self.e_a23())/48/pow(self.u1(spot,texp),3)
-    
+        return (
+            pow(z, 6)
+            * (4 * self.e_a1_a2_a3() + self.e_a23())
+            / 48
+            / pow(self.u1(spot, texp), 3)
+        )
+
     def func_c4(self, z):
-        return self.func_a1(z)*self.func_a2(z)-2*pow(self.func_a1(z),3)/3-self.func_a3(z)/6
-    
+        return (
+            self.func_a1(z) * self.func_a2(z)
+            - 2 * pow(self.func_a1(z), 3) / 3
+            - self.func_a3(z) / 6
+        )
+
     def func_d1(self, spot, texp, z):
-        return 0.5*(6*pow(self.func_a1(z),2)+self.func_a2(z)-4*self.func_b1(spot, texp, z)+2*self.func_b2(z))-1/6*(120*pow(self.func_a1(z),3)-self.func_a3(z)+6*(24*self.func_c1(spot, texp, z)-6*self.func_c2(spot, texp, z)+2*self.func_c3(spot, texp, z)-self.func_c4(z)))
-    
+        return 0.5 * (
+            6 * pow(self.func_a1(z), 2)
+            + self.func_a2(z)
+            - 4 * self.func_b1(spot, texp, z)
+            + 2 * self.func_b2(z)
+        ) - 1 / 6 * (
+            120 * pow(self.func_a1(z), 3)
+            - self.func_a3(z)
+            + 6
+            * (
+                24 * self.func_c1(spot, texp, z)
+                - 6 * self.func_c2(spot, texp, z)
+                + 2 * self.func_c3(spot, texp, z)
+                - self.func_c4(z)
+            )
+        )
+
     def func_d2(self, spot, texp, z):
-        return 0.5*(10*pow(self.func_a1(z),2)+self.func_a2(z)-6*self.func_b1(spot, texp, z)+2*self.func_b2(z))-(128*pow(self.func_a1(z),3)/3-self.func_a3(z)/6+2*self.func_a1(z)*self.func_b1(spot, texp, z)-self.func_a1(z)*self.func_b2(z)+50*self.func_c1(spot, texp, z)-11*self.func_c2(spot, texp, z)+3*self.func_c3(spot, texp, z)-self.func_c4(z))
-    
+        return 0.5 * (
+            10 * pow(self.func_a1(z), 2)
+            + self.func_a2(z)
+            - 6 * self.func_b1(spot, texp, z)
+            + 2 * self.func_b2(z)
+        ) - (
+            128 * pow(self.func_a1(z), 3) / 3
+            - self.func_a3(z) / 6
+            + 2 * self.func_a1(z) * self.func_b1(spot, texp, z)
+            - self.func_a1(z) * self.func_b2(z)
+            + 50 * self.func_c1(spot, texp, z)
+            - 11 * self.func_c2(spot, texp, z)
+            + 3 * self.func_c3(spot, texp, z)
+            - self.func_c4(z)
+        )
+
     def func_d3(self, spot, texp, z):
-        return 2*pow(self.func_a1(z),2)-self.func_b1(spot, texp, z)-1/3*(88*pow(self.func_a1(z),3)+3*self.func_a1(z)*(5*self.func_b1(spot, texp, z)-2*self.func_b2(z))+3*(35*self.func_c1(spot, texp, z)-6*self.func_c2(spot, texp, z)+self.func_c3(spot, texp, z)))
-                     
+        return (
+            2 * pow(self.func_a1(z), 2)
+            - self.func_b1(spot, texp, z)
+            - 1
+            / 3
+            * (
+                88 * pow(self.func_a1(z), 3)
+                + 3
+                * self.func_a1(z)
+                * (5 * self.func_b1(spot, texp, z) - 2 * self.func_b2(z))
+                + 3
+                * (
+                    35 * self.func_c1(spot, texp, z)
+                    - 6 * self.func_c2(spot, texp, z)
+                    + self.func_c3(spot, texp, z)
+                )
+            )
+        )
+
     def func_d4(self, spot, texp, z):
-        return -20*pow(self.func_a1(z),3)/3+self.func_a1(z)*(-4*self.func_b1(spot, texp, z)+self.func_b2(z))-10*self.func_c1(spot, texp, z)+self.func_c2(spot, texp, z)
+        return (
+            -20 * pow(self.func_a1(z), 3) / 3
+            + self.func_a1(z) * (-4 * self.func_b1(spot, texp, z) + self.func_b2(z))
+            - 10 * self.func_c1(spot, texp, z)
+            + self.func_c2(spot, texp, z)
+        )
 
-    def price(self, strike, spot, texp, cp=1, basket = True):
+    def price(self, strike, spot, texp, cp=1, basket=True):
         if np.isscalar(spot):
-             spot = np.full(num_asset,spot)
+            spot = np.full(num_asset, spot)
         if np.isscalar(self.divr):
-             self.divr = np.full(num_asset,self.divr)
-        if(basket):
+            self.divr = np.full(num_asset, self.divr)
+        if basket:
             self.average_s(spot, texp)
             self.average_rho(texp)
         else:
             self.average_s(spot, texp, False)
             self.average_rho(texp, False)
         self.ak_bar()
-        m1=2*np.log(self.u1(spot,texp))-0.5*np.log(self.u2(1))
-        v1=np.log(self.u2(1))-2*np.log(self.u1(spot,texp))
-        sqrtv1=np.sqrt(v1)
-        y=np.log(strike)
-        y1=(m1-y)/np.sqrt(v1)+sqrtv1
-        y2=y1-sqrtv1
-        z1=self.func_d2(spot,texp,1)-self.func_d3(spot,texp,1)+self.func_d4(spot,texp,1)
-        z2=self.func_d3(spot,texp,1)-self.func_d4(spot,texp,1)
-        z3=self.func_d4(spot,texp,1)
-        bc= self.u1(spot,texp) * np.exp(-self.intr*texp) * spst.norm.cdf(y1, loc=0, scale=1) - strike * np.exp(-self.intr * texp) * spst.norm.cdf(y2, loc=0, scale=1) + np.exp(-self.intr * texp) * strike * (z1 * spst.norm.pdf(y, loc=m1, scale=sqrtv1) + z2 * spst.norm.pdf(y, loc=m1, scale=sqrtv1) * (m1 - y) / v1 + z3 * ((y - m1) * (y - m1) / v1 / v1 - 1 / v1) * spst.norm.pdf(y, loc=m1, scale=sqrtv1))
+        m1 = 2 * np.log(self.u1(spot, texp)) - 0.5 * np.log(self.u2(1))
+        v1 = np.log(self.u2(1)) - 2 * np.log(self.u1(spot, texp))
+        sqrtv1 = np.sqrt(v1)
+        y = np.log(strike)
+        y1 = (m1 - y) / np.sqrt(v1) + sqrtv1
+        y2 = y1 - sqrtv1
+        z1 = (
+            self.func_d2(spot, texp, 1)
+            - self.func_d3(spot, texp, 1)
+            + self.func_d4(spot, texp, 1)
+        )
+        z2 = self.func_d3(spot, texp, 1) - self.func_d4(spot, texp, 1)
+        z3 = self.func_d4(spot, texp, 1)
+        bc = (
+            self.u1(spot, texp)
+            * np.exp(-self.intr * texp)
+            * spst.norm.cdf(y1, loc=0, scale=1)
+            - strike * np.exp(-self.intr * texp) * spst.norm.cdf(y2, loc=0, scale=1)
+            + np.exp(-self.intr * texp)
+            * strike
+            * (
+                z1 * spst.norm.pdf(y, loc=m1, scale=sqrtv1)
+                + z2 * spst.norm.pdf(y, loc=m1, scale=sqrtv1) * (m1 - y) / v1
+                + z3
+                * ((y - m1) * (y - m1) / v1 / v1 - 1 / v1)
+                * spst.norm.pdf(y, loc=m1, scale=sqrtv1)
+            )
+        )
         if cp == 1:
             return bc
         elif cp == -1:
-            return np.exp(-self.intr*texp)* (strike-self.u1(spot, texp))+bc
+            return np.exp(-self.intr * texp) * (strike - self.u1(spot, texp)) + bc
         else:
             return -1
+
+
 class BsmContinuousAsianJu2002(opt.OptABC):
     def price(self, strike, spot, texp, cp=1):
 
-        if (np.isscalar(spot) == False):
+        if np.isscalar(spot) == False:
             print("spot should not be array")
             return 0
-        elif (np.isscalar(self.divr) == False):
+        elif np.isscalar(self.divr) == False:
             print("dividend should not be array")
             return 0
         else:
-            g = self.intr-self.divr
-            gt = g*texp
-            u1 = spot*(np.exp(gt)-1)/g/texp
-            u2 = 2*spot**2*((np.exp((2*g+self.sigma**2)*texp)-1)/(2*g+self.sigma**2)-(np.exp(gt)-1)/g)/texp/texp/(g+self.sigma**2)
-            z1 = -pow(self.sigma,4)*texp**2*(1/45+gt/180-11*gt**2/15120-pow(gt,3)/2520+pow(gt,4)/113400)-pow(self.sigma,6)*pow(texp,3)*(1/11340-13*gt/30240-17*gt**2/226800+23*pow(gt,3)/453600+59*pow(gt,4)/5987520)
-            z2 = -pow(self.sigma,4)*texp**2*(1/90+gt/360-11*gt**2/30240-pow(gt,3)/5040+pow(gt,4)/226800)-pow(self.sigma,6)*pow(texp,3)*(31/22680-11*gt/60480-37*gt**2/151200-19*pow(gt,3)/302400+953*pow(gt,4)/59875200)
-            z3 = pow(self.sigma,6)*pow(texp,3)*(2/2835-gt/60480-2*gt**2/14175-17*pow(gt,3)/907200+13*pow(gt,4)/1247400)
-            m1 = 2*np.log(u1)-0.5*np.log(u2)
-            v1 = np.log(u2)-2*np.log(u1)
+            g = self.intr - self.divr
+            gt = g * texp
+            u1 = spot * (np.exp(gt) - 1) / g / texp
+            u2 = (
+                2
+                * spot ** 2
+                * (
+                    (np.exp((2 * g + self.sigma ** 2) * texp) - 1)
+                    / (2 * g + self.sigma ** 2)
+                    - (np.exp(gt) - 1) / g
+                )
+                / texp
+                / texp
+                / (g + self.sigma ** 2)
+            )
+            z1 = -pow(self.sigma, 4) * texp ** 2 * (
+                1 / 45
+                + gt / 180
+                - 11 * gt ** 2 / 15120
+                - pow(gt, 3) / 2520
+                + pow(gt, 4) / 113400
+            ) - pow(self.sigma, 6) * pow(texp, 3) * (
+                1 / 11340
+                - 13 * gt / 30240
+                - 17 * gt ** 2 / 226800
+                + 23 * pow(gt, 3) / 453600
+                + 59 * pow(gt, 4) / 5987520
+            )
+            z2 = -pow(self.sigma, 4) * texp ** 2 * (
+                1 / 90
+                + gt / 360
+                - 11 * gt ** 2 / 30240
+                - pow(gt, 3) / 5040
+                + pow(gt, 4) / 226800
+            ) - pow(self.sigma, 6) * pow(texp, 3) * (
+                31 / 22680
+                - 11 * gt / 60480
+                - 37 * gt ** 2 / 151200
+                - 19 * pow(gt, 3) / 302400
+                + 953 * pow(gt, 4) / 59875200
+            )
+            z3 = (
+                pow(self.sigma, 6)
+                * pow(texp, 3)
+                * (
+                    2 / 2835
+                    - gt / 60480
+                    - 2 * gt ** 2 / 14175
+                    - 17 * pow(gt, 3) / 907200
+                    + 13 * pow(gt, 4) / 1247400
+                )
+            )
+            m1 = 2 * np.log(u1) - 0.5 * np.log(u2)
+            v1 = np.log(u2) - 2 * np.log(u1)
             sqrtv1 = np.sqrt(v1)
             y = np.log(strike)
-            y1 = (m1-y)/np.sqrt(v1)+sqrtv1
-            y2 = y1-sqrtv1
-            bc = u1 * np.exp(-self.intr*texp) * spst.norm.cdf(y1, loc=0, scale=1) - strike * np.exp(-self.intr * texp) * spst.norm.cdf(y2, loc=0, scale=1) + np.exp(-self.intr * texp) * strike * (z1 * spst.norm.pdf(y, loc=m1, scale=sqrtv1) + z2 * spst.norm.pdf(y, loc=m1, scale=sqrtv1) * (m1 - y) / v1 + z3 * ((y - m1) * (y - m1) / v1 / v1 - 1 / v1) * spst.norm.pdf(y, loc=m1, scale=sqrtv1))
+            y1 = (m1 - y) / np.sqrt(v1) + sqrtv1
+            y2 = y1 - sqrtv1
+            bc = (
+                u1 * np.exp(-self.intr * texp) * spst.norm.cdf(y1, loc=0, scale=1)
+                - strike * np.exp(-self.intr * texp) * spst.norm.cdf(y2, loc=0, scale=1)
+                + np.exp(-self.intr * texp)
+                * strike
+                * (
+                    z1 * spst.norm.pdf(y, loc=m1, scale=sqrtv1)
+                    + z2 * spst.norm.pdf(y, loc=m1, scale=sqrtv1) * (m1 - y) / v1
+                    + z3
+                    * ((y - m1) * (y - m1) / v1 / v1 - 1 / v1)
+                    * spst.norm.pdf(y, loc=m1, scale=sqrtv1)
+                )
+            )
         if cp == 1:
             return bc
         elif cp == -1:
-            return np.exp(-self.intr*texp)* (strike-u1)+bc
+            return np.exp(-self.intr * texp) * (strike - u1) + bc
         else:
             return -1
-        
-        
-        
-        
-        
-        
diff --git a/pyfeng/multiasset_mc.py b/pyfeng/multiasset_mc.py
index 3d2be8e..cdd24c7 100644
--- a/pyfeng/multiasset_mc.py
+++ b/pyfeng/multiasset_mc.py
@@ -23,7 +23,7 @@ class BsmNdMc(opt.OptMaABC):
     """
 
     spot = np.ones(2)
-    sigma = np.ones(2)*0.1
+    sigma = np.ones(2) * 0.1
 
     # MC params
     n_path = 0
@@ -71,12 +71,18 @@ def _bm_incr(self, tobs, n_path=None):
         if self.antithetic:
             # generate random number in the order of path, time, asset and transposed
             # in this way, the same paths are generated when increasing n_path
-            bm_incr = self.rng.normal(size=(n_path//2, n_t, self.n_asset)).transpose((1, 0, 2))
+            bm_incr = self.rng.normal(size=(n_path // 2, n_t, self.n_asset)).transpose(
+                (1, 0, 2)
+            )
             bm_incr *= np.sqrt(dt[:, None, None])
             bm_incr = np.dot(bm_incr, self.chol_m.T)
-            bm_incr = np.stack([bm_incr, -bm_incr], axis=2).reshape((n_t, n_path, self.n_asset))
+            bm_incr = np.stack([bm_incr, -bm_incr], axis=2).reshape(
+                (n_t, n_path, self.n_asset)
+            )
         else:
-            bm_incr = np.random.randn(n_path, n_t, self.n_asset).transpose((1, 0, 2)) * np.sqrt(dt[:, None, None])
+            bm_incr = np.random.randn(n_path, n_t, self.n_asset).transpose(
+                (1, 0, 2)
+            ) * np.sqrt(dt[:, None, None])
             bm_incr = np.dot(bm_incr, self.chol_m.T)
 
         return bm_incr
@@ -100,12 +106,12 @@ def simulate(self, n_path=None, tobs=None, store=1):
         elif store == 2 and tobs is not None:
             # make sure that tobs == self.tobs
             if not np.all(np.isclose(self.tobs, tobs)):
-                raise ValueError('tobs is different from the saved value.')
+                raise ValueError("tobs is different from the saved value.")
 
         path = self._bm_incr(self.tobs, n_path)
         # Add drift and convexity
         dt = np.diff(self.tobs, prepend=0)
-        path += (self.intr - self.divr - 0.5*self.sigma**2)*dt[:, None, None]
+        path += (self.intr - self.divr - 0.5 * self.sigma ** 2) * dt[:, None, None]
         np.cumsum(path, axis=0, out=path)
         np.exp(path, out=path)
 
@@ -131,14 +137,19 @@ def price_european(self, spot, texp, payoff):
             The MC price of the payoff
         """
         if self.n_path == 0:
-            raise ValueError('Simulated paths are not available. Run simulate() first.')
+            raise ValueError("Simulated paths are not available. Run simulate() first.")
 
         # check if texp is in tobs
         ind, *_ = np.where(np.isclose(self.tobs, texp))
         if len(ind) == 0:
-            raise ValueError(f'Stored tobs does not contain t={texp}')
-
-        path = self.path[ind[0], ] * spot
+            raise ValueError(f"Stored tobs does not contain t={texp}")
+
+        path = (
+            self.path[
+                ind[0],
+            ]
+            * spot
+        )
         price = np.exp(-self.intr * texp) * np.mean(payoff(path), axis=0)
         return price
 
@@ -176,13 +187,18 @@ def simulate(self, tobs, n_path=None, store=True):
 
     def price_european(self, spot, texp, payoff):
         if self.path is None:
-            raise ValueError('Simulated paths are not available. Run simulate() first.')
+            raise ValueError("Simulated paths are not available. Run simulate() first.")
 
         # check if texp is in tobs
         ind, *_ = np.where(np.isclose(self.tobs, texp))
         if len(ind) == 0:
-            raise ValueError(f'Stored path does not contain t = {texp}')
-
-        path = self.path[ind[0], ] + spot
+            raise ValueError(f"Stored path does not contain t = {texp}")
+
+        path = (
+            self.path[
+                ind[0],
+            ]
+            + spot
+        )
         price = np.exp(-self.intr * texp) * np.mean(payoff(path), axis=0)
-        return price
\ No newline at end of file
+        return price
diff --git a/pyfeng/norm.py b/pyfeng/norm.py
index a937f3f..965fdd3 100644
--- a/pyfeng/norm.py
+++ b/pyfeng/norm.py
@@ -32,20 +32,37 @@ class Norm(opt.OptAnalyticABC):
 
     # Coefficients for _impvol_Choi2009
     _POLY_NU = [
-        1.266458051348246e4, 2.493415285349361e4, 6.106322407867059e3, 1.848489695437094e3,
-        5.988761102690991e2, 4.980340217855084e1, 2.100960795068497e1, 3.994961687345134e-1]
+        1.266458051348246e4,
+        2.493415285349361e4,
+        6.106322407867059e3,
+        1.848489695437094e3,
+        5.988761102690991e2,
+        4.980340217855084e1,
+        2.100960795068497e1,
+        3.994961687345134e-1,
+    ]
 
     _POLY_DE = [
-        1.174240599306013e1, -2.067719486400926e2, 3.608817108375034e3, 2.392008891720782e4,
-        1.598919697679745e4, 1.323614537899738e3, 1.495105008310999e3, 3.093573936743112e1,
-        4.990534153589422e1, 1.0]
+        1.174240599306013e1,
+        -2.067719486400926e2,
+        3.608817108375034e3,
+        2.392008891720782e4,
+        1.598919697679745e4,
+        1.323614537899738e3,
+        1.495105008310999e3,
+        3.093573936743112e1,
+        4.990534153589422e1,
+        1.0,
+    ]
 
     # option value below the intrinsic value by IMPVOL_TOL is considered to be numerical error.
     # volatility will be set to 0
-    IMPVOL_TOL = 1000*np.finfo(float).eps
+    IMPVOL_TOL = 1000 * np.finfo(float).eps
 
     @staticmethod
-    def price_formula(strike, spot, sigma, texp, cp=1, intr=0.0, divr=0.0, is_fwd=False):
+    def price_formula(
+        strike, spot, sigma, texp, cp=1, intr=0.0, divr=0.0, is_fwd=False
+    ):
         """
         Bachelier model call/put option pricing formula (static method)
 
@@ -70,7 +87,9 @@ def price_formula(strike, spot, sigma, texp, cp=1, intr=0.0, divr=0.0, is_fwd=Fa
         d = (fwd - strike) / sigma_std
 
         cp = np.array(cp)
-        price = df * (cp*(fwd - strike)*spst.norm.cdf(cp * d) + sigma_std * spst.norm.pdf(d))
+        price = df * (
+            cp * (fwd - strike) * spst.norm.cdf(cp * d) + sigma_std * spst.norm.pdf(d)
+        )
         return price
 
     def _impvol_Choi2009(self, price, strike, spot, texp, cp=1, setval=False):
@@ -95,24 +114,38 @@ def _impvol_Choi2009(self, price, strike, spot, texp, cp=1, setval=False):
         """
         fwd, df, _ = self._fwd_factor(spot, texp)
         price_fwd = np.array(price) / df
-        strike_std = np.array(cp)*(fwd - strike)
+        strike_std = np.array(cp) * (fwd - strike)
 
         time_val = price_fwd - np.maximum(0, strike_std)  # option time value
-        strd = 2*price_fwd - strike_std  # straddle value (=call + put)
+        strd = 2 * price_fwd - strike_std  # straddle value (=call + put)
         # Note: time_val > 0  => strd >= time_val > 0
         # v = |fwd - strike|/(call + put) = np.fabs(strike_std) / strd, ATM when v=0
         # v1 = 1 - v, ATM when v1=1
         # Use v1 instead of v to preserve the option value < machine epsilon
-        v1 = np.clip(2*time_val/strd, 0, 1)  # bound between 0 and 1 for now. Out-of-bound value will be handleded later
-        v_sq = (1 - v1)**2
+        v1 = np.clip(
+            2 * time_val / strd, 0, 1
+        )  # bound between 0 and 1 for now. Out-of-bound value will be handleded later
+        v_sq = (1 - v1) ** 2
 
         # ignore 'divide by 0' error when v1 = 0 for now.
         # eta = v / atanh(v) = 2v / log((1+v)/(1-v)) = 2v / log((2-v1)/v1)
-        with np.errstate(divide='ignore', invalid='ignore'):
-            eta = np.where(v1 < 0.999, 2*(1-v1)/(np.log((2.0-v1)/v1)), 1/(1 + v_sq*(1/3 + v_sq/5)))
-        h_a = np.sqrt(eta) * np.polyval(Norm._POLY_NU, eta) / np.polyval(Norm._POLY_DE, eta)
+        with np.errstate(divide="ignore", invalid="ignore"):
+            eta = np.where(
+                v1 < 0.999,
+                2 * (1 - v1) / (np.log((2.0 - v1) / v1)),
+                1 / (1 + v_sq * (1 / 3 + v_sq / 5)),
+            )
+        h_a = (
+            np.sqrt(eta)
+            * np.polyval(Norm._POLY_NU, eta)
+            / np.polyval(Norm._POLY_DE, eta)
+        )
         # sigma = sqrt(pi/2T) * (call + put) * h_a
-        _sigma = np.where(time_val >= -self.IMPVOL_TOL, np.sqrt(np.pi/(2*texp)) * strd * h_a, np.nan)
+        _sigma = np.where(
+            time_val >= -self.IMPVOL_TOL,
+            np.sqrt(np.pi / (2 * texp)) * strd * h_a,
+            np.nan,
+        )
 
         if setval:
             self.sigma = _sigma
@@ -126,7 +159,9 @@ def vega(self, strike, spot, texp, cp=1):
         sigma_std = np.maximum(self.sigma * np.sqrt(texp), np.finfo(float).eps)
         d = (fwd - strike) / sigma_std
 
-        vega = df * spst.norm.pdf(d) * np.sqrt(texp)  # formula according to lecture notes
+        vega = (
+            df * spst.norm.pdf(d) * np.sqrt(texp)
+        )  # formula according to lecture notes
         return vega
 
     def delta(self, strike, spot, texp, cp=1):
@@ -156,9 +191,9 @@ def gamma(self, strike, spot, texp, cp=1):
         sigma_std = np.maximum(self.sigma * np.sqrt(texp), np.finfo(float).eps)
         d = (fwd - strike) / sigma_std
 
-        gamma = df * spst.norm.pdf(d)/sigma_std  # formula according to wikipedia
+        gamma = df * spst.norm.pdf(d) / sigma_std  # formula according to wikipedia
         if not self.is_fwd:
-            gamma *= (divf/df)**2
+            gamma *= (divf / df) ** 2
         return gamma
 
     def theta(self, strike, spot, texp, cp=1):
@@ -169,18 +204,22 @@ def theta(self, strike, spot, texp, cp=1):
         d = (fwd - strike) / sigma_std
 
         # still not perfect; need to consider the derivative w.r.t. divr and is_fwd = True
-        theta = sigma_std*spst.norm.pdf(d)
-        theta = -0.5*theta/texp + self.intr*(theta - cp*strike*spst.norm.cdf(cp*d))
+        theta = sigma_std * spst.norm.pdf(d)
+        theta = -0.5 * theta / texp + self.intr * (
+            theta - cp * strike * spst.norm.cdf(cp * d)
+        )
         return df * theta
 
-    def price_binary(self, strike, spot, texp, cp=1, opt_type='cash'):
+    def price_binary(self, strike, spot, texp, cp=1, opt_type="cash"):
         fwd, df, divf = self._fwd_factor(spot, texp)
 
         sigma_std = np.maximum(self.sigma * np.sqrt(texp), np.finfo(float).eps)
         d = (fwd - strike) / sigma_std
 
-        if opt_type.lower() == 'asset':
-            price = df * (cp * fwd * spst.norm.cdf(cp * d) + sigma_std * spst.norm.pdf(d))
+        if opt_type.lower() == "asset":
+            price = df * (
+                cp * fwd * spst.norm.cdf(cp * d) + sigma_std * spst.norm.pdf(d)
+            )
         else:
             price = df * spst.norm.cdf(cp * d)
 
@@ -189,7 +228,7 @@ def price_binary(self, strike, spot, texp, cp=1, opt_type='cash'):
     ####
     impvol = _impvol_Choi2009
 
-    def vol_smile(self, strike, spot, texp, cp=1, model='bsm'):
+    def vol_smile(self, strike, spot, texp, cp=1, model="bsm"):
         """
         Equivalent volatility smile for a given model
 
@@ -203,21 +242,21 @@ def vol_smile(self, strike, spot, texp, cp=1, model='bsm'):
         Returns:
             volatility smile under the specified model
         """
-        if model.lower() == 'norm':
+        if model.lower() == "norm":
             return self.sigma * np.ones_like(strike + spot + texp + cp)
-        if model.lower() == 'bsm':
+        if model.lower() == "bsm":
             price = self.price(strike, spot, texp, cp=cp)
             return bsm.Bsm(None).impvol(price, strike, spot, texp, cp=cp)
-        elif model.lower() == 'bsm-approx':
+        elif model.lower() == "bsm-approx":
             fwd, _, _ = self._fwd_factor(spot, texp)
             sigma_std = self.sigma / fwd
             kk = strike / fwd
             lnk = np.log(kk)
             vol = sigma_std / np.sqrt(kk)
-            vol *= (1 + vol**2 * texp/24) / (1 + lnk**2/24)
+            vol *= (1 + vol ** 2 * texp / 24) / (1 + lnk ** 2 / 24)
             return vol
         else:
-            raise ValueError(f'Unknown model: {model}')
+            raise ValueError(f"Unknown model: {model}")
 
     def _price_suboptimal(self, strike, spot, texp, cp=1, strike2=None):
         fwd, df, _ = self._fwd_factor(spot, texp)
@@ -226,7 +265,9 @@ def _price_suboptimal(self, strike, spot, texp, cp=1, strike2=None):
         sigma_std = np.maximum(self.sigma * np.sqrt(texp), np.finfo(float).eps)
         d2 = (fwd - strike2) / sigma_std
 
-        price = cp*(fwd - strike)*spst.norm.cdf(cp * d2) + sigma_std*spst.norm.pdf(d2)
+        price = cp * (fwd - strike) * spst.norm.cdf(
+            cp * d2
+        ) + sigma_std * spst.norm.pdf(d2)
         price *= df
         return price
 
@@ -242,7 +283,7 @@ def _barrier_params(self, barrier, spot):
             barrier option pricing parameters (psi, spot_mirror)
         """
         psi = 1  # need to fix this
-        spot_mirror = 2*barrier - spot
+        spot_mirror = 2 * barrier - spot
         return psi, spot_mirror
 
     """
diff --git a/pyfeng/nsvh.py b/pyfeng/nsvh.py
index 0320a6c..a9fb995 100644
--- a/pyfeng/nsvh.py
+++ b/pyfeng/nsvh.py
@@ -23,7 +23,17 @@ class Nsvh1(sabr.SabrABC):
     beta = 0.0  # beta is already defined in the parent class, but the default value set as 0
     is_atmvol = False
 
-    def __init__(self, sigma, vov=0.0, rho=0.0, beta=None, intr=0.0, divr=0.0, is_fwd=False, is_atmvol=False):
+    def __init__(
+        self,
+        sigma,
+        vov=0.0,
+        rho=0.0,
+        beta=None,
+        intr=0.0,
+        divr=0.0,
+        is_fwd=False,
+        is_atmvol=False,
+    ):
         """
         Args:
             sigma: model volatility at t=0
@@ -37,22 +47,27 @@ def __init__(self, sigma, vov=0.0, rho=0.0, beta=None, intr=0.0, divr=0.0, is_fw
         """
         # Make sure beta = 0
         if beta is not None and not np.isclose(beta, 0.0):
-            print(f'Ignoring beta = {beta}...')
+            print(f"Ignoring beta = {beta}...")
         self.is_atmvol = is_atmvol
         super().__init__(sigma, vov, rho, beta=0, intr=intr, divr=divr, is_fwd=is_fwd)
 
     def _sig0_from_atmvol(self, texp):
         s_sqrt = self.vov * np.sqrt(texp)
-        vov_var = np.exp(0.5 * s_sqrt**2)
-        rhoc = np.sqrt(1 - self.rho**2)
+        vov_var = np.exp(0.5 * s_sqrt ** 2)
+        rhoc = np.sqrt(1 - self.rho ** 2)
 
-        d = (np.arctanh(self.rho) - np.arcsinh(self.rho*vov_var/rhoc)) / s_sqrt
+        d = (np.arctanh(self.rho) - np.arcsinh(self.rho * vov_var / rhoc)) / s_sqrt
         ncdf_p = spst.norm.cdf(d + s_sqrt)
         ncdf_m = spst.norm.cdf(d - s_sqrt)
         ncdf = spst.norm.cdf(d)
 
-        price = 0.5/self.vov*vov_var * ((1+self.rho)*ncdf_p - (1-self.rho)*ncdf_m - 2*self.rho*ncdf)
-        sig0 = self.sigma*np.sqrt(texp/2/np.pi)/price
+        price = (
+            0.5
+            / self.vov
+            * vov_var
+            * ((1 + self.rho) * ncdf_p - (1 - self.rho) * ncdf_m - 2 * self.rho * ncdf)
+        )
+        sig0 = self.sigma * np.sqrt(texp / 2 / np.pi) / price
         return sig0
 
     def price(self, strike, spot, texp, cp=1):
@@ -66,17 +81,27 @@ def price(self, strike, spot, texp, cp=1):
 
         sig_sqrt = sig0 * np.sqrt(texp)
 
-        vov_var = np.exp(0.5 * s_sqrt**2)
-        rhoc = np.sqrt(1 - self.rho**2)
-
-        d = (np.arctanh(self.rho) + np.arcsinh(((fwd-strike)*s_sqrt/sig_sqrt - self.rho*vov_var)/rhoc)) / s_sqrt
-        ncdf_p = spst.norm.cdf(cp*(d + s_sqrt))
-        ncdf_m = spst.norm.cdf(cp*(d - s_sqrt))
-        ncdf = spst.norm.cdf(cp*d)
-
-        price = 0.5*sig_sqrt/s_sqrt*vov_var\
-            * ((1+self.rho)*ncdf_p - (1-self.rho)*ncdf_m - 2*self.rho*ncdf)\
-            + (fwd-strike) * ncdf
+        vov_var = np.exp(0.5 * s_sqrt ** 2)
+        rhoc = np.sqrt(1 - self.rho ** 2)
+
+        d = (
+            np.arctanh(self.rho)
+            + np.arcsinh(
+                ((fwd - strike) * s_sqrt / sig_sqrt - self.rho * vov_var) / rhoc
+            )
+        ) / s_sqrt
+        ncdf_p = spst.norm.cdf(cp * (d + s_sqrt))
+        ncdf_m = spst.norm.cdf(cp * (d - s_sqrt))
+        ncdf = spst.norm.cdf(cp * d)
+
+        price = (
+            0.5
+            * sig_sqrt
+            / s_sqrt
+            * vov_var
+            * ((1 + self.rho) * ncdf_p - (1 - self.rho) * ncdf_m - 2 * self.rho * ncdf)
+            + (fwd - strike) * ncdf
+        )
         price *= cp * df
         return price
 
@@ -85,12 +110,16 @@ def cdf(self, strike, spot, texp, cp=-1):
 
         s_sqrt = self.vov * np.sqrt(texp)
         sig_sqrt = self.sigma * np.sqrt(texp)
-        vov_var = np.exp(0.5 * s_sqrt**2)
-        rhoc = np.sqrt(1 - self.rho**2)
+        vov_var = np.exp(0.5 * s_sqrt ** 2)
+        rhoc = np.sqrt(1 - self.rho ** 2)
 
-        d = (np.arctanh(self.rho) +
-             np.arcsinh(((fwd - strike) * s_sqrt / sig_sqrt - self.rho * vov_var) / rhoc)) / s_sqrt
-        return spst.norm.cdf(cp*d)
+        d = (
+            np.arctanh(self.rho)
+            + np.arcsinh(
+                ((fwd - strike) * s_sqrt / sig_sqrt - self.rho * vov_var) / rhoc
+            )
+        ) / s_sqrt
+        return spst.norm.cdf(cp * d)
 
     def moments_vsk(self, texp=1):
         """
@@ -102,34 +131,40 @@ def moments_vsk(self, texp=1):
         Returns:
             (variance, skewness, and ex-kurtosis)
         """
-        vol_std = self.vov*np.sqrt(texp)
-        ww = np.exp(vol_std**2)
-        rho2 = self.rho**2
-        rho3 = self.rho*rho2
-        rho4 = rho2**2
+        vol_std = self.vov * np.sqrt(texp)
+        ww = np.exp(vol_std ** 2)
+        rho2 = self.rho ** 2
+        rho3 = self.rho * rho2
+        rho4 = rho2 ** 2
 
         c20 = ww + 1
         c22 = ww - 1
 
-        c31 = 3*(ww + 1)**2
-        c33 = (ww - 1)*(ww + 3)
+        c31 = 3 * (ww + 1) ** 2
+        c33 = (ww - 1) * (ww + 3)
 
-        #C40 = (ww + 1)**2*(ww**4 + 2*ww**2 + 3)
-        #C42 = 6*(ww**2 - 1)*(ww**4 + 2*ww**3 + 4*ww**2 + 2*ww + 1)
-        #C44 = (ww - 1)**2*(ww**4 + 4*ww**3 + 10*ww**2 + 12*ww + 3)
-        #M3 = (1 / 4)*np.sqrt(ww)*(ww - 1)**2*(self.rho*C_31 + rho3*C_33)
-        #M4 = (1 / 8)*(ww - 1)**2*(C_40 + rho2*C_42 + rho4*C_44)
+        # C40 = (ww + 1)**2*(ww**4 + 2*ww**2 + 3)
+        # C42 = 6*(ww**2 - 1)*(ww**4 + 2*ww**3 + 4*ww**2 + 2*ww + 1)
+        # C44 = (ww - 1)**2*(ww**4 + 4*ww**3 + 10*ww**2 + 12*ww + 3)
+        # M3 = (1 / 4)*np.sqrt(ww)*(ww - 1)**2*(self.rho*C_31 + rho3*C_33)
+        # M4 = (1 / 8)*(ww - 1)**2*(C_40 + rho2*C_42 + rho4*C_44)
 
-        m2 = (1/2)*(ww - 1)*(c20 + rho2*c22)
+        m2 = (1 / 2) * (ww - 1) * (c20 + rho2 * c22)
 
-        k0 = (ww + 1)**3*(ww**2 + 3)
-        k2 = 6*(ww + 1)**2*(ww**3 + ww**2 + 3*ww - 1)
-        k4 = (ww - 1)*(ww**4 + 4*ww**3 + 10*ww**2 + 12*ww - 3)
+        k0 = (ww + 1) ** 3 * (ww ** 2 + 3)
+        k2 = 6 * (ww + 1) ** 2 * (ww ** 3 + ww ** 2 + 3 * ww - 1)
+        k4 = (ww - 1) * (ww ** 4 + 4 * ww ** 3 + 10 * ww ** 2 + 12 * ww - 3)
 
-        skew = np.sqrt(ww*(ww - 1) / 2)*(self.rho*c31 + rho3*c33) / np.power(c20 + rho2*c22, 1.5)
-        exkurt = (1/2)*(ww - 1)*(k0 + rho2*k2 + rho4*k4) / (c20 + rho2*c22)**2
+        skew = (
+            np.sqrt(ww * (ww - 1) / 2)
+            * (self.rho * c31 + rho3 * c33)
+            / np.power(c20 + rho2 * c22, 1.5)
+        )
+        exkurt = (
+            (1 / 2) * (ww - 1) * (k0 + rho2 * k2 + rho4 * k4) / (c20 + rho2 * c22) ** 2
+        )
 
-        return m2*(self.sigma/self.vov)**2, skew, exkurt
+        return m2 * (self.sigma / self.vov) ** 2, skew, exkurt
 
     def calibrate_vsk(self, var, skew, exkurt, texp=1, setval=False):
         """
@@ -147,36 +182,40 @@ def calibrate_vsk(self, var, skew, exkurt, texp=1, setval=False):
             Tuenter, H. J. H. (2001). An algorithm to determine the parameters of SU-curves in the johnson system of probabillity distributions by moment matching. Journal of Statistical Computation and Simulation, 70(4), 325–347. https://doi.org/10.1080/00949650108812126
         """
         assert exkurt > 0
-        beta1 = skew**2
+        beta1 = skew ** 2
         beta2 = exkurt + 3
 
         # min of w search
-        roots = np.roots(np.array([1, 2, 3, 0, -3-beta2]))
+        roots = np.roots(np.array([1, 2, 3, 0, -3 - beta2]))
         roots = roots[(roots.real > 0) & np.isclose(roots.imag, 0)]
         assert len(roots) == 1
         w_min = roots.real[0]
-        w_max = np.sqrt(-1 + np.sqrt(2*(beta2 - 1)))
+        w_max = np.sqrt(-1 + np.sqrt(2 * (beta2 - 1)))
 
         def f_beta1(w):
-            term1 = np.sqrt(4 + 2*(w*w - (beta2 + 3) / (w*w + 2*w + 3)))
-            return (w + 1 - term1)*(w + 1 + 0.5*term1)**2 - beta1
+            term1 = np.sqrt(4 + 2 * (w * w - (beta2 + 3) / (w * w + 2 * w + 3)))
+            return (w + 1 - term1) * (w + 1 + 0.5 * term1) ** 2 - beta1
 
         assert f_beta1(w_min) >= 0
-        #print(w_min, f_beta1(w_min), w_max, f_beta1(w_max))
+        # print(w_min, f_beta1(w_min), w_max, f_beta1(w_max))
 
         # root finding for w = exp(S) = exp(vov^2 texp)
         w_root = spop.brentq(f_beta1, w_min, w_max)
-        m = -2 + np.sqrt(4 + 2*(w_root**2 - (beta2 + 3)/(w_root**2 + 2*w_root + 3)))
-        term = (w_root + 1)/(2*w_root)*((w_root - 1)/m - 1)  # - sinh(Omega) = rho / rho_*
+        m = -2 + np.sqrt(
+            4 + 2 * (w_root ** 2 - (beta2 + 3) / (w_root ** 2 + 2 * w_root + 3))
+        )
+        term = (
+            (w_root + 1) / (2 * w_root) * ((w_root - 1) / m - 1)
+        )  # - sinh(Omega) = rho / rho_*
 
         # if term is slightly negative, next line error in sqrt
-        if abs(term) < np.finfo(float).eps*100:
+        if abs(term) < np.finfo(float).eps * 100:
             term = 0.0
 
-        rho = np.sign(skew)*np.sqrt(1 - 1/(1 + term))
-        vov = np.sqrt(np.log(w_root)/texp)
-        m2 = 0.5*(w_root - 1)*((w_root + 1) + rho**2*(w_root - 1))
-        sig0 = np.sqrt(var/m2)*vov
+        rho = np.sign(skew) * np.sqrt(1 - 1 / (1 + term))
+        vov = np.sqrt(np.log(w_root) / texp)
+        m2 = 0.5 * (w_root - 1) * ((w_root + 1) + rho ** 2 * (w_root - 1))
+        sig0 = np.sqrt(var / m2) * vov
 
         if setval:
             self.sigma = sig0
@@ -201,7 +240,17 @@ class NsvhMc(sabr.SabrABC):
     rng = np.random.default_rng(None)
     antithetic = True
 
-    def __init__(self, sigma, vov=0.0, rho=0.0, lam=0, beta=None, intr=0.0, divr=0.0, is_fwd=False):
+    def __init__(
+        self,
+        sigma,
+        vov=0.0,
+        rho=0.0,
+        lam=0,
+        beta=None,
+        intr=0.0,
+        divr=0.0,
+        is_fwd=False,
+    ):
         """
         Args:
             sigma: model volatility at t=0
@@ -215,7 +264,7 @@ def __init__(self, sigma, vov=0.0, rho=0.0, lam=0, beta=None, intr=0.0, divr=0.0
         """
         # Make sure beta = 0
         if beta is not None and not np.isclose(beta, 0.0):
-            print(f'Ignoring beta = {beta}...')
+            print(f"Ignoring beta = {beta}...")
         self.lam = lam
         super().__init__(sigma, vov, rho, beta=0, intr=intr, divr=divr, is_fwd=is_fwd)
 
@@ -239,25 +288,34 @@ def mc_vol_price(self, texp):
         # forward path starting from zero
         # returns both sigma and price
 
-        rhoc = np.sqrt(1-self.rho**2)
-        vol_var = self.vov**2 * texp
+        rhoc = np.sqrt(1 - self.rho ** 2)
+        vol_var = self.vov ** 2 * texp
         vol_std = np.sqrt(vol_var)
 
-        z_rn = self.rng.normal(size=(int(self.n_path/2), 3))
+        z_rn = self.rng.normal(size=(int(self.n_path / 2), 3))
         z_rn = np.stack([z_rn, -z_rn], axis=1).reshape((-1, 3))
-        z_rn[:, 2] += 0.5*(self.lam-1)*vol_std  # add shift
+        z_rn[:, 2] += 0.5 * (self.lam - 1) * vol_std  # add shift
 
-        r2 = np.sum(z_rn[:, 0:2]**2, axis=1)
+        r2 = np.sum(z_rn[:, 0:2] ** 2, axis=1)
         exp_plus = np.exp(0.5 * vol_std * z_rn[:, 2])
 
-        phi_r1 = np.sqrt(2/r2) * np.sqrt(np.cosh(np.sqrt(r2 + z_rn[:, 2]**2)*vol_std) - np.cosh(z_rn[:, 2]*vol_std))
+        phi_r1 = np.sqrt(2 / r2) * np.sqrt(
+            np.cosh(np.sqrt(r2 + z_rn[:, 2] ** 2) * vol_std)
+            - np.cosh(z_rn[:, 2] * vol_std)
+        )
 
-        df_z = exp_plus**2
-        df_w = exp_plus[:, None] * z_rn[:, 0:2] * phi_r1[:, None]  # use both X and Y components
+        df_z = exp_plus ** 2
+        df_w = (
+            exp_plus[:, None] * z_rn[:, 0:2] * phi_r1[:, None]
+        )  # use both X and Y components
 
         path = (
             df_z,
-            (self.sigma/self.vov) * (self.rho*(df_z[:, None] - np.exp(0.5*self.lam*vol_var)) + rhoc*df_w)
+            (self.sigma / self.vov)
+            * (
+                self.rho * (df_z[:, None] - np.exp(0.5 * self.lam * vol_var))
+                + rhoc * df_w
+            ),
         )
         return path
 
@@ -281,7 +339,12 @@ def price(self, strike, spot, texp, cp=1):
         scalar_output = np.isscalar(strike_std)
         cp *= np.ones_like(strike_std)
 
-        price = np.array([np.mean(np.fmax(cp[k] * (mc_path[1] - strike_std[k]), 0)) for k in range(len(strike_std))])
+        price = np.array(
+            [
+                np.mean(np.fmax(cp[k] * (mc_path[1] - strike_std[k]), 0))
+                for k in range(len(strike_std))
+            ]
+        )
         if scalar_output:
             price = price[0]
-        return df*price
+        return df * price
diff --git a/pyfeng/opt_abc.py b/pyfeng/opt_abc.py
index b60f745..88552c7 100644
--- a/pyfeng/opt_abc.py
+++ b/pyfeng/opt_abc.py
@@ -10,7 +10,7 @@ class OptABC(abc.ABC):
 
     IMPVOL_TOL = 1e-10
     IMPVOL_MAXVOL = 99.99
-    
+
     def __init__(self, sigma, intr=0.0, divr=0.0, is_fwd=False):
         """
         Args:
@@ -28,7 +28,12 @@ def params_kw(self):
         """
         Model parameters in dictionary
         """
-        params = {"sigma": self.sigma, "intr": self.intr, "divr": self.divr, "is_fwd": self.is_fwd}
+        params = {
+            "sigma": self.sigma,
+            "intr": self.intr,
+            "divr": self.divr,
+            "is_fwd": self.is_fwd,
+        }
         return params
 
     def forward(self, spot, texp):
@@ -45,7 +50,7 @@ def forward(self, spot, texp):
         if self.is_fwd:
             return np.array(spot)
         else:
-            return np.array(spot) * np.exp((self.intr - self.divr)*np.array(texp))
+            return np.array(spot) * np.exp((self.intr - self.divr) * np.array(texp))
 
     def _fwd_factor(self, spot, texp):
         """
@@ -171,7 +176,7 @@ def _delta_shock(self, strike=100, spot=100, texp=1, cp=1):
         Returns:
             shock size
         """
-        return spot*0.001  # 10 bps of spot price
+        return spot * 0.001  # 10 bps of spot price
 
     def delta_numeric(self, strike, spot, texp, cp=1):
         """
@@ -187,7 +192,10 @@ def delta_numeric(self, strike, spot, texp, cp=1):
             delta value
         """
         h = self._delta_shock(strike, spot, texp, cp)
-        delta = (self.price(strike, spot+h, texp, cp)-self.price(strike, spot-h, texp, cp))/(2*h)
+        delta = (
+            self.price(strike, spot + h, texp, cp)
+            - self.price(strike, spot - h, texp, cp)
+        ) / (2 * h)
         return delta
 
     def gamma_numeric(self, strike, spot, texp, cp=1):
@@ -204,8 +212,11 @@ def gamma_numeric(self, strike, spot, texp, cp=1):
             Delta with numerical derivative
         """
         h = self._delta_shock(strike, spot, texp, cp)
-        gamma = (self.price(strike, spot+h, texp, cp) - 2*self.price(strike, spot, texp, cp)
-                 + self.price(strike, spot-h, texp, cp))/(h*h)
+        gamma = (
+            self.price(strike, spot + h, texp, cp)
+            - 2 * self.price(strike, spot, texp, cp)
+            + self.price(strike, spot - h, texp, cp)
+        ) / (h * h)
         return gamma
 
     def _vega_shock(self, strike=100, spot=100, texp=1, cp=1):
@@ -240,10 +251,10 @@ def vega_numeric(self, strike, spot, texp, cp=1):
         model = copy.copy(self)
         model.sigma += h
         p_up = model.price(strike, spot, texp, cp)
-        model.sigma -= 2*h
+        model.sigma -= 2 * h
         p_dn = model.price(strike, spot, texp, cp)
 
-        vega = (p_up - p_dn)/(2*h)
+        vega = (p_up - p_dn) / (2 * h)
         return vega
 
     def volga_numeric(self, strike, spot, texp, cp=1):
@@ -265,10 +276,10 @@ def volga_numeric(self, strike, spot, texp, cp=1):
         p_0 = model.price(strike, spot, texp, cp)
         model.sigma += h
         p_up = model.price(strike, spot, texp, cp)
-        model.sigma -= 2*h
+        model.sigma -= 2 * h
         p_dn = model.price(strike, spot, texp, cp)
 
-        volga = (p_up+p_dn-2*p_0)/(h*h)
+        volga = (p_up + p_dn - 2 * p_0) / (h * h)
         return volga
 
     def vanna_numeric(self, strike, spot, texp, cp=1):
@@ -285,10 +296,10 @@ def vanna_numeric(self, strike, spot, texp, cp=1):
             vanna value
         """
         h = self._delta_shock(strike, spot, texp, cp)
-        vega_up = self.vega_numeric(strike, spot+h, texp, cp)
-        vega_dn = self.vega_numeric(strike, spot-h, texp, cp)
+        vega_up = self.vega_numeric(strike, spot + h, texp, cp)
+        vega_dn = self.vega_numeric(strike, spot - h, texp, cp)
 
-        vanna = (vega_up - vega_dn)/(2*h)
+        vanna = (vega_up - vega_dn) / (2 * h)
         return vanna
 
     def _theta_shock(self, strike=100, spot=100, texp=1, cp=1):
@@ -304,7 +315,7 @@ def _theta_shock(self, strike=100, spot=100, texp=1, cp=1):
         Returns:
             shock size
         """
-        return np.minimum(1/365.25, texp)  # one day
+        return np.minimum(1 / 365.25, texp)  # one day
 
     def theta_numeric(self, strike, spot, texp, cp=1):
         """
@@ -320,7 +331,9 @@ def theta_numeric(self, strike, spot, texp, cp=1):
             theta value
         """
         dt = self._delta_shock(strike, spot, texp, cp)
-        theta = self.price(strike, spot, texp-dt, cp)-self.price(strike, spot, texp, cp)
+        theta = self.price(strike, spot, texp - dt, cp) - self.price(
+            strike, spot, texp, cp
+        )
         theta /= dt
         return theta
 
@@ -337,12 +350,12 @@ def pdf_numeric(self, strike, spot, texp, cp=-1, h=0.001):
         Returns:
             probability densitiy
         """
-        fwd = spot*(1.0 if self.is_fwd else np.exp(texp*(self.intr - self.divr)))
+        fwd = spot * (1.0 if self.is_fwd else np.exp(texp * (self.intr - self.divr)))
         kk = strike / fwd
         kk_arr = np.array([kk - h, kk, kk + h]).flatten()
         price = self.price(kk_arr, 1, texp, cp=cp)
         price = price.reshape(3, -1)
-        pdf = (price[2] + price[0] - 2.0*price[1])/(h*h)
+        pdf = (price[2] + price[0] - 2.0 * price[1]) / (h * h)
         return pdf
 
     # create aliases
@@ -383,9 +396,11 @@ def price_formula(strike, spot, sigma, texp, cp=1, *args, **kwargs):
 
     def price(self, strike, spot, texp, cp=1):
         if self.THROW_NEGATIVE_TEXP:
-            assert(~np.any(texp < 0))
+            assert ~np.any(texp < 0)
 
-        return self.price_formula(strike=strike, spot=spot, texp=texp, cp=cp, **self.params_kw())
+        return self.price_formula(
+            strike=strike, spot=spot, texp=texp, cp=cp, **self.params_kw()
+        )
 
     @abc.abstractmethod
     def delta(self, strike, spot, texp, cp=1):
@@ -494,11 +509,13 @@ def __init__(self, sigma, cor=None, intr=0.0, divr=0.0, is_fwd=False):
 
         if self.n_asset == 1:
             if cor is not None:
-                print(f'Ignoring cor={cor} for a single asset')
+                print(f"Ignoring cor={cor} for a single asset")
             self.rho = None
             self.cor_m = np.array([[1.0]])
         elif np.isscalar(cor):
-            self.cor_m = cor * np.ones((self.n_asset, self.n_asset)) + (1 - cor) * np.eye(self.n_asset)
+            self.cor_m = cor * np.ones((self.n_asset, self.n_asset)) + (
+                1 - cor
+            ) * np.eye(self.n_asset)
             self.rho = cor
         else:
             assert cor.shape == (self.n_asset, self.n_asset)
diff --git a/pyfeng/opt_smile_abc.py b/pyfeng/opt_smile_abc.py
index 4ccefc0..895c349 100644
--- a/pyfeng/opt_smile_abc.py
+++ b/pyfeng/opt_smile_abc.py
@@ -10,18 +10,19 @@ class OptSmileABC(opt.OptABC, abc.ABC):
     """
     Abstract class to model with volatility smile
     """
-    def _m_smile(self, model='bsm', is_fwd=None):
+
+    def _m_smile(self, model="bsm", is_fwd=None):
         if is_fwd is None:
             is_fwd = self.is_fwd
-        if model.lower() == 'bsm':
+        if model.lower() == "bsm":
             base_model = bsm.Bsm(None, intr=self.intr, divr=self.divr, is_fwd=is_fwd)
-        elif model.lower() == 'norm':
+        elif model.lower() == "norm":
             base_model = norm.Norm(None, intr=self.intr, divr=self.divr, is_fwd=is_fwd)
         else:
             base_model = None
         return base_model
 
-    def vol_smile(self, strike, spot, texp, cp=1, model='bsm'):
+    def vol_smile(self, strike, spot, texp, cp=1, model="bsm"):
         """
         Equivalent volatility smile for a given model
 
@@ -94,7 +95,7 @@ def vol_from_mass_zero(self, strike, spot, texp, mass=None):
         leading = tmp / np.sqrt(texp)
 
         qq = spst.norm.ppf(mass)
-        vol = 1 + (qq + 0.5*((2 + qq**2) + qq/tmp)/tmp)/tmp
+        vol = 1 + (qq + 0.5 * ((2 + qq ** 2) + qq / tmp) / tmp) / tmp
         vol *= leading
         return vol
 
diff --git a/pyfeng/ousv.py b/pyfeng/ousv.py
index c442f49..0bcd7ca 100644
--- a/pyfeng/ousv.py
+++ b/pyfeng/ousv.py
@@ -35,10 +35,17 @@ def D_B_C(self, s1, s2, s3, texp):
         s2g3 = vov ** 2 * gamma1 ** 3
 
         D = (mr - gamma1 * sincos / cossin) / vov ** 2
-        B = ((ktg3 + gamma3 * sincos) / cossin - mr * theta * gamma1) / (vov ** 2 * gamma1)
-        C = -0.5 * np.log(cossin) + 0.5 * mr * texp + ((mr * theta * gamma1) ** 2 - gamma3 ** 2) / (2 * s2g3) * (
-                sinh / cossin - gamma1 * texp) \
+        B = ((ktg3 + gamma3 * sincos) / cossin - mr * theta * gamma1) / (
+            vov ** 2 * gamma1
+        )
+        C = (
+            -0.5 * np.log(cossin)
+            + 0.5 * mr * texp
+            + ((mr * theta * gamma1) ** 2 - gamma3 ** 2)
+            / (2 * s2g3)
+            * (sinh / cossin - gamma1 * texp)
             + ktg3 * gamma3 / s2g3 * ((cosh - 1) / cossin)
+        )
 
         return D, B, C
 
@@ -50,7 +57,9 @@ def f_1(self, phi, fwd, texp):
         s2 = tmp * mr * theta * rho / vov
         s3 = 0.5 * tmp * rho / vov
 
-        res = 1j * phi * np.log(fwd) - 0.5 * rho * (1 + 1j * phi) * (self.sigma ** 2 / vov + vov * texp)
+        res = 1j * phi * np.log(fwd) - 0.5 * rho * (1 + 1j * phi) * (
+            self.sigma ** 2 / vov + vov * texp
+        )
         D, B, C = self.D_B_C(s1, s2, s3, texp)
         res += 0.5 * D * self.sigma ** 2 + B * self.sigma + C
         return np.exp(res)
@@ -63,7 +72,9 @@ def f_2(self, phi, fwd, texp):
         s2 = 1j * phi * mr * theta * rho / vov
         s3 = 0.5 * 1j * phi * rho / vov
 
-        res = 1j * phi * np.log(fwd) - 0.5 * 1j * phi * rho * (self.sigma ** 2 / vov + vov * texp)
+        res = 1j * phi * np.log(fwd) - 0.5 * 1j * phi * rho * (
+            self.sigma ** 2 / vov + vov * texp
+        )
         D, B, C = self.D_B_C(s1, s2, s3, texp)
         res += 0.5 * D * self.sigma ** 2 + B * self.sigma + C
         return np.exp(res)
@@ -75,12 +86,20 @@ def price(self, strike, spot, texp, cp=1):
         log_k = np.log(strike)
         J, h = 100001, 0.001
         phi = (np.arange(J)[:, None] + 1) * h  # shape=(J,1)
-        ff1 = 0.5 + 1 / np.pi * scint.simps((self.f_1(phi, fwd, texp) * np.exp(-1j * phi * log_k) / (1j * phi)).real,
-                                            dx=h, axis=0)
-        ff2 = 0.5 + 1 / np.pi * scint.simps((self.f_2(phi, fwd, texp) * np.exp(-1j * phi * log_k) / (1j * phi)).real,
-                                            dx=h, axis=0)
-
-        price = np.where(cp > 0, fwd * ff1 - strike * ff2, strike * (1 - ff2) - fwd * (1 - ff1))
+        ff1 = 0.5 + 1 / np.pi * scint.simps(
+            (self.f_1(phi, fwd, texp) * np.exp(-1j * phi * log_k) / (1j * phi)).real,
+            dx=h,
+            axis=0,
+        )
+        ff2 = 0.5 + 1 / np.pi * scint.simps(
+            (self.f_2(phi, fwd, texp) * np.exp(-1j * phi * log_k) / (1j * phi)).real,
+            dx=h,
+            axis=0,
+        )
+
+        price = np.where(
+            cp > 0, fwd * ff1 - strike * ff2, strike * (1 - ff2) - fwd * (1 - ff1)
+        )
         if len(price) == 1:
             price = price[0]
 
@@ -89,9 +108,9 @@ def price(self, strike, spot, texp, cp=1):
 
 class OusvCondMC(sv.SvABC, sv.CondMcBsmABC):
     """
-        OUSV model with conditional Monte-Carlo simulation
-        The SDE of SV is: dsigma_t = mr (theta - sigma_t) dt + vov dB_T
-        """
+    OUSV model with conditional Monte-Carlo simulation
+    The SDE of SV is: dsigma_t = mr (theta - sigma_t) dt + vov dB_T
+    """
 
     def _bm_incr(self, tobs, cum=False, n_path=None):
         """
@@ -116,7 +135,9 @@ def _bm_incr(self, tobs, cum=False, n_path=None):
         if self.antithetic:
             # generate random number in the order of path, time, asset and transposed
             # in this way, the same paths are generated when increasing n_path
-            bm_incr = self.rng.normal(size=(int(n_path / 2), n_dt)).T * np.sqrt(bm_var[:, None])
+            bm_incr = self.rng.normal(size=(int(n_path / 2), n_dt)).T * np.sqrt(
+                bm_var[:, None]
+            )
             bm_incr = np.stack([bm_incr, -bm_incr], axis=-1).reshape((-1, n_path))
         else:
             # bm_incr = np.random.randn(n_path, n_dt).T * np.sqrt(bm_var[:, None])
@@ -139,35 +160,44 @@ def vol_paths(self, tobs):
         Returns: volatility path (time, path) including the value at t=0
         """
         bm_path = self._bm_incr(tobs, cum=True)  # B_s (0 <= s <= 1)
-        sigma_t = self.theta + (self.sigma - self.theta + self.vov / np.sqrt(2 * self.mr) * bm_path) * np.exp(
-            -self.mr * tobs[:, None])
-        sigma_t = np.insert(sigma_t, 0, np.array([self.sigma] * sigma_t.shape[1]), axis=0)
+        sigma_t = self.theta + (
+            self.sigma - self.theta + self.vov / np.sqrt(2 * self.mr) * bm_path
+        ) * np.exp(-self.mr * tobs[:, None])
+        sigma_t = np.insert(
+            sigma_t, 0, np.array([self.sigma] * sigma_t.shape[1]), axis=0
+        )
         return sigma_t
 
     def cond_fwd_vol(self, texp):
         """
-            Returns new forward and volatility conditional on volatility path (e.g., sigma_T, integrated variance)
-            The forward and volatility are standardized in the sense that F_0 = 1 and sigma_0 = 1
-            Therefore, they should be scaled by the original F_0 and sigma_0 values
+        Returns new forward and volatility conditional on volatility path (e.g., sigma_T, integrated variance)
+        The forward and volatility are standardized in the sense that F_0 = 1 and sigma_0 = 1
+        Therefore, they should be scaled by the original F_0 and sigma_0 values
 
-            Args:
-                theta: the long term average
-                mr: coefficient of dt
-                texp: time-to-expiry
+        Args:
+            theta: the long term average
+            mr: coefficient of dt
+            texp: time-to-expiry
 
-            Returns: (forward, volatility)
+        Returns: (forward, volatility)
         """
         rhoc = np.sqrt(1.0 - self.rho ** 2)
         tobs = self.tobs(texp)
         n_dt = len(tobs)
         sigma_paths = self.vol_paths(tobs)
         sigma_final = sigma_paths[-1, :]
-        int_sigma = scint.simps(sigma_paths, dx=texp/n_dt, axis=0)
-        int_var = scint.simps(sigma_paths ** 2, dx=texp/n_dt, axis=0)
-
-        fwd_cond = np.exp(self.rho * ((sigma_final ** 2 - self.sigma ** 2) / (2 * self.vov) - self.vov * 0.5 * texp -
-                                      self.mr * self.theta / self.vov * int_sigma +
-                                      (self.mr / self.vov - self.rho * 0.5) * int_var))  # scaled by initial value
+        int_sigma = scint.simps(sigma_paths, dx=texp / n_dt, axis=0)
+        int_var = scint.simps(sigma_paths ** 2, dx=texp / n_dt, axis=0)
+
+        fwd_cond = np.exp(
+            self.rho
+            * (
+                (sigma_final ** 2 - self.sigma ** 2) / (2 * self.vov)
+                - self.vov * 0.5 * texp
+                - self.mr * self.theta / self.vov * int_sigma
+                + (self.mr / self.vov - self.rho * 0.5) * int_var
+            )
+        )  # scaled by initial value
 
         vol_cond = rhoc * np.sqrt(int_var / texp)
 
@@ -195,7 +225,7 @@ def price(self, strike, spot, texp, cp=1):
             base_model = self.base_model(vol_cond)
             price_grid = base_model.price(kk[:, None], fwd_cond, texp=t, cp=cp)
 
-            #np.set_printoptions(suppress=True, precision=6)
+            # np.set_printoptions(suppress=True, precision=6)
             price.append(spot * np.mean(price_grid, axis=1))  # in cond_fwd_vol, S_0 = 1
 
         return np.array(price).T
diff --git a/pyfeng/sabr.py b/pyfeng/sabr.py
index f360ea0..b24fa0f 100644
--- a/pyfeng/sabr.py
+++ b/pyfeng/sabr.py
@@ -22,7 +22,9 @@ class SabrABC(smile.OptSmileABC, abc.ABC):
     vov, beta, rho = 0.0, 1.0, 0.0
     _base_beta = None
 
-    def __init__(self, sigma, vov=0.0, rho=0.0, beta=1.0, intr=0.0, divr=0.0, is_fwd=False):
+    def __init__(
+        self, sigma, vov=0.0, rho=0.0, beta=1.0, intr=0.0, divr=0.0, is_fwd=False
+    ):
         """
         Args:
             sigma: model volatility at t=0
@@ -45,7 +47,7 @@ def params_kw(self):
 
     def _variables(self, fwd, texp):
         betac = 1.0 - self.beta
-        alpha = self.sigma / np.power(fwd, betac) #if self.beta > 0.0 else self.sigma
+        alpha = self.sigma / np.power(fwd, betac)  # if self.beta > 0.0 else self.sigma
         rho2 = self.rho * self.rho
         rhoc = np.sqrt(1.0 - rho2)
         vovn = self.vov * np.sqrt(np.maximum(texp, 1e-64))
@@ -69,11 +71,19 @@ def _m_base(self, vol, is_fwd=None):
         elif np.isclose(base_beta, 0):
             return norm.Norm(vol, intr=self.intr, divr=self.divr, is_fwd=is_fwd)
         else:
-            return cev.Cev(vol, beta=base_beta, intr=self.intr, divr=self.divr, is_fwd=is_fwd)
+            return cev.Cev(
+                vol, beta=base_beta, intr=self.intr, divr=self.divr, is_fwd=is_fwd
+            )
 
     def vol_smile(self, strike, spot, texp, cp=1, model=None):
         if model is None:
-            model = 'norm' if np.isclose(self.beta, 0) else 'bsm' if self.beta > 0.0 else None
+            model = (
+                "norm"
+                if np.isclose(self.beta, 0)
+                else "bsm"
+                if self.beta > 0.0
+                else None
+            )
 
         return super().vol_smile(strike, spot, texp, cp=1, model=model)
 
@@ -99,30 +109,32 @@ def init_benchmark(cls, set_no=None):
             - von Sydow, L., ..., Haentjens, T., & Waldén, J. (2018). BENCHOP - SLV: The BENCHmarking project in Option Pricing – Stochastic and Local Volatility problems. International Journal of Computer Mathematics, 1–14. https://doi.org/10.1080/00207160.2018.1544368
         """
         this_dir, _ = os.path.split(__file__)
-        file = os.path.join(this_dir, 'data/sabr_benchmark.xlsx')
-        df_param = pd.read_excel(file, sheet_name='Param', index_col='Sheet')
+        file = os.path.join(this_dir, "data/sabr_benchmark.xlsx")
+        df_param = pd.read_excel(file, sheet_name="Param", index_col="Sheet")
 
         if set_no is None:
             return df_param
         else:
             df_val = pd.read_excel(file, sheet_name=str(set_no))
             param = df_param.loc[set_no]
-            args_model = {k: param[k] for k in ('sigma', 'vov', 'rho', 'beta')}
-            args_pricing = {k: param[k] for k in ('texp', 'spot')}
+            args_model = {k: param[k] for k in ("sigma", "vov", "rho", "beta")}
+            args_pricing = {k: param[k] for k in ("texp", "spot")}
 
-            assert(df_val.columns[0] == 'k' or df_val.columns[0] == 'K')
-            args_pricing['strike'] = df_val.values[:, 0]
-            if df_val.columns[0] == 'k':
-                args_pricing['strike'] *= param['spot']
+            assert df_val.columns[0] == "k" or df_val.columns[0] == "K"
+            args_pricing["strike"] = df_val.values[:, 0]
+            if df_val.columns[0] == "k":
+                args_pricing["strike"] *= param["spot"]
 
-            val = df_val[param['col_name']].values
-            is_iv = param['col_name'].startswith('IV')
+            val = df_val[param["col_name"]].values
+            is_iv = param["col_name"].startswith("IV")
 
             m = cls(**args_model)
 
             param_dict = {
-                'args_pricing': args_pricing,
-                'ref': param['Reference'], 'val': val, 'is_iv': is_iv
+                "args_pricing": args_pricing,
+                "ref": param["Reference"],
+                "val": val,
+                "is_iv": is_iv,
             }
 
             return m, df_val, param_dict
@@ -132,7 +144,8 @@ class SabrVolApproxABC(SabrABC):
     """
     Abstract class for SABR models with volatility approximation (asymptotic expansion)
     """
-    approx_order = 1  # order in texp: 0: leading order, 1: first order, -1: reserved for 2/(1+exp(-2*eps)) smoothing 
+
+    approx_order = 1  # order in texp: 0: leading order, 1: first order, -1: reserved for 2/(1+exp(-2*eps)) smoothing
 
     @staticmethod
     def _vv(zz, rho):
@@ -162,11 +175,13 @@ def _int_inv_locvol(strike, beta, fwd=1.0):
         # fall-back indices and values first
         ind_atm = np.fabs(kk - 1.0) < 1e-6
         val = 1 - beta / 2 * (kk - 1) * (1 - (1 + beta) / 3 * (kk - 1))
-        with np.errstate(divide='ignore', invalid='ignore'):
+        with np.errstate(divide="ignore", invalid="ignore"):
             if abs(betac) < 1e-4:
                 val = np.where(ind_atm, val, np.log(kk) / (kk - 1))
             else:
-                val = np.where(ind_atm, val, (np.power(kk, betac) - 1) / betac / (kk - 1))
+                val = np.where(
+                    ind_atm, val, (np.power(kk, betac) - 1) / betac / (kk - 1)
+                )
 
         return val
 
@@ -184,16 +199,22 @@ def _hh(zz, rho):
         """
         rho2 = rho * rho
         # initalization with expansion for for small |zz|
-        xx_zz = 1 - (zz/2)*(rho - zz*(rho2 - 1/3 - (5*rho2 - 3)/4*rho*zz))
+        xx_zz = 1 - (zz / 2) * (
+            rho - zz * (rho2 - 1 / 3 - (5 * rho2 - 3) / 4 * rho * zz)
+        )
 
         yy = SabrVolApproxABC._vv(zz, rho)
         eps = 1e-5
 
-        with np.errstate(divide='ignore', invalid='ignore'):  # suppress error for zz=0
+        with np.errstate(divide="ignore", invalid="ignore"):  # suppress error for zz=0
             # replace negative zz
-            xx_zz = np.where(zz > -eps, xx_zz, np.log((1 - rho) / (yy - (zz + rho))) / zz)
+            xx_zz = np.where(
+                zz > -eps, xx_zz, np.log((1 - rho) / (yy - (zz + rho))) / zz
+            )
             # replace positive zz
-            xx_zz = np.where(zz < eps, xx_zz, np.log((yy + (zz + rho)) / (1 + rho)) / zz)
+            xx_zz = np.where(
+                zz < eps, xx_zz, np.log((yy + (zz + rho)) / (1 + rho)) / zz
+            )
 
         return 1.0 / xx_zz
 
@@ -211,7 +232,7 @@ def vol_for_price(self, strike, spot, texp):
             equivalent volatility
         """
         pass
-    
+
     def price(self, strike, spot, texp, cp=1):
         vol = self.vol_for_price(strike, spot, texp)
         m_vol = self._m_base(vol)
@@ -235,11 +256,18 @@ def iv_func(_sigma):
 
     def vol_smile(self, strike, spot, texp, cp=1, model=None):
         if model is None:
-            model = 'norm' if np.isclose(self.beta, 0) else 'bsm' if self.beta > 0.0 else None
+            model = (
+                "norm"
+                if np.isclose(self.beta, 0)
+                else "bsm"
+                if self.beta > 0.0
+                else None
+            )
 
         vol_beta = self._base_beta or self.beta
-        if (model.lower() == 'bsm' and np.isclose(vol_beta, 1)) or \
-                (model.lower() == 'norm' and np.isclose(vol_beta, 0)):
+        if (model.lower() == "bsm" and np.isclose(vol_beta, 1)) or (
+            model.lower() == "norm" and np.isclose(vol_beta, 0)
+        ):
             vol = self.vol_for_price(strike, spot, texp)
         else:
             vol = super().vol_smile(strike, spot, texp, cp=cp, model=model)
@@ -263,6 +291,7 @@ class SabrHagan2002(SabrVolApproxABC):
         >>> m.price(np.arange(80, 121, 10), 100, 1.2)
         array([22.04862858, 14.56226187,  8.74170415,  4.72352155,  2.28891776])
     """
+
     _base_beta = 1.0  # should not be changed
 
     def vol_for_price(self, strike, spot, texp):
@@ -273,29 +302,33 @@ def vol_for_price(self, strike, spot, texp):
             return 0.0
         fwd, _, _ = self._fwd_factor(spot, texp)
         alpha, betac, rhoc, rho2, vovn = self._variables(spot, texp)
-        betac2 = betac**2
+        betac2 = betac ** 2
 
-        log_kk = np.log(fwd/strike)
-        log_kk2 = log_kk*log_kk
-        pow_kk = np.power(strike/fwd, betac/2)
+        log_kk = np.log(fwd / strike)
+        log_kk2 = log_kk * log_kk
+        pow_kk = np.power(strike / fwd, betac / 2)
 
-        pre1 = pow_kk*(1 + betac2/24*log_kk2*(1 + betac2/80*log_kk2))
+        pre1 = pow_kk * (1 + betac2 / 24 * log_kk2 * (1 + betac2 / 80 * log_kk2))
 
-        term02 = (2 - 3*rho2)/24 * self.vov**2
-        term11 = alpha*self.vov*self.rho*self.beta/4/pow_kk
-        term20 = (betac*alpha/pow_kk)**2 / 24
+        term02 = (2 - 3 * rho2) / 24 * self.vov ** 2
+        term11 = alpha * self.vov * self.rho * self.beta / 4 / pow_kk
+        term20 = (betac * alpha / pow_kk) ** 2 / 24
 
         if self.approx_order == 0:
             vol = 1.0
         else:
-            vol = 1.0 + texp*(term02 + term11 + term20)
+            vol = 1.0 + texp * (term02 + term11 + term20)
 
-        zz = pow_kk*log_kk*self.vov/np.maximum(alpha, np.finfo(float).eps)
-        hh = self._hh(-zz, self.rho)  # note we pass -zz becaues hh(zz) definition is different
-        vol *= alpha*hh/pre1   # bsm vol
+        zz = pow_kk * log_kk * self.vov / np.maximum(alpha, np.finfo(float).eps)
+        hh = self._hh(
+            -zz, self.rho
+        )  # note we pass -zz becaues hh(zz) definition is different
+        vol *= alpha * hh / pre1  # bsm vol
         return vol
 
-    def calibrate3(self, price_or_vol3, strike3, spot, texp, cp=1, setval=False, is_vol=True):
+    def calibrate3(
+        self, price_or_vol3, strike3, spot, texp, cp=1, setval=False, is_vol=True
+    ):
         """
         Given option prices or implied vols at 3 strikes, compute the sigma, vov, rho to fit the data using `scipy.optimize.root`.
         If prices are given (is_vol=False) convert the prices to vol first.
@@ -327,10 +360,18 @@ def iv_func(x):
             return err
 
         sol = spop.root(iv_func, np.array([np.log(vol3[1]), -1, 0.0]))
-        params = {"sigma": np.exp(sol.x[0]), "vov": np.exp(sol.x[1]), "rho": np.tanh(sol.x[2])}
+        params = {
+            "sigma": np.exp(sol.x[0]),
+            "vov": np.exp(sol.x[1]),
+            "rho": np.tanh(sol.x[2]),
+        }
 
         if setval:
-            self.sigma, self.vov, self.rho = params["sigma"], params["vov"], params["rho"]
+            self.sigma, self.vov, self.rho = (
+                params["sigma"],
+                params["vov"],
+                params["rho"],
+            )
 
         return params
 
@@ -348,10 +389,21 @@ class SabrNorm(SabrVolApproxABC):
         >>> m.price(np.arange(80, 121, 10), 100, 1.2)
         array([22.04862858, 14.56226187,  8.74170415,  4.72352155,  2.28891776])
     """
+
     _base_beta = 0  # should not be changed
     is_atmvol = False
 
-    def __init__(self, sigma, vov=0.0, rho=0.0, beta=None, intr=0.0, divr=0.0, is_fwd=False, is_atmvol=False):
+    def __init__(
+        self,
+        sigma,
+        vov=0.0,
+        rho=0.0,
+        beta=None,
+        intr=0.0,
+        divr=0.0,
+        is_fwd=False,
+        is_atmvol=False,
+    ):
         """
         Args:
             sigma: model volatility at t=0
@@ -365,7 +417,7 @@ def __init__(self, sigma, vov=0.0, rho=0.0, beta=None, intr=0.0, divr=0.0, is_fw
         """
         # Make sure beta = 0
         if beta is not None and not np.isclose(beta, 0.0):
-            print(f'Ignoring beta = {beta}...')
+            print(f"Ignoring beta = {beta}...")
         self.is_atmvol = is_atmvol
         super().__init__(sigma, vov, rho, beta=0, intr=intr, divr=divr, is_fwd=is_fwd)
 
@@ -375,9 +427,9 @@ def vol_for_price(self, strike, spot, texp):
         if self.approx_order == 0 or self.is_atmvol:
             vol = 1.0
         else:
-            vol = 1.0 + texp*(2 - 3*self.rho**2)/24 * self.vov**2
+            vol = 1.0 + texp * (2 - 3 * self.rho ** 2) / 24 * self.vov ** 2
 
-        zz = self.vov/self.sigma*(strike-fwd)
+        zz = self.vov / self.sigma * (strike - fwd)
         hh = self._hh(zz, self.rho)
         vol *= self.sigma * hh
         return vol
@@ -399,7 +451,17 @@ class SabrChoiWu2021H(SabrVolApproxABC, smile.MassZeroABC):
         array([22.04897988, 14.56240351,  8.74169054,  4.72340753,  2.28876105])
     """
 
-    def __init__(self, sigma, vov=0.0, rho=0.0, beta=1.0, intr=0.0, divr=0.0, is_fwd=False, vol_beta=None):
+    def __init__(
+        self,
+        sigma,
+        vov=0.0,
+        rho=0.0,
+        beta=1.0,
+        intr=0.0,
+        divr=0.0,
+        is_fwd=False,
+        vol_beta=None,
+    ):
         """
         Args:
             sigma: model volatility at t=0
@@ -426,46 +488,60 @@ def vol_for_price(self, strike, spot, texp):
 
         kk = strike / fwd  # standardized strike
 
-        vov_over_alpha_safe = self.vov/np.maximum(alpha, np.finfo(float).eps)
+        vov_over_alpha_safe = self.vov / np.maximum(alpha, np.finfo(float).eps)
         tmp = self._int_inv_locvol(kk, self.beta)
-        qq_ratio = 1.0 if self._base_beta is None else self._int_inv_locvol(kk, vol_beta) / tmp
+        qq_ratio = (
+            1.0 if self._base_beta is None else self._int_inv_locvol(kk, vol_beta) / tmp
+        )
 
-        qq = tmp * (kk-1.0)
+        qq = tmp * (kk - 1.0)
         zz = vov_over_alpha_safe * qq  # zeta = (vov/sigma0) q
         hh = self._hh(zz, self.rho)
         # term02: O(vov^2)
-        term02 = (2 - 3*rho2)/24 * self.vov**2
+        term02 = (2 - 3 * rho2) / 24 * self.vov ** 2
 
         # term11: O(alpha*vov)
         # C(k)-C(1)/(k-1). Notice that 1/beta comes from int_inv_locvol
-        term11 = self.rho*self.beta/4 * self.vov*alpha * self._int_inv_locvol(kk, betac)
+        term11 = (
+            self.rho
+            * self.beta
+            / 4
+            * self.vov
+            * alpha
+            * self._int_inv_locvol(kk, betac)
+        )
 
         # term20: O(alpha^2)
         if np.isclose(self.beta, vol_beta):
             term20 = 0.0
         else:
-            with np.errstate(divide='ignore', invalid='ignore'):
+            with np.errstate(divide="ignore", invalid="ignore"):
                 ## Override ATM (qq=0)
                 term20 = np.where(
                     np.fabs(qq) < 1e-6,
                     (np.square(betac) - np.square(vol_betac)) / 24 * np.square(alpha),
-                    (0.5*(self.beta-self._base_beta)*np.log(kk) - np.log(qq_ratio)) * np.square(alpha)/qq
+                    (
+                        0.5 * (self.beta - self._base_beta) * np.log(kk)
+                        - np.log(qq_ratio)
+                    )
+                    * np.square(alpha)
+                    / qq,
                 )
-        #else:
+        # else:
         #    raise ValueError('Cannot handle this vol_beta different from beta')
 
         ## Return
         if self.approx_order == 0:
             vol = 1.0
         else:
-            vol = 1.0 + texp*(term20 + term11 + term02)
+            vol = 1.0 + texp * (term20 + term11 + term02)
 
         pre_fac = alpha * np.power(fwd, vol_betac) * qq_ratio
-        vol *= pre_fac*hh
+        vol *= pre_fac * hh
         return vol
 
     def mass_zero(self, spot, texp, log=False):
-        assert(self._base_beta is None)
+        assert self._base_beta is None
         vol_cev = self.vol_for_price(0.0, spot, texp)
         cev_m = cev.Cev(sigma=vol_cev, beta=self.beta)
         mass = cev_m.mass_zero(spot, texp, log=log)
@@ -485,7 +561,7 @@ def mass_zero_t0(self, spot, texp):
         fwd = self.forward(spot, texp)
         alpha, betac, rhoc, rho2, _ = self._variables(fwd, texp)
         hh = self._hh(-self.vov / (alpha * betac), self.rho)
-        t0 = 0.5/(betac*alpha*hh)**2
+        t0 = 0.5 / (betac * alpha * hh) ** 2
         return t0
 
 
@@ -523,60 +599,94 @@ def vol_for_price(self, strike, spot, texp):
             kk = np.array([kk])
 
         ## Eq 32 (leading order)
-        vov_over_alpha_safe = self.vov/np.maximum(alpha, np.finfo(float).eps)
+        vov_over_alpha_safe = self.vov / np.maximum(alpha, np.finfo(float).eps)
 
         # qq_ratio = qq_vol_beta / qq_beta
         tmp = self._int_inv_locvol(kk, self.beta)
-        qq_ratio = 1.0 if self._base_beta is None else self._int_inv_locvol(kk, vol_beta) / tmp
+        qq_ratio = (
+            1.0 if self._base_beta is None else self._int_inv_locvol(kk, vol_beta) / tmp
+        )
 
         zz = vov_over_alpha_safe * tmp * (kk - 1.0)  # zeta = (vov/sigma0) q
         hh = self._hh(zz, self.rho)
         v_m = self._vv(zz, self.rho)
 
         if abs(betac) < 0.001:
-            gg_diff = 0.5*self.beta*self.rho/(1.0-rho2) * (v_m - 1.0 - self.rho*zz)
+            gg_diff = (
+                0.5 * self.beta * self.rho / (1.0 - rho2) * (v_m - 1.0 - self.rho * zz)
+            )
         else:
-            t1 = (v_m + self.rho + zz)/rhoc  # array
-            t2 = (1 + self.rho)/rhoc  # scalar
+            t1 = (v_m + self.rho + zz) / rhoc  # array
+            t2 = (1 + self.rho) / rhoc  # scalar
 
-            eta = self.vov/(betac*alpha)*np.power(kk, betac) * (rhoc/v_m)
-            eta2_m_1 = eta*eta - 1.0
+            eta = self.vov / (betac * alpha) * np.power(kk, betac) * (rhoc / v_m)
+            eta2_m_1 = eta * eta - 1.0
             sqrt_eta2_m_1 = np.sqrt(abs(eta2_m_1))
 
-            num_t1 = self.rho + (eta-rhoc)*t1
-            num_t2 = self.rho + (eta-rhoc)*t2
+            num_t1 = self.rho + (eta - rhoc) * t1
+            num_t2 = self.rho + (eta - rhoc) * t2
 
             gg1 = np.arctan(t1)
             gg2 = np.arctan(t2) * np.ones_like(t1)  # t2 is a scalar, so vectorize
 
             eps = 1e-6
 
-            ind = (eta2_m_1 > eps)
-            gg2[ind] += eta[ind] / sqrt_eta2_m_1[ind] * np.arctan2(sqrt_eta2_m_1[ind], num_t2[ind])
-            gg1[ind] += eta[ind] / sqrt_eta2_m_1[ind] * np.arctan2(sqrt_eta2_m_1[ind], num_t1[ind])
+            ind = eta2_m_1 > eps
+            gg2[ind] += (
+                eta[ind]
+                / sqrt_eta2_m_1[ind]
+                * np.arctan2(sqrt_eta2_m_1[ind], num_t2[ind])
+            )
+            gg1[ind] += (
+                eta[ind]
+                / sqrt_eta2_m_1[ind]
+                * np.arctan2(sqrt_eta2_m_1[ind], num_t1[ind])
+            )
 
-            ind = (eta2_m_1 < -eps)
-            gg1[ind] += eta[ind] / sqrt_eta2_m_1[ind] \
-                        * 0.5*np.log(abs( (num_t1[ind]+sqrt_eta2_m_1[ind])/(num_t1[ind]-sqrt_eta2_m_1[ind]) ))
+            ind = eta2_m_1 < -eps
+            gg1[ind] += (
+                eta[ind]
+                / sqrt_eta2_m_1[ind]
+                * 0.5
+                * np.log(
+                    abs(
+                        (num_t1[ind] + sqrt_eta2_m_1[ind])
+                        / (num_t1[ind] - sqrt_eta2_m_1[ind])
+                    )
+                )
+            )
 
             # as eta->0 (k->0) the numerator diverge. but the gg2 value should be 0.
-            ind = (eta2_m_1 < -eps)*(eta > eps)
-            gg2[ind] += eta[ind] / sqrt_eta2_m_1[ind] \
-                        * 0.5*np.log(abs( (num_t2[ind]+sqrt_eta2_m_1[ind])/(num_t2[ind]-sqrt_eta2_m_1[ind]) ))
+            ind = (eta2_m_1 < -eps) * (eta > eps)
+            gg2[ind] += (
+                eta[ind]
+                / sqrt_eta2_m_1[ind]
+                * 0.5
+                * np.log(
+                    abs(
+                        (num_t2[ind] + sqrt_eta2_m_1[ind])
+                        / (num_t2[ind] - sqrt_eta2_m_1[ind])
+                    )
+                )
+            )
             # when eta is very small, the term above is zero, so do nothing.
 
-            ind = (abs(eta2_m_1) <= eps)
-            gg1[ind] += eta[ind]/num_t1[ind] * (1 - eta2_m_1[ind]/num_t1[ind]**2/3)
-            gg2[ind] += eta[ind]/num_t2[ind] * (1 - eta2_m_1[ind]/num_t2[ind]**2/3)
+            ind = abs(eta2_m_1) <= eps
+            gg1[ind] += (
+                eta[ind] / num_t1[ind] * (1 - eta2_m_1[ind] / num_t1[ind] ** 2 / 3)
+            )
+            gg2[ind] += (
+                eta[ind] / num_t2[ind] * (1 - eta2_m_1[ind] / num_t2[ind] ** 2 / 3)
+            )
 
-            gg_diff = self.rho*self.beta/(rhoc*betac) * (gg2 - gg1)
+            gg_diff = self.rho * self.beta / (rhoc * betac) * (gg2 - gg1)
 
         zz2_safe = np.maximum(np.square(zz), np.finfo(float).eps)
 
         if np.isclose(self.beta, vol_beta):
             tmp = 0.0
         else:
-            tmp = 0.5*(self.beta-self._base_beta)*np.log(kk) - np.log(qq_ratio)
+            tmp = 0.5 * (self.beta - self._base_beta) * np.log(kk) - np.log(qq_ratio)
 
         """
         with np.errstate(divide='ignore', invalid='ignore'):
@@ -587,19 +697,25 @@ def vol_for_price(self, strike, spot, texp):
                 (0.5 * (self.beta - self._base_beta) * np.log(kk) - np.log(qq_ratio)) * np.square(alpha) / qq
             )
         """
-        order1 = np.square(self.vov*hh)/zz2_safe *\
-                 (tmp + 0.5*np.log(v_m/np.square(hh)) + gg_diff)
+        order1 = (
+            np.square(self.vov * hh)
+            / zz2_safe
+            * (tmp + 0.5 * np.log(v_m / np.square(hh)) + gg_diff)
+        )
 
         ## Override ATM (z=0)
-        ind_atm = (abs(zz) < 1e-6)
-        order1[ind_atm] = (np.square(betac)-np.square(vol_betac))/24*alpha*alpha \
-            + ((self.rho*self.beta/4)*alpha + (2 - 3*rho2)/24*self.vov)*self.vov  # RHS scalar
+        ind_atm = abs(zz) < 1e-6
+        order1[ind_atm] = (
+            np.square(betac) - np.square(vol_betac)
+        ) / 24 * alpha * alpha + (
+            (self.rho * self.beta / 4) * alpha + (2 - 3 * rho2) / 24 * self.vov
+        ) * self.vov  # RHS scalar
 
         # return value
         if self.approx_order == 0:
             vol = 1.0
         else:
-            vol = 1.0 + order1*texp
+            vol = 1.0 + order1 * texp
 
         # order0 = (z_base' / z_base) * hh
         vol *= alpha * np.power(fwd, vol_betac) * qq_ratio * hh
@@ -630,58 +746,104 @@ def vol_for_price(self, strike, spot, texp):
         kk = strike / fwd  # standardized strike
 
         kmx = np.log(kk)
-        kmx_2 = kmx*kmx
-        kmx_3 = kmx*kmx_2
+        kmx_2 = kmx * kmx
+        kmx_3 = kmx * kmx_2
 
-        t_2 = texp*texp if self.approx_order > 1 else 0.0
-        t_3 = texp*t_2 if self.approx_order > 2 else 0.0
+        t_2 = texp * texp if self.approx_order > 1 else 0.0
+        t_3 = texp * t_2 if self.approx_order > 2 else 0.0
 
-        betac_2 = betac*betac
-        betac_3 = betac*betac_2
+        betac_2 = betac * betac
+        betac_3 = betac * betac_2
 
         rho = self.rho
         vov = self.vov
-        vov_2 = vov*vov
-        vov_3 = vov*vov_2
+        vov_2 = vov * vov
+        vov_3 = vov * vov_2
 
         s0 = alpha
-        s0_2 = s0*s0
-        s0_3 = s0*s0_2
-        s0_4 = s0*s0_3
-        s0_5 = s0*s0_4
-
-        s10 = s0/2*(-betac)*kmx
-        s01 = vov/4 * ( 2*kmx*rho + texp*s0*(-vov+rho*s0) )
-
-        s20 = s0*betac_2 * (texp/24*s0_2 - t_2/96*s0_4 + kmx_2/12) * self.lv_factor
-
-        s11 = -betac*vov
-        s11 *= texp/4*rho*s0_2 - t_2/48*rho*s0_4 - texp/24*s0*(3*vov-5*rho*s0)*kmx
+        s0_2 = s0 * s0
+        s0_3 = s0 * s0_2
+        s0_4 = s0 * s0_3
+        s0_5 = s0 * s0_4
+
+        s10 = s0 / 2 * (-betac) * kmx
+        s01 = vov / 4 * (2 * kmx * rho + texp * s0 * (-vov + rho * s0))
+
+        s20 = (
+            s0
+            * betac_2
+            * (texp / 24 * s0_2 - t_2 / 96 * s0_4 + kmx_2 / 12)
+            * self.lv_factor
+        )
+
+        s11 = -betac * vov
+        s11 *= (
+            texp / 4 * rho * s0_2
+            - t_2 / 48 * rho * s0_4
+            - texp / 24 * s0 * (3 * vov - 5 * rho * s0) * kmx
+        )
 
         s02 = vov_2
-        s02 *= texp/24*(8-3*rho2)*s0 + t_2/96*s0*(5*vov_2 + 2*s0*((6*rho2-2)*s0 - 7*vov*rho)) \
-            - texp/24*rho*(vov-3*rho*s0)*kmx + (2-3*rho2)/(12*s0)*kmx_2
-
-        s30 = -betac_3*s0_3*kmx
-        s30 *= texp/16 - 5/192*t_2*s0_2
-
-        s21 = betac_2*vov
-        s21 *= t_2/96*s0_3*(7*rho*s0 - 3*vov) + t_3/384*s0_5*(5*vov-7*rho*s0) \
-            + 13/48*texp*rho*s0_2*kmx - 13/192*t_2*rho*s0_4*kmx \
-            - texp/48*s0*(vov-3*rho*s0)*kmx_2
-
-        s12 = -betac*vov_2
-        s12 *= t_2/48*rho*s0_2*(13*rho*s0-7*vov) + t_3/192*rho*s0_4*(5*vov-7*rho*s0) \
-            + texp/48*(3*rho2+8)*s0*kmx \
-            + t_2/192*s0*(5*vov_2 - 22*vov*rho*s0 + 4*(5*rho2-3)*s0_2)*kmx \
-            + texp/24*rho2*s0*kmx_2 + (3*rho2-2)/(24*s0)*kmx_3
+        s02 *= (
+            texp / 24 * (8 - 3 * rho2) * s0
+            + t_2
+            / 96
+            * s0
+            * (5 * vov_2 + 2 * s0 * ((6 * rho2 - 2) * s0 - 7 * vov * rho))
+            - texp / 24 * rho * (vov - 3 * rho * s0) * kmx
+            + (2 - 3 * rho2) / (12 * s0) * kmx_2
+        )
+
+        s30 = -betac_3 * s0_3 * kmx
+        s30 *= texp / 16 - 5 / 192 * t_2 * s0_2
+
+        s21 = betac_2 * vov
+        s21 *= (
+            t_2 / 96 * s0_3 * (7 * rho * s0 - 3 * vov)
+            + t_3 / 384 * s0_5 * (5 * vov - 7 * rho * s0)
+            + 13 / 48 * texp * rho * s0_2 * kmx
+            - 13 / 192 * t_2 * rho * s0_4 * kmx
+            - texp / 48 * s0 * (vov - 3 * rho * s0) * kmx_2
+        )
+
+        s12 = -betac * vov_2
+        s12 *= (
+            t_2 / 48 * rho * s0_2 * (13 * rho * s0 - 7 * vov)
+            + t_3 / 192 * rho * s0_4 * (5 * vov - 7 * rho * s0)
+            + texp / 48 * (3 * rho2 + 8) * s0 * kmx
+            + t_2
+            / 192
+            * s0
+            * (5 * vov_2 - 22 * vov * rho * s0 + 4 * (5 * rho2 - 3) * s0_2)
+            * kmx
+            + texp / 24 * rho2 * s0 * kmx_2
+            + (3 * rho2 - 2) / (24 * s0) * kmx_3
+        )
 
         s03 = vov_3
-        s03 *= t_2/96*s0*(3*vov*(rho2-4) + rho*(26-9*rho2)*s0) \
-            + t_3/384*s0*(s0*(19*vov_2*rho + 2*s0*(vov*(8-21*rho2) + rho*(15*rho2-11)*s0)) - 3*vov_3) \
-            + texp/48*rho*(3*rho2-2)*kmx + t_2/192*rho*(vov_2+6*s0*(vov*rho+(1-2*rho2)*s0))*kmx \
-            - texp/16*rho*(rho2-1)*kmx_2 + rho*(6*rho2-5)/(24*s0_2)*kmx_3
-
-        vol = s0 + (s10+s01) + (s20+s11+s02) + (s30+s21+s12+s03)
+        s03 *= (
+            t_2 / 96 * s0 * (3 * vov * (rho2 - 4) + rho * (26 - 9 * rho2) * s0)
+            + t_3
+            / 384
+            * s0
+            * (
+                s0
+                * (
+                    19 * vov_2 * rho
+                    + 2 * s0 * (vov * (8 - 21 * rho2) + rho * (15 * rho2 - 11) * s0)
+                )
+                - 3 * vov_3
+            )
+            + texp / 48 * rho * (3 * rho2 - 2) * kmx
+            + t_2
+            / 192
+            * rho
+            * (vov_2 + 6 * s0 * (vov * rho + (1 - 2 * rho2) * s0))
+            * kmx
+            - texp / 16 * rho * (rho2 - 1) * kmx_2
+            + rho * (6 * rho2 - 5) / (24 * s0_2) * kmx_3
+        )
+
+        vol = s0 + (s10 + s01) + (s20 + s11 + s02) + (s30 + s21 + s12 + s03)
 
         return vol
diff --git a/pyfeng/sabr_int.py b/pyfeng/sabr_int.py
index 0af4518..f68f191 100644
--- a/pyfeng/sabr_int.py
+++ b/pyfeng/sabr_int.py
@@ -30,7 +30,17 @@ class SabrUncorrChoiWu2021(sabr.SabrABC, smile.MassZeroABC):
     _base_beta = None
     n_quad = 9
 
-    def __init__(self, sigma, vov=0.0, rho=0.0, beta=1.0, intr=0.0, divr=0.0, is_fwd=False, n_quad=9):
+    def __init__(
+        self,
+        sigma,
+        vov=0.0,
+        rho=0.0,
+        beta=1.0,
+        intr=0.0,
+        divr=0.0,
+        is_fwd=False,
+        n_quad=9,
+    ):
         """
         Args:
             sigma: model volatility at t=0
@@ -59,15 +69,15 @@ def int_var_lndist(vovn):
             (m1, sig)
             True distribution should be multiplied by sigma^2*t
         """
-        v2 = vovn**2
+        v2 = vovn ** 2
         w = np.exp(v2)
-        m1 = np.where(v2 > 1e-6, (w-1)/v2, 1+v2/2*(1+v2/3))
-        m2m1ratio = (5 + w*(4 + w*(3 + w*(2 + w))))/15
-        sig = np.sqrt(np.where(v2 > 1e-8, np.log(m2m1ratio), 4/3*v2))
+        m1 = np.where(v2 > 1e-6, (w - 1) / v2, 1 + v2 / 2 * (1 + v2 / 3))
+        m2m1ratio = (5 + w * (4 + w * (3 + w * (2 + w)))) / 15
+        sig = np.sqrt(np.where(v2 > 1e-8, np.log(m2m1ratio), 4 / 3 * v2))
         return m1, sig
 
     def price(self, strike, spot, texp, cp=1):
-        assert (self._base_beta is None)
+        assert self._base_beta is None
         m1, fac = self.int_var_lndist(self.vov * np.sqrt(texp))
 
         zz, ww = spsp.roots_hermitenorm(self.n_quad)
@@ -90,7 +100,11 @@ def mass_zero(self, spot, texp, log=False, mu=0):
         vol = self.sigma * np.sqrt(m1) * np.exp(0.5 * (zz - 0.5 * fac) * fac)
 
         if log:
-            log_mass = np.log(ww) + log_rn_deriv + self._m_base(vol).mass_zero(spot, texp, log=True)
+            log_mass = (
+                np.log(ww)
+                + log_rn_deriv
+                + self._m_base(vol).mass_zero(spot, texp, log=True)
+            )
             log_max = np.amax(log_mass)
             log_mass -= log_max
             log_mass = log_max + np.log(np.sum(np.exp(log_mass)))
diff --git a/pyfeng/sabr_mc.py b/pyfeng/sabr_mc.py
index f393d08..e381c5f 100644
--- a/pyfeng/sabr_mc.py
+++ b/pyfeng/sabr_mc.py
@@ -27,26 +27,32 @@ def vol_paths(self, tobs, mu=0):
         vov_std = self.vov * np.sqrt(texp)
 
         log_sig_s = self._bm_incr(tobs01, cum=True)  # B_s (0 <= s <= 1)
-        log_rn_deriv = 0.0 if mu == 0 else -mu*(log_sig_s[-1, :] + 0.5*mu)
+        log_rn_deriv = 0.0 if mu == 0 else -mu * (log_sig_s[-1, :] + 0.5 * mu)
 
-        log_sig_s = vov_std*(log_sig_s + (mu - 0.5*vov_std) * tobs01[:, None])
+        log_sig_s = vov_std * (log_sig_s + (mu - 0.5 * vov_std) * tobs01[:, None])
         log_sig_s = np.insert(log_sig_s, 0, np.zeros(log_sig_s.shape[1]), axis=0)
         return np.exp(log_sig_s), log_rn_deriv
 
     def cond_fwd_vol(self, texp, mu=0):
-        rhoc = np.sqrt(1.0 - self.rho**2)
-        rho_sigma = self.rho*self.sigma
+        rhoc = np.sqrt(1.0 - self.rho ** 2)
+        rho_sigma = self.rho * self.sigma
 
         tobs = self.tobs(texp)
         n_dt = len(tobs)
         sigma_paths, log_rn_deriv = self.vol_paths(tobs, mu=mu)
         sigma_final = sigma_paths[-1, :]
-        int_var = scint.simps(sigma_paths**2, dx=1, axis=0) / n_dt
-        vol_cond = rhoc*np.sqrt(int_var)
+        int_var = scint.simps(sigma_paths ** 2, dx=1, axis=0) / n_dt
+        vol_cond = rhoc * np.sqrt(int_var)
         if np.isclose(self.beta, 0):
-            fwd_cond = rho_sigma/self.vov*(sigma_final - 1)
+            fwd_cond = rho_sigma / self.vov * (sigma_final - 1)
         else:
-            fwd_cond = np.exp(rho_sigma*(1.0/self.vov*(sigma_final - 1) - 0.5*rho_sigma*int_var*texp))
+            fwd_cond = np.exp(
+                rho_sigma
+                * (
+                    1.0 / self.vov * (sigma_final - 1)
+                    - 0.5 * rho_sigma * int_var * texp
+                )
+            )
 
         return fwd_cond, vol_cond, log_rn_deriv
 
@@ -54,16 +60,16 @@ def price(self, strike, spot, texp, cp=1):
         fwd = self.forward(spot, texp)
         fwd_cond, vol_cond, log_rn_deriv = self.cond_fwd_vol(texp)
         if np.isclose(self.beta, 0):
-            base_model = self._m_base(self.sigma*vol_cond, is_fwd=True)
+            base_model = self._m_base(self.sigma * vol_cond, is_fwd=True)
             price_grid = base_model.price(strike[:, None], fwd + fwd_cond, texp, cp=cp)
-            price = np.mean(price_grid*np.exp(log_rn_deriv), axis=1)
+            price = np.mean(price_grid * np.exp(log_rn_deriv), axis=1)
         else:
             alpha = self.sigma / np.power(spot, 1.0 - self.beta)
             kk = strike / fwd
 
-            base_model = self._m_base(alpha*vol_cond, is_fwd=True)
+            base_model = self._m_base(alpha * vol_cond, is_fwd=True)
             price_grid = base_model.price(kk[:, None], fwd_cond, texp, cp=cp)
-            price = fwd * np.mean(price_grid*np.exp(log_rn_deriv), axis=1)
+            price = fwd * np.mean(price_grid * np.exp(log_rn_deriv), axis=1)
 
         return price
 
@@ -71,15 +77,19 @@ def mass_zero(self, spot, texp, log=False, mu=0):
         assert 0 < self.beta < 1
         assert self.rho == 0
 
-        eta = self.vov*np.power(spot, 1.0-self.beta)/(self.sigma*(1.0-self.beta))
+        eta = (
+            self.vov
+            * np.power(spot, 1.0 - self.beta)
+            / (self.sigma * (1.0 - self.beta))
+        )
         vov_std = self.vov * np.sqrt(texp)
 
         if mu is None:
-            mu = 0.5*(vov_std + np.log(1+eta**2)/vov_std)
-            #print(f'mu = {mu}')
+            mu = 0.5 * (vov_std + np.log(1 + eta ** 2) / vov_std)
+            # print(f'mu = {mu}')
 
         fwd_cond, vol_cond, log_rn_deriv = self.cond_fwd_vol(texp, mu=mu)
-        base_model = cev.Cev(sigma=self.sigma*vol_cond, beta=self.beta)
+        base_model = cev.Cev(sigma=self.sigma * vol_cond, beta=self.beta)
         if log:
             log_mass_grid = base_model.mass_zero(spot, texp, log=True) + log_rn_deriv
             log_mass_max = np.amax(log_mass_grid)
@@ -87,6 +97,8 @@ def mass_zero(self, spot, texp, log=False, mu=0):
             log_mass = log_mass_max + np.log(np.mean(np.exp(log_mass_grid)))
             return log_mass
         else:
-            mass_grid = base_model.mass_zero(spot, texp, log=False)*np.exp(log_rn_deriv)
+            mass_grid = base_model.mass_zero(spot, texp, log=False) * np.exp(
+                log_rn_deriv
+            )
             mass = np.mean(mass_grid)
             return mass
diff --git a/pyfeng/sv32_mc.py b/pyfeng/sv32_mc.py
index fb83c1e..5662161 100644
--- a/pyfeng/sv32_mc.py
+++ b/pyfeng/sv32_mc.py
@@ -17,7 +17,7 @@
 
 
 class Sv32CondMcQE:
-    '''
+    """
     Conditional MC for 3/2 model based on QE discretization scheme by Andersen(2008)
 
     Underlying price is assumed to follow a geometric Brownian motion.
@@ -41,10 +41,10 @@ class Sv32CondMcQE:
                [22.93527518, 10.4758762 , 38.87971674],
                [22.9298084 , 10.47613694, 38.88556212],
                [23.12806844, 10.56484306, 39.16893668]])
-    '''
+    """
 
     def __init__(self, vov=1, kappa=0.5, rho=-0.9, theta=0.04):
-        '''
+        """
         Initiate a 3/2 model
 
         Args:
@@ -52,7 +52,7 @@ def __init__(self, vov=1, kappa=0.5, rho=-0.9, theta=0.04):
             kappa: speed of variance's mean-reversion, strictly positive
             rho: correlation between BMs of price and vol
             theta: long-term mean (equilibirum level) of the variance, strictly positive
-        '''
+        """
         self.vov = vov
         self.kappa = kappa
         self.rho = rho
@@ -62,8 +62,21 @@ def __init__(self, vov=1, kappa=0.5, rho=-0.9, theta=0.04):
         self.rx_results = None
         self.dis = 1e-3
 
-    def price(self, strike, spot, texp, sigma, delta, intr=0, divr=0, psi_c=1.5, path=10000, scheme='QE', seed=None):
-        '''
+    def price(
+        self,
+        strike,
+        spot,
+        texp,
+        sigma,
+        delta,
+        intr=0,
+        divr=0,
+        psi_c=1.5,
+        path=10000,
+        scheme="QE",
+        seed=None,
+    ):
+        """
         Conditional MC routine for 3/2 model
         Generate paths for vol only using QE discretization scheme.
         Compute integrated variance and get BSM prices vector for all strikes.
@@ -83,7 +96,7 @@ def price(self, strike, spot, texp, sigma, delta, intr=0, divr=0, psi_c=1.5, pat
 
         Return:
             BSM price vector for all strikes
-        '''
+        """
         self.sigma = sigma
         self.bsm_model = pf.Bsm(self.sigma, intr=intr, divr=divr)
         self.delta = delta
@@ -91,22 +104,25 @@ def price(self, strike, spot, texp, sigma, delta, intr=0, divr=0, psi_c=1.5, pat
         self.step = int(texp / self.delta)
 
         # xt = 1 / vt
-        xt = 1 / self.sigma**2 * np.ones([self.path, self.step + 1])
+        xt = 1 / self.sigma ** 2 * np.ones([self.path, self.step + 1])
         np.random.seed(seed)
 
         # equivalent kappa and theta for xt to follow a Heston model
         kappa_new = self.kappa * self.theta
         theta_new = (self.kappa + self.vov ** 2) / (self.kappa * self.theta)
-        vov_new = - self.vov
-        if scheme == 'QE':
+        vov_new = -self.vov
+        if scheme == "QE":
             u = np.random.uniform(size=(self.path, self.step))
 
             expo = np.exp(-kappa_new * self.delta)
             for i in range(self.step):
                 # compute m, s_square, psi given xt(i)
                 m = theta_new + (xt[:, i] - theta_new) * expo
-                s2 = xt[:, i] * (vov_new ** 2) * expo * (1 - expo) / kappa_new + theta_new * (vov_new ** 2) * \
-                     ((1 - expo) ** 2) / (2 * kappa_new)
+                s2 = xt[:, i] * (vov_new ** 2) * expo * (
+                    1 - expo
+                ) / kappa_new + theta_new * (vov_new ** 2) * ((1 - expo) ** 2) / (
+                    2 * kappa_new
+                )
                 psi = s2 / m ** 2
 
                 # compute xt(i+1) given psi
@@ -116,18 +132,20 @@ def price(self, strike, spot, texp, sigma, delta, intr=0, divr=0, psi_c=1.5, pat
                 b = np.sqrt(b2)
                 a = m[below] / (1 + b2)
                 z = spst.norm.ppf(u[below, i])
-                xt[below, i+1] = a * (b + z) ** 2
+                xt[below, i + 1] = a * (b + z) ** 2
 
                 above = np.where(psi > psi_c)[0]
                 p = (psi[above] - 1) / (psi[above] + 1)
                 beta = (1 - p) / m[above]
                 for k in range(len(above)):
                     if u[above[k], i] > p[k]:
-                        xt[above[k], i+1] = beta[k] ** -1 * np.log((1 - p[k]) / (1 - u[above[k], i]))
+                        xt[above[k], i + 1] = beta[k] ** -1 * np.log(
+                            (1 - p[k]) / (1 - u[above[k], i])
+                        )
                     else:
-                        xt[above[k], i+1] = 0
+                        xt[above[k], i + 1] = 0
 
-        elif scheme == 'TG':
+        elif scheme == "TG":
             if np.all(self.rx_results) == None:
                 self.psi_points, self.rx_results = self.prepare_rx()
 
@@ -135,8 +153,11 @@ def price(self, strike, spot, texp, sigma, delta, intr=0, divr=0, psi_c=1.5, pat
             for i in range(self.step):
                 # compute m, s_square, psi given vt(i)
                 m = theta_new + (xt[:, i] - theta_new) * expo
-                s2 = xt[:, i] * (vov_new ** 2) * expo * (1 - expo) / kappa_new + theta_new * (vov_new ** 2) * \
-                     ((1 - expo) ** 2) / (2 * kappa_new)
+                s2 = xt[:, i] * (vov_new ** 2) * expo * (
+                    1 - expo
+                ) / kappa_new + theta_new * (vov_new ** 2) * ((1 - expo) ** 2) / (
+                    2 * kappa_new
+                )
                 psi = s2 / m ** 2
 
                 rx = np.array([self.find_rx(j) for j in psi])
@@ -145,29 +166,42 @@ def price(self, strike, spot, texp, sigma, delta, intr=0, divr=0, psi_c=1.5, pat
                 mu_v = np.zeros_like(z)
                 sigma_v = np.zeros_like(z)
                 mu_v[:, i] = rx * m / (spst.norm.pdf(rx) + rx * spst.norm.cdf(rx))
-                sigma_v[:, i] = np.sqrt(s2) * psi ** (-0.5) / (spst.norm.pdf(rx) + rx * spst.norm.cdf(rx))
+                sigma_v[:, i] = (
+                    np.sqrt(s2)
+                    * psi ** (-0.5)
+                    / (spst.norm.pdf(rx) + rx * spst.norm.cdf(rx))
+                )
 
-                xt[:, i+1] = np.fmax(mu_v[:, i] + sigma_v[:, i] * z[:, i], 0)
+                xt[:, i + 1] = np.fmax(mu_v[:, i] + sigma_v[:, i] * z[:, i], 0)
 
-        elif scheme == 'Euler':
+        elif scheme == "Euler":
             z = np.random.normal(size=(self.path, self.step))
             for i in range(self.step):
-                xt[:, i+1] = xt[:, i] + kappa_new * (theta_new - np.max(xt[:, i], 0)) * self.delta + \
-                             vov_new * np.sqrt(np.max(xt[:, i], 0) * self.delta) * z[:, i]
+                xt[:, i + 1] = (
+                    xt[:, i]
+                    + kappa_new * (theta_new - np.max(xt[:, i], 0)) * self.delta
+                    + vov_new * np.sqrt(np.max(xt[:, i], 0) * self.delta) * z[:, i]
+                )
 
-        elif scheme == 'Milstein':
+        elif scheme == "Milstein":
             z = np.random.normal(size=(self.path, self.step))
             for i in range(self.step):
-                xt[:, i+1] = xt[:, i] + kappa_new * (theta_new - np.max(xt[:, i], 0)) * self.delta + vov_new * \
-                             np.sqrt(np.max(xt[:, i], 0) * self.delta) * z[:, i] + \
-                             vov_new**2 * 0.25 * (z[:, i]**2 - 1) * self.delta
-
-        elif scheme == 'KJ':
+                xt[:, i + 1] = (
+                    xt[:, i]
+                    + kappa_new * (theta_new - np.max(xt[:, i], 0)) * self.delta
+                    + vov_new * np.sqrt(np.max(xt[:, i], 0) * self.delta) * z[:, i]
+                    + vov_new ** 2 * 0.25 * (z[:, i] ** 2 - 1) * self.delta
+                )
+
+        elif scheme == "KJ":
             z = np.random.normal(size=(self.path, self.step))
             for i in range(self.step):
-                xt[:, i+1] = (xt[:, i] + kappa_new * theta_new * self.delta + vov_new * \
-                             np.sqrt(np.max(xt[:, i], 0) * self.delta) * z[:, i] + \
-                             vov_new**2 * 0.25 * (z[:, i]**2 - 1) * self.delta) / (1 + kappa_new * self.delta)
+                xt[:, i + 1] = (
+                    xt[:, i]
+                    + kappa_new * theta_new * self.delta
+                    + vov_new * np.sqrt(np.max(xt[:, i], 0) * self.delta) * z[:, i]
+                    + vov_new ** 2 * 0.25 * (z[:, i] ** 2 - 1) * self.delta
+                ) / (1 + kappa_new * self.delta)
 
         # compute integral of vt, equivalent spot and vol
         vt = 1 / xt
@@ -175,33 +209,50 @@ def price(self, strike, spot, texp, sigma, delta, intr=0, divr=0, psi_c=1.5, pat
         vt[below_0] = 0
         vt_int = spint.simps(vt, dx=self.delta)
 
-        spot_cmc = spot * np.exp(self.rho / self.vov * (np.log(vt[:, -1] / vt[:, 0]) - self.kappa *
-                    (self.theta * texp - vt_int * (1 + self.vov**2 * 0.5 / self.kappa))) - self.rho ** 2 * vt_int / 2)
+        spot_cmc = spot * np.exp(
+            self.rho
+            / self.vov
+            * (
+                np.log(vt[:, -1] / vt[:, 0])
+                - self.kappa
+                * (self.theta * texp - vt_int * (1 + self.vov ** 2 * 0.5 / self.kappa))
+            )
+            - self.rho ** 2 * vt_int / 2
+        )
         vol_cmc = np.sqrt((1 - self.rho ** 2) * vt_int / texp)
 
         # compute bsm price vector for the given strike vector
         price_cmc = np.zeros_like(strike)
         for j in range(len(strike)):
-            price_cmc[j] = np.mean(self.bsm_model.price_formula(strike[j], spot_cmc, vol_cmc, texp, intr=intr, divr=divr))
+            price_cmc[j] = np.mean(
+                self.bsm_model.price_formula(
+                    strike[j], spot_cmc, vol_cmc, texp, intr=intr, divr=divr
+                )
+            )
 
         return price_cmc
 
     def prepare_rx(self):
-        '''
+        """
         Pre-calculate r(x) and store the result
         for TG scheme only
-        '''
-        fx = lambda rx: rx * spst.norm.pdf(rx) + spst.norm.cdf(rx) * (1 + rx ** 2) / \
-                        ((spst.norm.pdf(rx) + rx * spst.norm.cdf(rx)) ** 2) - 1
+        """
+        fx = (
+            lambda rx: rx * spst.norm.pdf(rx)
+            + spst.norm.cdf(rx)
+            * (1 + rx ** 2)
+            / ((spst.norm.pdf(rx) + rx * spst.norm.cdf(rx)) ** 2)
+            - 1
+        )
         rx_results = np.linspace(-2, 100, 10 ** 5)
         psi_points = fx(rx_results)
 
         return psi_points, rx_results
 
     def find_rx(self, psi):
-        '''
+        """
         Return r(psi) according to the pre_calculated results
-        '''
+        """
 
         if self.rx_results[self.psi_points >= psi].size == 0:
             print("Caution: input psi too large")
@@ -210,18 +261,35 @@ def find_rx(self, psi):
             print("Caution: input psi too small")
             return self.rx_results[0]
         else:
-            return (self.rx_results[self.psi_points >= psi][0] + self.rx_results[self.psi_points <= psi][-1])/2
+            return (
+                self.rx_results[self.psi_points >= psi][0]
+                + self.rx_results[self.psi_points <= psi][-1]
+            ) / 2
 
 
 class Sv32McAe:
-    def __init__(self, S0, Ks, T, r=0.05, sigma_0=1, beta=1, rho=-0.5, theta=1.5, kappa=2, vov=0.2, path_num = 100, cp=1):
-        assert(beta==1 or beta==0), 'Beta must be 0 or 1.'
+    def __init__(
+        self,
+        S0,
+        Ks,
+        T,
+        r=0.05,
+        sigma_0=1,
+        beta=1,
+        rho=-0.5,
+        theta=1.5,
+        kappa=2,
+        vov=0.2,
+        path_num=100,
+        cp=1,
+    ):
+        assert beta == 1 or beta == 0, "Beta must be 0 or 1."
         self.S0 = S0
         if not isinstance(Ks, np.ndarray):
             self.Ks = np.array([Ks])
         else:
             self.Ks = Ks
-        self.Ks = self.Ks.reshape(1,-1)
+        self.Ks = self.Ks.reshape(1, -1)
         self.T = T
         self.r = r
         self.sigma_0 = sigma_0
@@ -234,7 +302,7 @@ def __init__(self, S0, Ks, T, r=0.05, sigma_0=1, beta=1, rho=-0.5, theta=1.5, ka
         self.cp = cp
 
     def optionPrice_version1(self):
-        '''
+        """
         version-1: This one follows Professor Choi's idea exactly.
         First, we will simulate sigma_T with Heston model.
         Second, we will use every simulated sigma_T with characteristic function to simulate U_T from conditional distribution,
@@ -249,84 +317,99 @@ def optionPrice_version1(self):
         That is to say, the shape of self.Ks is the same as the shape of the first element,
         and our calculation supports for vectorization on strike price, not for time to maturity.
         Likewise, the second element is a np.ndarray, too, which contains all the stds of prices under every strike price.
-        '''
+        """
         self.simulate_sigma_2_T()
         self.simulate_U_T_version1()
         self.calForwardAndVolatility_version1()
         return self.calOutput()
 
     def simulate_sigma_2_T(self):
-        '''
+        """
         Generate path_num of sigma_T^2 with Heston Model,
         where self.X_T follow Heston Process, in which the parameters have been changed for 3/2 Model.
         And V_T are the inverses of X_T.
         Returns:
         None
-        '''
+        """
         self.X_T = self.simulate_Heston_outlayer()
-        self.V_T = 1/self.X_T
+        self.V_T = 1 / self.X_T
 
     def simulate_U_T_version1(self):
-        '''
+        """
         With characteristic function and path_num of sigma_T^2, we can get conditional distribution of U_T,
         where U_T is the integral of sigma_t^2 from t=0 to t=T.
         With conditional distribution, we will do one RNG on U_T with each sigma_T^2
         Returns:
             None
-        '''
+        """
         chfs = self.charFuncs()
         M1, M2 = self.calTwoMoments(chfs)
         mu, sigma = self.calLogNormalParas_version1(M1, M2)
-        self.U_T = spst.lognorm.rvs(mu, sigma) # Generate random numbers from lognormal distribution with mu and sigma
+        self.U_T = spst.lognorm.rvs(
+            mu, sigma
+        )  # Generate random numbers from lognormal distribution with mu and sigma
 
     def calForwardAndVolatility_version1(self):
-        '''
+        """
         Calculate E_F_T, expected forward price, and sigma_N or sigma_BS.
         Returns:
             None
-        '''
-        internal_term = self.rho / self.vov * (
-                    np.log(self.V_T / self.sigma_0 ** 2) - self.kappa * self.T * self.theta + (
-                    (self.kappa + self.vov ** 2 / 2) * self.U_T
-            ))
-        if self.beta==0:
-            self.E_F_T = self.S0+internal_term
-            self.sigma_N = np.sqrt((1-self.rho**2)*self.U_T/self.T)
+        """
+        internal_term = (
+            self.rho
+            / self.vov
+            * (
+                np.log(self.V_T / self.sigma_0 ** 2)
+                - self.kappa * self.T * self.theta
+                + ((self.kappa + self.vov ** 2 / 2) * self.U_T)
+            )
+        )
+        if self.beta == 0:
+            self.E_F_T = self.S0 + internal_term
+            self.sigma_N = np.sqrt((1 - self.rho ** 2) * self.U_T / self.T)
             self.sigma_N = self.sigma_N.reshape(-1, 1)
         else:
-            outside_term = self.rho**2*self.U_T/2
-            self.E_F_T = self.S0 * np.exp(internal_term-outside_term)
+            outside_term = self.rho ** 2 * self.U_T / 2
+            self.E_F_T = self.S0 * np.exp(internal_term - outside_term)
             self.sigma_BS = np.sqrt((1 - self.rho ** 2) * self.U_T / self.T)
             self.sigma_BS = self.sigma_BS.reshape(-1, 1)
-        self.E_F_T = self.E_F_T.reshape(-1,1)
+        self.E_F_T = self.E_F_T.reshape(-1, 1)
 
     def get_prices(self):
-        '''
+        """
         After running the simulation, we can get prices with this function, without run the simulation again.
         Returns:
             np.ndarray: Mean and Std of price output from AEMC
-        '''
+        """
         try:
             return self.option_prices, self.option_prices_std
         except:
-            print('Please run object function optionPrice_version1()/optionPrice_version2()/optionPrice_version3() first.')
+            print(
+                "Please run object function optionPrice_version1()/optionPrice_version2()/optionPrice_version3() first."
+            )
 
     def calOutput(self):
-        '''
+        """
         Calculate the mean and std of price output from AEMC
         Returns:
             np.ndarray: Mean and Std of price output from AEMC
-        '''
+        """
         if self.beta == 0:
-            prices = self.Bachelier(self.E_F_T, self.Ks, self.r, self.T, self.sigma_N, self.cp)
+            prices = self.Bachelier(
+                self.E_F_T, self.Ks, self.r, self.T, self.sigma_N, self.cp
+            )
         else:
-            prices = self.BSM(self.E_F_T, self.Ks, self.r, self.T, self.sigma_BS, self.cp)
-        self.option_prices, self.option_prices_std = prices.mean(axis=0), prices.std(axis=0)
+            prices = self.BSM(
+                self.E_F_T, self.Ks, self.r, self.T, self.sigma_BS, self.cp
+            )
+        self.option_prices, self.option_prices_std = prices.mean(axis=0), prices.std(
+            axis=0
+        )
         return self.option_prices, self.option_prices_std
 
     @staticmethod
     def calLogNormalParas_version1(M1, M2):
-        '''
+        """
         Calculate mu and sigma of lognormal distribution for first moment and second moment.
         Args:
             M1: np.ndarray, first moment
@@ -335,15 +418,15 @@ def calLogNormalParas_version1(M1, M2):
         Returns:
             M1: np.ndarray, mu of lognormal distribution
             np.sqrt(np.log(M2/M1**2)): np.ndarray, sigma of lognormal distribution
-        '''
+        """
         M1 = np.array(np.abs(M1).tolist(), dtype=float)
         M2 = np.array(np.abs(M2).tolist(), dtype=float)
         M2[np.isnan(M2)] = np.mean(M2)
-        return M1, np.sqrt(np.log(M2/M1**2))
+        return M1, np.sqrt(np.log(M2 / M1 ** 2))
 
     @staticmethod
     def calTwoMoments(chfs):
-        '''
+        """
         Calculate two moments from the characteristic function.
         Args:
             chfs: function with only one arg_in
@@ -351,32 +434,34 @@ def calTwoMoments(chfs):
         Returns:
             M1: np.ndarray, first moment
             M2: np.ndarray, second moment
-        '''
+        """
         M1 = derivative(chfs, x0=0, dx=0.00001, n=1)
         M2 = -derivative(chfs, x0=0, dx=0.00001, n=2)
         return M1, M2
 
     def charFuncs(self):
-        '''
+        """
         This is the characteristic function for the integration of sigma_T^2 from t=0 to t=T
         Returns:
             It return the characteristic function with only one independent variable a.
-        '''
+        """
         n = 4 * (self.kappa + self.vov ** 2) / self.vov ** 2
         j = -2 * self.kappa * self.theta / self.vov ** 2
         vega = n / 2 - 1
         delta = self.vov ** 2 * self.T / 4
         X_0 = 1 / self.sigma_0 ** 2
-        denominator = np.sinh(j*delta)
-        arg_in_Iv = j*np.sqrt(self.X_T*X_0)/denominator
+        denominator = np.sinh(j * delta)
+        arg_in_Iv = j * np.sqrt(self.X_T * X_0) / denominator
         besseli_ufun = np.frompyfunc(besseli, 2, 1)
+
         def char_func(a):
-            order_1 = np.sqrt(vega ** 2 + 8 * a *(-1j)/ self.vov ** 2)
-            return (besseli_ufun(order_1, arg_in_Iv)/besseli_ufun(vega, arg_in_Iv))
+            order_1 = np.sqrt(vega ** 2 + 8 * a * (-1j) / self.vov ** 2)
+            return besseli_ufun(order_1, arg_in_Iv) / besseli_ufun(vega, arg_in_Iv)
+
         return char_func
 
     def optionPrice_version2(self):
-        '''
+        """
         version-2: This one follows Professor Choi's idea exactly.
         First, we will simulate sigma_T with Heston model.
         Second, we direct use the first moment calculated from characteristic function as the proxy for U_T,
@@ -393,28 +478,33 @@ def optionPrice_version2(self):
         That is to say, the shape of self.Ks is the same as the shape of the first element,
         and our calculation supports for vectorization on strike price, not for time to maturity.
         Likewise, the second element is a np.ndarray, too, which contains all the stds of prices under every strike price.
-        '''
+        """
         self.simulate_sigma_2_T()
         M1 = self.simulate_M1()
         self.calForwardAndVolatility_version2(M1)
         return self.calOutput()
 
     def calForwardAndVolatility_version2(self, M1):
-        '''
+        """
         Calculate E_F_T, expected forward price, and sigma_N or sigma_BS.
         The slight difference with version 1 is that we substitute U_T with the mean of distribtion of U_T,
         which is M1 in this function.
         Returns:
             None
-        '''
-        internal_term = self.rho / self.vov * (
-                np.log(self.V_T / self.sigma_0 ** 2) - self.kappa * self.T * self.theta + (
-                (self.kappa + self.vov ** 2 / 2) * M1
-        ))
+        """
+        internal_term = (
+            self.rho
+            / self.vov
+            * (
+                np.log(self.V_T / self.sigma_0 ** 2)
+                - self.kappa * self.T * self.theta
+                + ((self.kappa + self.vov ** 2 / 2) * M1)
+            )
+        )
         if self.beta == 0:
             self.E_F_T = self.S0 + internal_term
             self.sigma_N = np.sqrt((1 - self.rho ** 2) * M1 / self.T)
-            self.sigma_N = self.sigma_N.reshape(-1,1)
+            self.sigma_N = self.sigma_N.reshape(-1, 1)
         else:
             outside_term = self.rho ** 2 * M1 / 2
             self.E_F_T = self.S0 * np.exp(internal_term - outside_term)
@@ -423,17 +513,17 @@ def calForwardAndVolatility_version2(self, M1):
         self.E_F_T = self.E_F_T.reshape(-1, 1)
 
     def simulate_M1(self):
-        '''
+        """
         First, we need the characteristic function.
         Second, we direct use the first moment calculated from characteristic function as the proxy for U_T,
         as the RNG from conditional distribution have oscillating output,
         which is not good for the standard deviation of option pricing in AEMC approach.
         Returns:
             M1: first moment as the proxy for U_T
-        '''
+        """
         chfs = self.charFuncs()
         M1 = self.calOneMoment(chfs)
-        M1 = np.array(np.abs(M1).tolist(),dtype=float)
+        M1 = np.array(np.abs(M1).tolist(), dtype=float)
         return M1
 
     @staticmethod
@@ -464,21 +554,25 @@ def calOneMoment(chfs):
     #     return mu, scale
 
     def simulate_Heston_outlayer(self):
-        '''
+        """
         Generate X_T for the preparation for V_T, which is sigma_T^2
         Returns:
             X_T: np.ndarray
                 the shape is (1,path_num)
-        '''
-        X_0 = 1/self.sigma_0**2
+        """
+        X_0 = 1 / self.sigma_0 ** 2
         delta = 4 * (self.kappa + self.vov ** 2) / self.vov ** 2
         c_T = self.c_function(self.T)
         alpha = X_0 / c_T
-        exp_term = np.exp(self.kappa*self.theta*self.T)
-        return np.random.noncentral_chisquare(delta, alpha, size = self.path_num)/exp_term*c_T
-
-    def c_function(self,t):
-        '''
+        exp_term = np.exp(self.kappa * self.theta * self.T)
+        return (
+            np.random.noncentral_chisquare(delta, alpha, size=self.path_num)
+            / exp_term
+            * c_T
+        )
+
+    def c_function(self, t):
+        """
         c function in the paper: Exact Simulation of the 3/2 Model. The formula is in page-6
         Args:
             t: float
@@ -487,12 +581,16 @@ def c_function(self,t):
         Returns:
             float
                 c(t)
-        '''
-        return self.vov**2*(np.exp(self.kappa*self.theta*t)-1)/(4*self.kappa*self.theta)
+        """
+        return (
+            self.vov ** 2
+            * (np.exp(self.kappa * self.theta * t) - 1)
+            / (4 * self.kappa * self.theta)
+        )
 
     @staticmethod
-    def BSM(S0,Ks,r,T,sigma,cp):
-        '''
+    def BSM(S0, Ks, r, T, sigma, cp):
+        """
         import from pyfeng package, to ues BSM model pricing formula
         Args:
             S0: Stock initial price
@@ -504,13 +602,13 @@ def BSM(S0,Ks,r,T,sigma,cp):
 
         Returns:
             prices: np,ndarray
-        '''
-        prices = Bsm.price_formula(Ks,S0,sigma,T,cp,r)
+        """
+        prices = Bsm.price_formula(Ks, S0, sigma, T, cp, r)
         return prices
 
     @staticmethod
     def Bachelier(S0, Ks, r, T, sigma, cp):
-        '''
+        """
         import from pyfeng package, to ues Normal model pricing formula
         Args:
             S0: Stock initial price
@@ -522,55 +620,60 @@ def Bachelier(S0, Ks, r, T, sigma, cp):
 
         Returns:
             prices: np,ndarray
-        '''
-        prices = Norm.price_formula(Ks,S0,sigma,T,cp,r)
+        """
+        prices = Norm.price_formula(Ks, S0, sigma, T, cp, r)
         return prices
 
     def imp_vol_for_bsm(self, K, price):
-        '''
+        """
         This function is used to calculate the implied volatility in BSM model, given strike prices and option prices.
         spot price in given in self.S0
-        '''
-        iv_func = lambda _vol: \
-            self.BSM(self.S0, K, self.r, self.T, _vol, self.cp) - price
+        """
+        iv_func = (
+            lambda _vol: self.BSM(self.S0, K, self.r, self.T, _vol, self.cp) - price
+        )
         vol = spop.brentq(iv_func, 0, 10000)
         return vol
 
     def imp_vol_for_normal(self, K, price):
-        '''
+        """
         This function is used to calculate the implied volatility in Normal model, given strike prices and option prices.
         spot price in given in self.S0
-        '''
-        iv_func = lambda _vol: \
-            self.Bachelier(self.S0, K, self.r, self.T, _vol, self.cp) - price
+        """
+        iv_func = (
+            lambda _vol: self.Bachelier(self.S0, K, self.r, self.T, _vol, self.cp)
+            - price
+        )
         vol = spop.brentq(iv_func, 0, 10000)
         return vol
 
     def impliedVolatility1(self):
-        '''
+        """
         Calculate the implied volatility with default strike prices, spot prices, and option price.
         Returns:
             output: np.ndarray, the shape is the same as the shape of strike prices.
-        '''
+        """
         output = []
         if self.beta == 1:
-            for k,p in zip(self.Ks.ravel(), self.option_prices.ravel()):
+            for k, p in zip(self.Ks.ravel(), self.option_prices.ravel()):
                 output.append(self.imp_vol_for_bsm(k, p))
         else:
             for k, p in zip(self.Ks.ravel(), self.option_prices.ravel()):
                 output.append(self.imp_vol_for_normal(k, p))
         return np.array(output)
 
-    def impliedVolatility2(self, model_type='bsm'):
-        '''
+    def impliedVolatility2(self, model_type="bsm"):
+        """
         Calculate the implied volatility with default strike prices, spot prices, and option price.
         Compared with impliedVolatility1, you can assign the type of implied volatility that you want.
         Returns:
             output: np.ndarray, the shape is the same as the shape of strike prices.
-        '''
-        assert model_type=='bsm' or  model_type=='normal', 'The model_type arg should be bsm or normal'
+        """
+        assert (
+            model_type == "bsm" or model_type == "normal"
+        ), "The model_type arg should be bsm or normal"
         output = []
-        if model_type=='bsm':
+        if model_type == "bsm":
             for k, p in zip(self.Ks.ravel(), self.option_prices.ravel()):
                 output.append(self.imp_vol_for_bsm(k, p))
         else:
diff --git a/pyfeng/sv_abc.py b/pyfeng/sv_abc.py
index 11de2ba..8f9b55d 100644
--- a/pyfeng/sv_abc.py
+++ b/pyfeng/sv_abc.py
@@ -8,7 +8,17 @@ class SvABC(smile.OptSmileABC, abc.ABC):
 
     vov, rho, mr, theta = 0.01, 0.0, 0.01, 1.0
 
-    def __init__(self, sigma, vov=0.01, rho=0.0, mr=0.01, theta=None, intr=0.0, divr=0.0, is_fwd=False):
+    def __init__(
+        self,
+        sigma,
+        vov=0.01,
+        rho=0.0,
+        mr=0.01,
+        theta=None,
+        intr=0.0,
+        divr=0.0,
+        is_fwd=False,
+    ):
         """
         Args:
             sigma: model volatility at t=0. variance = sigma**2
@@ -30,7 +40,12 @@ def __init__(self, sigma, vov=0.01, rho=0.0, mr=0.01, theta=None, intr=0.0, divr
 
     def params_kw(self):
         params1 = super().params_kw()
-        params2 = {"vov": self.vov, "rho": self.rho, "mr": self.mr, "sig_inf": self.theta}
+        params2 = {
+            "vov": self.vov,
+            "rho": self.rho,
+            "mr": self.mr,
+            "sig_inf": self.theta,
+        }
         return {**params1, **params2}
 
 
@@ -45,7 +60,6 @@ class CondMcBsmABC(smile.OptSmileABC, abc.ABC):
     rng = np.random.default_rng(None)
     antithetic = True
 
-    
     def set_mc_params(self, n_path=10000, dt=0.05, rn_seed=None, antithetic=True):
         """
         Set MC parameters
@@ -76,7 +90,7 @@ def tobs(self, texp):
         Returns:
             array of observation time
         """
-        n_steps = (texp//(2*self.dt)+1)*2
+        n_steps = (texp // (2 * self.dt) + 1) * 2
         tobs = np.arange(1, n_steps + 0.1) / n_steps * texp
         return tobs
 
@@ -100,7 +114,9 @@ def _bm_incr(self, tobs, cum=False, n_path=None):
         if self.antithetic:
             # generate random number in the order of path, time, asset and transposed
             # in this way, the same paths are generated when increasing n_path
-            bm_incr = self.rng.normal(size=(int(n_path/2), n_dt)).T * np.sqrt(dt[:, None])
+            bm_incr = self.rng.normal(size=(int(n_path / 2), n_dt)).T * np.sqrt(
+                dt[:, None]
+            )
             bm_incr = np.stack([bm_incr, -bm_incr], axis=-1).reshape((-1, n_path))
         else:
             bm_incr = np.random.randn(n_path, n_dt).T * np.sqrt(dt[:, None])
@@ -147,7 +163,7 @@ def price(self, strike, spot, texp, cp=1):
 
         fwd_cond, vol_cond = self.cond_fwd_vol(texp)
 
-        base_model = self.base_model(self.sigma*vol_cond)
+        base_model = self.base_model(self.sigma * vol_cond)
         price_grid = base_model.price(kk[:, None], fwd_cond, texp=texp, cp=cp)
 
         price = spot * np.mean(price_grid, axis=1)
diff --git a/setup.py b/setup.py
index 0b6cd6f..eb705d9 100644
--- a/setup.py
+++ b/setup.py
@@ -4,21 +4,19 @@
     long_description = fh.read()
 
 setup(
-    name='pyfeng',
-    version='0.1.6',
-    description='Python Financial Engineering',
+    name="pyfeng",
+    version="0.1.6",
+    description="Python Financial Engineering",
     long_description=long_description,
     long_description_content_type="text/markdown",
-    author='Jaehyuk Choi',
-    author_email='pyfe@eml.cc',
-    url='https://github.com/PyFE/PyFENG',
+    author="Jaehyuk Choi",
+    author_email="pyfe@eml.cc",
+    url="https://github.com/PyFE/PyFENG",
     project_urls={
-        'Documentation': 'https://pyfeng.readthedocs.io/',
-        'Bug Tracker': 'https://github.com/PyFE/PyFENG/issues',
+        "Documentation": "https://pyfeng.readthedocs.io/",
+        "Bug Tracker": "https://github.com/PyFE/PyFENG/issues",
     },
-    packages=['pyfeng'],
+    packages=["pyfeng"],
     test_suite="tests",
-    install_requires=[
-        'numpy', 'scipy'
-    ]
+    install_requires=["numpy", "scipy"],
 )


Test patch:
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/test_bsm_norm.py b/tests/test_bsm_norm.py
index 9336a9c..7ab9a8c 100644
--- a/tests/test_bsm_norm.py
+++ b/tests/test_bsm_norm.py
@@ -3,28 +3,48 @@
 import numpy as np
 import sys
 import os
+
 sys.path.insert(0, os.getcwd())
 import pyfeng as pf
 
 
 class TestBsmMethods(unittest.TestCase):
-
     def test_bsm_price(self):
         bsm = pf.Bsm(sigma=0.2, intr=0.05, divr=0.1)
         result = bsm.price(strike=np.arange(80, 126, 5), spot=100, texp=1.2)
-        expect_result = np.array([
-            15.71361973, 12.46799006, 9.69250803, 7.38869609,  5.52948546,
-            4.06773495, 2.94558338, 2.10255008, 1.48139131, 1.03159130
-        ])
+        expect_result = np.array(
+            [
+                15.71361973,
+                12.46799006,
+                9.69250803,
+                7.38869609,
+                5.52948546,
+                4.06773495,
+                2.94558338,
+                2.10255008,
+                1.48139131,
+                1.03159130,
+            ]
+        )
         np.testing.assert_almost_equal(result, expect_result)
 
     def test_norm_price(self):
         bsm = pf.Norm(sigma=20, intr=0.05, divr=0.1)
         result = bsm.price(strike=np.arange(80, 126, 5), spot=100, texp=1.2)
-        expect_result = np.array([
-            16.572334463, 13.264066668, 10.347114013, 7.849408092, 5.778270261,
-            4.119308420, 2.838573670, 1.887447196, 1.209104773, 0.745157542
-        ])
+        expect_result = np.array(
+            [
+                16.572334463,
+                13.264066668,
+                10.347114013,
+                7.849408092,
+                5.778270261,
+                4.119308420,
+                2.838573670,
+                1.887447196,
+                1.209104773,
+                0.745157542,
+            ]
+        )
         np.testing.assert_almost_equal(result, expect_result)
 
     def test_bsm_iv_greeks(self):
@@ -36,11 +56,11 @@ def test_bsm_iv_greeks(self):
             intr = np.random.uniform(0, 0.3)
             divr = np.random.uniform(0, 0.3)
             cp = 1 if np.random.rand() > 0.5 else -1
-            is_fwd = (np.random.rand() > 0.5)
+            is_fwd = np.random.rand() > 0.5
 
             # print( spot, strike, vol, texp, intr, divr, cp)
             m_bsm = pf.Bsm(sigma, intr=intr, divr=divr, is_fwd=is_fwd)
-            price = m_bsm.price(strike, spot, texp, cp),
+            price = (m_bsm.price(strike, spot, texp, cp),)
 
             # get implied vol
             iv = m_bsm.impvol(price, strike, spot, texp=texp, cp=cp)
@@ -74,7 +94,7 @@ def test_norm_iv_greeks(self):
             intr = np.random.uniform(0, 0.3)
             divr = np.random.uniform(0, 0.3)
             cp = 1 if np.random.rand() > 0.5 else -1
-            is_fwd = (np.random.rand() > 0.5)
+            is_fwd = np.random.rand() > 0.5
 
             # print( spot, strike, vol, texp, intr, divr, cp)
             m_norm = pf.Norm(sigma, intr=intr, divr=divr, is_fwd=is_fwd)
@@ -116,17 +136,19 @@ def test_BsmDisp(self):
         # DBS = Norm if beta=0
         dbs.beta = 0.0001
         dbs.is_fwd = True
-        norm = pf.Norm(sigma=dbs.sigma_disp*dbs.pivot, intr=0.05, divr=0.1, is_fwd=True)
+        norm = pf.Norm(
+            sigma=dbs.sigma_disp * dbs.pivot, intr=0.05, divr=0.1, is_fwd=True
+        )
         r1 = norm.price(strike, 100, 2.5, cp=-1)
         r2 = dbs.price(strike, 100, 2.5, cp=-1)
-        np.testing.assert_almost_equal(r1/r2, 1, decimal=4)
+        np.testing.assert_almost_equal(r1 / r2, 1, decimal=4)
         dbs.is_fwd = False
 
         # Approximate BSM vol
         dbs.beta = 0.2
-        v1 = dbs.vol_smile(strike, 100, 2.5, model='bsm')
-        v2 = dbs.vol_smile(strike, 100, 2.5, model='bsm-approx')
-        np.testing.assert_almost_equal(v1/v2, 1, decimal=4)
+        v1 = dbs.vol_smile(strike, 100, 2.5, model="bsm")
+        v2 = dbs.vol_smile(strike, 100, 2.5, model="bsm-approx")
+        np.testing.assert_almost_equal(v1 / v2, 1, decimal=4)
 
         p1 = dbs.price(strike, 100, 2.5)
         p2 = pf.Bsm(v1, intr=0.05, divr=0.1).price(strike, 100, 2.5)
@@ -134,15 +156,15 @@ def test_BsmDisp(self):
 
         # Approximate Bachelier vol
         dbs.beta = 0.8
-        v1 = dbs.vol_smile(strike, 100, 2.5, model='norm')
-        v2 = dbs.vol_smile(strike, 100, 2.5, model='norm-approx')
-        np.testing.assert_almost_equal(v1/v2, 1, decimal=4)
+        v1 = dbs.vol_smile(strike, 100, 2.5, model="norm")
+        v2 = dbs.vol_smile(strike, 100, 2.5, model="norm-approx")
+        np.testing.assert_almost_equal(v1 / v2, 1, decimal=4)
 
         p1 = dbs.price(strike, 100, 2.5)
         p2 = pf.Norm(v1, intr=0.05, divr=0.1).price(strike, 100, 2.5)
         np.testing.assert_almost_equal(p1, p2)
 
 
-if __name__ == '__main__':
-    print(f'Pyfeng loaded from {pf.__path__}')
+if __name__ == "__main__":
+    print(f"Pyfeng loaded from {pf.__path__}")
     unittest.main()
diff --git a/tests/test_multiasset.py b/tests/test_multiasset.py
index b28037d..7896c3b 100644
--- a/tests/test_multiasset.py
+++ b/tests/test_multiasset.py
@@ -2,6 +2,7 @@
 import numpy as np
 import sys
 import os
+
 sys.path.insert(0, os.getcwd())
 import pyfeng as pf
 
@@ -40,13 +41,17 @@ def test_BsmBasketLevy1992(self):
         # Table 2
         m = pf.BsmBasketLevy1992(sigma, rho)
         result = np.round(m.price(p_grid, fwd, texp), 2)
-        result2 = np.array([54.34, 47.52, 41.57, 36.4, 31.92, 28.05, 24.7, 21.8, 19.28, 17.1, 15.19])
+        result2 = np.array(
+            [54.34, 47.52, 41.57, 36.4, 31.92, 28.05, 24.7, 21.8, 19.28, 17.1, 15.19]
+        )
         np.testing.assert_almost_equal(result, result2)
 
         # Table 3
         m = pf.BsmBasketLevy1992(sigma, rho)
-        result = np.round(m.price(100, p_grid[:, None]*o4, texp), 2)
-        result2 = np.array([4.34, 7.52, 11.57, 16.4, 21.92, 28.05, 34.7, 41.8, 49.28, 57.1, 65.19])
+        result = np.round(m.price(100, p_grid[:, None] * o4, texp), 2)
+        result2 = np.array(
+            [4.34, 7.52, 11.57, 16.4, 21.92, 28.05, 34.7, 41.8, 49.28, 57.1, 65.19]
+        )
         np.testing.assert_almost_equal(result, result2)
 
         # Table 1
@@ -74,13 +79,17 @@ def test_BsmBasketMilevsky1998(self):
         m = pf.BsmBasketMilevsky1998(sigma, rho)
         result = np.round(m.price(p_grid, fwd, texp), 2)
         # Replaced 38.01 (Krekel et al., 2004) with 38.03
-        result2 = np.array([51.93, 44.41, 38.03, 32.68, 28.22, 24.5, 21.39, 18.77, 16.57, 14.7, 13.1])
+        result2 = np.array(
+            [51.93, 44.41, 38.03, 32.68, 28.22, 24.5, 21.39, 18.77, 16.57, 14.7, 13.1]
+        )
         np.testing.assert_almost_equal(result, result2)
 
         # Table 3
         m = pf.BsmBasketMilevsky1998(sigma, rho)
-        result = np.round(m.price(100, p_grid[:, None]*o4, texp), 2)
-        result2 = np.array([3.93, 6.56, 9.95, 14.1, 18.97, 24.5, 30.63, 37.32, 44.49, 52.08, 60.05])
+        result = np.round(m.price(100, p_grid[:, None] * o4, texp), 2)
+        result2 = np.array(
+            [3.93, 6.56, 9.95, 14.1, 18.97, 24.5, 30.63, 37.32, 44.49, 52.08, 60.05]
+        )
         np.testing.assert_almost_equal(result, result2)
 
         # Table 1
@@ -99,14 +108,14 @@ def test_BsmRainbow2(self):
         result = np.zeros_like(result2)
         fwds = [90, 100, 110]
         for k in range(len(fwds)):
-            result[k] = m.price(100, fwds[k]*o2, 3)
+            result[k] = m.price(100, fwds[k] * o2, 3)
 
     def test_BsmNormNdMc(self):
-        spot = np.ones(4)*100
-        sigma = np.ones(4)*0.4
+        spot = np.ones(4) * 100
+        sigma = np.ones(4) * 0.4
         texp = 5
         # Basket Option with equal weight
-        payoff = lambda x: np.fmax(np.mean(x,axis=1) - strike, 0) # Basket option
+        payoff = lambda x: np.fmax(np.mean(x, axis=1) - strike, 0)  # Basket option
         strikes = np.arange(80, 121, 10)
 
         # Test BsmNd
@@ -114,17 +123,17 @@ def test_BsmNormNdMc(self):
         m.simulate(tobs=[texp], n_path=20000)
         p = []
         for strike in strikes:
-           p.append(m.price_european(spot, texp, payoff))
+            p.append(m.price_european(spot, texp, payoff))
         p = np.array(p)
         p2 = np.array([36.31612946, 31.80861014, 27.91269315, 24.55319506, 21.62677625])
         np.testing.assert_almost_equal(p, p2)
 
         # Test NormNd
-        m = pf.NormNdMc(sigma*spot, cor=0.5, rn_seed=1234)
+        m = pf.NormNdMc(sigma * spot, cor=0.5, rn_seed=1234)
         m.simulate(tobs=[texp], n_path=20000)
         p = []
         for strike in strikes:
-           p.append(m.price_european(spot, texp, payoff))
+            p.append(m.price_european(spot, texp, payoff))
         p = np.array(p)
         p2 = np.array([39.42304794, 33.60383167, 28.32667559, 23.60383167, 19.42304794])
         np.testing.assert_almost_equal(p, p2)
@@ -135,7 +144,7 @@ def test_BsmBasket1Bm(self):
             n = np.random.randint(1, 8)
             spot = np.random.uniform(80, 120, size=n)
             strike = np.random.uniform(80, 120, size=10)
-            sigma = np.random.uniform(0.01, 1)*np.ones(n)
+            sigma = np.random.uniform(0.01, 1) * np.ones(n)
             texp = np.random.uniform(0.1, 10)
             intr = np.random.uniform(0, 0.1)
             divr = np.random.uniform(0, 0.1)
@@ -143,17 +152,18 @@ def test_BsmBasket1Bm(self):
             weight /= np.sum(weight)
 
             cp = np.where(np.random.rand(10) > 0.5, 1, -1)
-            is_fwd = (np.random.rand() > 0.5)
+            is_fwd = np.random.rand() > 0.5
 
-            m = pf.BsmBasket1Bm(sigma, weight=weight, intr=intr, divr=divr, is_fwd=is_fwd)
+            m = pf.BsmBasket1Bm(
+                sigma, weight=weight, intr=intr, divr=divr, is_fwd=is_fwd
+            )
             p = m.price(strike, spot, texp, cp)
 
             m2 = pf.Bsm(sigma[0], intr=intr, divr=divr, is_fwd=is_fwd)
-            p2 = m2.price(strike, np.sum(spot*weight), texp, cp)
+            p2 = m2.price(strike, np.sum(spot * weight), texp, cp)
             np.testing.assert_almost_equal(p, p2)
 
 
-
-if __name__ == '__main__':
-    print(f'Pyfeng loaded from {pf.__path__}')
+if __name__ == "__main__":
+    print(f"Pyfeng loaded from {pf.__path__}")
     unittest.main()
diff --git a/tests/test_sabr.py b/tests/test_sabr.py
index 7d1d922..fb96ab0 100644
--- a/tests/test_sabr.py
+++ b/tests/test_sabr.py
@@ -3,27 +3,27 @@
 import numpy as np
 import sys
 import os
+
 sys.path.insert(0, os.getcwd())
 import pyfeng as pf
 
 
 class TestSabr(unittest.TestCase):
-
     def test_Hagan2002(self):
         for k in list(range(1, 19)) + [22, 23]:
             m, df, rv = pf.SabrHagan2002.init_benchmark(k)
-            #ref = rv['ref']
-            #print(f'Sheet {k:02d}: {ref}')
-            v1 = np.round(m.vol_for_price(**rv['args_pricing']), 4)
-            v2 = df['IV Hagan'].values
+            # ref = rv['ref']
+            # print(f'Sheet {k:02d}: {ref}')
+            v1 = np.round(m.vol_for_price(**rv["args_pricing"]), 4)
+            v2 = df["IV Hagan"].values
             np.testing.assert_almost_equal(v1, v2)
 
     def test_SabrNorm(self):
         for k in [22, 23]:
             m, df, rv = pf.SabrNorm.init_benchmark(k)
-            v1 = m.price(**rv['args_pricing'])
+            v1 = m.price(**rv["args_pricing"])
             m, df, rv = pf.SabrChoiWu2021H.init_benchmark(k)
-            v2 = m.price(**rv['args_pricing'])
+            v2 = m.price(**rv["args_pricing"])
             np.testing.assert_almost_equal(v1, v2)
 
     def test_SabrNormATM(self):
@@ -42,21 +42,21 @@ def test_PaulotBsm(self):
         for k in list(range(1, 19)):
             m, df, rv = pf.SabrChoiWu2021P.init_benchmark(k)
             m._base_beta = 1.0  # For Paulot's BS volatility approximation
-            #print(f'Sheet {k:02d}: {ref}')
-            v1 = np.round(m.vol_for_price(**rv['args_pricing']), 4)
-            v2 = df['IV HL-P'].values
+            # print(f'Sheet {k:02d}: {ref}')
+            v1 = np.round(m.vol_for_price(**rv["args_pricing"]), 4)
+            v2 = df["IV HL-P"].values
             np.testing.assert_almost_equal(v1, v2)
 
     def test_UnCorrChoiWu2021(self):
         # Param Set 19: Table 7 (Case III.C) in Cai et al. (2017). https://doi.org/10.1287/opre.2017.1617
         m, df, rv = pf.SabrUncorrChoiWu2021.init_benchmark(19)
-        mass = m.mass_zero(rv['args_pricing']['spot'], rv['args_pricing']['texp'])
-        p = m.price(**rv['args_pricing'])
+        mass = m.mass_zero(rv["args_pricing"]["spot"], rv["args_pricing"]["texp"])
+        p = m.price(**rv["args_pricing"])
 
         mass2 = 0.7623543217183134
-        p2 = np.array([
-            0.04533777, 0.04095806, 0.03889591, 0.03692339, 0.03324944, 0.02992918
-        ])
+        p2 = np.array(
+            [0.04533777, 0.04095806, 0.03889591, 0.03692339, 0.03324944, 0.02992918]
+        )
 
         np.testing.assert_almost_equal(mass, mass2)
         np.testing.assert_almost_equal(p, p2)
@@ -65,10 +65,10 @@ def test_CondMc(self):
         for k in [19, 20]:  # can test 22 (Korn&Tang) also, but difficult to pass
             m, df, rv = pf.SabrCondMc.init_benchmark(k)
             m.set_mc_params(n_path=5e4, dt=0.05, rn_seed=1234)
-            p = m.price(**rv['args_pricing'])
-            np.testing.assert_almost_equal(p, rv['val'], decimal=4)
+            p = m.price(**rv["args_pricing"])
+            np.testing.assert_almost_equal(p, rv["val"], decimal=4)
 
 
-if __name__ == '__main__':
-    print(f'Pyfeng loaded from {pf.__path__}')
+if __name__ == "__main__":
+    print(f"Pyfeng loaded from {pf.__path__}")
     unittest.main()