Merge pull request #135 from PyFE/heston-mc

Ball & Roma (1994) implemented

Author: jaehyukchoi

MANIFEST.in

@@ -1 +1,3 @@
+include pyfeng/data/ousv_benchmark.xlsx
+include pyfeng/data/heston_benchmark.xlsx
 include pyfeng/data/sabr_benchmark.xlsx

pyfeng/__init__.py

@@ -13,6 +13,7 @@
     SabrChoiWu2021P,
 )
 from .garch import GarchMcTimeStep, GarchUncorrBaroneAdesi2004
+from .heston import HestonUncorrBallRoma1994
 from .ousv import OusvUncorrBallRoma1994
 from .sabr_int import SabrUncorrChoiWu2021
 from .sabr_mc import SabrMcCond

pyfeng/bsm.py

@@ -70,11 +70,38 @@ def vega(self, strike, spot, texp, cp=1):
         sigma_std = np.maximum(self.sigma * np.sqrt(texp), np.finfo(float).eps)
         d1 = np.log(fwd / strike) / sigma_std + 0.5 * sigma_std
 
-        vega = (
-            df * fwd * spst.norm.pdf(d1) * np.sqrt(texp)
-        )  # formula according to wikipedia
+        # formula according to wikipedia
+        vega = df * fwd * spst.norm.pdf(d1) * np.sqrt(texp)
         return vega
 
+    def d2_var(self, strike, spot, texp, cp=1):
+        """
+        2nd derivative w.r.t. variance (=sigma^2)
+        Eq. (9) in Hull & White (1987)
+
+        Args:
+            strike: strike price
+            spot: spot (or forward)
+            texp: time to expiry
+            cp: 1/-1 for call/put option
+
+        Returns: d^2 price / d var^2
+
+        References:
+            - Hull J, White A (1987) The Pricing of Options on Assets with Stochastic Volatilities. The Journal of Finance 42:281–300. https://doi.org/10.1111/j.1540-6261.1987.tb02568.x
+        """
+        fwd, df, _ = self._fwd_factor(spot, texp)
+
+        sigma_std = np.maximum(self.sigma * np.sqrt(texp), np.finfo(float).eps)
+        d1 = np.log(fwd / strike) / sigma_std
+        d2 = d1 - 0.5 * sigma_std
+        d1 += 0.5 * sigma_std
+
+        d_sig2 = df * spot * np.sqrt(texp) * spst.norm.pdf(d1) * (d1*d2-1) / (4*self.sigma**3)
+
+        return d_sig2
+
+
     def delta(self, strike, spot, texp, cp=1):
 
         fwd, df, divf = self._fwd_factor(spot, texp)

pyfeng/data/heston_benchmark.xlsx

@@ -6,69 +6,75 @@
 
 class GarchUncorrBaroneAdesi2004(sv.SvABC):
     """
-    The implementation of 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:
         - Barone-Adesi G, Rasmussen H, Ravanelli C (2005) An option pricing formula for the GARCH diffusion model. Computational Statistics & Data Analysis 49:287–310. https://doi.org/10.1016/j.csda.2004.05.014
 
-    This method is only used to compare with the method GarchCondMC.
+    See Also: OusvUncorrBallRoma1994, HestonUncorrBallRoma1994
     """
+
     model_type = "GarchDiff"
+    order = 2
 
-    def price(self, strike, spot, texp, cp=1):
+    def avgvar_mv(self, var0, texp):
+        """
+        Mean and variance of the average variance given V(0) = var0.
+        Eqs. (12)-(13) in Barone-Adesi et al. (2005)
 
-        if not np.isclose(self.rho, 0.0):
-            print(f"Pricing ignores rho = {self.rho}.")
+        Args:
+            var0: initial variance
+            texp: time step
 
-        var0, mr, vov, theta = self.sigma, self.mr, self.vov, self.theta
+        Returns:
+            mean, variance
+        """
+
+        mr, vov, theta = self.mr, self.vov, self.theta
 
         mr2 = mr * mr
         vov2 = vov * vov
         theta2 = theta*theta
-        decay = np.exp(-mr * texp)
+        mr_t = self.mr * texp
+        e_mr = np.exp(-mr_t)
 
         # Eq (12) of Barone-Adesi et al. (2005)
-        M1 = theta + (var0 - theta) * (1 - decay) / (mr * texp)
+        M1 = theta + (var0 - theta) * (1 - e_mr) / mr_t
 
         term1 = vov2 - mr
         term2 = vov2 - 2*mr
 
         # Eq (13)
-        M2c_1 = - (decay * (var0 - theta))**2
+        M2c_1 = -(e_mr * (var0 - theta))**2
         M2c_2 = 2*np.exp(term2 * texp) * (2*mr * theta * (mr * theta + term2 * var0) + term1 * term2 * var0**2)
-        M2c_3 = -vov2 * (theta2 * (4*mr * (3 - texp * mr) + (2*texp * mr - 5) * vov2) + term2 * var0 * (2*theta + var0))
+        M2c_3 = -vov2 * (theta2 * (4*mr * (3 - mr_t) + (2*mr_t - 5) * vov2) + term2 * var0 * (2*theta + var0))
 
-        M2c_4 = 2*decay * vov2
-        M2c_4 *= 2*theta2 * (texp * mr2 - (1 + texp * mr) * vov2) \
+        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 = M2c_1 / mr2 + M2c_2 / (term1 * term2)**2 + M2c_3 / mr2 / term2**2 + M2c_4 / mr2 / term1**2
         M2c /= texp**2
-        # M3c=None
-        # M4c=None
 
-        # Eq. (11)
-        logk = np.log(spot / strike)
-        sigma_std = np.sqrt(self.sigma * texp)
+        return M1, M2c
 
-        m = logk + (self.intr - self.divr) * texp
-        d1 = logk / sigma_std + sigma_std / 2
+    def price(self, strike, spot, texp, cp=1):
 
-        m_bs = bsm.Bsm(np.sqrt(M1), intr=self.intr, divr=self.divr)
-        c_bs = m_bs.price(strike, spot, texp, cp)
+        if not np.isclose(self.rho, 0.0):
+            print(f"Pricing ignores rho = {self.rho}.")
 
-        c_bs_p1 = np.exp(-self.divr * texp) * strike * np.sqrt(texp / 4 / M1) * spst.norm.pdf(d1)
-        c_bs_p2 = c_bs_p1 * texp / 2 * ((m / M1 / texp)**2 - 1 / M1 / texp - 1 / 4)
+        avgvar, var = self.avgvar_mv(self.sigma, texp)
 
-        # C_bs_p3=C_bs_p1*(m**4/(4*(M1*texp)**4)-m**2*(12+M1*texp)/(8*(M1*texp)**3)+(48+8*M1*texp+(M1*texp)**2)/(64*(M1*texp)**2))*texp**2
-        # C_bs_p4=C_bs_p1*(m**6/(8*(M1*texp)**6)-3*m**4*(20+M1*texp)/(32*(M1*texp)**5)+3*m**2*(240+24*M1*texp+(M1*texp)**2)/(128*(M1*texp)**4)-(960+144*M1*texp+12*(M1*texp)**2+(M1*texp)**3)/(512*(M1*texp)**3))*texp**3
+        m_bs = bsm.Bsm(np.sqrt(avgvar), intr=self.intr, divr=self.divr)
+        price = m_bs.price(strike, spot, texp, cp)
 
-        c_ga_2 = c_bs + (M2c / 2) * c_bs_p2
-        # C_ga_3=C_ga_2+(M3c/6)*C_bs_p3
-        # C_ga_4=C_ga_3+(M4c/24)*C_bs_p4
+        if self.order == 2:
+            price += 0.5 * var * m_bs.d2_var(strike, spot, texp, cp)
+        elif self.order > 2:
+            raise ValueError(f"Not implemented for approx order: {self.order}")
 
-        return c_ga_2
+        return price
 
 
 class GarchMcTimeStep(sv.SvABC, sv.CondMcBsmABC):

pyfeng/garch.py

@@ -0,0 +1,81 @@
+import abc
+import numpy as np
+import scipy.stats as spst
+from . import sv_abc as sv
+from . import bsm
+
+
+class HestonABC(sv.SvABC, abc.ABC):
+    model_type = "Heston"
+
+    def var_mv(self, var0, dt):
+        """
+        Mean and variance of the variance V(t+dt) given V(0) = var_0
+
+        Args:
+            var0: initial variance
+            dt: time step
+
+        Returns:
+            mean, variance
+        """
+
+        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
+        return m, s2
+
+    def avgvar_mv(self, var0, texp):
+        """
+        Mean and variance of the average variance given V(0) = var_0.
+        Appnedix B in Ball & Roma (1994)
+
+        Args:
+            var0: initial variance
+            texp: time step
+
+        Returns:
+            mean, variance
+        """
+
+        mr_t = self.mr * texp
+        e_mr = np.exp(-mr_t)
+        x0 = var0 - self.theta
+        vovn = self.vov * np.sqrt(texp)  # normalized vov
+
+        m = self.theta + x0 * (1 - e_mr)/mr_t
+
+        var = (self.theta - 2*x0*e_mr) + \
+              (var0 - 2.5*self.theta + (2*self.theta + (0.5*self.theta - var0)*e_mr)*e_mr)/mr_t
+        var *= (vovn/mr_t)**2
+
+        return m, var
+
+
+class HestonUncorrBallRoma1994(HestonABC):
+    """
+    Ball & Roma (1994)'s approximation pricing formula for European options under uncorrelated (rho=0) Heston model.
+    Up to 2nd order is implemented.
+
+    See Also: OusvUncorrBallRoma1994, GarchUncorrBaroneAdesi2004
+    """
+
+    order = 2
+
+    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)
+
+        m_bs = bsm.Bsm(np.sqrt(avgvar), intr=self.intr, divr=self.divr)
+        price = m_bs.price(strike, spot, texp, cp)
+
+        if self.order == 2:
+            price += 0.5 * var * m_bs.d2_var(strike, spot, texp, cp)
+        elif self.order > 2:
+            raise ValueError(f"Not implemented for approx order: {self.order}")
+
+        return price

pyfeng/heston.py

@@ -8,11 +8,11 @@
 from scipy.misc import derivative
 import functools
 from . import sv_abc as sv
+from . import heston
 
 
-class HestonMcABC(sv.SvABC, sv.CondMcBsmABC, abc.ABC):
+class HestonMcABC(heston.HestonABC, sv.CondMcBsmABC, abc.ABC):
     var_process = True
-    model_type = "Heston"
     scheme = None
 
     def chi_dim(self):
@@ -41,24 +41,6 @@ def phi_exp(self, texp):
         phi = 4*self.mr / self.vov**2 / (1/exp - exp)
         return phi, exp
 
-    def var_mv(self, var_0, dt):
-        """
-        Mean and variance of the variance V(t+dt) given V(0) = var_0
-
-        Args:
-            var_0: initial variance
-            dt: time step
-
-        Returns:
-            mean, variance
-        """
-
-        expo = np.exp(-self.mr * dt)
-        m = self.theta + (var_0 - self.theta) * expo
-        s2 = var_0 * expo + self.theta * (1 - expo) / 2
-        s2 *= self.vov**2 * (1 - expo) / self.mr
-        return m, s2
-
     def var_step_euler(self, var_0, dt, milstein=False):
         """
         Simulate final variance with Euler/Milstein schemes (scheme = 0, 1)

pyfeng/heston_mc.py

@@ -6,24 +6,25 @@
 
 
 class OusvABC(sv.SvABC, abc.ABC):
+
     model_type = "OUSV"
 
-    def avgvar_mv(self, var_0, texp):
+    def avgvar_mv(self, var0, 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
+            var0: initial variance
             texp: time step
 
         Returns:
-            mean, variance(=NULL)
+            mean, variance(=None)
         """
 
         mr_t = self.mr * texp
         e_mr = np.exp(-mr_t)
-        x0 = var_0 - self.theta
+        x0 = var0 - 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, None
@@ -129,16 +130,28 @@ def price(self, strike, spot, texp, cp=1):
 
 
 class OusvUncorrBallRoma1994(OusvABC):
+    """
+    Ball & Roma (1994)'s approximation pricing formula for European options under uncorrelated (rho=0) OUSV model.
+    Just 0th order (average variance) is implemented because higher order terms depends on numerical derivative of MGF
+
+    See Also: HestonUncorrBallRoma1994, GarchUncorrBaroneAdesi2004
+    """
+
+    order = 0
 
     def price(self, strike, spot, texp, cp=1):
 
         if not np.isclose(self.rho, 0.0):
             print(f"Pricing ignores rho = {self.rho}.")
 
         avgvar, _ = self.avgvar_mv(self.sigma, texp)
+
         m_bs = bsm.Bsm(np.sqrt(avgvar), intr=self.intr, divr=self.divr)
         price = m_bs.price(strike, spot, texp, cp)
 
+        if self.order > 0:
+            raise ValueError(f"Not implemented for approx order: {self.order}")
+
         return price
 
 
@@ -458,12 +471,12 @@ def cond_states(self, vol_0, texp):
             var_mean /= texp
         return vol_t, var_mean, vol_mean
 
-    def cond_states_step(self, vol_0, dt, zn=None):
+    def cond_states_step(self, vol0, dt, zn=None):
         """
         Incremental conditional states
 
         Args:
-            vol_0: initial volatility
+            vol0: initial volatility
             dt: time step
             zn: specified normal rvs to use. (n_sin + 1, n_path)
 
@@ -479,11 +492,11 @@ 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 = vol_0 - self.theta
+        x_0 = vol0 - self.theta
         if zn is None:
-            x_t = self.vol_step(vol_0, dt) - self.theta
+            x_t = self.vol_step(vol0, dt) - self.theta
         else:
-            x_t = self.vol_step(vol_0, dt, zn=zn[0, :]) - self.theta
+            x_t = self.vol_step(vol0, dt, zn=zn[0, :]) - self.theta
         sighat = x_t - x_0 * e_mr
 
         if zn is None: