Repository: PyFE/PyFENG (language: python)

Issue:
Experimental (ex) refactoring



Discussion / hints:

Fix patch:
diff --git a/pyfeng/__init__.py b/pyfeng/__init__.py
index 9e78288..b7e04cf 100644
--- a/pyfeng/__init__.py
+++ b/pyfeng/__init__.py
@@ -4,8 +4,12 @@
 from .gamma import InvGam, InvGauss
 from .sabr import SabrHagan2002, SabrNorm, SabrLorig2017, SabrChoiWu2021H, SabrChoiWu2021P
 from .sabr_int import SabrUncorrChoiWu2021
-from .nsvh import Nsvh1
+from .sabr_mc import SabrCondMc
+from .nsvh import Nsvh1, NsvhMc
 from .multiasset import BsmSpreadKirk, BsmSpreadBjerksund2014, NormBasket, NormSpread, BsmBasketLevy1992, BsmMax2, \
-    BsmBasketMilevsky1998, BsmBasket1Bm
+    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/_ex.py b/pyfeng/_ex.py
new file mode 100644
index 0000000..52e7155
--- /dev/null
+++ b/pyfeng/_ex.py
@@ -0,0 +1,5 @@
+# SV models (CMC, AE) from ASP 2021
+from .heston_mc import HestonCondMc, HestonCondMcQE, HestonMcAe
+from .sv32_mc import Sv32CondMcQE, Sv32McAe
+from .ousv import OusvIft, OusvIft, OusvCondMC
+from .garch import GarchCondMC, GarchApproxUncor
diff --git a/pyfeng/asian.py b/pyfeng/asian.py
new file mode 100644
index 0000000..289a8a4
--- /dev/null
+++ b/pyfeng/asian.py
@@ -0,0 +1,220 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue May  4 18:29:04 2021
+
+@author: cy-wang15 / Zay
+"""
+
+import math
+import numpy as np
+import mpmath as m
+import sympy
+from scipy.misc import derivative
+from . import opt_abc as opt
+from . import nsvh
+
+
+class BsmAsianJsu(opt.OptMaABC):
+    """
+
+    Johnson's SU distribution approximation for Asian option pricing under the BSM model.
+
+    Note: Johnson's SU distribution is the solution of NSVh with NSVh with lambda = 1.
+
+    References:
+        [1] Posner, S. E., & Milevsky, M. A. (1998). Valuing exotic options by approximating the SPD
+        with higher moments. The Journal of Financial Engineering, 7(2). https://ssrn.com/abstract=108539
+
+        [2] Choi, J., Liu, C., & Seo, B. K. (2019). Hyperbolic normal stochastic volatility model.
+        Journal of Futures Markets, 39(2), 186–204. https://doi.org/10.1002/fut.21967
+
+    """
+
+    def moments(self, spot, texp, n=1):
+        """
+
+        The moments of the Asian option with a log-normal underlying asset.
+
+        Args:
+            spot: spot price
+            texp: time to expiry
+            n: moment order
+
+        References:
+            Geman, H., & Yor, M. (1993). Bessel processes, Asian options, and perpetuities.
+            Mathematical finance, 3(4), 349-375.
+
+        Returns: the nth moment
+
+        """
+        lam = self.sigma[0]
+        v = (self.intr - self.divr - lam**2 / 2) / lam
+        beta = v / lam
+
+        ds = []
+        for j in range(n+1):
+            item0 = 2 ** n
+            for i in range(n+1):
+                if i != j:
+                    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
+
+        return moment
+
+    def moment_mvsk(self, spot, texp):
+        """
+
+        Return mean, variance, skewness, kurtosis for Asian options.
+
+        Args:
+            spot: spot price
+            texp: time to expiry
+
+        Returns: mean, variance, skewness, kurtosis of Asian options
+
+        """
+        moments = []
+        for i in range(1, 5):
+            moments.append(self.moments(spot, texp, i))
+
+        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
+
+        return mu, var, skew, kurt
+
+    def price(self, strike, spot, texp, cp=1):
+        """
+
+        Asian options price.
+        Args:
+            strike: strike price
+            spot: spot price
+            texp: time to expiry
+            cp: 1/-1 for call/put option
+        Returns: Basket options price
+
+        """
+        df = np.exp(-texp * self.intr)
+
+        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)
+        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
+        self.vol = vol
+        self.texp = texp
+        self.strike = strike
+        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):
+        eigenval = []
+        if nu <= -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):
+        eigval_p = []
+        for i in range(len(eigval)):
+            real_i = np.array(eigval[i], dtype=complex).real[0]
+            if real_i > 0:
+                eigval_p.append(real_i)
+            else:
+                pass
+        return eigval_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)
+            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))
+        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))
+        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))
+        return p1 * p2 * const * w1 * m1
+        
+    def exact_asian_price(self,n_eig):
+        intr = self.intr
+        divr = self.divr
+        vol = self.vol
+        texp = self.texp
+        strike = self.strike
+        spot = self.spot
+        b = self.b
+        call = self.call
+        nu = self.nu
+        tau = self.tau
+        k = self.k
+        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_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_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)) 
+        
+        if call == False:
+            return(put_price.real)
+        elif call == True:
+            return(call_price.real)
+        else:
+            return("Wrong")
diff --git a/pyfeng/bsm.py b/pyfeng/bsm.py
index 6477eff..4af7940 100644
--- a/pyfeng/bsm.py
+++ b/pyfeng/bsm.py
@@ -307,6 +307,25 @@ def price_barrier(self, strike, barrier, spot, texp, cp=1, io=-1):
 
         return p
 
+    def moments_vsk(self, texp=1):
+        """
+        Variance, skewness, and ex-kurtosis. Assume mean=1.
+
+        Args:
+            texp: time-to-expiry
+
+        Returns:
+            (variance, skewness, and ex-kurtosis)
+
+        References:
+            https://en.wikipedia.org/wiki/Log-normal_distribution
+        """
+        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
+        return var, skew, exkurt
+
 
 class BsmDisp(smile.OptSmileABC, Bsm):
     """
@@ -460,7 +479,12 @@ def vol_smile(self, strike, spot, texp, cp=1, model='bsm'):
 
         return vol
 
-    def price_barrier(self, strike, barrier, spot, texp, cp=1, io=-1):
-        fwd = self.forward(spot, texp)
-        return (1/self.beta)*super().price_barrier(
-            self.disp_spot(strike), self.disp_strike(barrier, strike), self.disp_spot(fwd), texp, cp=cp, io=io)
+    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)
+
+    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
diff --git a/pyfeng/garch.py b/pyfeng/garch.py
new file mode 100644
index 0000000..b6c9147
--- /dev/null
+++ b/pyfeng/garch.py
@@ -0,0 +1,106 @@
+import numpy as np
+import scipy.integrate as scint
+from . import sv_abc as sv
+
+
+class GarchApproxUncor(sv.SvABC):
+    """
+    The implementation of Barone-Adesi et al. (2004)'s approximation pricing formula for European
+    options under uncorrelated GARCH diffusion model.
+
+    References: Barone-Adesi, G., Rasmussen, H., Ravanelli, C., 2005. An option pricing formula for the GARCH diffusion model. Computational Statistics & Data Analysis, 2nd CSDA Special Issue on Computational Econometrics 49, 287–310. https://doi.org/10.1016/j.csda.2004.05.014
+
+    This method is only used to compare with the method GarchCondMC.
+    """
+
+
+class GarchCondMC(sv.SvABC, sv.CondMcBsmABC):
+    """
+    Garch model with conditional Monte-Carlo simulation
+    The SDE of SV is: dv_t = mr * (theta - v_t) dt + vov * v_t dB_T
+    """
+
+    def vol_paths(self, tobs, n_path=None):
+        """
+        Milstein Schemes:
+        w_(t+dt) = w_t + (mr * theta * exp(-w_t) - mr - vov^2 / 2) * dt + vov * Z * sqrt(dt)
+        v_t = exp(w_t)
+        Args:
+            mr: coefficient of dt
+            theta: the long term average
+            Z : std normal distributed RN
+            dt : delta t, time step
+
+        Returns: volatility path (time, path) including the value at t=0
+        """
+        n_dt = len(tobs)
+        n_path = n_path or self.n_path
+        rn_norm = self.rng.normal(size=(n_dt, int(n_path)))
+
+        w_t = np.zeros((n_dt + 1, int(n_path)))
+        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, :]
+
+        v_t = np.exp(w_t)
+
+        return v_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
+
+            Args:
+                theta: the long term average
+                mr: coefficient of dt
+                texp: time-to-expiry
+
+            Returns: (forward, volatility)
+        """
+        rhoc = np.sqrt(1.0 - self.rho ** 2)
+        tobs = self.tobs(texp)
+        n_dt = len(tobs)
+        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
+
+        vol_cond = rhoc * np.sqrt(int_var / texp)
+
+        return fwd_cond, vol_cond
+
+    def price(self, strike, spot, texp, cp=1):
+        """
+        Calculate option price based on BSM
+        Args:
+            strike: strike price
+            spot: spot price
+            texp: time to maturity
+            cp: cp=1 if call option else put option
+
+        Returns: price
+        """
+        price = []
+        texp = [texp] if isinstance(texp, (int, float)) else texp
+        for t in texp:
+            kk = strike / spot
+            kk = np.atleast_1d(kk)
+
+            fwd_cond, vol_cond = self.cond_fwd_vol(t)
+
+            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)
+            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/heston_mc.py b/pyfeng/heston_mc.py
new file mode 100644
index 0000000..a88d756
--- /dev/null
+++ b/pyfeng/heston_mc.py
@@ -0,0 +1,397 @@
+import numpy as np
+import scipy.integrate as scint
+import scipy.stats as spst
+import scipy.special as spsp
+import scipy.integrate as spint
+from scipy.misc import derivative
+from . import bsm
+from . import sv_abc as sv
+
+
+class HestonCondMc(sv.SvABC, sv.CondMcBsmABC):
+    """
+    Heston model with conditional Monte-Carlo simulation
+    """
+
+    def vol_paths(self, tobs):
+        dt = np.diff(tobs, prepend=0)
+        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[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 < 0] = 0  # variance should be larger than zero
+            vv_path[i+1, :] = vv
+
+        # return normalized sigma, e.g., sigma(0) = 1
+        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
+        vv_ratio = sigma_paths[-1, :]
+        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)
+
+        # 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.
+    Volatility (variance) of the price is assumed to follow a CIR process.
+
+    Example:
+        >>> import numpy as np
+        >>> import pyfeng as pf
+        >>> strike = [100.0, 140.0, 70.0]
+        >>> forward = 100
+        >>> delta = [1, 1/2, 1/4, 1/8, 1/16, 1/32]
+        >>> vov, kappa, rho, texp, theta, sigma = [1, 0.5, -0.9, 10, 0.04, 0.2]
+        >>> heston_cmc_qe = pf.HestonCondMcQE(vov=vov, kappa=kappa, rho=rho, theta=theta)
+        >>> price_cmc = np.zeros([len(delta), len(strike)])
+        >>> for d in range(len(delta)):
+        >>>     price_cmc[d, :] = heston_cmc_qe.price(strike, forward, texp, sigma=sigma, delta=delta[d], path=1e5, seed=123456)
+        >>> price_cmc
+        array([[14.52722285,  0.19584722, 37.20591415],
+               [13.56691261,  0.26568546, 36.12295964],
+               [13.22061601,  0.29003533, 35.9154245 ],
+               [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:
+            vov: volatility of variance, strictly positive
+            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
+        self.theta = theta
+
+        self.psi_points = None  # for TG scheme only
+        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):
+        '''
+        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.
+
+        Args:
+            strike: strike price, in vector form
+            spot: spot (or forward)
+            texp: time to expiry
+            sigma: initial volatility
+            delta: length of each time step
+            intr: interest rate (domestic interest rate)
+            divr: dividend/convenience yield (foreign interest rate)
+            psi_c: critical value for psi, lying in [1, 2]
+            path: number of vol paths generated
+            scheme: discretization scheme for vt, {'QE', 'TG', 'Euler', 'Milstein', 'KJ'}
+            seed: random seed for rv generation
+
+        Return:
+            BSM price vector for all strikes
+        '''
+        self.sigma = sigma
+        self.bsm_model = bsm.Bsm(self.sigma, intr=intr, divr=divr)
+        self.delta = delta
+        self.path = int(path)
+        self.step = int(texp / self.delta)
+
+        vt = self.sigma ** 2 * np.ones([self.path, self.step + 1])
+        np.random.seed(seed)
+
+        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)
+                psi = s2 / m ** 2
+
+                # compute vt(i+1) given psi
+                below = np.where(psi <= psi_c)[0]
+                ins = 2 * psi[below] ** -1
+                b2 = ins - 1 + np.sqrt(ins * (ins - 1))
+                b = np.sqrt(b2)
+                a = m[below] / (1 + b2)
+                z = spst.norm.ppf(u[below, i])
+                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]))
+                    else:
+                        vt[above[k], i+1] = 0
+
+        elif scheme == 'TG':
+            if np.all(self.rx_results) == None:
+                self.psi_points, self.rx_results = self.prepare_rx()
+
+            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)
+                psi = s2 / m ** 2
+
+                rx = np.array([self.find_rx(j) for j in psi])
+
+                z = np.random.normal(size=(self.path, self.step))
+                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))
+
+                vt[:, i+1] = np.fmax(mu_v[:, i] + sigma_v[:, i] * z[:, i], 0)
+
+        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]
+            below_0 = np.where(vt < 0)
+            vt[below_0] = 0
+
+        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
+            below_0 = np.where(vt < 0)
+            vt[below_0] = 0
+
+        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)
+            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)
+        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))
+
+        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
+        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")
+            return self.rx_results[-1]
+        elif self.rx_results[self.psi_points <= psi].size == 0:
+            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
+
+
+class HestonMcAe:
+    """
+    Almost exact MC for Heston model.
+
+    Underlying price is assumed to follow a geometric Brownian motion.
+    Volatility (variance) of the price is assumed to follow a CIR process.
+    Example:
+        >>> strike = 100
+        >>> spot = 100
+        >>> vov, kappa, rho, texp, theta, sigma = 0.61, 6.21, -0.7, 1, 0.019, 0.10201
+        >>> heston_ae = HestonMcAe(vov, kappa, rho, theta, r)
+        >>> price_ae = heston_ae.price(strike, spot, texp, sigma_0, intr=0, divr=0)
+        >>> price_ae
+        8.946951375550809
+    """
+    def __init__(self, vov=1, kappa=0.5, rho=-0.9, theta=0.04, r=0):
+        """
+        Initiate a Heston model
+
+        Args:
+            vov: volatility of variance, strictly positive
+            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
+            r：the drift item
+        """
+        self.vov = vov
+        self.theta = theta
+        self.kappa = kappa
+        self.rho = rho
+        self.r = r
+
+    def ch_f(self, texp, sigma_0, sigma_t, chi_dim):
+        """
+
+        Args:
+            texp: time to expiry
+            sigma_0: initial volatility
+            sigma_t: volatility at time T
+            chi_dim: dimensions of chisquare distribution
+
+        Returns:
+            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)
+
+        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 = lambda a: ch_f_part_1(a) * ch_f_part_2(a) * ch_f_part_3(a)
+        return ch_f
+
+    def gen_vov_t(self, chi_dim, chi_lambda, texp, n_paths):
+        """
+
+        Args:
+            chi_dim: dimensions of chisquare distribution
+            chi_lambda: the skewing item of chisquare distribution
+            texp: time to expiry
+            n_paths: number of vol paths generated
+
+        Returns:
+            sigma_t: volatility at time T
+
+        """
+        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))
+        return sigma_t
+
+    def gen_s_t(self, spot, sigma_t, sigma_0, texp, integral_sigma_t, n_paths):
+        """
+
+        Args:
+            spot: spot (or forward)
+            sigma_t: volatility at time T
+            sigma_0: initial volatility
+            texp: time to expiry
+            integral_sigma_t: samples from the distribution of integral sigma_t
+            n_paths: number of vol paths generated
+
+        Returns:
+            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)
+        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):
+        """
+        Args:
+            strike: strike price
+            spot: spot (or forward)
+            texp: time to expiry
+            sigma_0: initial volatility
+            intr: interest rate (domestic interest rate)
+            divr: dividend/convenience yield (foreign interest rate)
+            n_paths: number of vol paths generated
+            seed: random seed for rv generation
+        Returns:
+            price_ae: option price
+        """
+        if 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
+
+        sigma_t = self.gen_vov_t(chi_dim, chi_lambda, texp, n_paths)
+
+        ch_f = self.ch_f(texp, sigma_0, sigma_t, chi_dim)
+
+        moment_1st = (derivative(ch_f, 0, n=1, dx=1e-5) / 1j).real
+        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)
+            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)
+
+        elif dis_can == "Log-normal":
+            scale_ln = np.sqrt(np.log(moment_2st) - 2 * np.log(moment_1st))
+            miu_ln = np.log(moment_1st) - 0.5 * scale_ln ** 2
+            integral_sigma_t = np.random.lognormal(miu_ln, scale_ln)
+            s_t = self.gen_s_t(spot, sigma_t, sigma_0, texp, integral_sigma_t, n_paths)
+        else:
+            print("This function is not currently a candidate function!")
+            return -1
+
+        if call:
+            price_ae = np.fmax(s_t - strike, 0).mean()
+        else:
+            price_ae = np.fmax(strike - s_t, 0).mean()
+
+        return np.exp(- self.r * texp) * price_ae
diff --git a/pyfeng/multiasset.py b/pyfeng/multiasset.py
index 5ae3e4a..5407e99 100644
--- a/pyfeng/multiasset.py
+++ b/pyfeng/multiasset.py
@@ -1,9 +1,11 @@
 import scipy.stats as spst
 import warnings
 import numpy as np
+from itertools import product, combinations
 from . import bsm
 from . import norm
 from . import gamma
+from . import nsvh
 from . import opt_abc as opt
 
 
@@ -388,4 +390,122 @@ def price(self, strike, spot, texp, cp=1):
         price = np.sum(fwd_basket*spst.norm.cdf(d1), axis=-1)
         price -= strike * spst.norm.cdf(d2)
         price *= cp*df
-        return price
\ No newline at end of file
+        return price
+
+
+class BsmBasketLowerBound(NormBasket):
+
+    def V_func(self, fwd):
+        """
+        Generate factor matrix V
+        Args:
+            fwd: forward rates of basket
+        Returns:
+            V matrix
+        """
+        gg = self.weight * 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
+        V1 = self.chol_m @ Q1
+        V1 = V1[:, None]
+
+        # obtain full V
+        e1 = np.zeros_like(self.sigma)
+        e1[0] = 1
+        v = (Q1 - e1) / np.linalg.norm(Q1 - e1)
+        v = v[:, None]
+        R = np.eye(self.n_asset) - 2 * v @ v.T
+
+        # singular value decomposition
+        U, D, Q = np.linalg.svd(self.chol_m @ R[:, 1:], full_matrices=False)
+        V = np.hstack((V1, U @ np.diag(D)))
+
+        return V
+
+    def price(self, strike, spot, texp, cp=1):
+        fwd, df, _ = self._fwd_factor(spot, texp)
+        V = self.V_func(fwd)
+        sigma1 = V[:, 0]
+        m = BsmBasket1Bm(sigma1, is_fwd=True)
+        price = m.price(strike, fwd, texp)
+        return price
+
+
+class BsmBasketJsu(NormBasket):
+    """
+
+    Johnson's SU distribution approximation for Basket option pricing under the multiasset BSM model.
+
+    Note: Johnson's SU distribution is the solution of NSVh with NSVh with lambda = 1.
+
+    References:
+        [1] Posner, S. E., & Milevsky, M. A. (1998). Valuing exotic options by approximating the SPD
+        with higher moments. The Journal of Financial Engineering, 7(2). https://ssrn.com/abstract=108539
+
+        [2] Choi, J., Liu, C., & Seo, B. K. (2019). Hyperbolic normal stochastic volatility model.
+        Journal of Futures Markets, 39(2), 186–204. https://doi.org/10.1002/fut.21967
+
+    """
+
+    def moment_vsk(self, fwd, texp):
+        """
+
+        Return variance, skewness, kurtosis for Basket options.
+
+        Args:
+            fwd: forward price
+            texp: time to expiry
+
+        Returns: variance, skewness, kurtosis of Basket options
+
+        """
+        n = len(self.weight)
+
+        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)
+
+        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)
+
+        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)
+
+        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
+
+        return var, skew, kurt
+
+    def price(self, strike, spot, texp, cp=1):
+        """
+
+        Basket options price.
+        Args:
+            strike: strike price
+            spot: spot price
+            texp: time to expiry
+            cp: 1/-1 for call/put option
+        Returns: Basket options price
+
+        """
+        fwd, df, _ = self._fwd_factor(spot, texp)
+        assert fwd.shape[-1] == self.n_asset
+
+        fwd_basket = fwd @ self.weight
+
+        var, skew, kurt = self.moment_vsk(fwd, texp)
+
+        m = nsvh.Nsvh1(sigma=self.sigma)
+        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
new file mode 100644
index 0000000..e867c0c
--- /dev/null
+++ b/pyfeng/multiasset_Ju2002.py
@@ -0,0 +1,244 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Apr 28 09:08:29 2021
+
+@author: Yuze, Lantian
+"""
+from . import multiasset as ma
+from . import opt_abc as opt
+import numpy as np
+import scipy.stats as spst
+
+
+class BsmBasketAsianJu2002(ma.NormBasket):
+    def __init__(self, sigma, cor=None, weight=None, intr=0.0, divr=0.0, is_fwd=False):
+        """
+        Args:
+            sigma: model volatilities of `n_asset` assets. (n_asset, ) array
+            cor: correlation. If matrix, used as it is. (n_asset, n_asset)
+                If scalar, correlation matrix is constructed with all same off-diagonal values.
+            weight: asset weights, If None, equally weighted as 1/n_asset
+                If scalar, equal weights of the value
+                If 1-D array, uses as it is. (n_asset, )
+            intr: interest rate (domestic interest rate)
+            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)
+        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):
+            if np.isscalar(spot):
+                 spot = np.full(num_asset,spot)
+            if np.isscalar(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))
+        else:
+            if np.isscalar(spot):
+                 spot = np.full(num_asset,spot)
+            if np.isscalar(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))
+        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))
+            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
+        else:
+            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]
+        self.av_rho = av_rho
+    
+    def u1(self, spot, texp):
+        #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#
+        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])
+        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]
+        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)
+        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)
+        return u2_3rd_value
+    
+    def ak_bar(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)
+    
+    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()
+        return value
+    
+    def e_a13_a3(self):
+        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 *= 6
+        return value
+    
+    def e_a23(self):
+        value = 0
+        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)
+        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 *= 8
+        return value
+    
+    def func_a1(self, z):
+        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)
+    
+    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)
+    
+    def func_b1(self, spot, texp, z):
+        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
+    
+    def func_c1(self, 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)
+    
+    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)
+    
+    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
+    
+    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)))
+    
+    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))
+    
+    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)))
+                     
+    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)
+
+    def price(self, strike, spot, texp, cp=1, basket = True):
+        if np.isscalar(spot):
+             spot = np.full(num_asset,spot)
+        if np.isscalar(self.divr):
+             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))
+        if cp == 1:
+            return bc
+        elif cp == -1:
+            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):
+            print("spot should not be array")
+            return 0
+        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)
+            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))
+        if cp == 1:
+            return bc
+        elif cp == -1:
+            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 f0ee7aa..3d2be8e 100644
--- a/pyfeng/multiasset_mc.py
+++ b/pyfeng/multiasset_mc.py
@@ -40,6 +40,14 @@ def __init__(self, sigma, cor=None, intr=0.0, divr=0.0, rn_seed=None):
         super().__init__(sigma, cor=cor, intr=intr, divr=divr, is_fwd=False)
 
     def set_mc_params(self, n_path, rn_seed=None, antithetic=True):
+        """
+        Set MC parameters
+
+        Args:
+            n_path: number of paths
+            rn_seed: random number seed for rn generator.
+            antithetic: antithetic
+        """
         self.n_path = n_path
         self.rn_seed = rn_seed
         self.antithetic = antithetic
@@ -51,7 +59,6 @@ def _bm_incr(self, tobs, n_path=None):
         Args:
             tobs: array of observation times
             n_path: number of paths. If None (default), use the stored one.
-            store: if True (default), save the result to self.path_stored
 
         Returns:
             price path (time, path, asset)
diff --git a/pyfeng/nsvh.py b/pyfeng/nsvh.py
index 29c647c..0320a6c 100644
--- a/pyfeng/nsvh.py
+++ b/pyfeng/nsvh.py
@@ -1,7 +1,7 @@
 import numpy as np
 import scipy.stats as spst
-#import numpy.polynomial as nppoly
-from. import sabr
+import scipy.optimize as spop
+from . import sabr
 
 
 class Nsvh1(sabr.SabrABC):
@@ -91,3 +91,197 @@ def cdf(self, strike, spot, texp, cp=-1):
         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):
+        """
+        Variance, skewness, and ex-kurtosis
+
+        Args:
+            texp: time-to-expiry
+
+        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
+
+        c20 = ww + 1
+        c22 = ww - 1
+
+        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)
+
+        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)
+
+        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
+
+    def calibrate_vsk(self, var, skew, exkurt, texp=1, setval=False):
+        """
+        Calibrate parameters to the moments: variance, skewness, ex-kurtosis.
+
+        Args:
+            texp: time-to-expiry
+            var: variance
+            skew: skewness
+            exkurt: ex-kurtosis. should be > 0.
+
+        Returns: (sigma, vov, rho)
+
+        References:
+            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
+        beta2 = exkurt + 3
+
+        # min of w search
+        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)))
+
+        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
+
+        assert f_beta1(w_min) >= 0
+        #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_*
+
+        # if term is slightly negative, next line error in sqrt
+        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
+
+        if setval:
+            self.sigma = sig0
+            self.vov = vov
+            self.rho = rho
+
+        return sig0, vov, rho
+
+
+class NsvhMc(sabr.SabrABC):
+    """
+    Monte-Carlo model of Hyperbolic Normal Stochastic Volatility (NSVh) model
+
+    References:
+        Choi, J., Liu, C., & Seo, B. K. (2019). Hyperbolic normal stochastic volatility model.
+        Journal of Futures Markets, 39(2), 186–204. https://doi.org/10.1002/fut.21967
+    """
+
+    lam = 0
+    n_path = int(1e6)
+    rn_seed = None
+    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):
+        """
+        Args:
+            sigma: model volatility at t=0
+            vov: volatility of volatility
+            rho: correlation between price and volatility
+            lam: lambda. Norma SABR if 0, Johnson's SU if 1 (same as `Nsvh1`)
+            beta: elasticity parameter. should be 0 or None.
+            intr: interest rate (domestic interest rate)
+            divr: dividend/convenience yield (foreign interest rate)
+            is_fwd: if True, treat `spot` as forward price. False by default.
+        """
+        # Make sure beta = 0
+        if beta is not None and not np.isclose(beta, 0.0):
+            print(f'Ignoring beta = {beta}...')
+        self.lam = lam
+        super().__init__(sigma, vov, rho, beta=0, intr=intr, divr=divr, is_fwd=is_fwd)
+
+    def set_mc_params(self, n_path=1e6, rn_seed=None, antithetic=True):
+        self.n_path = int(n_path)
+        self.rn_seed = rn_seed
+        self.antithetic = antithetic
+        self.rn_seed = rn_seed
+        self.rng = np.random.default_rng(rn_seed)
+
+    def mc_vol_price(self, texp):
+        """
+        Simulate volatility and price pair
+
+        Args:
+            texp: time-to-expiry
+
+        Returns: (vol, price). vol: (n_path, ), price: (n_path, 2)
+
+        """
+        # forward path starting from zero
+        # returns both sigma and price
+
+        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 = np.stack([z_rn, -z_rn], axis=1).reshape((-1, 3))
+        z_rn[:, 2] += 0.5*(self.lam-1)*vol_std  # add shift
+
+        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))
+
+        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)
+        )
+        return path
+
+    def price(self, strike, spot, texp, cp=1):
+        """
+        Vanilla option price from MC simulation of NSVh model.
+
+        Args:
+            strike: strike price
+            spot: spot price
+            texp: time to expiry
+            cp: 1 or call, -1 for put option
+
+        Returns:
+            vanilla option price
+        """
+
+        fwd, df, _ = self._fwd_factor(spot, texp)
+        mc_path = self.mc_vol_price(texp)
+        strike_std = strike - fwd
+        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))])
+        if scalar_output:
+            price = price[0]
+        return df*price
diff --git a/pyfeng/opt_smile_abc.py b/pyfeng/opt_smile_abc.py
index cfaf234..4ccefc0 100644
--- a/pyfeng/opt_smile_abc.py
+++ b/pyfeng/opt_smile_abc.py
@@ -10,11 +10,13 @@ class OptSmileABC(opt.OptABC, abc.ABC):
     """
     Abstract class to model with volatility smile
     """
-    def _m_smile(self, model='bsm'):
+    def _m_smile(self, model='bsm', is_fwd=None):
+        if is_fwd is None:
+            is_fwd = self.is_fwd
         if model.lower() == 'bsm':
-            base_model = bsm.Bsm(None, intr=self.intr, divr=self.divr, is_fwd=self.is_fwd)
+            base_model = bsm.Bsm(None, intr=self.intr, divr=self.divr, is_fwd=is_fwd)
         elif model.lower() == 'norm':
-            base_model = norm.Norm(None, intr=self.intr, divr=self.divr, is_fwd=self.is_fwd)
+            base_model = norm.Norm(None, intr=self.intr, divr=self.divr, is_fwd=is_fwd)
         else:
             base_model = None
         return base_model
diff --git a/pyfeng/ousv.py b/pyfeng/ousv.py
new file mode 100644
index 0000000..c442f49
--- /dev/null
+++ b/pyfeng/ousv.py
@@ -0,0 +1,201 @@
+import numpy as np
+import scipy.integrate as scint
+from . import sv_abc as sv
+
+
+class OusvIft(sv.SvABC):
+    """
+    The implementation of Schobel & Zhu (1998)'s inverse FT pricing formula for European
+    options the Ornstein-Uhlenbeck driven stochastic volatility process.
+
+    References: Schöbel, R., & Zhu, J. (1999). Stochastic Volatility With an Ornstein–Uhlenbeck Process: An Extension. Review of Finance, 3(1), 23–46. https://doi.org/10.1023/A:1009803506170
+
+    Examples:
+        >>> import pyfeng as pf
+        >>> model = pf.OusvIft(0.2, mr=4, vov=0.1, rho=-0.7, intr=0.09531)
+        >>> model.price(100, 100, texp=np.array([1, 5, 10]))
+        array([13.21493, 40.79773, 62.76312])
+        >>> model = pf.OusvIft(0.25, mr=8, vov=0.3, rho=-0.6, intr=0.09531)
+        >>> model.price(np.array([90, 100, 110]), 100, texp=1)
+        array([21.41873, 15.16798, 10.17448])
+    """
+
+    def D_B_C(self, s1, s2, s3, texp):
+        # implement the formula for D(t,T), B(t,T), C(t,T) in paper appendix
+        mr, theta, vov = self.mr, self.theta, self.vov
+
+        gamma1 = np.sqrt(2 * vov ** 2 * s1 + mr ** 2)
+        gamma2 = (mr - 2 * vov ** 2 * s3) / gamma1
+        gamma3 = mr ** 2 * theta - s2 * vov ** 2
+        sinh = np.sinh(gamma1 * texp)
+        cosh = np.cosh(gamma1 * texp)
+        sincos = sinh + gamma2 * cosh
+        cossin = cosh + gamma2 * sinh
+        ktg3 = mr * theta * gamma1 - gamma2 * gamma3
+        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) \
+            + ktg3 * gamma3 / s2g3 * ((cosh - 1) / cossin)
+
+        return D, B, C
+
+    def f_1(self, phi, fwd, texp):
+        # implement the formula (12)
+        mr, theta, vov, rho = self.mr, self.theta, self.vov, self.rho
+        tmp = 1 + 1j * phi
+        s1 = 0.5 * tmp * (-tmp * (1 - rho ** 2) + (1 - 2 * mr * rho / vov))
+        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)
+        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)
+
+    def f_2(self, phi, fwd, texp):
+        # implement the formula (13)
+        mr, theta, vov, rho = self.mr, self.theta, self.vov, self.rho
+
+        s1 = 0.5 * (phi ** 2 * (1 - rho ** 2) + 1j * phi * (1 - 2 * mr * rho / vov))
+        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)
+        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)
+
+    def price(self, strike, spot, texp, cp=1):
+        # implement the formula (14) and (15)
+        fwd, df, _ = self._fwd_factor(spot, texp)
+
+        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))
+        if len(price) == 1:
+            price = price[0]
+
+        return df * price
+
+
+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
+        """
+
+    def _bm_incr(self, tobs, cum=False, n_path=None):
+        """
+        Calculate incremental Brownian Motions
+
+        Args:
+            tobs: observation times (array). 0 is not included.
+            cum: return cumulative values if True
+            n_path: number of paths. If None (default), use the stored one.
+
+        Returns:
+            price path (time, path)
+        """
+        # dt = np.diff(np.atleast_1d(tobs), prepend=0)
+        n_dt = len(tobs)
+
+        tobs_lag = tobs[:-1]
+        tobs_lag = np.insert(tobs_lag, 0, 0)
+        bm_var = np.exp(2 * self.mr * tobs) - np.exp(2 * self.mr * tobs_lag)
+        n_path = n_path or self.n_path
+
+        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 = 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])
+            bm_incr = self.rng.normal(size=(n_path, n_dt)).T * np.sqrt(bm_var[:, None])
+
+        if cum:
+            np.cumsum(bm_incr, axis=0, out=bm_incr)
+
+        return bm_incr
+
+    def vol_paths(self, tobs):
+        """
+        sigma_t = np.exp(-mr * tobs) * (sigma0 - theta * mr + vov / np.sqrt(2 * mr) * bm) + theta * mr
+        Args:
+            tobs: observation time (array)
+            mr: coefficient of dt
+            theta: the long term average
+            mu: rn-derivative
+
+        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)
+        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
+
+            Args:
+                theta: the long term average
+                mr: coefficient of dt
+                texp: time-to-expiry
+
+            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
+
+        vol_cond = rhoc * np.sqrt(int_var / texp)
+
+        return fwd_cond, vol_cond
+
+    def price(self, strike, spot, texp, cp=1):
+        """
+        Calculate option price based on BSM
+        Args:
+            strike: strike price
+            spot: spot price
+            texp: time to maturity
+            cp: cp=1 if call option else put option
+
+        Returns: price
+        """
+        price = []
+        texp = [texp] if isinstance(texp, (int, float)) else texp
+        for t in texp:
+            kk = strike / spot
+            kk = np.atleast_1d(kk)
+
+            fwd_cond, vol_cond = self.cond_fwd_vol(t)
+
+            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)
+            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 11c0275..7f14edc 100644
--- a/pyfeng/sabr.py
+++ b/pyfeng/sabr.py
@@ -51,7 +51,7 @@ def _variables(self, fwd, texp):
         vovn = self.vov * np.sqrt(np.maximum(texp, 1e-64))
         return alpha, betac, rhoc, rho2, vovn
 
-    def _m_base(self, vol):
+    def _m_base(self, vol, is_fwd=None):
         """
         Create base model based on _base_beta value: `Norm` for 0, Cev for (0,1), and `Bsm` for 1
         If `_base_beta` is None, use `base` instead.
@@ -62,13 +62,14 @@ def _m_base(self, vol):
         Returns: model
         """
         base_beta = self._base_beta or self.beta
-
+        if is_fwd is None:
+            is_fwd = self.is_fwd
         if np.isclose(base_beta, 1):
-            return bsm.Bsm(vol, intr=self.intr, divr=self.divr)
+            return bsm.Bsm(vol, intr=self.intr, divr=self.divr, is_fwd=is_fwd)
         elif np.isclose(base_beta, 0):
-            return norm.Norm(vol, intr=self.intr, divr=self.divr)
+            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)
+            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:
diff --git a/pyfeng/sabr_mc.py b/pyfeng/sabr_mc.py
new file mode 100644
index 0000000..f393d08
--- /dev/null
+++ b/pyfeng/sabr_mc.py
@@ -0,0 +1,92 @@
+import numpy as np
+import scipy.integrate as scint
+from . import sabr
+from . import sv_abc as sv
+from . import cev
+
+
+class SabrCondMc(sabr.SabrABC, sv.CondMcBsmABC):
+    """
+    Conditional MC for SABR model (beta=0,1 or rho=0) with conditional Monte-Carlo simulation
+
+    """
+
+    def vol_paths(self, tobs, mu=0):
+        """
+        exp(vov_std B_s - 0.5*vov_std^2 * s)  where s = 0, ..., 1, vov_std = vov*sqrt(T)
+
+        Args:
+            tobs: observation time (array)
+            mu: rn-derivative
+
+        Returns: volatility path (time, path) including the value at t=0
+        """
+
+        texp = tobs[-1]
+        tobs01 = tobs / texp  # normalized time: 0<s<1
+        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_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
+
+        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)
+        if np.isclose(self.beta, 0):
+            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))
+
+        return fwd_cond, vol_cond, log_rn_deriv
+
+    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)
+            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)
+        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)
+            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)
+
+        return price
+
+    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))
+        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}')
+
+        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)
+        if log:
+            log_mass_grid = base_model.mass_zero(spot, texp, log=True) + log_rn_deriv
+            log_mass_max = np.amax(log_mass_grid)
+            log_mass_grid -= log_mass_max
+            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 = np.mean(mass_grid)
+            return mass
diff --git a/pyfeng/sv32_mc.py b/pyfeng/sv32_mc.py
new file mode 100644
index 0000000..fb83c1e
--- /dev/null
+++ b/pyfeng/sv32_mc.py
@@ -0,0 +1,579 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon, May 3, 2021
+Last modified on Fri, May 7, 2021
+Conditional MC for 3/2 model based on QE discretization scheme by Andersen(2008)
+@author: Xueyang & Xiaoyin
+"""
+import numpy as np
+import pyfeng as pf
+import scipy.stats as spst
+import scipy.integrate as spint
+import scipy.optimize as spop
+from scipy.misc import derivative
+from mpmath import besseli
+from .bsm import Bsm
+from .norm import Norm
+
+
+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.
+    Volatility (variance) of the price is assumed to follow 3/2 model by Heston (1997) and Lewis (2000).
+
+    Example:
+        >>> import numpy as np
+        >>> import pyfeng as pf
+        >>> strike = [100.0, 140.0, 70.0]
+        >>> forward = 100
+        >>> delta = [1, 1/2, 1/4, 1/8, 1/16, 1/32]
+        >>> vov, kappa, rho, texp, theta, sigma = [1, 0.5, -0.9, 10, 0.04, np.sqrt(0.04)]
+        >>> sv32_cmc_qe = pf.Sv32CondMcQE(vov=vov, kappa=kappa, rho=rho, theta=theta)
+        >>> price_cmc = np.zeros([len(delta), len(strike)])
+        >>> for d in range(len(delta)):
+        >>>     price_cmc[d, :] = sv32_cmc_qe.price(strike, forward, texp, sigma=sigma, delta=delta[d], path=1e5, seed=123456)
+        >>> price_cmc
+        array([[22.95314785, 10.44093783, 38.98364955],
+               [23.2425543 , 10.67162543, 39.26731165],
+               [23.20965635, 10.64143576, 39.21865023],
+               [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:
+            vov: volatility of variance, strictly positive
+            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
+        self.theta = theta
+
+        self.psi_points = None  # for TG scheme only
+        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):
+        '''
+        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.
+
+        Args:
+            strike: strike price, in vector form
+            spot: spot (or forward)
+            texp: time to expiry
+            sigma: initial volatility
+            delta: length of each time step
+            intr: interest rate (domestic interest rate)
+            divr: dividend/convenience yield (foreign interest rate)
+            psi_c: critical value for psi, lying in [1, 2]
+            path: number of vol paths generated
+            scheme: discretization scheme for vt, {'QE', 'TG', 'Euler', 'Milstein', 'KJ'}
+            seed: random seed for rv generation
+
+        Return:
+            BSM price vector for all strikes
+        '''
+        self.sigma = sigma
+        self.bsm_model = pf.Bsm(self.sigma, intr=intr, divr=divr)
+        self.delta = delta
+        self.path = int(path)
+        self.step = int(texp / self.delta)
+
+        # xt = 1 / vt
+        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':
+            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)
+                psi = s2 / m ** 2
+
+                # compute xt(i+1) given psi
+                below = np.where(psi <= psi_c)[0]
+                ins = 2 * psi[below] ** -1
+                b2 = ins - 1 + np.sqrt(ins * (ins - 1))
+                b = np.sqrt(b2)
+                a = m[below] / (1 + b2)
+                z = spst.norm.ppf(u[below, i])
+                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]))
+                    else:
+                        xt[above[k], i+1] = 0
+
+        elif scheme == 'TG':
+            if np.all(self.rx_results) == None:
+                self.psi_points, self.rx_results = self.prepare_rx()
+
+            expo = np.exp(-self.kappa * self.delta)
+            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)
+                psi = s2 / m ** 2
+
+                rx = np.array([self.find_rx(j) for j in psi])
+
+                z = np.random.normal(size=(self.path, self.step))
+                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))
+
+                xt[:, i+1] = np.fmax(mu_v[:, i] + sigma_v[:, i] * z[:, i], 0)
+
+        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]
+
+        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':
+            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)
+
+        # compute integral of vt, equivalent spot and vol
+        vt = 1 / xt
+        below_0 = np.where(vt < 0)
+        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)
+        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))
+
+        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
+        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")
+            return self.rx_results[-1]
+        elif self.rx_results[self.psi_points <= psi].size == 0:
+            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
+
+
+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.'
+        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.T = T
+        self.r = r
+        self.sigma_0 = sigma_0
+        self.beta = beta
+        self.rho = rho
+        self.theta = theta
+        self.kappa = kappa
+        self.vov = vov
+        self.path_num = path_num
+        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,
+        where U_T is the integral of sigma_t^2 from t=0 to t=T.
+        Third, we use 100 (by default, we can change this parameter when we construct the Sv32McAe object) pairs of sigma_T and U_T
+        to calculate the corresponding E_F_T, expected forward price, and sigma_N or sigma_BS
+        Finally, we will use these pairs of (E_F_T, sigma_N or sigma_BS) to calculate corresponding option prices.
+
+        Returns:
+        A tuple with two elements.
+        The first element is a np.ndarray, which includes all the means of prices under every strike price.
+        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
+
+    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
+
+    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)
+            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)
+            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)
+
+    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.')
+
+    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)
+        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)
+        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
+            M2: np.ndarray, second moment
+
+        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))
+
+    @staticmethod
+    def calTwoMoments(chfs):
+        '''
+        Calculate two moments from the characteristic function.
+        Args:
+            chfs: function with only one arg_in
+
+        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
+        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))
+        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,
+        as the RNG from conditional distribution have oscillating output,
+        which is not good for the standard deviation of option pricing in AEMC approach.
+        This is the innovation on insider layer simulation.
+        Third, we use 100 (by default, we can change this parameter when we construct the Sv32McAe object) pairs of sigma_T and U_T
+        to calculate the corresponding E_F_T, expected forward price, and sigma_N or sigma_BS
+        Finally, we will use these pairs of (E_F_T, sigma_N or sigma_BS) to calculate corresponding option prices.
+
+        Returns:
+        A tuple with two elements.
+        The first element is a np.ndarray, which includes all the means of prices under every strike price.
+        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
+        ))
+        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)
+        else:
+            outside_term = self.rho ** 2 * M1 / 2
+            self.E_F_T = self.S0 * np.exp(internal_term - outside_term)
+            self.sigma_BS = np.sqrt((1 - self.rho ** 2) * M1 / self.T)
+            self.sigma_BS = self.sigma_BS.reshape(-1, 1)
+        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)
+        return M1
+
+    @staticmethod
+    def calOneMoment(chfs):
+        M1 = derivative(chfs, x0=0, dx=0.00001, n=1)
+        return M1
+
+    # def optionPrice_version3(self):
+    #     self.simulate_sigma_2_T()
+    #     self.simulate_U_T_version3()
+    #     self.calForwardAndVolatility_version1()
+    #     return self.calOutput()
+
+    # def simulate_U_T_version3(self):
+    #     chfs = self.charFuncs()
+    #     M1, M2 = self.calTwoMoments(chfs)
+    #     mu, scale = self.calInverseGaussianParas(M1, M2)
+    #     self.U_T = st.lognorm.rvs(mu, 0, scale)
+
+    # @staticmethod
+    # def calInverseGaussianParas(M1, M2):
+    #     M1 = np.array(np.abs(M1).tolist(), dtype=float)
+    #     M2 = np.array(np.abs(M2).tolist(), dtype=float)
+    #     M2 = np.sqrt(np.log(np.array(M2/M1**2, dtype=float)))
+    #     M2[np.isnan(M2)] = 0
+    #     scale = np.power(M1,3)/(M2-M1**2)
+    #     mu = M1/scale
+    #     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
+        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):
+        '''
+        c function in the paper: Exact Simulation of the 3/2 Model. The formula is in page-6
+        Args:
+            t: float
+                t is the independent variable of c function
+
+        Returns:
+            float
+                c(t)
+        '''
+        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):
+        '''
+        import from pyfeng package, to ues BSM model pricing formula
+        Args:
+            S0: Stock initial price
+            Ks: Array of strike prices
+            r: interest rate
+            T: time to maturity of options
+            sigma: volatility of stock
+            cp: call or put
+
+        Returns:
+            prices: np,ndarray
+        '''
+        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
+            Ks: Array of strike prices
+            r: interest rate
+            T: time to maturity of options
+            sigma: volatility of stock
+            cp: call or put
+
+        Returns:
+            prices: np,ndarray
+        '''
+        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
+        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
+        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()):
+                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'):
+        '''
+        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'
+        output = []
+        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:
+            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)
diff --git a/pyfeng/sv_abc.py b/pyfeng/sv_abc.py
new file mode 100644
index 0000000..11de2ba
--- /dev/null
+++ b/pyfeng/sv_abc.py
@@ -0,0 +1,155 @@
+import numpy as np
+import abc
+from . import bsm
+from . import opt_smile_abc as smile
+
+
+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):
+        """
+        Args:
+            sigma: model volatility at t=0. variance = sigma**2
+            vov: volatility of volatility
+            rho: correlation between price and volatility
+            mr: mean-reversion speed (kappa)
+            theta: long-term mean of volatility. For variance process, use theta**2. If None, same as sigma
+            intr: interest rate (domestic interest rate)
+            divr: dividend/convenience yield (foreign interest rate)
+            is_fwd: if True, treat `spot` as forward price. False by default.
+        """
+
+        super().__init__(sigma, intr=intr, divr=divr, is_fwd=is_fwd)
+
+        self.vov = vov
+        self.rho = rho
+        self.mr = mr
+        self.theta = sigma if theta is None else theta
+
+    def params_kw(self):
+        params1 = super().params_kw()
+        params2 = {"vov": self.vov, "rho": self.rho, "mr": self.mr, "sig_inf": self.theta}
+        return {**params1, **params2}
+
+
+class CondMcBsmABC(smile.OptSmileABC, abc.ABC):
+    """
+    Abstract Class for conditional Monte-Carlo method for BSM-based stochastic volatility models
+    """
+
+    dt = 0.05
+    n_path = 10000
+    rn_seed = None
+    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
+
+        Args:
+            n_path: number of paths
+            dt: time step for Euler/Milstein steps
+            rn_seed: random number seed
+            antithetic: antithetic
+        """
+        self.n_path = int(n_path)
+        self.dt = dt
+        self.rn_seed = rn_seed
+        self.antithetic = antithetic
+        self.rn_seed = rn_seed
+        self.rng = np.random.default_rng(rn_seed)
+
+    def base_model(self, vol):
+        return bsm.Bsm(vol, intr=self.intr, divr=self.divr, is_fwd=self.is_fwd)
+
+    def tobs(self, texp):
+        """
+        Return array of observation time in even size. 0 is not included.
+
+        Args:
+            texp: time-to-expiry
+
+        Returns:
+            array of observation time
+        """
+        n_steps = (texp//(2*self.dt)+1)*2
+        tobs = np.arange(1, n_steps + 0.1) / n_steps * texp
+        return tobs
+
+    def _bm_incr(self, tobs, cum=False, n_path=None):
+        """
+        Calculate incremental Brownian Motions
+
+        Args:
+            tobs: observation times (array). 0 is not included.
+            cum: return cumulative values if True
+            n_path: number of paths. If None (default), use the stored one.
+
+        Returns:
+            price path (time, path)
+        """
+        dt = np.diff(np.atleast_1d(tobs), prepend=0)
+        n_dt = len(dt)
+
+        n_path = n_path or self.n_path
+
+        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 = 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])
+
+        if cum:
+            np.cumsum(bm_incr, axis=0, out=bm_incr)
+
+        return bm_incr
+
+    @abc.abstractmethod
+    def vol_paths(self, tobs):
+        """
+        Volatility or variance paths at 0 and tobs.
+
+        Args:
+            tobs: observation time (array)
+
+        Returns:
+            2d array of (time, path)
+        """
+
+        return np.ones(size=(len(tobs), self.n_path))
+
+    @abc.abstractmethod
+    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
+
+        Args:
+            texp: time-to-expiry
+
+        Returns: (forward, volatility)
+        """
+
+        return np.ones(self.n_path), self.sigma * np.ones(self.n_path)
+
+    def price(self, strike, spot, texp, cp=1):
+
+        kk = strike / spot
+        scalar_output = np.isscalar(kk)
+        kk = np.atleast_1d(kk)
+
+        fwd_cond, vol_cond = self.cond_fwd_vol(texp)
+
+        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)
+
+        return price[0] if scalar_output else price


Test patch:
diff --git a/tests/test_sabr.py b/tests/test_sabr.py
index 184306d..7d1d922 100644
--- a/tests/test_sabr.py
+++ b/tests/test_sabr.py
@@ -61,6 +61,13 @@ def test_UnCorrChoiWu2021(self):
         np.testing.assert_almost_equal(mass, mass2)
         np.testing.assert_almost_equal(p, p2)
 
+    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)
+
 
 if __name__ == '__main__':
     print(f'Pyfeng loaded from {pf.__path__}')