Merge pull request #133 from PyFE/heston-mc

Heston mc

Author: jaehyukchoi

pyfeng/__init__.py

@@ -12,6 +12,7 @@
     SabrChoiWu2021H,
     SabrChoiWu2021P,
 )
+from .garch import GarchMcTimeStep, GarchUncorrBaroneAdesi2004
 from .sabr_int import SabrUncorrChoiWu2021
 from .sabr_mc import SabrMcCond
 from .nsvh import Nsvh1, NsvhMc

pyfeng/data/heston_benchmark.xlsx

@@ -4,7 +4,6 @@
 from .heston_mc import HestonMcAndersen2008, HestonMcGlassermanKim2011, HestonMcTseWan2013, HestonMcChoiKwok2023
 from .sv32_mc import Sv32McCondQE, Sv32McAe2
 from .sv32_mc2 import Sv32McTimeStep, Sv32McExactBaldeaux2012, Sv32McExactChoiKwok2023
-from .garch import GarchMcTimeStep, GarchUncorrBaroneAdesi2004
 from .subord_bm import VarGammaQuad, ExpNigQuad
 
 # SABR / OUSV models for research

pyfeng/ex.py

@@ -81,7 +81,7 @@ class GarchMcTimeStep(sv.SvABC, sv.CondMcBsmABC):
     var_process = True
     scheme = 1  #
 
-    def set_mc_params(self, n_path=10000, dt=0.05, rn_seed=None, antithetic=True, scheme=1):
+    def set_num_params(self, n_path=10000, dt=0.05, rn_seed=None, antithetic=True, scheme=1):
         """
         Set MC parameters
 
@@ -95,7 +95,7 @@ def set_mc_params(self, n_path=10000, dt=0.05, rn_seed=None, antithetic=True, sc
         References:
             - Andersen L (2008) Simple and efficient simulation of the Heston stochastic volatility model. Journal of Computational Finance 11:1–42. https://doi.org/10.21314/JCF.2008.189
         """
-        super().set_mc_params(n_path, dt, rn_seed, antithetic)
+        super().set_num_params(n_path, dt, rn_seed, antithetic)
         self.scheme = scheme
 
     def var_step_euler(self, var_0, dt, milstein=True):

pyfeng/garch.py

@@ -162,15 +162,15 @@ class HestonMcAndersen2008(HestonMcABC):
         >>> spot = 100
         >>> sigma, vov, mr, rho, texp = 0.04, 1, 0.5, -0.9, 10
         >>> m = pfex.HestonMcAndersen2008(sigma, vov=vov, mr=mr, rho=rho)
-        >>> m.set_mc_params(n_path=1e5, dt=1/8, rn_seed=123456)
+        >>> m.set_num_params(n_path=1e5, dt=1/8, rn_seed=123456)
         >>> m.price(strike, spot, texp)
         >>> # true price: 44.330, 13.085, 0.296
         array([44.31943535, 13.09371251,  0.29580431])
     """
     psi_c = 1.5  # parameter used by the Andersen QE scheme
     scheme = 4
 
-    def set_mc_params(self, n_path=10000, dt=0.05, rn_seed=None, antithetic=True, scheme=4):
+    def set_num_params(self, n_path=10000, dt=0.05, rn_seed=None, antithetic=True, scheme=4):
         """
         Set MC parameters
 
@@ -184,7 +184,7 @@ def set_mc_params(self, n_path=10000, dt=0.05, rn_seed=None, antithetic=True, sc
         References:
             - Andersen L (2008) Simple and efficient simulation of the Heston stochastic volatility model. Journal of Computational Finance 11:1–42. https://doi.org/10.21314/JCF.2008.189
         """
-        super().set_mc_params(n_path, dt, rn_seed, antithetic)
+        super().set_num_params(n_path, dt, rn_seed, antithetic)
         self.scheme = scheme
 
     def var_step_qe(self, var_0, dt):
@@ -305,7 +305,7 @@ class HestonMcGlassermanKim2011(HestonMcABC):
     kk = 1  # K for series truncation.
     tabulate_x2_z = False
 
-    def set_mc_params(self, n_path=10000, dt=None, rn_seed=None, scheme=3, kk=1):
+    def set_num_params(self, n_path=10000, dt=None, rn_seed=None, scheme=3, kk=1):
         """
         Set MC parameters
 
@@ -317,7 +317,7 @@ def set_mc_params(self, n_path=10000, dt=None, rn_seed=None, scheme=3, kk=1):
             kk: truncation index
 
         """
-        super().set_mc_params(n_path, dt, rn_seed, antithetic=False)
+        super().set_num_params(n_path, dt, rn_seed, antithetic=False)
         self.scheme = scheme
         self.kk = kk
 
@@ -828,14 +828,14 @@ class HestonMcTseWan2013(HestonMcGlassermanKim2011):
         >>> spot = 100
         >>> sigma, vov, mr, rho, texp = 0.04, 1, 0.5, -0.9, 10
         >>> m = pfex.HestonMcTseWan2013(sigma, vov=vov, mr=mr, rho=rho)
-        >>> m.set_mc_params(n_path=1e4, rn_seed=123456)
+        >>> m.set_num_params(n_path=1e4, rn_seed=123456)
         >>> m.price(strike, spot, texp)
         >>> # true price: 44.330, 13.085, 0.296
         array([12.08981758,  0.33379748, 42.28798189])  # not close so far
     """
     dist = 'ig'
 
-    def set_mc_params(self, n_path=10000, dt=None, rn_seed=None, scheme=3, dist=None):
+    def set_num_params(self, n_path=10000, dt=None, rn_seed=None, scheme=3, dist=None):
         """
         Set MC parameters
 
@@ -847,7 +847,7 @@ def set_mc_params(self, n_path=10000, dt=None, rn_seed=None, scheme=3, dist=None
             dist: distribution to use for approximation.
                 'ig' for inverse Gaussian (default), 'ga' for Gamma, 'ln' for LN
         """
-        super().set_mc_params(n_path, dt, rn_seed, scheme=scheme)
+        super().set_num_params(n_path, dt, rn_seed, scheme=scheme)
         if dist is not None:
             self.dist = dist
 
@@ -891,7 +891,7 @@ class HestonMcChoiKwok2023(HestonMcGlassermanKim2011):
 
     dist = 'ig'
 
-    def set_mc_params(self, n_path=10000, dt=None, rn_seed=None, scheme=3, kk=0, dist=None):
+    def set_num_params(self, n_path=10000, dt=None, rn_seed=None, scheme=3, kk=0, dist=None):
         """
         Set MC parameters
 
@@ -903,7 +903,7 @@ def set_mc_params(self, n_path=10000, dt=None, rn_seed=None, scheme=3, kk=0, dis
             dist: distribution to use for approximation.
                 'ig' for inverse Gaussian (default), 'ga' for Gamma, 'ln' for LN
         """
-        super().set_mc_params(n_path, dt, rn_seed, scheme=scheme, kk=kk)
+        super().set_num_params(n_path, dt, rn_seed, scheme=scheme, kk=kk)
         if dist is not None:
             self.dist = dist
 

pyfeng/heston_mc.py

@@ -156,6 +156,26 @@ def cond_spot_sigma(self, vol_0, texp):
         sigma_cond = np.sqrt((1 - self.rho**2) * var_mean) / vol_0
         return spot_cond, sigma_cond
 
+    def avgvar_mv(self, var_0, texp):
+        """
+        Mean and variance of the variance V(t+dt) given V(0) = var_0
+        (variance is not implemented yet)
+
+        Args:
+            var_0: initial variance
+            texp: time step
+
+        Returns:
+            mean, variance
+        """
+
+        mr_t = self.mr * texp
+        e_mr = np.exp(-mr_t)
+        x0 = var_0 - self.theta
+        vv = self.vov**2/2/self.mr + self.theta**2 + \
+             ((x0**2 - self.vov**2/2/self.mr)*(1 + e_mr) + 4*self.theta * x0)*(1 - e_mr)/(2*self.mr*texp)
+        return vv
+
 
 class OusvMcTimeStep(OusvMcABC):
     """
@@ -211,7 +231,7 @@ class OusvMcChoi2023(OusvMcABC):
 
     n_sin = 2
 
-    def set_mc_params(self, n_path=10000, dt=None, rn_seed=None, antithetic=True, n_sin=2):
+    def set_num_params(self, n_path=10000, dt=None, rn_seed=None, antithetic=True, n_sin=2):
         """
         Set MC parameters
 
@@ -224,7 +244,7 @@ def set_mc_params(self, n_path=10000, dt=None, rn_seed=None, antithetic=True, n_
         assert n_sin % 2 == 0
         self.n_sin = n_sin
 
-        super().set_mc_params(n_path, dt, rn_seed, antithetic)
+        super().set_num_params(n_path, dt, rn_seed, antithetic)
 
     @classmethod
     def _a2sum(cls, mr_t, ns=0, odd=None):
@@ -444,7 +464,7 @@ def cond_states_step(self, vol_0, dt, zn=None):
         cosh = np.cosh(mr_t)
         vovn = self.vov * np.sqrt(dt)  # normalized vov
 
-        x_0 = self.sigma - self.theta
+        x_0 = vol_0 - self.theta
         if zn is None:
             x_t = self.vol_step(vol_0, dt) - self.theta
         else:
@@ -512,6 +532,11 @@ def cond_states_step(self, vol_0, dt, zn=None):
 
         return x_t, vv_t, uu_t
 
+    def unexplained_var_ratio(self, mr_t, ns):
+        ns = ns or self.n_sin
+        rv = self._a4sum(mr_t, ns=ns) / self._a4sum(mr_t)
+        return rv
+
     def vol_path_sin(self, tobs, zn=None):
         """
         vol path composed of sin terms
@@ -548,3 +573,21 @@ def vol_path_sin(self, tobs, zn=None):
                      + self.vov * np.sqrt(dt) * (an*sin) @ zn[1:,:]
 
         return sigma_path
+
+    def price_var_option(self, strike, texp, cp=1):
+        """
+        Price of variance option
+
+        Args:
+            strike:
+            texp:
+            cp:
+
+        Returns:
+
+        """
+        df = np.exp(-self.intr * texp)
+        vol_t, vv_t, uu_t = self.cond_states(self.sigma, texp)
+        # vv_t is the average variance
+        price = df * np.fmax(np.sign(cp)*(vv_t[:, None] - strike), 0).mean(axis=0)
+        return price

pyfeng/ousv.py

@@ -117,7 +117,7 @@ class SabrMcExactCai2017(sabr.SabrABC, sv.CondMcBsmABC):
     comb_coef = None
     nn = None
 
-    def set_mc_params(self, n_path=10000, m_inv=20, m_euler=20, n_euler=35, rn_seed=None, antithetic=True):
+    def set_num_params(self, n_path=10000, m_inv=20, m_euler=20, n_euler=35, rn_seed=None, antithetic=True):
         """
         Set MC parameters
 

pyfeng/sabr_mc.py

@@ -57,7 +57,7 @@ class Sv32McTimeStep(Sv32McABC):
     """
     scheme = 1  # Milstein
 
-    def set_mc_params(self, n_path=10000, dt=0.05, rn_seed=None, antithetic=True, scheme=1):
+    def set_num_params(self, n_path=10000, dt=0.05, rn_seed=None, antithetic=True, scheme=1):
         """
         Set MC parameters
 
@@ -71,13 +71,13 @@ def set_mc_params(self, n_path=10000, dt=0.05, rn_seed=None, antithetic=True, sc
         References:
             - Andersen L (2008) Simple and efficient simulation of the Heston stochastic volatility model. Journal of Computational Finance 11:1–42. https://doi.org/10.21314/JCF.2008.189
         """
-        super().set_mc_params(n_path, dt, rn_seed, antithetic)
+        super().set_num_params(n_path, dt, rn_seed, antithetic)
 
         self.scheme = scheme
         mr = self.mr * self.theta
         theta = (self.mr + self.vov**2)/mr
         self._m_heston = heston_mc.HestonMcAndersen2008(1/self.sigma, self.vov, self.rho, mr, theta)
-        self._m_heston.set_mc_params(n_path, dt, rn_seed, antithetic, scheme=scheme)
+        self._m_heston.set_num_params(n_path, dt, rn_seed, antithetic, scheme=scheme)
 
     def var_step_euler(self, var_0, dt, milstein=True):
         """
@@ -155,7 +155,7 @@ class Sv32McExactBaldeaux2012(Sv32McABC):
         is_fwd: Bool, true if asset price is forward
     """
 
-    def set_mc_params(self, n_path=10000, dt=None, rn_seed=None, antithetic=True, scheme=1):
+    def set_num_params(self, n_path=10000, dt=None, rn_seed=None, antithetic=True, scheme=1):
         """
         Set MC parameters
     
@@ -169,13 +169,13 @@ def set_mc_params(self, n_path=10000, dt=None, rn_seed=None, antithetic=True, sc
         References:
             - Andersen L (2008) Simple and efficient simulation of the Heston stochastic volatility model. Journal of Computational Finance 11:1–42. https://doi.org/10.21314/JCF.2008.189
         """
-        super().set_mc_params(n_path, dt, rn_seed, antithetic)
+        super().set_num_params(n_path, dt, rn_seed, antithetic)
 
         self.scheme = scheme
         mr = self.mr * self.theta
         theta = (self.mr + self.vov**2)/mr
         self._m_heston = heston_mc.HestonMcAndersen2008(1/self.sigma, self.vov, self.rho, mr, theta)
-        self._m_heston.set_mc_params(n_path, dt, rn_seed, antithetic, scheme=scheme)
+        self._m_heston.set_num_params(n_path, dt, rn_seed, antithetic, scheme=scheme)
 
     def laplace(self, bb, var_0, var_t, dt):
         phi, _ = self._m_heston.phi_exp(dt)

pyfeng/sv32_mc2.py

@@ -116,7 +116,7 @@ class CondMcBsmABC(smile.OptSmileABC, abc.ABC):
     correct_fwd = True
     result = {}
 
-    def set_mc_params(self, n_path=10000, dt=0.05, rn_seed=None, antithetic=True):
+    def set_num_params(self, n_path=10000, dt=0.05, rn_seed=None, antithetic=True):
         """
         Set MC parameters
 
@@ -152,7 +152,7 @@ def tobs(self, texp):
         if self.dt is None:
             return np.array([texp])
         else:
-            n_dt = np.ceil(texp / (2 * self.dt)) * 2
+            n_dt = np.ceil(texp / self.dt)
             tobs = np.arange(1, n_dt + 1) / n_dt * texp
         return tobs
 

pyfeng/sv_abc.py

@@ -55,7 +55,7 @@ def test_price_mc(self):
         """
         for no in [1, 2, 3]:
             m, p, rv = pfex.HestonMcAndersen2008.init_benchmark(no)
-            m.set_mc_params(n_path=1e5, dt=1/8, rn_seed=123456)
+            m.set_num_params(n_path=1e5, dt=1/8, rn_seed=123456)
             m.correct_fwd = False
 
             vol0 = pf.Bsm(None, intr=m.intr, divr=m.divr).impvol(rv['val'], **rv['args_pricing'])
@@ -65,27 +65,27 @@ def test_price_mc(self):
             np.testing.assert_allclose(m.result['spot error'], 0, atol=2e-3)
 
             m, *_ = pfex.HestonMcGlassermanKim2011.init_benchmark(no)
-            m.set_mc_params(n_path=1e5, rn_seed=123456, kk=10)
+            m.set_num_params(n_path=1e5, rn_seed=123456, kk=10)
             m.correct_fwd = False
             vol1 = m.vol_smile(**rv['args_pricing'])
             np.testing.assert_allclose(vol0, vol1, atol=5e-3)
             np.testing.assert_allclose(m.result['spot error'], 0, atol=2e-3)
 
             m, *_ = pfex.HestonMcTseWan2013.init_benchmark(no)
-            m.set_mc_params(n_path=1e5, rn_seed=123456, dt=1)
+            m.set_num_params(n_path=1e5, rn_seed=123456, dt=1)
             m.correct_fwd = False
             vol1 = m.vol_smile(**rv['args_pricing'])
             np.testing.assert_allclose(vol0, vol1, atol=5e-3)
             np.testing.assert_allclose(m.result['spot error'], 0, atol=2e-3)
 
             m, *_ = pfex.HestonMcChoiKwok2023.init_benchmark(no)
             m.correct_fwd = False
-            m.set_mc_params(n_path=1e5, rn_seed=123456, kk=10, dt=None)
+            m.set_num_params(n_path=1e5, rn_seed=123456, kk=10, dt=None)
             vol1 = m.vol_smile(**rv['args_pricing'])
             np.testing.assert_allclose(vol0, vol1, atol=5e-3)
             np.testing.assert_allclose(m.result['spot error'], 0, atol=2e-3)
 
-            m.set_mc_params(n_path=1e5, rn_seed=123456, kk=1, dt=1/4)
+            m.set_num_params(n_path=1e5, rn_seed=123456, kk=1, dt=1/4)
             vol1 = m.vol_smile(**rv['args_pricing'])
             np.testing.assert_allclose(vol0, vol1, atol=5e-3)
             np.testing.assert_allclose(m.result['spot error'], 0, atol=2e-3)