Merge pull request PyFE#2 from PyFE/main

Upstream

pyfeng/__init__.py

@@ -1,7 +1,11 @@
 from .norm import Norm  # the order is sensitive because of `price_barrier` method. Put it before .bsm
 from .bsm import Bsm, BsmDisp
 from .cev import Cev
-from .sabr import SabrHagan2002, SabrLorig2017, SabrChoiWu2021H, SabrChoiWu2021P
+from .gamma import Invgam
+from .sabr import SabrHagan2002, SabrNorm, SabrLorig2017, SabrChoiWu2021H, SabrChoiWu2021P
 from .sabr_int import SabrUncorrChoiWu2021
 from .nsvh import Nsvh1
-from .multiasset import BsmSpreadKirk, BsmSpreadBjerksund2014, NormBasket, NormSpread, BsmBasketLevy1992, BsmMax2
+from .multiasset import BsmSpreadKirk, BsmSpreadBjerksund2014, NormBasket, NormSpread, BsmBasketLevy1992, BsmMax2, \
+    BsmBasketMilevsky1998
+from .multiasset_mc import BsmNdMc, NormNdMc
+

pyfeng/bsm.py

@@ -81,6 +81,15 @@ def delta(self, strike, spot, texp, cp=1):
         delta *= df if self.is_fwd else divf
         return delta
 
+    def cdf(self, strike, spot, texp, cp=1):
+
+        fwd, df, divf = self._fwd_factor(spot, texp)
+
+        sigma_std = np.maximum(self.sigma * np.sqrt(texp), np.finfo(float).eps)
+        d2 = np.log(fwd / strike) / sigma_std - 0.5 * sigma_std
+        cdf = spst.norm.cdf(cp * d2)  # formula according to wikipedia
+        return cdf
+
     def gamma(self, strike, spot, texp, cp=1):
 
         fwd, df, divf = self._fwd_factor(spot, texp)
@@ -357,6 +366,10 @@ def delta(self, strike, spot, *args, **kwargs):
         return self.bsm_model.delta(
             self.disp(strike), self.disp(spot), *args, **kwargs)
 
+    def cdf(self, strike, spot, *args, **kwargs):
+        return self.bsm_model.cdf(
+            self.disp(strike), self.disp(spot), *args, **kwargs)
+
     def vega(self, strike, spot, *args, **kwargs):
         return self.bsm_model.vega(
             self.disp(strike), self.disp(spot), *args, **kwargs)

pyfeng/cev.py

@@ -74,7 +74,7 @@ def mass_zero_t0(self, spot, texp):
         return t0
 
     @staticmethod
-    def price_formula(strike, spot, texp, cp=1, sigma=None, beta=0.5, intr=0.0, divr=0.0, is_fwd=False):
+    def price_formula(strike, spot, texp, sigma=None, cp=1, beta=0.5, intr=0.0, divr=0.0, is_fwd=False):
         """
 
         Args:
@@ -118,7 +118,7 @@ def price_formula(strike, spot, texp, cp=1, sigma=None, beta=0.5, intr=0.0, divr
         )
         return disc_fac * price
 
-    def delta(self, strike, spot, texp, cp_sign=1):
+    def delta(self, strike, spot, texp, cp=1):
         fwd, df, divf = self._fwd_factor(spot, texp)
         betac_inv = 1 / (1 - self.beta)
 
@@ -127,16 +127,30 @@ def delta(self, strike, spot, texp, cp_sign=1):
         y = k_star * np.power(strike, 2 / betac_inv)
 
         if self.beta < 1.0:
-            delta = 0.5 * (cp_sign - 1) + spst.ncx2.sf(y, 2 + betac_inv, x) + 2 * x / betac_inv * \
+            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))
         else:
-            delta = 0.5 * (cp_sign - 1) + spst.ncx2.sf(x, -betac_inv, y) - 2 * x / betac_inv * \
+            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 *= df if self.is_fwd else divf
         return delta
 
-    def gamma(self, strike, spot, texp, cp_sign=1):
+    def cdf(self, strike, spot, texp, cp=1):
+        fwd = self.forward(spot, texp)
+
+        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
+
+        cdf = np.where(cp > 0, spst.ncx2.cdf(x, betac_inv, y), spst.ncx2.sf(x, betac_inv, y))
+        return cdf
+
+    def gamma(self, strike, spot, texp, cp=1):
         fwd, df, divf = self._fwd_factor(spot, texp)
         betac_inv = 1 / (1 - self.beta)
 
@@ -159,7 +173,7 @@ def gamma(self, strike, spot, texp, cp_sign=1):
 
         return gamma
 
-    def vega(self, strike, spot, texp, cp_sign=1):
+    def vega(self, strike, spot, texp, cp=1):
         fwd, df, divf = self._fwd_factor(spot, texp)
         spot = fwd * df / divf
 

pyfeng/gamma.py

@@ -0,0 +1,66 @@
+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
+
+
+class Invgam(smile.OptSmileABC, opt.OptABC):
+    """
+    Option pricing model with the inverse gamma (reciprocal gamma) distribution.
+
+    The parameters (alpha, beta) is from Wikipedia. https://en.wikipedia.org/wiki/Inverse-gamma_distribution
+    Note that the n-th moment of the inverse gamma RV is beta^n / (alpha-1)*...*(alpha-n).
+    Alpha and beta is calibrated to match the first two moments of the lognormal distribution with volatility sigma
+    so that the option price is similar to that of the BSM model with volatility sigma.
+
+    Examples:
+        >>> import numpy as np
+        >>> import pyfeng as pf
+        >>> m = pf.Invgam(sigma=0.2, intr=0.05, divr=0.1)
+        >>> m.price(np.arange(80, 121, 10), 100, 1.2)
+        array([21.34327542, 13.99490086,  8.60288219,  5.02287171,  2.82349989])
+    """
+    sigma = None
+
+    @staticmethod
+    def price_formula(strike, spot, texp, alpha, beta, cp=1, intr=0.0, divr=0.0, is_fwd=False):
+        disc_fac = np.exp(-texp * intr)
+        fwd_scale = spot * (1.0 if is_fwd else np.exp(-texp * divr) / disc_fac) / beta
+        kk = strike/beta
+
+        price = np.where(
+            cp > 0,
+            fwd_scale*spst.gamma.cdf(x=1/kk, a=alpha-1) - kk*spst.gamma.cdf(x=1/kk, a=alpha),
+            kk*spst.gamma.sf(x=1/kk, a=alpha) - fwd_scale*spst.gamma.sf(x=1/kk, a=alpha-1)
+        )
+        return disc_fac*beta * price
+
+    def alpha_beta(self, spot, texp):
+        """
+        Computes the inverse gamma distribution parameters (alpha, beta) from sigma, spot, texp.
+
+        m1 = beta/(alpha-1)
+        m2/m1^2 = exp(sigma^2 T) = (alpha-1)/(alpha-2)
+
+        Args:
+            spot: spot (or forward) price
+            texp: time to expiry
+
+        Returns: (alpha, beta)
+        """
+
+        fwd = self.forward(spot, texp)
+        alpha = 1/(np.exp(self.sigma**2*texp)-1) + 2
+        beta = (alpha-1)*fwd
+        return alpha, beta
+
+    def price(self, strike, spot, texp, cp=1):
+        alpha, beta = self.alpha_beta(spot, texp)
+        return self.price_formula(strike, spot, texp, alpha, beta, cp=cp, intr=self.intr, divr=self.divr)
+
+    def cdf(self, strike, spot, texp, cp=1):
+        alpha, beta = self.alpha_beta(spot, texp)
+        x = strike/beta
+        cdf = np.where(cp > 0, spst.gamma.cdf(1/x, a=alpha), spst.gamma.sf(1/x, a=alpha))
+        return cdf

pyfeng/multiasset.py

@@ -2,6 +2,7 @@
 import numpy as np
 from . import bsm
 from . import norm
+from . import gamma
 from . import opt_abc as opt
 
 
@@ -178,6 +179,43 @@ def price(self, strike, spot, texp, cp=1):
         return df * price
 
 
+class BsmBasketMilevsky1998(NormBasket):
+    """
+    Basket option pricing with the inverse gamma distribution of Milevsky & Posner (1998)
+
+    References:
+        Milevsky, M. A., & Posner, S. E. (1998). A Closed-Form Approximation for Valuing Basket Options.
+        The Journal of Derivatives, 5(4), 54–61. https://doi.org/10.3905/jod.1998.408005
+
+        Krekel, M., de Kock, J., Korn, R., & Man, T.-K. (2004). An analysis of pricing methods for basket options.
+        Wilmott Magazine, 2004(7), 82–89.
+
+    Examples:
+        >>> import numpy as np
+        >>> import pyfeng as pf
+        >>> strike = np.arange(50, 151, 10)
+        >>> m = pf.BsmBasketMilevsky1998(sigma=0.4*np.ones(4), cor=0.5)
+        >>> m.price(strike, spot=100*np.ones(4), texp=5)
+        array([51.93069524, 44.40986   , 38.02596564, 32.67653542, 28.21560931,
+               24.49577509, 21.38543199, 18.77356434, 16.56909804, 14.69831445,
+               13.10186928])
+    """
+    def price(self, strike, spot, texp, cp=1):
+        df = np.exp(-texp * self.intr)
+        fwd = np.array(spot) * (1.0 if self.is_fwd else np.exp(-texp * np.array(self.divr)) / df)
+        assert fwd.shape[-1] == self.n_asset
+
+        fwd_basket = fwd * self.weight
+        m1 = np.sum(fwd_basket, axis=-1)
+        m2 = np.sum(fwd_basket @ np.exp(self.cov_m * texp) * fwd_basket, axis=-1)
+
+        alpha = 1/(m2/m1**2-1) + 2
+        beta = (alpha-1)*m1
+
+        price = gamma.Invgam.price_formula(strike, m1, texp, alpha, beta, cp=cp, is_fwd=True)
+        return df * price
+
+
 class BsmMax2(opt.OptMaABC):
     """
     Option on the max of two assets.

pyfeng/multiasset_mc.py

@@ -0,0 +1,172 @@
+import numpy as np
+from . import opt_abc as opt
+
+
+class BsmNdMc(opt.OptMaABC):
+    """
+    Monte-Carlo simulation of multiasset (N-d) BSM (geometric Brownian Motion)
+
+    Examples:
+        >>> import pyfeng as pf
+        >>> spot = np.ones(4)*100
+        >>> sigma = np.ones(4)*0.4
+        >>> texp = 5
+        >>> payoff = lambda x: np.fmax(np.mean(x,axis=1) - strike, 0) # Basket option
+        >>> strikes = np.arange(80, 121, 10)
+        >>> m = pf.BsmNdMc(sigma, cor=0.5, rn_seed=1234)
+        >>> m.simulate(tobs=[texp], n_path=20000)
+        >>> p = []
+        >>> for strike in strikes:
+        >>>    p.append(m.price_european(spot, texp, payoff))
+        >>> np.array(p)
+        array([36.31612946, 31.80861014, 27.91269315, 24.55319506, 21.62677625])
+    """
+
+    spot = np.ones(2)
+    sigma = np.ones(2)*0.1
+
+    # MC params
+    n_path = 100
+    rn_seed = None
+    rng = None
+    antithetic = True
+
+    # path
+    path, tobs = None, None
+
+    def __init__(self, sigma, cor=None, intr=0.0, divr=0.0, rn_seed=None):
+        self.rn_seed = rn_seed
+        self.rng = np.random.default_rng(rn_seed)
+        super().__init__(sigma, cor=cor, intr=intr, divr=divr, is_fwd=False)
+
+    def set_mc_params(self, n_path, rn_seed=None, antithetic=True):
+        self.n_path = n_path
+        self.rn_seed = rn_seed
+        self.antithetic = antithetic
+
+    def _bm_incr(self, tobs, n_path=None):
+        """
+        Calculate incremental Brownian Motions
+
+        Args:
+            tobs: array of observation times
+            n_path: number of paths. If None (default), use the stored one.
+            store: if True (default), save the result to self.path_stored
+
+        Returns:
+            price path (time, path, asset)
+        """
+        dt = np.diff(np.atleast_1d(tobs), prepend=0)
+        n_t = len(dt)
+
+        n_path = self.n_path if n_path is None else n_path
+
+        if self.antithetic:
+            # generate random number in the order of path, time, asset and transposed
+            # in this way, the same paths are generated when increasing n_path
+            bm_incr = self.rng.normal(size=(n_path//2, n_t, self.n_asset)).transpose((1, 0, 2))
+            bm_incr *= np.sqrt(dt[:, None, None])
+            bm_incr = np.dot(bm_incr, self.chol_m.T)
+            bm_incr = np.stack([bm_incr, -bm_incr], axis=2).reshape((n_t, n_path, self.n_asset))
+        else:
+            bm_incr = np.random.randn(n_path, n_t, self.n_asset).transpose((1, 0, 2)) * np.sqrt(dt[:, None, None])
+            bm_incr = np.dot(bm_incr, self.chol_m.T)
+
+        return bm_incr
+
+    def simulate(self, tobs, n_path=None, store=True):
+        """
+        Simulate the price paths and store in the class.
+        The initial prices are normalized to 1.
+
+        Args:
+            tobs: array of observation times
+            n_path: number of paths. If None (default), use the stored one.
+            store: if True (default), save the result to self.path_stored
+
+        Returns:
+            price path (time, path, asset)
+        """
+        # (n_t, n_path, n_asset) * (n_asset, n_asset)
+        path = self._bm_incr(tobs, n_path)
+        # Add drift and convexity
+        dt = np.diff(np.atleast_1d(tobs), prepend=0)
+        path += (self.intr - self.divr - 0.5*self.sigma**2)*dt[:, None, None]
+        np.cumsum(path, axis=0, out=path)
+        np.exp(path, out=path)
+
+        if store:
+            self.n_path = n_path
+            self.path = path
+            self.tobs = tobs
+
+        return path
+
+    def price_european(self, spot, texp, payoff):
+        """
+        The European price of that payoff at the expiry.
+
+        Args:
+            spot: array of spot prices
+            texp: time-to-expiry
+            payoff: payoff function applicable to the time-slice of price path
+
+        Returns:
+            The MC price of the payoff
+        """
+        if self.path is None:
+            raise ValueError('Simulated paths are not available. Run simulate() first.')
+
+        # check if texp is in tobs
+        ind, *_ = np.where(np.isclose(self.tobs, texp))
+        if len(ind) == 0:
+            raise ValueError(f'Stored path does not contain t = {texp}')
+
+        path = self.path[ind[0], ] * spot
+        price = np.exp(-self.intr * texp) * np.mean(payoff(path), axis=0)
+        return price
+
+
+class NormNdMc(BsmNdMc):
+    """
+    Monte-Carlo simulation of multiasset (N-d) Normal/Bachelier model (arithmetic Brownian Motion)
+
+    Examples:
+        >>> import pyfeng as pf
+        >>> spot = np.ones(4)*100
+        >>> sigma = np.ones(4)*0.4
+        >>> texp = 5
+        >>> payoff = lambda x: np.fmax(np.mean(x,axis=1) - strike, 0) # Basket option
+        >>> strikes = np.arange(80, 121, 10)
+        >>> m = pf.NormNdMc(sigma*spot, cor=0.5, rn_seed=1234)
+        >>> m.simulate(tobs=[texp], n_path=20000)
+        >>> p = []
+        >>> for strike in strikes:
+        >>>    p.append(m.price_european(spot, texp, payoff))
+        >>> np.array(p)
+        array([39.42304794, 33.60383167, 28.32667559, 23.60383167, 19.42304794])
+    """
+
+    def simulate(self, tobs, n_path=None, store=True):
+        path = self._bm_incr(tobs, n_path)
+        np.cumsum(path, axis=0, out=path)
+
+        if store:
+            self.n_path = n_path
+            self.path = path
+            self.tobs = tobs
+
+        return path
+
+    def price_european(self, spot, texp, payoff):
+        if self.path is None:
+            raise ValueError('Simulated paths are not available. Run simulate() first.')
+
+        # check if texp is in tobs
+        ind, *_ = np.where(np.isclose(self.tobs, texp))
+        if len(ind) == 0:
+            raise ValueError(f'Stored path does not contain t = {texp}')
+
+        path = self.path[ind[0], ] + spot
+        price = np.exp(-self.intr * texp) * np.mean(payoff(path), axis=0)
+        return price