Repository: PyFE/PyFENG (language: python)

Issue:
Refactoring model names
* TimeStep -> TimeDisc
* pfex -> pf


Discussion / hints:

Fix patch:
diff --git a/.idea/PyFENG.iml b/.idea/PyFENG.iml
index 06b0085..0aea259 100644
--- a/.idea/PyFENG.iml
+++ b/.idea/PyFENG.iml
@@ -2,7 +2,7 @@
 <module type="PYTHON_MODULE" version="4">
   <component name="NewModuleRootManager">
     <content url="file://$MODULE_DIR$" />
-    <orderEntry type="inheritedJdk" />
+    <orderEntry type="jdk" jdkName="Python 3.9" jdkType="Python SDK" />
     <orderEntry type="sourceFolder" forTests="false" />
   </component>
   <component name="PackageRequirementsSettings">
diff --git a/.idea/misc.xml b/.idea/misc.xml
index d1e22ec..d56657a 100644
--- a/.idea/misc.xml
+++ b/.idea/misc.xml
@@ -1,4 +1,4 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <project version="4">
-  <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.8" project-jdk-type="Python SDK" />
+  <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.9" project-jdk-type="Python SDK" />
 </project>
\ No newline at end of file
diff --git a/pyfeng/__init__.py b/pyfeng/__init__.py
index d4beebd..a30a734 100644
--- a/pyfeng/__init__.py
+++ b/pyfeng/__init__.py
@@ -1,9 +1,23 @@
-from .norm import (
-    Norm,
-)  # the order is sensitive because of `price_barrier` method. Put it before .bsm
+# the order is sensitive because of `price_barrier` method. Put it before .bsm
+from .norm import Norm
+
 from .bsm import Bsm, BsmDisp
-from .cev import Cev
+from .cev import Cev, CevMc
 from .gamma import InvGam, InvGauss
+<<<<<<< HEAD
+
+# FFT related models
+from .sv_fft import HestonFft, BsmFft, OusvFft, VarGammaFft, ExpNigFft, Sv32Fft
+
+# SABR/NSVh related models
+from .sabr import SabrHagan2002, SabrNormVolApprox, SabrLorig2017, SabrChoiWu2021H, SabrChoiWu2021P
+from .sabr_int import SabrUncorrChoiWu2021
+from .sabr_mc import SabrMcTimeDisc
+from .nsvh import Nsvh1, NsvhMc, NsvhGaussQuad
+
+# Other SV models
+from .garch import GarchMcTimeDisc, GarchUncorrBaroneAdesi2004
+=======
 from .sv_fft import HestonFft,RoughHestonFft, BsmFft, OusvFft, VarGammaFft, ExpNigFft, Sv32Fft
 from .sabr import (
     SabrHagan2002,
@@ -13,15 +27,15 @@
     SabrChoiWu2021P,
 )
 from .garch import GarchMcTimeStep, GarchUncorrBaroneAdesi2004
+>>>>>>> origin/rheston_sunan
 from .heston import HestonUncorrBallRoma1994
 from .heston_mc import (
     HestonMcAndersen2008, HestonMcGlassermanKim2011, HestonMcTseWan2013,
     HestonMcChoiKwok2023PoisGe, HestonMcChoiKwok2023PoisTd
 )
-from .ousv import OusvUncorrBallRoma1994
-from .sabr_int import SabrUncorrChoiWu2021
-from .sabr_mc import SabrMcTimeDisc
-from .nsvh import Nsvh1, NsvhMc, NsvhGaussQuad
+from .ousv import OusvUncorrBallRoma1994, OusvMcTimeDisc, OusvMcChoi2023KL
+from .svi import Svi
+
 from .multiasset import (
     BsmSpreadKirk,
     BsmSpreadBjerksund2014,
diff --git a/pyfeng/bsm.py b/pyfeng/bsm.py
index 2c1b1c2..03d99df 100644
--- a/pyfeng/bsm.py
+++ b/pyfeng/bsm.py
@@ -406,7 +406,7 @@ def moments_vsk(self, texp=1):
 
 class BsmDisp(smile.OptSmileABC, Bsm):
     """
-    Displaced Black-Scholes-Merton model for option pricing. Displace price,
+    Displaced Black-Scholes-Merton model for option pricing. Displaced price,
 
         D(F_t) = beta*F_t + (1-beta)*A
 
@@ -449,7 +449,7 @@ def __init__(self, sigma, beta=1, pivot=0, *args, **kwargs):
 
     @property
     def sigma(self):
-        return self.sigma_disp*self.beta
+        return self.sigma_disp * self.beta
 
     @sigma.setter
     def sigma(self, sigma):
@@ -463,7 +463,7 @@ def disp_spot(self, spot):
             spot: spot (or forward) price
 
         Returns:
-            Displaces spot
+            Displaced spot
         """
         return self.beta*spot + (1 - self.beta)*self.pivot
 
@@ -476,7 +476,7 @@ def disp_strike(self, strike, texp):
             texp: time to expiry
 
         Returns:
-            Displace strike
+            Displaced strike
         """
         return self.beta*strike + (1 - self.beta)*self.forward(self.pivot, texp)
 
@@ -540,13 +540,7 @@ def vol_smile(self, strike, spot, texp, cp=1, model="bsm"):
             kkd = self.disp_strike(strike, texp)/fwdd
             lnkd = np.log(kkd)
             # self.sigma actually means self.beta * self._sigma
-            vol = (
-                    self.sigma_disp
-                    *fwdd
-                    *np.sqrt(kkd)
-                    *(1 + lnkd**2/24)
-                    /(1 + self.sigma**2*texp/24)
-            )
+            vol = self.sigma_disp * fwdd * np.sqrt(kkd) * (1 + lnkd**2/24) / (1 + self.sigma**2*texp/24)
         elif model.lower() == "bsm-approx":
             fwd = self.forward(spot, texp)
             kk = strike/fwd
@@ -558,12 +552,7 @@ def vol_smile(self, strike, spot, texp, cp=1, model="bsm"):
 
             # self.sigma actually means self.beta * self.sigma_disp
             vol = self.sigma_disp*(fwdd/fwd)*np.sqrt(kkd/kk)
-            vol *= (
-                    (1 + lnkd**2/24)
-                    /(1 + lnk**2/24)
-                    *(1 + vol**2*texp/24)
-                    /(1 + self.sigma**2*texp/24)
-            )
+            vol *= (1 + lnkd**2/24) / (1 + lnk**2/24) * (1 + vol**2*texp/24) / (1 + self.sigma**2*texp/24)
         else:
             vol = super().vol_smile(strike, spot, texp, model=model, cp=cp)
 
diff --git a/pyfeng/cev.py b/pyfeng/cev.py
index 67b7187..774bd3c 100644
--- a/pyfeng/cev.py
+++ b/pyfeng/cev.py
@@ -1,26 +1,16 @@
+import abc
+
 import scipy.stats as spst
 import scipy.special as spsp
 import numpy as np
 from . import opt_abc as opt
 from . import opt_smile_abc as smile
+from . import sv_abc as sv
+from .util import avg_exp
 
 
-class Cev(opt.OptAnalyticABC, smile.OptSmileABC, smile.MassZeroABC):
-    """
-    Constant Elasticity of Variance (CEV) model.
-
-    Underlying price is assumed to follow CEV process:
-    dS_t = (r - q) S_t dt + sigma S_t^beta dW_t, where dW_t is a standard Brownian motion.
-
-    Examples:
-        >>> import numpy as np
-        >>> import pyfeng as pf
-        >>> m = pf.Cev(sigma=0.2, beta=0.5, intr=0.05, divr=0.1)
-        >>> m.price(np.arange(80, 121, 10), 100, 1.2)
-        array([16.11757214, 10.00786871,  5.64880408,  2.89028476,  1.34128656])
-    """
-
-    sigma = None
+class CevAbc(smile.OptSmileABC, abc.ABC):
+    model_type = "Cev"
     beta = 0.5
 
     def __init__(self, sigma, beta=0.5, intr=0.0, divr=0.0, is_fwd=False):
@@ -40,13 +30,29 @@ def params_kw(self):
         extra = {"beta": self.beta}
         return {**params, **extra}  # Py 3.9, params | extra
 
+
+class Cev(opt.OptAnalyticABC, CevAbc, smile.MassZeroABC):
+    """
+    Constant Elasticity of Variance (CEV) model.
+
+    Underlying price is assumed to follow CEV process:
+    dS_t = (r - q) S_t dt + sigma S_t^beta dW_t, where dW_t is a standard Brownian motion.
+
+    Examples:
+        >>> import numpy as np
+        >>> import pyfeng as pf
+        >>> m = pf.Cev(sigma=2, beta=0.5, intr=0.05, divr=0.1)
+        >>> m.price(np.arange(80, 121, 10), spot=100, texp=1.2)
+        array([16.0290539 ,  9.89374055,  5.53356595,  2.79581565,  1.27662838])
+    """
+
     def mass_zero(self, spot, texp, log=False):
         fwd = self.forward(spot, texp)
 
         betac = 1.0 - self.beta
-        a = 0.5/betac
-        sigma_std = np.maximum(self.sigma/np.power(fwd, betac)*np.sqrt(texp), np.finfo(float).eps)
-        x = 0.5/np.square(betac*sigma_std)
+        a = 0.5 / betac  # shape parameter of gamma
+        var = (betac * self.sigma)**2 * texp * avg_exp(2*(self.intr-self.divr)*betac*texp)
+        x = 0.5 * np.power(fwd, 2*betac) / var
 
         if log:
             log_mass = (a - 1)*np.log(x) - x - np.log(spsp.gamma(a))
@@ -70,7 +76,9 @@ def mass_zero_t0(self, spot, texp):
         fwd = self.forward(spot, texp)
         betac = 1.0 - self.beta
         alpha = self.sigma/np.power(fwd, betac)
-        t0 = 0.5/(betac*alpha)**2
+        var_t0 = (betac*alpha)**2 * avg_exp(2*(self.intr-self.divr)*betac*texp)
+        t0 = 0.5/var_t0
+
         return t0
 
     @staticmethod
@@ -93,16 +101,16 @@ def price_formula(strike, spot, texp, sigma=None, cp=1, beta=0.5, intr=0.0, divr
             Vanilla option price
         """
 
-        disc_fac = np.exp(-texp*intr)
-        fwd = spot*(1.0 if is_fwd else np.exp(-texp*divr)/disc_fac)
+        df = np.exp(-texp*intr)
+        fwd = spot * (1.0 if is_fwd else np.exp(-texp*divr)/df)
 
         betac = 1.0 - beta
         betac_inv = 1.0/betac
         alpha = sigma/np.power(fwd, betac)
-        sigma_std = np.maximum(alpha*np.sqrt(texp), np.finfo(float).eps)
-        kk = strike/fwd
-        x = 1.0/np.square(betac*sigma_std)
-        y = np.power(kk, 2*betac)*x
+        var = np.maximum(alpha**2 * texp, np.finfo(float).eps) * betac**2 * avg_exp(2*(intr-divr)*betac*texp)
+
+        xx = 1.0 / var
+        yy = np.power(strike/fwd, 2*betac) * xx
 
         # Need to clean up the case beta > 0
         if beta > 1.0:
@@ -114,25 +122,26 @@ def price_formula(strike, spot, texp, sigma=None, cp=1, beta=0.5, intr=0.0, divr
         # Computing call and put is a bit of computtion waste, but do this for vectorization.
         price = np.where(
             cp > 0,
-            fwd*ncx2_sf(y, 2 + betac_inv, x) - strike*ncx2_cdf(x, betac_inv, y),
-            strike*ncx2_sf(x, betac_inv, y) - fwd*ncx2_cdf(y, 2 + betac_inv, x),
+            fwd*ncx2_sf(yy, 2 + betac_inv, xx) - strike*ncx2_cdf(xx, betac_inv, yy),
+            strike*ncx2_sf(xx, betac_inv, yy) - fwd*ncx2_cdf(yy, 2 + betac_inv, xx),
         )
-        return disc_fac*price
+        return df*price
 
     def delta(self, strike, spot, texp, cp=1):
         fwd, df, divf = self._fwd_factor(spot, texp)
-        betac_inv = 1/(1 - self.beta)
+        betac = 1.0 - self.beta
+        betac_inv = 1.0/betac
 
-        k_star = 1.0/np.square(self.sigma/betac_inv)/texp
-        x = k_star*np.power(fwd, 2/betac_inv)
-        y = k_star*np.power(strike, 2/betac_inv)
+        var = (betac * self.sigma)**2 * texp * avg_exp(2*(self.intr-self.divr)*betac*texp)
+        xx = np.power(fwd, 2*betac) / var
+        yy = np.power(strike, 2*betac) / var
 
         if self.beta < 1.0:
-            delta = 0.5*(cp - 1) + spst.ncx2.sf(y, 2 + betac_inv, x)\
-                    + 2*x/betac_inv*(spst.ncx2.pdf(y, 4 + betac_inv, x) - strike/fwd*spst.ncx2.pdf(x, betac_inv, y))
+            delta = 0.5*(cp - 1) + spst.ncx2.sf(yy, 2 + betac_inv, xx)\
+                    + 2*xx*betac*(spst.ncx2.pdf(yy, 4 + betac_inv, xx) - (strike/fwd)*spst.ncx2.pdf(xx, betac_inv, yy))
         else:
-            delta = 0.5*(cp - 1) + spst.ncx2.sf(x, -betac_inv, y)\
-                    - 2*x/betac_inv*(spst.ncx2.pdf(x, -betac_inv, y) - strike/fwd*spst.ncx2.pdf(y, 4 - betac_inv, x))
+            delta = 0.5*(cp - 1) + spst.ncx2.sf(xx, -betac_inv, yy)\
+                    - 2*xx*betac*(spst.ncx2.pdf(xx, -betac_inv, yy) - (strike/fwd)*spst.ncx2.pdf(yy, 4 - betac_inv, xx))
 
         delta *= df if self.is_fwd else divf
         return delta
@@ -142,33 +151,35 @@ def cdf(self, strike, spot, texp, cp=1):
 
         betac = 1.0 - self.beta
         betac_inv = 1.0/betac
+
         alpha = self.sigma/np.power(fwd, betac)
-        sigma_std = np.maximum(alpha*np.sqrt(texp), np.finfo(float).eps)
-        kk = strike/fwd
-        x = 1.0/np.square(betac*sigma_std)
-        y = np.power(kk, 2*betac)*x
+        var = (betac * alpha)**2 * texp * avg_exp(2*(self.intr-self.divr)*betac*texp)
+        xx = 1.0/var
+        yy = np.power(strike/fwd, 2*betac) * xx
 
-        cdf = np.where(cp > 0, spst.ncx2.cdf(x, betac_inv, y), spst.ncx2.sf(x, betac_inv, y))
+        cdf = np.where(cp > 0, spst.ncx2.cdf(xx, betac_inv, yy), spst.ncx2.sf(xx, betac_inv, yy))
         return cdf
 
     def gamma(self, strike, spot, texp, cp=1):
         fwd, df, divf = self._fwd_factor(spot, texp)
-        betac_inv = 1/(1 - self.beta)
 
-        k_star = 1.0/np.square(self.sigma/betac_inv)/texp
-        x = k_star*np.power(fwd, 2/betac_inv)
-        y = k_star*np.power(strike, 2/betac_inv)
+        betac = 1.0 - self.beta
+        betac_inv = 1.0/betac
+
+        var = (betac * self.sigma)**2 * texp * avg_exp(2*(self.intr-self.divr)*betac*texp)
+        xx = np.power(fwd, 2*betac) / var
+        yy = np.power(strike, 2*betac) / var
 
         if self.beta < 1.0:
-            gamma = (2 + betac_inv - x)*spst.ncx2.pdf(y, 4 + betac_inv, x) \
-                    + x*spst.ncx2.pdf(y, 6 + betac_inv, x) \
-                    + strike/fwd*(x*spst.ncx2.pdf(x, betac_inv, y) - y*spst.ncx2.pdf(x, 2 + betac_inv, y))
+            gamma = (2 + betac_inv - xx)*spst.ncx2.pdf(yy, 4 + betac_inv, xx) \
+                    + xx*spst.ncx2.pdf(yy, 6 + betac_inv, xx) \
+                    + strike/fwd*(xx*spst.ncx2.pdf(xx, betac_inv, yy) - yy*spst.ncx2.pdf(xx, 2 + betac_inv, yy))
         else:
-            gamma = (x*spst.ncx2.pdf(x, -betac_inv, y) - y*spst.ncx2.pdf(x, 2 - betac_inv, y)) \
-                    + strike/fwd*((2 - betac_inv - x)*spst.ncx2.pdf(y, 4 - betac_inv, x)
-                                  + x*spst.ncx2.pdf(y, 6 - betac_inv, x))
+            gamma = (xx*spst.ncx2.pdf(xx, -betac_inv, yy) - yy*spst.ncx2.pdf(xx, 2 - betac_inv, yy)) \
+                    + strike/fwd*((2 - betac_inv - xx)*spst.ncx2.pdf(yy, 4 - betac_inv, xx)
+                                  + xx*spst.ncx2.pdf(yy, 6 - betac_inv, xx))
 
-        gamma *= 2*(divf/betac_inv)**2/df*x/fwd
+        gamma *= 2*(divf*betac)**2/df*xx/fwd
 
         if self.is_fwd:
             gamma *= (df/divf)**2
@@ -177,23 +188,134 @@ def gamma(self, strike, spot, texp, cp=1):
 
     def vega(self, strike, spot, texp, cp=1):
         fwd, df, divf = self._fwd_factor(spot, texp)
-        spot = fwd*df/divf
 
-        betac_inv = 1/(1 - self.beta)
+        betac = 1.0 - self.beta
+        betac_inv = 1.0/betac
 
-        k_star = 1.0/np.square(self.sigma/betac_inv)/texp
-        x = k_star*np.power(fwd, 2/betac_inv)
-        y = k_star*np.power(strike, 2/betac_inv)
+        var = (betac * self.sigma)**2 * texp * avg_exp(2*(self.intr-self.divr)*betac*texp)
+        xx = np.power(fwd, 2*betac) / var
+        yy = np.power(strike, 2*betac) / var
 
         if self.beta < 1.0:
-            vega = -fwd*spst.ncx2.pdf(y, 4 + betac_inv, x) + strike*spst.ncx2.pdf(x, betac_inv, y)
+            vega = -fwd*spst.ncx2.pdf(yy, 4 + betac_inv, xx) + strike*spst.ncx2.pdf(xx, betac_inv, yy)
         else:
-            vega = fwd*spst.ncx2.pdf(x, -betac_inv, y) - strike*spst.ncx2.pdf(y, 4 - betac_inv, x)
+            vega = fwd*spst.ncx2.pdf(xx, -betac_inv, yy) - strike*spst.ncx2.pdf(yy, 4 - betac_inv, xx)
 
-        sigma = self.sigma*np.power(spot, self.beta - 1)
-        vega *= df*2*x/sigma
+        sigma = self.sigma*np.power(spot, -betac)
+        vega *= df*2*xx/sigma
         return vega
 
     def theta(self, strike, spot, texp, cp=1):
         ### Need to implement this
         return self.theta_numeric(strike, spot, texp, cp=cp)
+
+
+class CevMc(CevAbc):
+    """
+    Constant Elasticity of Variance (CEV) model with exact Monte-Carlo method (Kang, 2014)
+
+    Underlying price is assumed to follow CEV process:
+    dS_t = (r - q) S_t dt + sigma S_t^beta dW_t, where dW_t is a standard Brownian motion.
+
+    Examples:
+        >>> import numpy as np
+        >>> import pyfeng as pf
+        >>> m = pf.CevMc(sigma=0.2*10, beta=0.5, intr=0.05, divr=0.1)
+        >>> m.set_num_params(n_path=160000, dt=None, rn_seed=123456)
+        >>> m.price(np.arange(80, 121, 10), spot=100, texp=1.2)
+        array([16.03769545,  9.88906886,  5.51705317,  2.78576033,  1.27297844])
+
+        References:
+            - Kang C (2014) Simulation of the Shifted Poisson Distribution with an Application to the CEV Model. Management Science and Financial Engineering 20:27–32. https://doi.org/10.7737/MSFE.2014.20.1.027
+    """
+
+    dt = 0.05
+    n_path = 10000
+    rn_seed = None
+    rng = np.random.default_rng(None)
+    correct_fwd = False
+
+    tobs = sv.CondMcBsmABC.tobs
+
+    def set_num_params(self, n_path=10000, dt=0.25, rn_seed=None):
+        """
+        Set MC parameters
+
+        Args:
+            n_path: number of paths
+            dt: time step for Euler/Milstein steps
+            rn_seed: random number seed
+        """
+        self.n_path = int(n_path)
+        self.dt = dt
+        self.rn_seed = rn_seed
+        self.rng = np.random.default_rng(rn_seed)
+
+    def price_step(self, spot, dt):
+        """
+        Simulated asset price after dt given the current price `spot` using the Algorithm 3 of Kang (2014)
+
+        Args:
+            spot: current price
+            dt: time step
+
+        Returns:
+            Simulated price after `dt`
+
+        References:
+            - Kang C (2014) Simulation of the Shifted Poisson Distribution with an Application to the CEV Model. Management Science and Financial Engineering 20:27–32. https://doi.org/10.7737/MSFE.2014.20.1.027
+        """
+
+        nz_idx = (spot > 0.0)
+        s_t = np.exp((self.intr - self.divr)*dt) * spot
+        betac = 1.0 - self.beta
+
+        var = (betac * self.sigma)**2 * dt * avg_exp(2*(self.intr-self.divr)*betac*dt)
+        z0 = np.power(s_t[nz_idx], 2*betac) / var
+        rv_gam = 2 * self.rng.standard_gamma(1/(2*betac), size=len(z0))
+        pois_lam = (z0 - rv_gam) / 2
+
+        # index for negative poisson value so that s_t becomes zero in this step
+        neg_idx = (pois_lam <= 0)
+
+        if np.any(neg_idx):
+            pois_lam = pois_lam[~neg_idx]  # positive only
+            tmp_idx = np.where(nz_idx)[0][neg_idx]
+            nz_idx[tmp_idx] = False
+            s_t[~nz_idx] = 0.0
+
+        zt = 2 * self.rng.standard_gamma(self.rng.poisson(pois_lam) + 1)
+        s_t[nz_idx] = np.power(var * zt, 1/(2*betac))
+        return s_t
+
+    def mass_zero(self, spot, texp, log=False):
+
+        tobs = self.tobs(texp)
+        dt = np.diff(tobs, prepend=0)
+        n_dt = len(dt)
+        s_t = np.full(self.n_path, spot)
+
+        for i in range(n_dt):
+            s_t = self.price_step(s_t, dt[i])
+
+        zm = (s_t <= 0).mean()
+        if log:
+            zm = np.log(zm)
+
+        return zm
+
+    def price(self, strike, spot, texp, cp=1):
+
+        tobs = self.tobs(texp)
+        dt = np.diff(tobs, prepend=0)
+        n_dt = len(dt)
+
+        s_t = np.full(self.n_path, spot)
+
+        for i in range(n_dt):
+            s_t = self.price_step(s_t, dt[i])
+
+        df = np.exp(-self.intr * texp)
+        p = np.mean(np.fmax(cp*(s_t - strike[:, None]), 0), axis=1)
+
+        return df * p
diff --git a/pyfeng/ex.py b/pyfeng/ex.py
index d38ca94..ad7648f 100644
--- a/pyfeng/ex.py
+++ b/pyfeng/ex.py
@@ -8,7 +8,6 @@
 # SABR / OUSV models for research
 from .sabr_int import SabrMixture
 from .sabr_mc import SabrMcCai2017Exact
-from .ousv import OusvSchobelZhu1998, OusvMcTimeStep, OusvMcChoi2023KL
 
 # Basket-Asian from ASP 2021
 from .multiasset_Ju2002 import BsmBasketAsianJu2002, BsmContinuousAsianJu2002
diff --git a/pyfeng/garch.py b/pyfeng/garch.py
index 42f61c2..d560158 100644
--- a/pyfeng/garch.py
+++ b/pyfeng/garch.py
@@ -7,7 +7,7 @@
 
 class GarchUncorrBaroneAdesi2004(sv.SvABC):
     """
-    Barone-Adesi et al (2004)'s approximation pricing formula for European options under uncorrelated (rho=0) GARCH diffusion model.
+    Barone-Adesi et al. (2004)'s approximation pricing formula for European options under uncorrelated (rho=0) GARCH diffusion model.
     Up to 2nd order is implemented.
 
     References:
@@ -26,8 +26,8 @@ def avgvar_mv(self, texp, var0):
         Eqs. (12)-(13) in Barone-Adesi et al. (2005)
 
         Args:
-            var0: initial variance
             texp: time step
+            var0: initial variance
 
         Returns:
             mean, variance
@@ -53,8 +53,7 @@ def avgvar_mv(self, texp, var0):
         M2c_3 = -vov2*(theta2*(4*mr*(3 - mr_t) + (2*mr_t - 5)*vov2) + term2*var0*(2*theta + var0))
 
         M2c_4 = 2*e_mr*vov2
-        M2c_4 *= 2*theta2*(mr_t*mr - (1 + mr_t)*vov2) \
-                 + var0*(2*mr*theta*(1 + texp*term1) + term1*var0)
+        M2c_4 *= 2*theta2*(mr_t*mr - (1 + mr_t)*vov2) + var0*(2*mr*theta*(1 + texp*term1) + term1*var0)
 
         M2c = M2c_1/mr2 + M2c_2/(term1*term2)**2 + M2c_3/mr2/term2**2 + M2c_4/mr2/term1**2
         M2c /= texp**2
@@ -79,7 +78,7 @@ def price(self, strike, spot, texp, cp=1):
         return price
 
 
-class GarchMcTimeStep(sv.SvABC, sv.CondMcBsmABC):
+class GarchMcTimeDisc(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
diff --git a/pyfeng/heston.py b/pyfeng/heston.py
index c14f03d..db5bbaa 100644
--- a/pyfeng/heston.py
+++ b/pyfeng/heston.py
@@ -34,7 +34,7 @@ def phi_exp(self, texp):
         phi = 4*self.mr / self.vov**2 / (1/exp - exp)
         return phi, exp
 
-    def var_mv(self, var0, dt):
+    def var_mv(self, dt, var0=None):
         """
         Mean and variance of the variance V(t+dt) given V(0) = var_0
 
@@ -45,11 +45,13 @@ def var_mv(self, var0, dt):
         Returns:
             mean, variance
         """
+        if var0 is None:
+            var0 = self.sigma
 
-        expo = np.exp(-self.mr*dt)
-        m = self.theta + (var0 - self.theta)*expo
-        s2 = var0*expo + self.theta*(1 - expo)/2
-        s2 *= self.vov**2*(1 - expo)/self.mr
+        e_mr = np.exp(-self.mr*dt)
+        m = self.theta + (var0 - self.theta)*e_mr
+        s2 = var0*e_mr + self.theta*(1 - e_mr)/2
+        s2 *= self.vov**2*(1 - e_mr)/self.mr
         return m, s2
 
     def avgvar_mv(self, texp, var0=None):
@@ -67,13 +69,16 @@ def avgvar_mv(self, texp, var0=None):
 
         ### We take var0 as argument instead of using self.sigma
         ### because it is used in var_step_qe method in HestonMcAndersen2008
-        var0 = var0 or self.sigma
+
+        if var0 is None:
+            var0 = self.sigma
 
         mr_t = self.mr*texp
-        e_mr_t = np.exp(-mr_t)
+        e_mr = np.exp(-mr_t)
+        phi = (1 - e_mr)/mr_t
         x0 = var0 - self.theta
-        mean = self.theta + x0*(1 - e_mr_t)/mr_t
-        var = (self.theta - 2*x0*e_mr_t) + (1 - e_mr_t)*(var0 - 2.5*self.theta + (var0 - self.theta/2)*e_mr_t)/mr_t
+        mean = self.theta + x0*phi
+        var = (self.theta - 2*x0*e_mr) + (var0 - 2.5*self.theta + (var0 - self.theta/2)*e_mr)*phi
         var *= (self.vov/mr_t)**2 * texp
         return mean, var
 
@@ -83,13 +88,14 @@ def strike_var_swap_analytic(self, texp, dt):
 
         Args:
             texp: time to expiry
-            dt: observation time step. If zero, continuous monitoring
+            dt: observation time step (e.g., dt=1/12 for monthly) For continuous monitoring, set dt=0
 
         Returns:
             Fair strike
 
         References:
             - Bernard C, Cui Z (2014) Prices and Asymptotics for Discrete Variance Swaps. Applied Mathematical Finance 21:140–173. https://doi.org/10.1080/1350486X.2013.820524
+
         """
 
         var0 = self.sigma
diff --git a/pyfeng/heston_mc.py b/pyfeng/heston_mc.py
index 115fab2..1de439e 100644
--- a/pyfeng/heston_mc.py
+++ b/pyfeng/heston_mc.py
@@ -114,15 +114,15 @@ def draw_from_mv(self, mean, var, dist):
     @abc.abstractmethod
     def cond_states_step(self, dt, var_0):
         """
-        Final variance and integrated variance over dt given var_0
-        The int_var is normalized by dt
+        Final variance after dt and average variance over (0, dt) given var_0.
+        `var_0` should be an array of (self.n_path, )
 
         Args:
             dt: time step
             var_0: initial variance
 
         Returns:
-            (var_t, avgvar)
+            (variance after dt, average variance during dt)
         """
         return NotImplementedError
 
@@ -194,11 +194,11 @@ def draw_log_return(self, dt, var_0, var_t, avgvar):
         Returns:
             log return
         """
-        mean_ln = self.rho / self.vov * ((var_t - var_0) + self.mr * dt * (avgvar - self.theta)) \
-                  + (self.intr - self.divr - 0.5 * avgvar) * dt
-        sigma_ln = np.sqrt((1.0 - self.rho**2) * dt * avgvar)
+        ln_m = self.rho/self.vov * ((var_t - var_0) + self.mr * dt * (avgvar - self.theta)) \
+            + (self.intr - self.divr - 0.5 * avgvar) * dt
+        ln_sig = np.sqrt((1.0 - self.rho**2) * dt * avgvar)
         zn = self.rv_normal(spawn=5)
-        return mean_ln + sigma_ln * zn
+        return ln_m + ln_sig * zn
 
     def return_var_realized(self, texp, cond=False):
         """
@@ -805,7 +805,7 @@ def var_step_qe(self, dt, var_0):
            variance after dt
         """
 
-        m, psi = self.var_mv(var_0, dt)  # put variance into psi
+        m, psi = self.var_mv(dt, var_0)  # put variance into psi
         psi /= m**2
 
         zz = self.rv_normal(spawn=0)
@@ -953,7 +953,9 @@ def avgvar_var_unexplained(self, texp, dt=None):
             ratio
         """
 
-        dt = dt or self.dt
+        if dt is None:
+            dt = self.dt
+
         mean, var = self.avgvar_mv(texp)
 
         m_x, v_x = self.x1star_avgvar_mv(dt, kk=0)
diff --git a/pyfeng/nsvh.py b/pyfeng/nsvh.py
index 95eb26d..64e089a 100644
--- a/pyfeng/nsvh.py
+++ b/pyfeng/nsvh.py
@@ -1,11 +1,132 @@
+import abc
+import os
+import pandas as pd
 import numpy as np
 import scipy.stats as spst
 import scipy.special as spsp
 import scipy.optimize as spop
-from . import sabr
 
+from . import opt_smile_abc as smile
 
-class Nsvh1(sabr.SabrABC):
+class NsvhABC(smile.OptSmileABC, abc.ABC):
+    beta = 0.0  ## should be fixed as 0
+    vov, rho, lam = 0.0, 0.0, 0.0
+    model_type = "Nsvh"
+    var_process = False
+
+    def __init__(self, sigma, vov=0.1, rho=0.0, lam=0.0, intr=0.0, divr=0.0, is_fwd=False, beta=None):
+        """
+        Args:
+            sigma: model volatility at t=0
+            vov: volatility of volatility
+            rho: correlation between price and volatility
+            lam: lambda. Normal 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
+        self.vov = vov
+        self.rho = rho
+        super().__init__(sigma, intr=intr, divr=divr, is_fwd=is_fwd)
+
+    def params_kw(self):
+        params = super().params_kw()
+        extra = {"vov": self.vov, "lam": self.lam, "rho": self.rho}
+        return {**params, **extra}  # Py 3.9, params | extra
+
+    def vol_smile(self, strike, spot, texp, cp=1, model=None):
+        return super().vol_smile(strike, spot, texp, cp=cp, model="norm")
+
+    @classmethod
+    def init_benchmark(cls, set_no=None):
+        """
+        Initiate a SABR model with stored benchmark parameter sets
+
+        Args:
+            set_no: set number
+
+        Returns:
+            Dataframe of all test cases if set_no = None
+            (model, Dataframe of result, params) if set_no is specified
+
+        References:
+            - Antonov, Alexander, Konikov, M., & Spector, M. (2013). SABR spreads its wings. Risk, 2013(Aug), 58–63.
+            - Antonov, Alexandre, Konikov, M., & Spector, M. (2019). Modern SABR Analytics. Springer International Publishing. https://doi.org/10.1007/978-3-030-10656-0
+            - Antonov, Alexandre, & Spector, M. (2012). Advanced analytics for the SABR model. Available at SSRN. https://ssrn.com/abstract=2026350
+            - Cai, N., Song, Y., & Chen, N. (2017). Exact Simulation of the SABR Model. Operations Research, 65(4), 931–951. https://doi.org/10.1287/opre.2017.1617
+            - Korn, R., & Tang, S. (2013). Exact analytical solution for the normal SABR model. Wilmott Magazine, 2013(7), 64–69. https://doi.org/10.1002/wilm.10235
+            - Lewis, A. L. (2016). Option valuation under stochastic volatility II: With Mathematica code. Finance Press.
+            - von Sydow, L., ..., Haentjens, T., & Waldén, J. (2018). BENCHOP - SLV: The BENCHmarking project in Option Pricing – Stochastic and Local Volatility problems. International Journal of Computer Mathematics, 1–14. https://doi.org/10.1080/00207160.2018.1544368
+        """
+        this_dir, _ = os.path.split(__file__)
+        file = os.path.join(this_dir, "data/sabr_benchmark.xlsx")
+        df_param = pd.read_excel(file, sheet_name="Param", index_col="Sheet")
+
+        if set_no is None:
+            return df_param
+        else:
+            df_val = pd.read_excel(file, sheet_name=str(set_no))
+            param = df_param.loc[set_no]
+            args_model = {k: param[k] for k in ("sigma", "vov", "rho", "beta")}
+            args_pricing = {k: param[k] for k in ("texp", "spot")}
+
+            assert df_val.columns[0] == "Strike"
+            args_pricing["strike"] = df_val.values[:, 0]
+
+            val = df_val[param["col_name"]].values
+            is_iv = param["col_name"].startswith("IV")
+
+            m = cls(**args_model)
+
+            param_dict = {
+                "args_pricing": args_pricing,
+                "ref": param["Reference"],
+                "val": val,
+                "is_iv": is_iv,
+            }
+
+            return m, df_val, param_dict
+
+    def price_vsk(self, texp=1):
+        """
+        Variance, skewness, and ex-kurtosis. Corollary 5 (Appendix B) in Choi et al. (2019)
+
+        Args:
+            texp: time to expiry
+
+        Returns:
+            (variance, skewness, and ex-kurtosis)
+
+        References:
+            - Choi J, Liu C, Seo BK (2019) Hyperbolic normal stochastic volatility model. J Futures Mark 39:186–204. https://doi.org/10.1002/fut.21967
+        """
+        vovn = self.vov * np.sqrt(texp)
+        ww = np.exp(vovn**2)
+        wlam = np.exp(self.lam * vovn**2)
+        rho2 = self.rho**2
+        rhoc2 = 1 - rho2
+
+        wf1 = (wlam * ww - 1) / (1 + self.lam)
+        wf3 = (wlam * ww**3 - 1) / (3 + self.lam)
+        wf5 = (wlam * ww**5 - 1) / (5 + self.lam)
+
+        m2 = rho2 * wlam*(ww-1) + rhoc2 * wf1
+        m3 = self.rho*np.sqrt(wlam) * (rho2*wlam*(ww-1)**2*(ww+2) + 3*rhoc2*(wf3 - wf1))
+        m4 = (rho2*wlam*(ww-1))**2 * (ww**2*(ww**2 + 2*ww + 3) - 3) + 6*rho2*rhoc2*wlam*(ww*wf5 - 2*wf3 + wf1) +\
+             1.5*rhoc2**2 * (-wlam*ww*wf5 + (wlam*ww**3 + 1)*wf3 - wf1)
+
+        skew = m3 / np.sqrt(m2**3)
+        exkurt = m4 / m2**2 - 3
+
+        return m2 * (self.sigma/self.vov)**2, skew, exkurt
+
+
+class Nsvh1(NsvhABC):
     """
     Hyperbolic Normal Stochastic Volatility (NSVh) model with lambda=1 by Choi et al. (2019)
 
@@ -20,8 +141,7 @@ class Nsvh1(sabr.SabrABC):
         array([25.51200027, 17.87539874, 11.47308947,  6.75128331,  3.79464422])
     """
 
-    model_type = "Nsvh"
-    beta = 0.0  # beta is already defined in the parent class, but the default value set as 0
+    lam = 1.0
     is_atmvol = False
 
     def _sig0_from_atmvol(self, texp):
@@ -38,7 +158,7 @@ def _sig0_from_atmvol(self, texp):
         sig0 = self.sigma * np.sqrt(texp/2/np.pi) / price
         return sig0
 
-    def __init__(self, sigma, vov=0.0, rho=0.0, beta=None, intr=0.0, divr=0.0, is_fwd=False, is_atmvol=False):
+    def __init__(self, sigma, vov=0.1, rho=0.0, intr=0.0, divr=0.0, is_fwd=False, is_atmvol=False, beta=None):
         """
         Args:
             sigma: model volatility at t=0
@@ -50,11 +170,8 @@ def __init__(self, sigma, vov=0.0, rho=0.0, beta=None, intr=0.0, divr=0.0, is_fw
             is_fwd: if True, treat `spot` as forward price. False by default.
             is_atmvol: If True, use `sigma` as the ATM normal vol. 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.is_atmvol = is_atmvol
-        super().__init__(sigma, vov, rho, beta=0, intr=intr, divr=divr, is_fwd=is_fwd)
+        super().__init__(sigma, vov, rho, lam=1.0, intr=intr, divr=divr, is_fwd=is_fwd)
 
     def price(self, strike, spot, texp, cp=1):
         fwd, df, _ = self._fwd_factor(spot, texp)
@@ -88,12 +205,13 @@ def cdf(self, strike, spot, texp, cp=-1):
         d = (np.arctanh(self.rho) + np.arcsinh(((fwd - strike)*vovn/sig_sqrt - self.rho*vov_var) / rhoc)) / vovn
         return spst.norm.cdf(cp * d)
 
-    def moments_vsk(self, texp=1):
+    def price_vsk(self, texp=1):
         """
-        Variance, skewness, and ex-kurtosis
+        Variance, skewness, and ex-kurtosis of forward price.
+        It is a special case (lambda=1) of NsvhABC.price_vsk()
 
         Args:
-            texp: time-to-expiry
+            texp: time to expiry
 
         Returns:
             (variance, skewness, and ex-kurtosis)
@@ -116,27 +234,23 @@ def moments_vsk(self, texp=1):
         # M3 = (1 / 4)*np.sqrt(ww)*(ww - 1)**2*(self.rho*C_31 + rho3*C_33)
         # M4 = (1 / 8)*(ww - 1)**2*(C_40 + rho2*C_42 + rho4*C_44)
 
-        m2 = (1 / 2) * (ww - 1) * (c20 + rho2 * c22)
+        m2 = 0.5 * (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)
-        )
+        skew = np.sqrt(ww * (ww - 1) / 2) * (self.rho * c31 + rho3 * c33) / np.power(c20 + rho2 * c22, 1.5)
         exkurt = 0.5 * (ww - 1) * (k0 + rho2 * k2 + rho4 * k4) / (c20 + rho2 * c22)**2
 
-        return m2 * (self.sigma / self.vov)**2, skew, exkurt
+        return m2 * (self.sigma/self.vov)**2, skew, exkurt
 
     def calibrate_vsk(self, var, skew, exkurt, texp=1, setval=False):
         """
         Calibrate parameters to the moments: variance, skewness, ex-kurtosis.
 
         Args:
-            texp: time-to-expiry
+            texp: time to expiry
             var: variance
             skew: skewness
             exkurt: ex-kurtosis. should be > 0.
@@ -186,7 +300,7 @@ def f_beta1(w):
         return sig0, vov, rho
 
 
-class NsvhMc(sabr.SabrABC):
+class NsvhMc(NsvhABC):
     """
     Monte-Carlo model of Hyperbolic Normal Stochastic Volatility (NSVh) model.
 
@@ -212,30 +326,11 @@ class NsvhMc(sabr.SabrABC):
         array([-0.00328887,  0.00523714,  0.00808885,  0.0069694 ,  0.00205566])
     """
 
-    lam = 0
     n_path = int(16e4)
     rn_seed = None
     rng = np.random.default_rng(None)
     antithetic = True
 
-    def __init__(self, sigma, vov=0.0, rho=0.0, lam=0.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. Normal 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_num_params(self, n_path=1e6, rn_seed=None, antithetic=True):
         self.n_path = int(n_path)
         self.rn_seed = rn_seed
@@ -248,7 +343,7 @@ def mc_vol_price(self, texp):
         Simulate volatility and price pair
 
         Args:
-            texp: time-to-expiry
+            texp: time to expiry
 
         Returns: (vol, price). vol: (n_path, ), price: (n_path, 2)
 
@@ -275,11 +370,7 @@ def mc_vol_price(self, texp):
 
         path = (
             df_z,
-            (self.sigma / self.vov)
-            * (
-                self.rho * (df_z[:, None] - np.exp(0.5 * self.lam * vol_var))
-                + rhoc * df_w
-            ),
+            (self.sigma/self.vov) * (self.rho * (df_z[:, None] - np.exp(0.5*self.lam*vol_var)) + rhoc * df_w),
         )
         return path
 
@@ -290,7 +381,7 @@ def price(self, strike, spot, texp, cp=1):
         Args:
             strike: strike price
             spot: spot price
-            texp: time to np.expiry
+            texp: time to expiry
             cp: 1/-1 for call/put
 
         Returns:
@@ -309,7 +400,7 @@ def price(self, strike, spot, texp, cp=1):
         return df * price
 
 
-class NsvhGaussQuad(sabr.SabrABC):
+class NsvhGaussQuad(NsvhABC):
     """
     Quadrature integration method of Hyperbolic Normal Stochastic Volatility (NSVh) model.
 
@@ -346,28 +437,12 @@ class NsvhGaussQuad(sabr.SabrABC):
 
     n_quad = (7, 7)
 
-    def __init__(self, sigma, vov=0.0, rho=0.0, lam=0.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. Normal 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 price(self, strike, spot, texp, cp=1):
 
         fwd, df, _ = self._fwd_factor(spot, texp)
-        _, _, rhoc, rho2, vovn = self._variables(1.0, texp)
+        rho2 = self.rho**2
+        rhoc = np.sqrt(1.0 - rho2)
+        vovn = self.vov * np.sqrt(np.maximum(texp, 1e-64))
 
         ### axis 1: nodes of x,y,z , get the weight of z,v
         z_value, z_weight = spsp.roots_hermitenorm(self.n_quad[0])
@@ -420,7 +495,9 @@ def price(self, strike, spot, texp, cp=1):
 
     def cdf(self, strike, spot, texp, cp=-1):
         fwd = self.forward(spot, texp)
-        _, _, rhoc, rho2, vovn = self._variables(1.0, texp)
+        rho2 = self.rho**2
+        rhoc = np.sqrt(1.0 - rho2)
+        vovn = self.vov * np.sqrt(np.maximum(texp, 1e-64))
 
         ### axis 1: nodes of x,y,z , get the weight of z,v
         z_value, z_weight = spsp.roots_hermitenorm(self.n_quad[0])
diff --git a/pyfeng/ousv.py b/pyfeng/ousv.py
index 4cbae4c..2fda536 100644
--- a/pyfeng/ousv.py
+++ b/pyfeng/ousv.py
@@ -4,34 +4,116 @@
 from . import sv_abc as sv
 from . import bsm
 
+#### Use of RN generation spawn:
+# 0: simulation of sigma: vol_step()
+# 1: truncated sine series: `z_gpqr`
+# 2: sine series: `z_sin`
+# 3: not used
+# 4: not used
+# 5: asset return
+
 
 class OusvABC(sv.SvABC, abc.ABC):
 
     model_type = "OUSV"
     var_process = False
 
-    def avgvar_mv(self, texp, vol0=None):
+    def avgvol_mv(self, texp, vol0=None, nz_theta=True):
         """
-        Mean and variance of the variance V(t+dt) given V(0) = var_0
+        Mean and variance of the volatility sigma(t+dt) given sigma(t) = var_0
         (variance is not implemented yet)
 
         Args:
             texp: time step
-            vol0: initial volatility
+            vol0: initial sigma
+            nz_theta: non-zero theta. True by default. If False, assume theta=0 making computation simpler.
 
         Returns:
             mean, variance(=None)
         """
 
-        vol0 = vol0 or self.sigma
+        if vol0 is None:
+            vol0 = self.sigma
+
+        if nz_theta:
+            vol0 = vol0 - self.theta
+
+        mr_t = self.mr * texp
+        e_mr = np.exp(-mr_t)
+        phi = (1 - e_mr)/mr_t
+        m = vol0 * phi
+
+        if nz_theta:
+            m += self.theta
+
+        return m, None
+
+    def avgvar_mv(self, texp, vol0=None, nz_theta=True):
+        """
+        Mean and variance of the variance sigma^2(t+dt) given sigma(t) = var_0
+        (variance is not implemented yet)
+
+        Args:
+            texp: time step
+            vol0: initial sigma
+            nz_theta: non-zero theta. True by default. If False, assume theta=0 making computation simpler.
+
+        Returns:
+            mean, variance(=None)
+        """
+
+        if vol0 is None:
+            vol0 = self.sigma
+
+        if nz_theta:
+            vol0 = vol0 - self.theta
+
+        mr_t = self.mr * texp
+        e_mr = np.exp(-mr_t)
+        phi = (1 - e_mr)/mr_t
+        phi2 = (1 + e_mr)/2 * phi
+        vv = vol0**2 * phi2 + 0.5*(self.vov**2/self.mr)*(1 - phi2)
+
+        if nz_theta:
+            vv += self.theta * (self.theta + 2*vol0 * phi)
 
-        mrt = self.mr * texp
-        e_mrt = np.exp(-mrt)
-        x0 = vol0 - self.theta
-        vv = self.vov**2/2/self.mr + self.theta**2 + \
-             ((x0**2 - self.vov**2/2/self.mr)*(1 + e_mrt) + 4*self.theta * x0)*(1 - e_mrt)/(2*mrt)
         return vv, None
 
+    def cond_avgvolvar_m(self, texp, vol_0, vol_t, nz_theta=True):
+        """
+        Conditional expectation of average vol and variance
+        Args:
+            texp: time step
+            vol_0: initial vol
+            vol_t: final vol
+            nz_theta: non-zero theta. True by default. If False, assume theta=0 making computation simpler.
+
+        Returns:
+            avgvol mean, avgvar mean
+        """
+        mr_t = self.mr * texp
+        e_mr = np.exp(-mr_t)
+        phi = (1 - e_mr)/mr_t
+        phi2 = (1 + e_mr)/2 * phi
+        sinh = np.sinh(mr_t)
+        cosh = np.cosh(mr_t)
+
+        if nz_theta:
+            vol_0 = vol_0 - self.theta  # Don't do `vol_0 -= self.theta`. Don't change `vol_0`
+            vol_t = vol_t - self.theta
+
+        volhat = vol_t - vol_0*e_mr
+
+        vol_m = (vol_0 + volhat/(1+e_mr))*phi
+        var_m = vol_0**2 * phi2 + 0.5*(self.vov**2/self.mr) * (cosh/sinh - 1/mr_t)
+        var_m += volhat * (volhat * (sinh*cosh - mr_t)/(2*mr_t * sinh**2) + vol_0 * (e_mr/mr_t) * (1/phi2 - 1))
+
+        if nz_theta:
+            var_m += self.theta*(self.theta + 2*vol_m)
+            vol_m += self.theta
+
+        return vol_m, var_m
+
     def strike_var_swap_analytic(self, texp, dt):
         """
         Analytic fair strike of variance swap. Eq (17), (24), (25) in Bernard & Cui (2014)
@@ -48,13 +130,12 @@ def strike_var_swap_analytic(self, texp, dt):
 
         """
 
-
         ### continuously monitored fair strike (same as mean of avgvar_mv)
         mrt = self.mr * texp
         e_mrt = np.exp(-mrt)
+        phi = (e_mrt - 1) / mrt  # D / T
         x0 = self.sigma - self.theta
-        strike = self.vov**2/2/self.mr + self.theta**2 + \
-             ((x0**2 - self.vov**2/2/self.mr)*(1 + e_mrt) + 4*self.theta * x0)*(1 - e_mrt)/(2*mrt)
+        strike, _ = self.avgvar_mv(texp, self.sigma)
 
         if not np.all(np.isclose(dt, 0.0)):
             mrt2 = mrt**2
@@ -65,16 +146,14 @@ def strike_var_swap_analytic(self, texp, dt):
             th = self.theta
             th2 = th**2
 
-            D = (e_mrt - 1) / mrt  # D / T
             E = 4*mrt2*(sig2**2 - th2**2) - 3*vov2t*(vov2t + 4*th2*mrt)  # E * T^2
-
-            d2 = (vov2t + 2*mrt*th2) + ((2*mrt*(th2-sig2) + vov2t) + mrt*(vov2t/2 - mrt*x0**2)*D)*D
+            d2 = (vov2t + 2*mrt*th2) + ((2*mrt*(th2-sig2) + vov2t) + mrt*(vov2t/2 - mrt*x0**2)*phi)*phi
 
             # d1 / T
-            d1 = sig2**2/4 - E*(1 + D)/(16*mrt2)
-            d1 += ((3*vov2t/4 - mrt*sig*x0/2)*sig2 + E/(32*mrt)) * D**2
-            d1 += ((2*th*sig/3 - sig2/6 - th2/2)*sig2*mrt2 - E/48 + (-mrt*sig*th + 3/4*sig2*mrt - vov2t/4)*vov2t) * D**3
-            d1 += (E/(8*mrt) - 3*vov2t*x0*th + 3*sig2*vov2t/2 - mrt*sig*x0*(2*th2 - th*sig + sig2)) * D**4 * mrt2 / 8
+            d1 = sig2**2/4 - E*(1 + phi)/(16*mrt2)
+            d1 += ((3*vov2t/4 - mrt*sig*x0/2)*sig2 + E/(32*mrt)) * phi**2
+            d1 += ((2*th*sig/3 - sig2/6 - th2/2)*sig2*mrt2 - E/48 + (-mrt*sig*th + 3/4*sig2*mrt - vov2t/4)*vov2t) * phi**3
+            d1 += (E/(8*mrt) - 3*vov2t*x0*th + 3*sig2*vov2t/2 - mrt*sig*x0*(2*th2 - th*sig + sig2)) * phi**4 * mrt2 / 8
 
             correction = self.intr * (self.intr - strike) + d1 - self.vov/(2*self.mr)*self.rho * d2/texp
             strike += correction * dt
@@ -122,7 +201,7 @@ def D_B_C(self, s1, s2, s3, texp):
             -0.5 * np.log(cossin)
             + 0.5 * mr * texp
             + ((mr * theta * gamma1)**2 - gamma3**2)
-            / (2 * s2g3)
+            / (2*s2g3)
             * (sinh / cossin - gamma1 * texp)
             + ktg3 * gamma3 / s2g3 * ((cosh - 1) / cossin)
         )
@@ -210,20 +289,21 @@ def price(self, strike, spot, texp, cp=1):
 class OusvMcABC(OusvABC, sv.CondMcBsmABC, abc.ABC):
 
     @abc.abstractmethod
-    def cond_states_step(self, dt, vol_0):
+    def cond_states_step(self, dt, vol_0, nz_theta=True):
         """
-        Final volatility (sigma), average variance and volatilityu over dt given vol_0
+        Final volatility (sigma), average variance and volatility over dt given vol_0
 
         Args:
             dt: time-to-expiry
             vol_0: initial volatility
+            nz_theta: non-zero theta. True by default. If False, assume theta=0 making computation simpler.
 
         Returns:
             (final vol, average var, average vol)
         """
         return NotImplementedError
 
-    def vol_step(self, dt, vol_0, zn=None):
+    def vol_step(self, dt, vol_0, zn=None, nz_theta=True):
         """
         Stepping volatility according to OU process dynamics
 
@@ -231,15 +311,23 @@ def vol_step(self, dt, vol_0, zn=None):
             vol_0: initial volatility
             dt: time step
             zn: specified normal rv to use (n_path, )
+            nz_theta: non-zero theta. True by default. If False, assume theta=0 making computation simpler.
 
         Returns:
             volatility after dt
         """
+        if nz_theta:
+            vol_0 = vol_0 - self.theta  # Don't do `vol_0 -= self.theta`. Don't change `vol_0`
+
         e_mr = np.exp(-self.mr * dt)
         if zn is None:
             zn = self.rv_normal(spawn=0)
 
-        vol_t = self.theta + (vol_0 - self.theta)*e_mr + self.vov*np.sqrt((1 - e_mr**2)/(2*self.mr))*zn
+        vol_t = vol_0*e_mr + self.vov*np.sqrt((1 - e_mr**2)/(2*self.mr))*zn
+
+        if nz_theta:
+            vol_t += self.theta
+
         return vol_t
 
     def cond_spot_sigma(self, texp, vol_0):
@@ -247,24 +335,27 @@ def cond_spot_sigma(self, texp, vol_0):
         dt = np.diff(tobs, prepend=0)
         n_dt = len(dt)
 
-        vol_t = np.full(self.n_path, vol_0)
+        vol_t = np.full(self.n_path, vol_0 - self.theta)
         avgvar = np.zeros(self.n_path)
         avgvol = np.zeros(self.n_path)
 
         for i in range(n_dt):
-            vol_t, avgvar_inc, avgvol_inc = self.cond_states_step(dt[i], vol_t)
+            vol_t, avgvar_inc, avgvol_inc = self.cond_states_step(dt[i], vol_t, nz_theta=False)
             avgvar += avgvar_inc * dt[i]
             avgvol += avgvol_inc * dt[i]
 
         avgvar /= texp
         avgvol /= texp
 
-        spot_cond = (vol_t**2 - vol_0**2) / (2 * self.vov) - self.vov * texp / 2 \
-            - (self.mr * self.theta / self.vov) * texp * avgvol \
-            + (self.mr / self.vov - self.rho / 2) * texp * avgvar
-        np.exp(self.rho * spot_cond, out=spot_cond)
+        avgvar += self.theta * (self.theta + 2*avgvol)
+        avgvol += self.theta
+        vol_t += self.theta
+
+        spot_cond = (vol_t**2 - vol_0**2 - self.vov**2*texp) +\
+            texp*(-2*(self.mr * self.theta) * avgvol + (2*self.mr - self.rho*self.vov) * avgvar)
+        np.exp(0.5*self.rho/self.vov * spot_cond, out=spot_cond)
 
-        sigma_cond = np.sqrt((1 - self.rho**2) * avgvar) / vol_0
+        sigma_cond = np.sqrt((1 - self.rho**2) * np.fmax(avgvar, 1e-64)) / vol_0
         return spot_cond, sigma_cond
 
     def strike_var_swap_analytic(self, texp, dt=None):
@@ -288,10 +379,11 @@ def cond_log_return_var(self, dt, vol_0, vol_t, avgvar, avgvol):
 
         """
         rho_vov = self.rho / self.vov
-        mean_ln = rho_vov * ((vol_t**2 - vol_0**2)/2 - self.mr*self.theta*dt*avgvol) - self.rho*self.vov*dt/2 \
-                  + (rho_vov*self.mr - 0.5)*dt*avgvar + (self.intr - self.divr)*dt
-        sigma_ln2 = (1.0 - self.rho**2) * dt * avgvar
-        return mean_ln**2 + sigma_ln2
+        ln_m = (self.intr - self.divr - self.rho*self.vov/2)*dt \
+               + rho_vov * ((vol_t**2 - vol_0**2)/2 - self.mr*self.theta*dt*avgvol) \
+               + (rho_vov*self.mr - 0.5)*dt*avgvar
+        ln_sig2 = (1.0 - self.rho**2) * dt * avgvar
+        return ln_m**2 + ln_sig2
 
     def draw_log_return(self, dt, vol_0, vol_t, avgvar, avgvol):
         """
@@ -308,11 +400,12 @@ def draw_log_return(self, dt, vol_0, vol_t, avgvar, avgvol):
             log return
         """
         rho_vov = self.rho / self.vov
-        mean_ln = rho_vov * ((vol_t**2 - vol_0**2)/2 - self.mr*self.theta*dt*avgvol) - self.rho*self.vov*dt/2 \
-                  + (rho_vov*self.mr - 0.5)*dt*avgvar + (self.intr - self.divr)*dt
-        sigma_ln = np.sqrt((1.0 - self.rho**2) * dt * avgvar)
+        ln_m = (self.intr - self.divr - self.rho*self.vov/2)*dt\
+               + rho_vov * ((vol_t**2 - vol_0**2)/2 - self.mr*self.theta*dt*avgvol)\
+               + (rho_vov*self.mr - 0.5)*dt*avgvar
+        ln_sig = np.sqrt((1.0 - self.rho**2) * dt * avgvar)
         zn = self.rv_normal(spawn=5)
-        return mean_ln + sigma_ln * zn
+        return ln_m + ln_sig * zn
 
     def return_var_realized(self, texp, cond=False):
         """
@@ -338,18 +431,19 @@ def return_var_realized(self, texp, cond=False):
             if cond:
                 var_r += self.cond_log_return_var(dt[i], vol_0, vol_t, avgvar_inc, avgvol_inc)
             else:
-                var_r += self.draw_log_return(dt[i], vol_0, vol_t, avgvar_inc, avgvol_inc) ** 2
+                var_r += self.draw_log_return(dt[i], vol_0, vol_t, avgvar_inc, avgvol_inc)**2
 
             vol_0 = vol_t
 
         return var_r / texp  # annualized
 
 
-class OusvMcTimeStep(OusvMcABC):
+class OusvMcTimeDisc(OusvMcABC):
     """
     OUSV model with conditional Monte-Carlo simulation
     The SDE of SV is: d sigma_t = mr (theta - sigma_t) dt + vov dB_T
     """
+    scheme = 0  ## 0 for trapezoidal, 1 for mean
 
     def vol_paths(self, tobs):
         # 2d array of (time, path) including t=0
@@ -372,11 +466,34 @@ def cond_states_full(self, texp, sig_0):
 
         return s_t, v_t_std, u_t_std
 
-    def cond_states_step(self, dt, vol_0):
+    def cond_states_step(self, dt, vol_0, nz_theta=True):
+        """
+        Final volatility (sigma), average variance and volatilityu over dt given vol_0
+
+        Args:
+            dt: time-to-expiry
+            vol_0: initial volatility
+            nz_theta: non-zero theta. True by default. If False, assume theta=0 making computation simpler.
+
+        Returns:
+            (final vol, average var, average vol)
+        """
+
+        if nz_theta:
+            vol_0 = vol_0 - self.theta  # Don't do `vol_0 -= self.theta`. Don't change `vol_0`
+
+        vol_t = self.vol_step(dt, vol_0, nz_theta=False)
 
-        vol_t = self.vol_step(dt, vol_0)
-        avgvol = (vol_0 + vol_t)/2
-        avgvar = (vol_0**2 + vol_t**2)/2
+        if self.scheme == 0:
+            avgvol = (vol_0 + vol_t) / 2
+            avgvar = (vol_0**2 + vol_t**2) / 2
+        elif self.scheme == 1:
+            avgvol, avgvar = self.cond_avgvolvar_m(dt, vol_0, vol_t, nz_theta=False)
+
+        if nz_theta:
+            avgvar += self.theta * (self.theta + 2*avgvol)
+            avgvol += self.theta
+            vol_t += self.theta
 
         return vol_t, avgvar, avgvol
 
@@ -521,11 +638,10 @@ def _a6sum(cls, mr_t, ns=0, odd=None):
         if odd == 2:  # even
             rv = cls._a6sum(mr_t / 2) / 2**6
         elif odd == 1:  # odd
-            rv = (3 * mr_t / np.tanh(mr_t) + (3 + 2 * mr_t / np.tanh(mr_t)) * mr_t**2 / np.sinh(mr_t)**2 - 8) / (
-                        2 * mr_t**6) - cls._a6sum(mr_t / 2) / 2**6
+            rv = (3*mr_t / np.tanh(mr_t) + (3 + 2*mr_t / np.tanh(mr_t))*mr_t**2 / np.sinh(mr_t)**2 - 8)/(2*mr_t**6) \
+                 - cls._a6sum(mr_t / 2) / 2**6
         else:  # all
-            rv = (3 * mr_t / np.tanh(mr_t) + (3 + 2 * mr_t / np.tanh(mr_t)) * mr_t**2 / np.sinh(mr_t)**2 - 8) / (
-                        2 * mr_t**6)
+            rv = (3*mr_t / np.tanh(mr_t) + (3 + 2*mr_t / np.tanh(mr_t))*mr_t**2 / np.sinh(mr_t)**2 - 8)/(2*mr_t**6)
 
         if ns == 0:
             return rv
@@ -576,7 +692,7 @@ def _a6n2sum(cls, mr_t, ns=0, odd=None):
             rv -= np.sum(a6n2)
         return rv
 
-    def cond_states_step(self, dt, vol_0, zn=None):
+    def cond_states_step(self, dt, vol_0, nz_theta=True, zn=None):
         """
         Final volatility (sigma), average variance and volatilityu over dt given vol_0
 
@@ -584,25 +700,22 @@ def cond_states_step(self, dt, vol_0, zn=None):
             dt: time-to-expiry
             vol_0: initial volatility
             zn: normal RVs to specify in the (1+n_sin, n_path) format. None by default.
+            nz_theta: non-zero theta. True by default. If False, assume theta=0 making computation simpler.
 
         Returns:
             (final vol, average var, average vol)
         """
 
-        mr, vov = self.mr, self.vov
-
-        mr_t = mr * dt
-        e_mr = np.exp(-mr_t)
-        sinh = np.sinh(mr_t)
-        cosh = np.cosh(mr_t)
+        mr_t = self.mr * dt
         vovn = self.vov * np.sqrt(dt)  # normalized vov
 
-        x_0 = vol_0 - self.theta
+        if nz_theta:
+            vol_0 = vol_0 - self.theta  # Don't do `vol_0 -= self.theta`. Don't change `vol_0`
+
         if zn is None:
-            x_t = self.vol_step(dt, vol_0) - self.theta
+            vol_t = self.vol_step(dt, vol_0, nz_theta=False)
         else:
-            x_t = self.vol_step(dt, vol_0, zn=zn[0, :]) - self.theta
-        sighat = x_t - x_0 * e_mr
+            vol_t = self.vol_step(dt, vol_0, zn=zn[0, :], nz_theta=False)
 
         if zn is None:
             n_sin = self.n_sin
@@ -629,52 +742,60 @@ def cond_states_step(self, dt, vol_0, zn=None):
             p_std = np.sqrt(self._a6n2sum(mr_t, ns=n_sin, odd=1))
             q_std = np.sqrt(self._a6n2sum(mr_t, ns=n_sin, odd=2))  # even
             corr = self._a4sum(mr_t, ns=n_sin, odd=1)/(g_std*p_std)
+            corr = np.clip(corr, -1.0, 1.0)
 
-            z_g[:] = (corr*z_p + np.sqrt(1 - corr**2)*z_g) * g_std
+            z_g[:] = (corr*z_p + np.sqrt(1.0 - corr**2)*z_g) * g_std
             z_p *= p_std
             z_q *= q_std
 
-            r_m = self._a2sum(mr_t, ns=n_sin)
             r_var = self._a4sum(mr_t, ns=n_sin)
-            z_r[:] = np.sqrt(r_var)*(z_r**2 - 1) + r_m
+            z_r[:] = np.sqrt(r_var)*(z_r**2 - 1)
         else:
             n_sin = zn.shape[0] - 1
+            r_m = self._a2sum(mr_t, ns=n_sin)
             z_sin = zn[1:, :]
-            z_g, z_p, z_q, z_r = 0.0, 0.0, 0.0, 0.0
+            z_g, z_p, z_q, z_r = 0.0, 0.0, 0.0, -r_m
+            #### -r_m is the correction to include only `an2 @ z_sin**2`
 
-        n_pi = np.pi * np.arange(1, n_sin + 1)
-        an2 = 2 / (mr_t**2 + n_pi**2)
-        an = np.sqrt(an2)
-        an3_n_pi = an2 * an * n_pi
+        if n_sin > 0:
+            n_pi = np.pi * np.arange(1, n_sin + 1)
+            an2 = 2 / (mr_t**2 + n_pi**2)
+            an = np.sqrt(an2)
+            an3_n_pi = an2 * an * n_pi
 
-        z_g += (an[::2] / n_pi[::2]) @ z_sin[::2, :]  # odd terms
-        z_p += an3_n_pi[::2] @ z_sin[::2, :]  # odd terms
-        z_q += an3_n_pi[1::2] @ z_sin[1::2, :]  # even terms
-        z_r += an2 @ z_sin**2
+            z_g += (an[::2] / n_pi[::2]) @ z_sin[::2, :]  # odd terms
+            z_p += an3_n_pi[::2] @ z_sin[::2, :]  # odd terms
+            z_q += an3_n_pi[1::2] @ z_sin[1::2, :]  # even terms
+            z_r += an2 @ (z_sin**2 - 1)
 
-        uu_t = x_0 * (1 - e_mr) / mr_t + (cosh - 1) / (mr_t * sinh) * sighat + 2 * vovn * z_g  # * dt
+        uu_t, vv_t = self.cond_avgvolvar_m(dt, vol_0, vol_t, nz_theta=False)
+        uu_t += 2 * vovn * z_g  # * dt
+        vv_t += vovn * ((vol_0*(z_p + z_q) + vol_t*(z_p - z_q)) + 0.5*vovn*z_r)
 
-        vv_t = (1 - e_mr**2) / (2 * mr_t) * x_0**2 + (sinh * cosh - mr_t) / (2 * mr_t * sinh**2) * sighat**2
-        vv_t += sighat * x_0 * (1/sinh - e_mr / mr_t) + vovn * (x_0 * (z_p + z_q) + x_t * (z_p - z_q))
-        vv_t += 0.5 * vovn**2 * z_r
+        if nz_theta:
+            vv_t += (2*uu_t + self.theta) * self.theta
+            uu_t += self.theta
+            vol_t += self.theta
 
-        ### sigma_t = x_t + theta
-        vv_t += (2*uu_t + self.theta) * self.theta
-        uu_t += self.theta
-        x_t += self.theta
+        return vol_t, vv_t, uu_t
 
-        return x_t, vv_t, uu_t
+    def unexplained_var_ratio(self, mr_t, ns=None):
 
-    def unexplained_var_ratio(self, mr_t, ns):
-        ns = ns or self.n_sin
+        if ns is None:
+            ns = self.n_sin
         rv = self._a4sum(mr_t, ns=ns) / self._a4sum(mr_t)
         return rv
 
+    def strike_var_swap_analytic(self, texp, dt=None):
+        if dt is None:
+            dt = self.dt
+        rv = super().strike_var_swap_analytic(texp, dt)
+        return rv
+
     def vol_path_sin(self, tobs, zn=None):
         """
         vol path composed of sin terms
         Args:
-            vol_t: terminal volatility (n_path, )
             tobs: observation time (n_time, )
             zn: specified normal rvs to use (n_sin + 1, n_path). None by default
 
@@ -686,23 +807,23 @@ def vol_path_sin(self, tobs, zn=None):
         e_mr = np.exp(-mr_t)
         e_mr_tobs = np.exp(-self.mr*tobs[:, None])
 
-        x_0 = self.sigma - self.theta
+        vol_0 = self.sigma - self.theta
 
         if zn is None:
-            vol_t = self.vol_step(dt, self.sigma)
-            zn = self.rng_spawn[1].standard_normal(size=(self.n_sin, self.n_path))
+            vol_t = self.vol_step(dt, vol_0, nz_theta=False)
+            zn = self.rng_spawn[2].standard_normal(size=(self.n_sin, self.n_path))
             n_sin, n_path = self.n_sin, self.n_path
         else:
-            vol_t = self.vol_step(dt, self.sigma, zn[0, :])
+            vol_t = self.vol_step(dt, vol_0, zn[0, :], nz_theta=False)
             n_sin = zn.shape[0] - 1
 
-        sighat = vol_t - self.theta - x_0 * e_mr
+        volhat = vol_t - vol_0 * e_mr
 
         n_pi = np.pi*np.arange(1, n_sin + 1)
         an = np.sqrt(2/(mr_t**2 + n_pi**2))
         sin = np.sin(n_pi*tobs[:, None]/dt)
-        sigma_path = self.theta + x_0 * e_mr_tobs \
-                     + 0.5*(1/e_mr_tobs - e_mr_tobs)/np.sinh(mr_t) * sighat \
-                     + self.vov * np.sqrt(dt) * (an*sin) @ zn[1:,:]
+        sigma_path = self.theta + vol_0 * e_mr_tobs \
+                     + 0.5*(1/e_mr_tobs - e_mr_tobs)/np.sinh(mr_t) * volhat \
+                     + self.vov * np.sqrt(dt) * (an*sin) @ zn[1:, :]
 
         return sigma_path
diff --git a/pyfeng/sabr.py b/pyfeng/sabr.py
index 03d98c9..03dfb28 100644
--- a/pyfeng/sabr.py
+++ b/pyfeng/sabr.py
@@ -11,6 +11,7 @@
 import pandas as pd
 import numpy as np
 import scipy.optimize as spop
+from scipy import stats as spst
 
 from . import opt_smile_abc as smile
 from . import bsm
@@ -22,9 +23,10 @@ class SabrABC(smile.OptSmileABC, abc.ABC):
     vov, beta, rho = 0.0, 1.0, 0.0
     model_type = "SABR"
     var_process = False
+    #### vol_beta: the beta for the volatility to choose _m_vol. If None (by default) vol_beta = beta
     _base_beta = None
 
-    def __init__(self, sigma, vov=0.0, rho=0.0, beta=1.0, intr=0.0, divr=0.0, is_fwd=False):
+    def __init__(self, sigma, vov=0.1, rho=0.0, beta=1.0, intr=0.0, divr=0.0, is_fwd=False):
         """
         Args:
             sigma: model volatility at t=0
@@ -53,26 +55,24 @@ def _variables(self, fwd, texp):
         vovn = self.vov * np.sqrt(np.maximum(texp, 1e-64))
         return alpha, betac, rhoc, rho2, vovn
 
-    def base_model(self, vol, is_fwd=None):
+    def base_model(self, vol):
         """
         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.
 
         Args:
             vol: base model volatility
-            is_fwd: if True, treat `spot` as forward price. False by default.
 
         Returns: model
         """
-        base_beta = self._base_beta or self.beta
-        if is_fwd is None:
-            is_fwd = self.is_fwd
+        base_beta = self.beta if self._base_beta is None else self._base_beta
+
         if np.isclose(base_beta, 1):
-            return bsm.Bsm(vol, intr=self.intr, divr=self.divr, is_fwd=is_fwd)
+            return bsm.Bsm(vol, intr=self.intr, divr=self.divr, is_fwd=self.is_fwd)
         elif np.isclose(base_beta, 0):
-            return norm.Norm(vol, intr=self.intr, divr=self.divr, is_fwd=is_fwd)
+            return norm.Norm(vol, intr=self.intr, divr=self.divr, is_fwd=self.is_fwd)
         else:
-            return cev.Cev(vol, beta=base_beta, intr=self.intr, divr=self.divr, is_fwd=is_fwd)
+            return cev.Cev(vol, beta=base_beta, intr=self.intr, divr=self.divr, is_fwd=self.is_fwd)
 
     def vol_smile(self, strike, spot, texp, cp=1, model=None):
         if model is None:
@@ -161,6 +161,64 @@ def _hh(zz, rho):
 
         return 1.0 / xx_zz
 
+    @staticmethod
+    def cond_avgvar_mv(vovn, z, remove_exp=False):
+        """
+        Mean and M2 (2nd moment) of the average variance conditional on z = log(sigma_T/sigma_0) / (vov*sqrt(T))
+        int_0^1 exp{2 vovn Z_s - vovn^2 s} ds | Z_1 - (vovn/2) = z
+        True mean and M2 should be multiplied by sigma_0 and sigma_0^2, respectively.
+        See Appendix B in Choi et al. (2019), Eqs. (B5) and (B6) in Kennedy et al. (2012)
+
+        Args:
+            vovn: vov * sqrt(T)
+            z: log(sigma_T/sigma_0) / (vov*sqrt(T)) = Z_1 - (vovn/2)
+            remove_exp: if True, return without multiplying exp(vovn * z). False by default.
+
+        Returns:
+            Mean, M2
+
+        References:
+            - Choi J, Liu C, Seo BK (2019) Hyperbolic normal stochastic volatility model. J Futures Mark 39:186–204. https://doi.org/10.1002/fut.21967
+            - Kennedy JE, Mitra S, Pham D (2012) On the Approximation of the SABR Model: A Probabilistic Approach. Applied Mathematical Finance 19:553–586. https://doi.org/10.1080/1350486X.2011.646523
+        """
+
+        def mean(vv):
+            ncdf_diff = spst.norm.cdf(-np.abs(z) + vv) - spst.norm.cdf(-np.abs(z) - vv)
+            m = np.sqrt(np.pi/2) / vv * ncdf_diff * np.exp((z**2 + vv**2)/2)
+            return m
+
+        m1 = mean(vovn)
+        m2 = (mean(2*vovn) - m1*np.cosh(z*vovn))/vovn**2
+
+        if not remove_exp:
+            exp = np.exp(vovn * z)
+            m1 *= exp
+            m2 *= exp**2
+
+        return m1, m2
+
+    @staticmethod
+    def avgvar_mv(vovn):
+        """
+        mean and variance of the normalized average variance:
+        (1/T) \int_0^T e^{2*vov Z_t - vov^2 Z_t} = \int_0^1 e^{2 vovn Z_s - vovn^2 s} ds
+        where vovn = vov*sqrt(T). See p.2 in Choi & Wu (2021).
+
+        Args:
+            vovn: vov * sqrt(texp)
+
+        Returns:
+            (mean, variance)
+
+        References
+            - Choi J, Wu L (2021) A note on the option price and ‘Mass at zero in the uncorrelated SABR model and implied volatility asymptotics.’ Quantitative Finance 21:1083–1086. https://doi.org/10.1080/14697688.2021.1876908
+        """
+        vovn2 = vovn**2
+        ww = np.exp(vovn2)
+        m = np.where(vovn2 > 1e-6, (ww - 1)/vovn2, 1 + vovn2/2 * (1 + vovn2/3))
+        v = (10 + ww*(6 + ww*(3 + ww))) / 15 * m**3 * vovn2
+        return m, v
+
 
 class SabrVolApproxABC(SabrABC):
     """
@@ -177,12 +235,12 @@ def _int_inv_locvol(strike, beta, fwd=1.0):
             = F^(-beta)/(1-beta) * (k^(1-beta) - 1)/(k-1) where k = K/F
 
         Args:
-            strike:
-            beta:
-            fwd:
+            strike: strike price
+            beta: beta
+            fwd: forward price
 
         Returns:
-
+            integrated inv local vol
         """
 
         assert np.isscalar(beta)
@@ -192,14 +250,12 @@ def _int_inv_locvol(strike, beta, fwd=1.0):
 
         # fall-back indices and values first
         ind_atm = np.fabs(kk - 1.0) < 1e-6
-        val = 1 - beta / 2 * (kk - 1) * (1 - (1 + beta) / 3 * (kk - 1))
+        val = 1 - beta/2 * (kk - 1) * (1 - (1 + beta)/3 * (kk - 1))
         with np.errstate(divide="ignore", invalid="ignore"):
             if abs(betac) < 1e-4:
                 val = np.where(ind_atm, val, np.log(kk) / (kk - 1))
             else:
-                val = np.where(
-                    ind_atm, val, (np.power(kk, betac) - 1) / betac / (kk - 1)
-                )
+                val = np.where(ind_atm, val, (np.power(kk, betac) - 1) / betac / (kk - 1))
 
         return val
 
@@ -243,10 +299,10 @@ def vol_smile(self, strike, spot, texp, cp=1, model=None):
         if model is None:
             model = "norm" if np.isclose(self.beta, 0) else "bsm"
 
-        vol_beta = self._base_beta or self.beta
-        if (model.lower() == "bsm" and np.isclose(vol_beta, 1)) or (
-            model.lower() == "norm" and np.isclose(vol_beta, 0)
-        ):
+        vol_beta = self.beta if self._base_beta is None else self._base_beta
+
+        if (model.lower() == "bsm" and np.isclose(vol_beta, 1)) \
+                or (model.lower() == "norm" and np.isclose(vol_beta, 0)):
             vol = self.vol_for_price(strike, spot, texp)
         else:
             vol = super().vol_smile(strike, spot, texp, cp=cp, model=model)
@@ -335,18 +391,10 @@ def iv_func(x):
             return err
 
         sol = spop.root(iv_func, np.array([np.log(vol3[1]), -1, 0.0]))
-        params = {
-            "sigma": np.exp(sol.x[0]),
-            "vov": np.exp(sol.x[1]),
-            "rho": np.tanh(sol.x[2]),
-        }
+        params = {"sigma": np.exp(sol.x[0]), "vov": np.exp(sol.x[1]), "rho": np.tanh(sol.x[2])}
 
         if setval:
-            self.sigma, self.vov, self.rho = (
-                params["sigma"],
-                params["vov"],
-                params["rho"],
-            )
+            self.sigma, self.vov, self.rho = params["sigma"], params["vov"], params["rho"]
 
         return params
 
@@ -368,7 +416,7 @@ class SabrNormVolApprox(SabrVolApproxABC):
     _base_beta = 0  # should not be changed
     is_atmvol = False
 
-    def __init__(self, sigma, vov=0.0, rho=0.0, beta=None, intr=0.0, divr=0.0, is_fwd=False, is_atmvol=False):
+    def __init__(self, sigma, vov=0.1, rho=0.0, beta=None, intr=0.0, divr=0.0, is_fwd=False, is_atmvol=False):
         """
         Args:
             sigma: model volatility at t=0
@@ -386,6 +434,29 @@ def __init__(self, sigma, vov=0.0, rho=0.0, beta=None, intr=0.0, divr=0.0, is_fw
         self.is_atmvol = is_atmvol
         super().__init__(sigma, vov, rho, beta=0, intr=intr, divr=divr, is_fwd=is_fwd)
 
+    def price_vsk(self, texp=1):
+        """
+        Variance, skewness, and ex-kurtosis of forward price. Eq. (21) in Choi et al. (2019)
+        It is a special case (lambda=0) of NsvhABC.price_vsk()
+
+        Args:
+            texp: time to expiry
+
+        Returns:
+            (variance, skewness, and ex-kurtosis)
+
+        References:
+            - Choi J, Liu C, Seo BK (2019) Hyperbolic normal stochastic volatility model. J Futures Mark 39:186–204. https://doi.org/10.1002/fut.21967
+        """
+        vovn = self.vov * np.sqrt(texp)
+        ww = np.exp(vovn**2)
+
+        m2 = ww - 1
+        skew = self.rho * (ww+2) * np.sqrt(m2)
+        exkurt = m2 * ((4*self.rho**2 + 1)/5 * (ww**3 + 3*ww**2 + 6*ww + 5) + 1)
+
+        return m2 * (self.sigma/self.vov)**2, skew, exkurt
+
     def vol_for_price(self, strike, spot, texp):
         fwd = self.forward(spot, texp)
 
@@ -417,26 +488,11 @@ class SabrChoiWu2021H(SabrVolApproxABC, smile.MassZeroABC):
         array([22.04897988, 14.56240351,  8.74169054,  4.72340753,  2.28876105])
     """
 
-    def __init__(self, sigma, vov=0.0, rho=0.0, beta=1.0, intr=0.0, divr=0.0, is_fwd=False, vol_beta=None):
-        """
-        Args:
-            sigma: model volatility at t=0
-            vov: volatility of volatility
-            rho: correlation between price and volatility
-            beta: elasticity parameter. 1.0 by default
-            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.
-            vol_beta: the beta for the volatility to choose _m_vol. If None (by default) vol_beta = beta
-        """
-        self._base_beta = vol_beta
-        super().__init__(sigma, vov, rho, beta, intr, divr, is_fwd)
-
     def vol_for_price(self, strike, spot, texp):
         # fwd, spot, sigma may be either scalar or np.array.
         # texp, vov, rho, beta should be scholar values
 
-        vol_beta = self._base_beta or self.beta
+        vol_beta = self.beta if self._base_beta is None else self._base_beta
         vol_betac = 1.0 - vol_beta
 
         fwd, _, _ = self._fwd_factor(spot, texp)
@@ -528,7 +584,7 @@ def vol_for_price(self, strike, spot, texp):
         # fwd, spot, sigma may be either scalar or np.array.
         # texp, vov, rho, beta should be scholar values
 
-        vol_beta = self._base_beta or self.beta
+        vol_beta = self.beta if self._base_beta is None else self._base_beta
         vol_betac = 1.0 - vol_beta
 
         fwd, _, _ = self._fwd_factor(spot, texp)
@@ -545,18 +601,14 @@ def vol_for_price(self, strike, spot, texp):
 
         # qq_ratio = qq_vol_beta / qq_beta
         tmp = self._int_inv_locvol(kk, self.beta)
-        qq_ratio = (
-            1.0 if self._base_beta is None else self._int_inv_locvol(kk, vol_beta) / tmp
-        )
+        qq_ratio = 1.0 if self._base_beta is None else self._int_inv_locvol(kk, vol_beta) / tmp
 
         zz = vov_over_alpha_safe * tmp * (kk - 1.0)  # zeta = (vov/sigma0) q
         hh = self._hh(zz, self.rho)
         v_m = self._vv(zz, self.rho)
 
         if abs(betac) < 0.001:
-            gg_diff = (
-                0.5 * self.beta * self.rho / (1.0 - rho2) * (v_m - 1.0 - self.rho * zz)
-            )
+            gg_diff = 0.5 * self.beta * self.rho / (1.0 - rho2) * (v_m - 1.0 - self.rho * zz)
         else:
             t1 = (v_m + self.rho + zz) / rhoc  # array
             t2 = (1 + self.rho) / rhoc  # scalar
diff --git a/pyfeng/sabr_int.py b/pyfeng/sabr_int.py
index 070c191..413da74 100644
--- a/pyfeng/sabr_int.py
+++ b/pyfeng/sabr_int.py
@@ -2,7 +2,6 @@
 import numpy as np
 from . import sabr
 import scipy.special as spsp
-import scipy.stats as spst
 from . import opt_smile_abc as smile
 
 
@@ -10,38 +9,22 @@ class SabrMixtureABC(sabr.SabrABC, smile.MassZeroABC, abc.ABC):
 
     correct_fwd = False
 
-    @staticmethod
-    def avgvar_mv(vovn):
-        """
-        mean and variance of the normalized average variance:
-        (1/T) \int_0^T e^{2*vov Z_t - vov^2 Z_t} = \int_0^1 e^{2 vovn Z_s - vovn^2 s} ds
-        where vovn = vov*sqrt(T)
-
-        Args:
-            vovn: vov * sqrt(texp)
-
-        Returns:
-            (mean, variance)
-        """
-        v2 = vovn**2
-        w = np.exp(v2)
-        m = np.where(v2 > 1e-6, (w - 1)/v2, 1 + v2/2 * (1 + v2/3))
-        v = (10 + w*(6 + w*(3 + w))) / 15 * m**3 * v2
-        return m, v
-
     @staticmethod
     def avgvar_lndist(vovn):
         """
         Lognormal distribution parameters (mean, sigma) of the normalized average variance:
         (1/T) \int_0^T e^{2*vov Z_t - vov^2 t} dt = \int_0^1 e^{2 vovn Z_s - vovn^2 s} ds
-        where vovn = vov*sqrt(T)
+        where vovn = vov*sqrt(T). See p.2 in Choi & Wu (2021).
 
         Args:
             vovn: vov * sqrt(texp)
 
         Returns:
             (m1, sig)
-            True distribution should be multiplied by sigma^2*t
+            True distribution should be multiplied by sigma^2 * texp
+
+        References
+            - Choi J, Wu L (2021) A note on the option price and ‘Mass at zero in the uncorrelated SABR model and implied volatility asymptotics.’ Quantitative Finance 21:1083–1086. https://doi.org/10.1080/14697688.2021.1876908
         """
         v2 = vovn**2
         w = np.exp(v2)
@@ -50,37 +33,6 @@ def avgvar_lndist(vovn):
         sig = np.sqrt(np.where(v2 > 1e-8, np.log(m2m1ratio), 4 / 3 * v2))
         return m1, sig
 
-    @staticmethod
-    def cond_avgvar_mv(vovn, z, remove_exp=False):
-        """
-        Mean and M2 (2nd moment) of the average variance conditional on z = log(sigma_T/sigma_0) / (vov*sqrt(T))
-        int_0^1 exp{2 vovn Z_s - vovn^2 s} ds | Z_1 - (vovn/2) = z
-        True mean and M2 should be multiplied by sigma_0 and sigma_0^2, respectively
-
-        Args:
-            vovn: vov * sqrt(T)
-            z: log(sigma_T/sigma_0) / (vov*sqrt(T)) = Z_1 - (vovn/2)
-            remove_exp: if True, return without multiplying exp(vovn * z). False by default.
-
-        Returns:
-            Mean, M2
-        """
-
-        def mean(vv):
-            ncdf_diff = spst.norm.cdf(-np.abs(z) + vv) - spst.norm.cdf(-np.abs(z) - vv)
-            m = np.sqrt(np.pi/2) / vv * ncdf_diff * np.exp((z**2 + vv**2)/2)
-            return m
-
-        m1 = mean(vovn)
-        m2 = (mean(2*vovn) - m1*np.cosh(z*vovn))/vovn**2
-
-        if not remove_exp:
-            exp = np.exp(vovn * z)
-            m1 *= exp
-            m2 *= exp**2
-
-        return m1, m2
-
     @abc.abstractmethod
     def cond_spot_sigma(self, texp, fwd):
         # return (fwd, vol, weight) each 1d array
@@ -112,7 +64,8 @@ def price(self, strike, spot, texp, cp=1):
         vol_ratio = np.expand_dims(vol_ratio[ind], -1)
         ww = np.expand_dims(ww[ind], -1)
 
-        base_model = self.base_model(alpha * vol_ratio, is_fwd=True)
+        base_model = self.base_model(alpha * vol_ratio)
+        base_model.is_fwd = True
         price_vec = base_model.price(kk, fwd_ratio, texp, cp=cp)
         price = fwd * np.sum(price_vec * ww, axis=0)
         return price
@@ -127,16 +80,17 @@ def mass_zero(self, spot, texp, log=False, mu=0):
             fwd_ratio /= np.sum(fwd_ratio*ww)
         assert np.isclose(np.sum(ww), 1)
 
-        base_m = self.base_model(alpha * vol_ratio, is_fwd=True)
+        base_model = self.base_model(alpha * vol_ratio)
+        base_model.is_fwd = True
 
         if log:
-            log_mass = np.log(ww) + base_m.mass_zero(fwd_ratio, texp, log=True)
+            log_mass = np.log(ww) + base_model.mass_zero(fwd_ratio, texp, log=True)
             log_max = np.amax(log_mass)
             log_mass -= log_max
             log_mass = log_max + np.log(np.sum(np.exp(log_mass)))
             return log_mass
         else:
-            mass = base_m.mass_zero(fwd_ratio, texp, log=False)
+            mass = base_model.mass_zero(fwd_ratio, texp, log=False)
             mass = np.sum(mass * ww)
             return mass
 
diff --git a/pyfeng/sabr_mc.py b/pyfeng/sabr_mc.py
index e2b42d3..b42fbb3 100644
--- a/pyfeng/sabr_mc.py
+++ b/pyfeng/sabr_mc.py
@@ -8,43 +8,64 @@
 import scipy.integrate as scint
 from . import sabr
 from . import sv_abc as sv
-from . import cev
 
 
 class SabrMcABC(sabr.SabrABC, sv.CondMcBsmABC, abc.ABC):
 
-    def vol_step(self, dt, sigma_0):
+    def vol_step(self, dt, log=False):
         """
-        Stepping volatility according to OU process dynamics
+        SABR sigma after dt according to the volatility dynamics (GBM).
+        Because of the multiplicative property of sigma, we assume sigma_0 = 1
 
         Args:
             dt: time step
-            sigma_0: initial volatility
-
+            log: if True, return log(sigma). False by default.
         Returns:
-            sigma after dt
+            sigma_dt
         """
         vovn = self.vov * np.sqrt(dt)
         zz = self.rv_normal(spawn=0)
-        sigma_t = sigma_0 * np.exp(vovn * (zz - vovn/2))
 
-        return sigma_t
+        if log:
+            return vovn * (zz - vovn/2)
+        else:
+            return np.exp(vovn * (zz - vovn/2))
 
     @abc.abstractmethod
     def cond_states_step(self, dt, sigma_0):
         """
-        Final variance and integrated variance over dt given var_0
-        The int_var is normalized by dt
+        Final variance after dt and average variance over (0, dt) given sigma_0.
+        `sigma_0` should be an array of (self.n_path, )
 
         Args:
             dt: time step
-            sigma_0: initial sigma
+            sigma_0: initial volatility
 
         Returns:
-            (sigma_t, avgvar)
+            (sigma after dt, average variance during dt)
         """
         return NotImplementedError
 
+    def draw_log_return(self, dt, sigma_0, sigma_t, avgvar):
+        """
+        Samples log return, log(S_t/S_0). Currently implemented only for beta=1
+
+        Args:
+            dt: time step
+            sigma_0: initial variance
+            sigma_t: final variance
+            avgvar: average variance
+
+        Returns:
+            log return (self.n_path, )
+        """
+
+        assert np.isclose(self.beta, 1.0)
+        ln_m = (self.intr - self.divr)*dt + self.rho/self.vov*(sigma_t - sigma_0) - 0.5*avgvar*dt
+        ln_sig = np.sqrt((1.0 - self.rho**2) * dt * avgvar)
+        zn = self.rv_normal(spawn=5)
+        return ln_m + ln_sig * zn
+
     def cond_spot_sigma(self, texp, fwd, mu=0):
         """
         Spot and sigma ratio.
@@ -57,7 +78,6 @@ def cond_spot_sigma(self, texp, fwd, mu=0):
         Returns:
             (spot ratio, sigma ratio)
         """
-
         tobs = self.tobs(texp)
         dt = np.diff(tobs, prepend=0)
         n_dt = len(dt)
@@ -102,7 +122,8 @@ def price(self, strike, spot, texp, cp=1):
         fwd_ratio = np.expand_dims(fwd_ratio[ind], -1)
         vol_ratio = np.expand_dims(vol_ratio[ind], -1)
 
-        base_model = self.base_model(alpha * vol_ratio, is_fwd=True)
+        base_model = self.base_model(alpha * vol_ratio)
+        base_model.is_fwd = True
         price_vec = base_model.price(kk, fwd_ratio, texp, cp=cp)
         price = fwd * np.sum(price_vec, axis=0) / self.n_path
 
@@ -116,8 +137,9 @@ class SabrMcTimeDisc(SabrMcABC):
     """
 
     def cond_states_step(self, dt, sigma_0):
-        sigma_t = self.vol_step(dt, sigma_0)
-        avgvar = (sigma_0**2 + sigma_t**2) / 2
+        sigma_t = self.vol_step(dt)
+        avgvar = sigma_0**2 * (1.0 + sigma_t**2) / 2
+        sigma_t *= sigma_0
 
         return sigma_t, avgvar
 
@@ -173,7 +195,6 @@ def cond_spot_sigma_volpath(self, texp, sigma_0, mu=0):
 
         return spot_cond, vol_cond, log_rn_deriv
 
-
     def mass_zero(self, spot, texp, log=False, mu=0):
 
         assert 0 < self.beta < 1
@@ -194,7 +215,9 @@ def mass_zero(self, spot, texp, log=False, mu=0):
         assert np.isclose(fwd_ratio, 1.0).all()
         log_rn_deriv = 0.0  ## currently not used
 
-        base_model = cev.Cev(sigma=alpha * vol_ratio, beta=self.beta)
+        base_model = self.base_model(alpha * vol_ratio)
+        base_model.is_fwd = True
+
         if log:
             log_mass_grid = base_model.mass_zero(1.0, texp, log=True) + log_rn_deriv
             log_mass_max = np.amax(log_mass_grid)
@@ -272,7 +295,7 @@ def inv_laplace(self, u, vovn, sigma_t):
         Return the original function from transform function
 
         Args:
-            u: original variable
+            u: dummy variable
             vovn: vov * sqrt(texp)
             sigma_t: final volatility
 
@@ -281,44 +304,63 @@ def inv_laplace(self, u, vovn, sigma_t):
         """
 
         ## index from 0 to m + n
-        ss_j = self.cond_laplace((self.m_inv - 2j * np.pi * self.nn) / (2*u), vovn, sigma_t).real
-        term1 = 0.5 * ss_j[0]
-        ss_j[1::2] *= -1
-        np.cumsum(ss_j, out=ss_j)
-        term2 = np.sum(self.comb_coef * ss_j[self.n_euler:])
+        ss_j = self.cond_laplace((self.m_inv - 2j * np.pi * self.nn[:, None]) / (2*u), vovn, sigma_t).real
+        term1 = 0.5 * ss_j[0, :]
+        ss_j[1::2, :] *= -1
+        np.cumsum(ss_j, axis=0, out=ss_j)
+        term2 = np.sum(self.comb_coef[:, None] * ss_j[self.n_euler:, :], axis=0)
 
         origin_L = np.exp(self.m_inv/2) / u * (-term1 + term2)
 
+        if np.isscalar(u):
+            origin_L = origin_L[0]
+
         return origin_L
 
-    def draw_cond_avgvar(self, dt, sigma_0, sigma_final):
+    def draw_cond_avgvar(self, dt, sigma_t):
         """
-        Normalized integraged variance samples.
+        Draw normalized average variance given sigma_t (and sigma_0 = 1)
+
+        int_0^1 exp{2 vovn Z_s - vovn^2 s} ds | exp(Z_1 - vovn/2) = sigma_t
 
         Args:
             dt: time step
-            sigma_final: final volatility
+            sigma_t: final sigma given sigma_0=1
 
         Returns:
-            (n_path, 1) array
+            (n_path, ) array
         """
         vovn = self.vov * np.sqrt(dt)
+        uu = self.rv_uniform(spawn=2)
+        avgvar = np.zeros_like(sigma_t)
 
-        u_rn = self.rv_uniform(spawn=2)
-        int_var = np.zeros(self.n_path)
-
-        for i in range(self.n_path):
-            obj_func = lambda x: self.inv_laplace(x, vovn, sigma_final[i]) - u_rn[i]
-
-            sol = spop.brentq(obj_func, 0.000001, 100)
-            int_var[i] = 1 / sol
+        for i in range(len(sigma_t)):
+            obj_func = lambda x: self.inv_laplace(x, vovn, sigma_t[i]) - uu[i]
+            sol = spop.brentq(obj_func, 1e-7, 1e5)
+            avgvar[i] = 1 / sol
 
-        return int_var
+        """
+        ## Vectorized newton method, but doesn't work well
+        
+        zz = self.rv_normal(spawn=2)
+        uu = spst.norm.cdf(zz)
+        ln_m, m2 = self.cond_avgvar_mv(vovn, np.log(sigma_t) / vovn)
+        ln_sig = np.sqrt(np.log(m2 / ln_m ** 2))
+        print(sigma_t, '\n', ln_sig, '\n', uu)
+
+        def obj_func(z):
+            x = ln_m * np.exp(ln_sig * (z - ln_sig / 2))
+            return self.inv_laplace(x, vovn, sigma_t) - uu
+
+        avgvar = 1 / spop.newton(obj_func, zz)
+        """
+        return avgvar
 
     def cond_states_step(self, dt, sigma_0):
 
-        sigma_t = self.vol_step(dt, sigma_0)
-        avgvar = self.draw_cond_avgvar(dt, sigma_0, sigma_t)
+        sigma_t = self.vol_step(dt)  ## ratio
+        avgvar = sigma_0**2 * self.draw_cond_avgvar(dt, sigma_t)
+        sigma_t *= sigma_0
 
         return sigma_t, avgvar
 
diff --git a/pyfeng/sv32_mc2.py b/pyfeng/sv32_mc2.py
index 8c8cbd9..42660df 100644
--- a/pyfeng/sv32_mc2.py
+++ b/pyfeng/sv32_mc2.py
@@ -234,7 +234,7 @@ def cond_avgvar_laplace(self, bb, dt, var_0, var_t, eta=None):
 
         return ret
 
-    def draw_avgvar(self, dt, var_0, var_t):
+    def draw_cond_avgvar(self, dt, var_0, var_t):
         def laplace_cond(bb):
             return self.cond_avgvar_laplace(bb, dt, var_0, var_t)
 
@@ -263,7 +263,7 @@ def laplace_cond(bb):
         phimat = laplace_cond(-1j * jj * h).real
 
         # Sample the conditional integrated variance by inverse transform sampling
-        zz = self.rv_normal()
+        zz = self.rv_normal(spawn=0)
         uu = spst.norm.cdf(zz)
 
         def root(xx):
@@ -289,7 +289,7 @@ def cond_states_step(self, dt, var_0):
         var_t = self._m_heston.var_step_ncx2(dt, 1 / var_0)
         np.divide(1.0, var_t, out=var_t)
 
-        avgvar = self.draw_avgvar(dt, var_0, var_t)
+        avgvar = self.draw_cond_avgvar(dt, var_0, var_t)
 
         return var_t, avgvar
 
diff --git a/pyfeng/sv_abc.py b/pyfeng/sv_abc.py
index 395c314..b8650c3 100644
--- a/pyfeng/sv_abc.py
+++ b/pyfeng/sv_abc.py
@@ -112,8 +112,7 @@ class CondMcBsmABC(smile.OptSmileABC, abc.ABC):
     rng = np.random.default_rng(None)
     rng_spawn = []
     antithetic = True
-
-    correct_fwd = True
+    correct_fwd = False
     result = {}
 
     def set_num_params(self, n_path=10000, dt=0.25, rn_seed=None, antithetic=True):
@@ -187,7 +186,8 @@ def _bm_incr(self, tobs, cum=False, n_path=None):
         dt = np.diff(np.atleast_1d(tobs), prepend=0)
         n_dt = len(dt)
 
-        n_path = n_path or self.n_path
+        if n_path is None:
+            n_path = self.n_path
 
         if self.antithetic:
             # generate random number in the order of (path, time) first and transposed
diff --git a/pyfeng/sv_fft.py b/pyfeng/sv_fft.py
index a3aeb9c..488871b 100644
--- a/pyfeng/sv_fft.py
+++ b/pyfeng/sv_fft.py
@@ -12,6 +12,7 @@
 from . import heston
 from . import rheston
 
+
 class FftABC(opt.OptABC, abc.ABC):
     n_x = 2**12  # number of grid. power of 2 for FFT
     x_lim = 200  # integratin limit
@@ -365,9 +366,7 @@ def mgf_logprice_schobelzhu1998(self, uu, texp):
         s2g3 = vov**2*gamma1**3
 
         D = (mr - gamma1*sincos/cossin)/vov**2
-        B = ((ktg3 + gamma3*sincos)/cossin - mr*theta*gamma1)/(
-                vov**2*gamma1
-        )
+        B = ((ktg3 + gamma3*sincos)/cossin - mr*theta*gamma1)/(vov**2*gamma1)
         C = (
                 -0.5*np.log(cossin)
                 + 0.5*mr*texp
diff --git a/pyfeng/svi.py b/pyfeng/svi.py
new file mode 100644
index 0000000..bd18b77
--- /dev/null
+++ b/pyfeng/svi.py
@@ -0,0 +1,84 @@
+import numpy as np
+from . import opt_smile_abc as smile
+from . import bsm
+
+
+class Svi(smile.OptSmileABC):
+    """
+    Stochastic Volatility-inspired (SVI) model by Gatheral.
+
+    References
+        - Gatheral J, Jacquier A (2013) Arbitrage-free SVI volatility surfaces. arXiv:12040646 [q-fin]
+    """
+
+    vov, rho, smooth, shift = 0.4, -0.4, 0.1, 0.0
+
+    def __init__(self, sigma=0.04, vov=0.4, rho=-0.4, smooth=0.1, shift=0.0, intr=0.0, divr=0.0, is_fwd=False):
+        """
+        Raw SVI parametrization
+
+        Args:
+            sigma: level (a)
+            vov: vol-of-vol (b)
+            rho: rotation (rho)
+            smooth: smoothness (sigma)
+            shift: translation (m)
+            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.
+        """
+
+        self.sigma, self.vov, self.rho, self.smooth, self.shift = sigma, vov, rho, smooth, shift
+        super().__init__(sigma, intr=intr, divr=divr, is_fwd=is_fwd)
+
+    def base_model(self, sigma=None):
+        base_model = bsm.Bsm(sigma, intr=self.intr, divr=self.divr, is_fwd=self.is_fwd)
+        return base_model
+
+    def vol_for_price(self, strike, spot, texp):
+        fwd = self.forward(spot, texp)
+        money = np.log(strike / fwd) - self.shift
+        vol = np.sqrt(self.sigma + self.vov * (self.rho * money + np.sqrt(money**2 + self.smooth**2)))
+        return vol
+
+    @classmethod
+    def init_from_heston(cls, sigma, vov=0.8, rho=-0.7, mr=0.5, theta=None, texp=1.0, intr=0.0, divr=0.0, is_fwd=False):
+        """
+        SVI initalization with equivalent Heston model by Gatheral & Jacquier (2011)
+
+        Args:
+            sigma: Heston sigma
+            vov: Heston vov
+            rho: Heston rho
+            mr: Heston mr
+            theta: Heston theta
+            texp: time to expiry
+            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.
+
+        Returns: model
+
+        References:
+            - Gatheral J, Jacquier A (2011) Convergence of Heston to SVI. Quantitative Finance 11:1129–1132. https://doi.org/10.1080/14697688.2010.550931
+        """
+
+        theta = sigma if theta is None else theta
+
+        rhoc2 = 1 - rho**2
+        w1 = np.sqrt((2*mr - rho*vov)**2 + vov**2 * rhoc2) - (2*mr - rho*vov)
+        w1 *= (4*mr*theta) / (vov**2 * rhoc2)
+        w2 = sigma / (mr * theta)
+        sigma_ = w1 * rhoc2 / 2
+        vov_ = (w1 * w2) / (2 * texp)
+        shift = - rho * texp / w2
+        smooth = np.sqrt(rhoc2) * texp / w2
+        m = cls(sigma_, vov=vov_, rho=rho, smooth=smooth, shift=shift, intr=intr, divr=divr, is_fwd=is_fwd)
+
+        return m
+
+    def price(self, strike, spot, texp, cp=1):
+        vol = self.vol_for_price(strike, spot, texp)
+        m_vol = self.base_model(vol)
+        price = m_vol.price(strike, spot, texp, cp=cp)
+        return price
diff --git a/pyfeng/util.py b/pyfeng/util.py
new file mode 100644
index 0000000..a3772e9
--- /dev/null
+++ b/pyfeng/util.py
@@ -0,0 +1,17 @@
+import numpy as np
+
+def avg_exp(x):
+    """
+    Integral_0^x exp(x) dx / x = ( exp(x) - 1 ) / x
+
+    Args:
+        x: argument
+
+    Returns:
+        value
+    """
+    with np.errstate(invalid="ignore"):
+        rv = np.where(np.abs(x) < 1e-5,
+                      1 + (x/2)*(1 + (x/3)*(1 + (x/4))),
+                      np.expm1(x)/x)
+    return rv
diff --git a/setup.py b/setup.py
index 1447206..96d03ba 100644
--- a/setup.py
+++ b/setup.py
@@ -5,7 +5,7 @@
 
 setup(
     name="pyfeng",
-    version="0.2.3",
+    version="0.2.4",
     description="Python Financial Engineering",
     long_description=long_description,
     long_description_content_type="text/markdown",


Test patch:
diff --git a/tests/test_heston.py b/tests/test_heston.py
index a963126..1bf5e76 100644
--- a/tests/test_heston.py
+++ b/tests/test_heston.py
@@ -89,6 +89,28 @@ def test_price_mc(self):
             np.testing.assert_allclose(vol0, vol1, atol=5e-3)
             np.testing.assert_allclose(m.result['spot error'], 0, atol=2e-3)
 
+    def test_avgvar_mv(self):
+        """
+        mean and variance of var_t and average variance
+        """
+        #sigma, vov, mr, rho, texp, spot = 0.04, 1, 0.5, -0.9, 10, 100
+        m, *_ = pf.HestonMcChoiKwok2023PoisGe.init_benchmark(1)
+        m.set_num_params(n_path=32e4, rn_seed=123456, kk=8)
+
+        for texp in (1.0, 3.0, 5.0):
+            # analytic mean and variance
+            var_m, var_v = m.var_mv(texp)
+            # analytic mean and average variance
+            avgvar_m, avgvar_v = m.avgvar_mv(texp)
+            # MC samples of variance and avgvar
+            var_t, avgvar, *_ = m.cond_states_step(texp, m.sigma)
+
+            np.testing.assert_allclose(np.mean(var_t), var_m, rtol=0.02)
+            np.testing.assert_allclose(np.var(var_t), var_v, rtol=0.02)
+
+            np.testing.assert_allclose(np.mean(avgvar), avgvar_m, rtol=0.02)
+            np.testing.assert_allclose(np.var(avgvar), avgvar_v, rtol=0.02)
+
 
 if __name__ == "__main__":
     print(f"Pyfeng loaded from {pf.__path__}")
diff --git a/tests/test_ousv.py b/tests/test_ousv.py
index ac199a9..8b13fbf 100644
--- a/tests/test_ousv.py
+++ b/tests/test_ousv.py
@@ -1,12 +1,11 @@
 import unittest
 import numpy as np
-import scipy.integrate as scin
+import scipy.special as spsp
 import sys
 import os
 
 sys.path.insert(0, os.getcwd())
 import pyfeng as pf
-import pyfeng.ex as pfex
 
 
 class TestOusvKL(unittest.TestCase):
@@ -29,7 +28,7 @@ def test_infinite_sum(self):
             An2n2 = An2 / n_pi_2
             An6n2 = n_pi_2 * An6
 
-            m = pfex.OusvMcChoi2023KL(1, 1, 0, 1, theta=mrt)
+            m = pf.OusvMcChoi2023KL(1, 1, 0, 1, theta=mrt)
 
             assert np.isclose(1, sum(An2[ns:])/m._a2sum(mrt, ns=ns), atol=1e-4)
             assert np.isclose(1, sum(An2[ns::2])/m._a2sum(mrt, ns=ns, odd=1), atol=1e-4)
@@ -51,6 +50,7 @@ def test_infinite_sum(self):
             assert np.isclose(1, sum(An6n2[ns::2])/m._a6n2sum(mrt, ns=ns, odd=1))
             assert np.isclose(1, sum(An6n2[ns+1::2])/m._a6n2sum(mrt, ns=ns, odd=2))
 
+
     def test_integarl_sin_path(self):
         """
         Vol/var integrated with Simpson method should equal to the analytic integral
@@ -59,7 +59,7 @@ def test_integarl_sin_path(self):
 
         """
         sheet_no = 1
-        m, p, rv = pfex.OusvMcChoi2023KL.init_benchmark(sheet_no)
+        m, p, rv = pf.OusvMcChoi2023KL.init_benchmark(sheet_no)
         n_sin = 6
         n_path = 50
         m.set_num_params(n_path=n_path, rn_seed=123456, n_sin=4, dt=None)
@@ -70,14 +70,44 @@ def test_integarl_sin_path(self):
             zn = np.random.normal(size=(n_sin + 1, n_path))
             vol_path = m.vol_path_sin(t_grid * texp, zn)
 
-            avgvol_simp = scin.simps(vol_path, dx=1, axis=0) / n_step
-            avgvar_simp = scin.simps(vol_path**2, dx=1, axis=0) / n_step
+            avgvol_simp = np.trapz(vol_path, dx=1, axis=0) / n_step
+            avgvar_simp = np.trapz(vol_path**2, dx=1, axis=0) / n_step
             sig_t, avgvar, avgvol = m.cond_states_step(texp, m.sigma, zn=zn)
 
             assert np.all(np.isclose(1, vol_path[-1, :]/sig_t))
             assert np.all(np.isclose(1, avgvol_simp/avgvol))
             assert np.all(np.isclose(1, avgvar_simp/avgvar))
 
+    def test_MomentsIntVariance(self):
+        """
+        Unconditional mean/var == E(conditional)
+        """
+        m = pf.OusvMcChoi2023KL(sigma=1, vov=0.75, mr=2.5)
+        zz, ww = spsp.roots_hermitenorm(31)
+        ww /= np.sqrt(2 * np.pi)
+
+        for texp in (1.0, 3.0, 5.0):
+            for sigma0 in m.sigma * np.array([0.5, 1, 2]):
+                ### Unconditional mean of avgvol and avgvar
+                mvol, _ = m.avgvol_mv(texp, sigma0)
+                mvar, _ = m.avgvar_mv(texp, sigma0)
+
+                ### Weighted average of conditional mean of avgvol and avgvar
+                sig_t = m.vol_step(texp, sigma0, zz)
+                cvol, cvar = m.cond_avgvolvar_m(texp, sigma0, sig_t)
+
+                np.testing.assert_allclose(mvol, np.sum(cvol * ww))
+                np.testing.assert_allclose(mvar, np.sum(cvar * ww))
+
+                mvol, _ = m.avgvol_mv(texp, sigma0 - m.theta, nz_theta=True)
+                mvar, _ = m.avgvar_mv(texp, sigma0 - m.theta, nz_theta=True)
+
+                sig_t = m.vol_step(texp, sigma0 - m.theta, zz, nz_theta=True)
+                cvol, cvar = m.cond_avgvolvar_m(texp, sigma0 - m.theta, sig_t, nz_theta=True)
+
+                np.testing.assert_allclose(mvol, np.sum(cvol * ww))
+                np.testing.assert_allclose(mvar, np.sum(cvar * ww))
+
 
 if __name__ == "__main__":
     print(f"Pyfeng loaded from {pf.__path__}")
diff --git a/tests/test_sabr.py b/tests/test_sabr.py
index 17c0d02..15b7252 100644
--- a/tests/test_sabr.py
+++ b/tests/test_sabr.py
@@ -3,6 +3,7 @@
 import numpy as np
 import sys
 import os
+import scipy.special as spsp
 
 sys.path.insert(0, os.getcwd())
 import pyfeng as pf
@@ -10,63 +11,118 @@
 
 class TestSabr(unittest.TestCase):
     def test_Hagan2002(self):
+        """
+        Hagan formula == Saved benchmark
+        """
         for k in list(range(1, 19)) + [22, 23]:
             m, df, rv = pf.SabrHagan2002.init_benchmark(k)
             # ref = rv['ref']
             # print(f'Sheet {k:02d}: {ref}')
             v1 = np.round(m.vol_for_price(**rv["args_pricing"]), 4)
             v2 = df["IV Hagan"].values
-            np.testing.assert_almost_equal(v1, v2)
+            np.testing.assert_allclose(v1, v2)
 
     def test_SabrNorm(self):
+        """
+        Choi & Wu (2021) Hagan == Normal Vol Approx
+        """
         for k in [22, 23]:
             m, df, rv = pf.SabrNormVolApprox.init_benchmark(k)
             v1 = m.price(**rv["args_pricing"])
             m, df, rv = pf.SabrChoiWu2021H.init_benchmark(k)
             v2 = m.price(**rv["args_pricing"])
-            np.testing.assert_almost_equal(v1, v2)
+            np.testing.assert_allclose(v1, v2)
 
     def test_SabrNormATM(self):
+        """
+        Test if is_atmvol works fine
+        """
         for k in [22, 23]:
             m, df, rv = pf.SabrNormVolApprox.init_benchmark(k)
             m.is_atmvol = True
-            np.testing.assert_almost_equal(m.vol_smile(0, 0, texp=0.1), m.sigma)
-            np.testing.assert_almost_equal(m.vol_smile(0, 0, texp=10), m.sigma)
+            np.testing.assert_allclose(m.vol_smile(0, 0, texp=0.1), m.sigma)
+            np.testing.assert_allclose(m.vol_smile(0, 0, texp=10), m.sigma)
 
             m, df, rv = pf.Nsvh1.init_benchmark(k)
             m.is_atmvol = True
-            np.testing.assert_almost_equal(m.vol_smile(0, 0, texp=0.1), m.sigma)
-            np.testing.assert_almost_equal(m.vol_smile(0, 0, texp=10), m.sigma)
+            np.testing.assert_allclose(m.vol_smile(0, 0, texp=0.1), m.sigma)
+            np.testing.assert_allclose(m.vol_smile(0, 0, texp=10), m.sigma)
 
     def test_PaulotBsm(self):
+        """
+        Paulot formula == Saved benchmark
+        """
         for k in list(range(1, 19)):
             m, df, rv = pf.SabrChoiWu2021P.init_benchmark(k)
             m._base_beta = 1.0  # For Paulot's BS volatility approximation
             # print(f'Sheet {k:02d}: {ref}')
             v1 = np.round(m.vol_for_price(**rv["args_pricing"]), 4)
             v2 = df["IV HL-P"].values
-            np.testing.assert_almost_equal(v1, v2)
+            np.testing.assert_allclose(v1, v2)
 
     def test_UnCorrChoiWu2021(self):
+        """
+        Uncorrelated SABR
+        """
         # Param Set 19: Table 7 (Case III.C) in Cai et al. (2017). https://doi.org/10.1287/opre.2017.1617
         m, df, rv = pf.SabrUncorrChoiWu2021.init_benchmark(19)
         mass = m.mass_zero(rv["args_pricing"]["spot"], rv["args_pricing"]["texp"])
         p = m.price(**rv["args_pricing"])
 
         mass2 = 0.7623543217183134
-        p2 = np.array(
-            [0.04533777, 0.04095806, 0.03889591, 0.03692339, 0.03324944, 0.02992918]
-        )
+        p2 = np.array([0.04533777, 0.04095806, 0.03889591, 0.03692339, 0.03324944, 0.02992918])
 
-        np.testing.assert_almost_equal(mass, mass2)
-        np.testing.assert_almost_equal(p, p2)
+        np.testing.assert_allclose(mass, mass2, atol=1e-8)
+        np.testing.assert_allclose(p, p2, atol=1e-8)
 
     def test_McTimeDisc(self):
+        """
+        Time discretization
+        """
         for k in [19, 20]:  # can test 22 (Korn&Tang) also, but difficult to pass
             m, df, rv = pf.SabrMcTimeDisc.init_benchmark(k)
             m.set_num_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)
+            np.testing.assert_allclose(p, rv["val"], rtol=5e-4)
+
+    def test_MomentsIntVariance(self):
+        """
+        Various test on momoents of SABR/NSVh
+        """
+
+        #### Unconditional mean/var == E(conditional)
+        m = pf.SabrNormVolApprox(1)
+        for vovn in [1.1, 1.2, 1.4, 1.6]:
+            zhat, ww = spsp.roots_hermitenorm(31)
+            ww /= np.sqrt(2*np.pi)
+            zhat -= 0.5*vovn
+
+            m1, v = m.avgvar_mv(vovn)
+            cond_m1, cond_m2 = m.cond_avgvar_mv(vovn, zhat, False)
+
+            np.testing.assert_allclose(np.sum(cond_m1 * ww), m1)
+            np.testing.assert_allclose(np.sum(cond_m2 * ww), (m1**2 + v))
+
+        #### Generic Nsvh (lambda = 0) == Normal SABR
+        #### Generic Nsvh (lambda = 1) == Nsvh1
+        for (rho, vov) in zip((-0.2, 0, 0.5), (0.1, 0.2, 0.5)):
+            m = pf.NsvhGaussQuad(1, rho=rho, vov=vov, lam=0)
+            m0 = pf.SabrNormVolApprox(1, rho=rho, vov=vov)
+            np.testing.assert_allclose(m0.price_vsk(texp=1.5), m.price_vsk(texp=1.5))
+
+            m = pf.NsvhGaussQuad(1, rho=rho, vov=vov, lam=1)
+            m1 = pf.Nsvh1(1, rho=rho, vov=vov)
+            np.testing.assert_allclose(m1.price_vsk(texp=1.5), m.price_vsk(texp=1.5))
+
+        for vovn in (0.1, 0.2, 0.5, 1.2):
+            m0 = pf.SabrNormVolApprox(sigma=vovn, rho=0, vov=vovn)
+            p_var, skew, kurt = m0.price_vsk(texp=1)
+            i_m, i_var = m0.avgvar_mv(vovn)
+
+            ### E(X^2_IntVar) = vovn^2 * E(I)
+            np.testing.assert_allclose(p_var, i_m*vovn**2)
+            ### E(X^4_IntVar) = 3 * vovn^4 * E(I^2)
+            np.testing.assert_allclose((kurt + 3)*p_var**2, 3*(i_var + i_m**2)*vovn**4)
 
 
 if __name__ == "__main__":