Merge pull request #141 from PyFE/heston-mc

Various changes in the SV/SABR models.

Author: jaehyukchoi

docs/source/index.rst

@@ -15,6 +15,7 @@
    nsvh
    sabr
    sabr_int
+   sv_fft
    gamma
    multiasset
    multiasset_mc

docs/source/nsvh.rst

@@ -1,5 +1,5 @@
 Hyperbolic Normal Stochastic Volatility (NSVh) Model
 =======================================
 .. automodule:: pyfeng.nsvh
-   :members: Nsvh1
+   :members: Nsvh1, NsvhQuadInt, NsvhMc
    :noindex:

docs/source/sabr.rst

@@ -1,7 +1,7 @@
 Stochastic-Alpha-Beta-Rho (SABR) Model
 ==========================
 .. automodule:: pyfeng.sabr
-   :members: SabrHagan2002, SabrLorig2017, SabrChoiWu2021H, SabrChoiWu2021P, SabrUncorrChoiWu2021
+   :members: SabrHagan2002, SabrNormVolApprox, SabrLorig2017, SabrChoiWu2021H, SabrChoiWu2021P, SabrUncorrChoiWu2021
    :inherited-members:
    :exclude-members: SabrABC, SabrVolApproxABC
    :noindex:

docs/source/sv_fft.rst

@@ -0,0 +1,7 @@
+Stochastic volatility with Fourier inversion
+==========================
+.. automodule:: pyfeng.sv_fft
+   :members: HestonFft, OusvFft
+   :inherited-members:
+   :exclude-members: FftABC
+   :noindex:

pyfeng/__init__.py

@@ -7,17 +7,21 @@
 from .sv_fft import HestonFft, BsmFft, OusvFft, VarGammaFft, ExpNigFft
 from .sabr import (
     SabrHagan2002,
-    SabrNorm,
+    SabrNormVolApprox,
     SabrLorig2017,
     SabrChoiWu2021H,
     SabrChoiWu2021P,
 )
 from .garch import GarchMcTimeStep, GarchUncorrBaroneAdesi2004
 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 SabrMcCond
-from .nsvh import Nsvh1, NsvhMc
+from .sabr_mc import SabrMcTimeDisc
+from .nsvh import Nsvh1, NsvhMc, NsvhQuadInt
 from .multiasset import (
     BsmSpreadKirk,
     BsmSpreadBjerksund2014,

pyfeng/bsm.py

@@ -204,9 +204,7 @@ def _impvol_newton(self, price, strike, spot, texp, cp=1, setval=False):
 
         # Exclude optoin price below intrinsic value or above max value (1 for call or k for put)
         # ind_solve can be scalar or array. scalar can be fine in np.abs(p_err[ind_solve])
-        ind_solve = (price_std - p_min > Bsm.IMPVOL_TOL) & (
-                p_max - price_std > Bsm.IMPVOL_TOL
-        )
+        ind_solve = (price_std - p_min > Bsm.IMPVOL_TOL) & (p_max - price_std > Bsm.IMPVOL_TOL)
 
         # initial guess = inflection point in sigma (volga=0)
         _sigma = np.ones_like(ind_solve)*np.sqrt(
@@ -287,21 +285,11 @@ def vol_smile(self, strike, spot, texp, cp=1, model="norm"):
             kk = strike/fwd
             lnk = np.log(kk)
             if model.lower() == "norm-approx":
-                return (
-                        self.sigma
-                        *fwd
-                        *np.sqrt(kk)
-                        *(1 + lnk**2/24)
-                        /(1 + self.sigma**2*texp/24)
-                )
+                return self.sigma * fwd * np.sqrt(kk) * (1 + lnk**2/24) / (1 + self.sigma**2*texp/24)
             else:
                 with np.errstate(divide="ignore", invalid="ignore"):
                     term1 = np.where(np.fabs(lnk) > 1e-8, (kk - 1)/lnk, 2/(3 - kk))
-                    term2 = np.where(
-                        np.fabs(lnk) > 1e-8,
-                        (np.log(term1) - lnk/2)/lnk**2,
-                        1/24,
-                    )
+                    term2 = np.where(np.fabs(lnk) > 1e-8, (np.log(term1) - lnk/2)/lnk**2, 1/24)
                 return self.sigma*fwd*term1*(1 - term2*self.sigma**2*texp)
         else:
             raise ValueError(f"Unknown model: {model}")

pyfeng/cev.py

@@ -45,9 +45,7 @@ def mass_zero(self, spot, texp, log=False):
 
         betac = 1.0 - self.beta
         a = 0.5/betac
-        sigma_std = np.maximum(
-            self.sigma/np.power(fwd, betac)*np.sqrt(texp), np.finfo(float).eps
-        )
+        sigma_std = np.maximum(self.sigma/np.power(fwd, betac)*np.sqrt(texp), np.finfo(float).eps)
         x = 0.5/np.square(betac*sigma_std)
 
         if log:

pyfeng/data/heston_benchmark.xlsx

@@ -1,14 +1,14 @@
 #from .heston_fft import HestonFft
 
 # SV models (CMC, AE) from ASP 2021
-from .heston_mc import HestonMcAndersen2008, HestonMcGlassermanKim2011, HestonMcTseWan2013, HestonMcChoiKwok2023
+from .heston_mixture import HestonMixture
 from .sv32_mc import Sv32McCondQE, Sv32McAe2
 from .sv32_mc2 import Sv32McTimeStep, Sv32McExactBaldeaux2012, Sv32McExactChoiKwok2023
 from .subord_bm import VarGammaQuad, ExpNigQuad
 
 # SABR / OUSV models for research
-from .sabr_int import SabrCondQuad
-from .sabr_mc import SabrMcExactCai2017
+from .sabr_int import SabrMixture
+from .sabr_mc import SabrMcCai2017Exact
 from .ousv import OusvSchobelZhu1998, OusvMcTimeStep, OusvMcChoi2023
 
 # Basket-Asian from ASP 2021

pyfeng/ex.py

@@ -2,6 +2,8 @@
 from . import sv_abc as sv
 from . import bsm
 
+#### Use of RN generation spawn:
+# 0: simulation of variance (gamma/ncx2/normal)
 
 class GarchUncorrBaroneAdesi2004(sv.SvABC):
     """
@@ -15,9 +17,10 @@ class GarchUncorrBaroneAdesi2004(sv.SvABC):
     """
 
     model_type = "GarchDiff"
+    var_process = True
     order = 2
 
-    def avgvar_mv(self, var0, texp):
+    def avgvar_mv(self, texp, var0):
         """
         Mean and variance of the average variance given V(0) = var0.
         Eqs. (12)-(13) in Barone-Adesi et al. (2005)
@@ -63,7 +66,7 @@ def price(self, strike, spot, texp, cp=1):
         if not np.isclose(self.rho, 0.0):
             print(f"Pricing ignores rho = {self.rho}.")
 
-        avgvar, var = self.avgvar_mv(self.sigma, texp)
+        avgvar, var = self.avgvar_mv(texp, self.sigma)
 
         m_bs = bsm.Bsm(np.sqrt(avgvar), intr=self.intr, divr=self.divr)
         price = m_bs.price(strike, spot, texp, cp)
@@ -103,7 +106,7 @@ def set_num_params(self, n_path=10000, dt=0.05, rn_seed=None, antithetic=True, s
         super().set_num_params(n_path, dt, rn_seed, antithetic)
         self.scheme = scheme
 
-    def var_step_euler(self, var_0, dt, milstein=True):
+    def vol_step_euler(self, dt, var_0, milstein=True):
         """
         Euler/Milstein Schemes:
         v_(t+dt) = v_t + mr * (theta - v_t) * dt + vov * v_t Z * sqrt(dt) + (vov^2/2) v_t (Z^2-1) dt
@@ -115,8 +118,7 @@ def var_step_euler(self, var_0, dt, milstein=True):
         Returns: Variance path (time, path) including the value at t=0
         """
 
-        zz = self.rv_normal()
-
+        zz = self.rv_normal(spawn=0)
         var_t = var_0 + self.mr*(self.theta - var_0)*dt + self.vov*var_0*np.sqrt(dt)*zz
         if milstein:
             var_t += (self.vov**2/2)*var_0*dt*(zz**2 - 1)
@@ -138,64 +140,69 @@ def var_step_log(self, log_var_0, dt):
         Returns: Variance path (time, path) including the value at t=0
         """
 
-        zz = self.rv_normal()
+        zz = self.rv_normal(spawn=0)
 
         log_var_t = log_var_0 + (self.mr*self.theta*np.exp(-log_var_0) - self.mr - self.vov**2/2)*dt \
                     + self.vov*np.sqrt(dt)*zz
 
         return log_var_t
 
-    def cond_states(self, var_0, texp):
-        tobs = self.tobs(texp)
-        n_dt = len(tobs)
-        dt = np.diff(tobs, prepend=0)
-
-        # precalculate the Simpson's rule weight
-        weight = np.ones(n_dt + 1)
-        weight[1:-1:2] = 4
-        weight[2:-1:2] = 2
-        weight /= weight.sum()
-
-        var_t = np.full(self.n_path, var_0)
-
-        mean_var = weight[0]*var_t
-        mean_vol = weight[0]*np.sqrt(var_t)
-        mean_inv_vol = weight[0]/np.sqrt(var_t)
+    def cond_states_step(self, dt, var_0):
 
         if self.scheme < 2:
             milstein = (self.scheme == 1)
+            var_t = self.vol_step_euler(dt, var_0, milstein=milstein)
 
-            for i in range(n_dt):
-                var_t = self.var_step_euler(var_t, dt[i], milstein=milstein)
-                mean_var += weight[i + 1]*var_t
-                vol_t = np.sqrt(var_t)
-                mean_vol += weight[i + 1]*vol_t
-                mean_inv_vol += weight[i + 1]/vol_t
+            mean_var = (var_0 + var_t)/2
+            vol_0 = np.sqrt(var_0)
+            vol_t = np.sqrt(var_t)
+            mean_vol = (vol_0 + vol_t)/2
+            mean_inv_vol = (1/vol_0 + 1/vol_t)/2
 
         elif self.scheme == 2:
-            log_var_t = np.full(self.n_path, np.log(var_0))
-            for i in range(n_dt):
-                # Euler scheme on Log(var)
-                log_var_t = self.var_step_log(log_var_t, dt[i])
-                vol_t = np.exp(log_var_t/2)
-                mean_var += weight[i + 1]*vol_t**2
-                mean_vol += weight[i + 1]*vol_t
-                mean_inv_vol += weight[i + 1]/vol_t
+
+            log_var_t = self.var_step_log(np.log(var_0), dt)
+
+            vol_0 = np.sqrt(var_0)
+            vol_t = np.exp(log_var_t/2)
+            var_t = vol_t**2
+
+            mean_var = (var_0 + var_t)/2
+            mean_vol = (vol_0 + vol_t)/2
+            mean_inv_vol = (1/vol_0 + 1/vol_t)/2
         else:
             raise ValueError(f'Invalid scheme: {self.scheme}')
 
-        return vol_t, mean_var, mean_vol, mean_inv_vol
-
-    def cond_spot_sigma(self, var_0, texp):
+        return var_t, mean_var, mean_vol, mean_inv_vol
 
-        vol_final, mean_var, mean_vol, mean_inv_vol = self.cond_states(var_0, texp)
+    def cond_spot_sigma(self, texp, var_0):
+        tobs = self.tobs(texp)
+        dt = np.diff(tobs, prepend=0)
+        n_dt = len(dt)
 
-        spot_cond = 2*(vol_final - np.sqrt(var_0))/self.vov \
-                    - self.mr*self.theta*mean_inv_vol*texp/self.vov \
-                    + (self.mr/self.vov + self.vov/4)*mean_vol*texp \
-                    - self.rho*mean_var*texp/2
+        var_t = np.full(self.n_path, var_0)
+        avgvar = np.zeros(self.n_path)
+        avgvol = np.zeros(self.n_path)
+        avgivol = np.zeros(self.n_path)
+
+        for i in range(n_dt):
+            var_t, avgvar_inc, avgvol_inc, avgivol_inc = self.cond_states_step(dt[i], var_t)
+            avgvar += avgvar_inc * dt[i]
+            avgvol += avgvol_inc * dt[i]
+            avgivol += avgivol_inc * dt[i]
+
+        avgvar /= texp
+        avgvol /= texp
+        avgivol /= texp
+
+        spot_cond = 2 * (np.sqrt(var_t) - np.sqrt(var_0)) / self.vov + \
+            (-self.mr * self.theta * avgivol / self.vov \
+            + (self.mr/self.vov + self.vov/4) * avgvol - self.rho * avgvar / 2) * texp
         np.exp(self.rho*spot_cond, out=spot_cond)
 
-        sigma_cond = np.sqrt((1.0 - self.rho**2)/var_0*mean_var)
+        cond_sigma = np.sqrt((1.0 - self.rho**2)/var_0*avgvar)
+
+        return spot_cond, cond_sigma
 
-        return spot_cond, sigma_cond
+    def return_var_realized(self, texp, cond):
+        return None