Merge pull request #26 from PyFE/v0.1.6

V0.1.6

Author: jaehyukchoi

MANIFEST.in

@@ -0,0 +1 @@
+include pyfeng/data/sabr_benchmark.xlsx

pyfeng/__init__.py

@@ -1,11 +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 .gamma import InvGam
+from .gamma import InvGam, InvGauss
 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, \
-    BsmBasketMilevsky1998
+    BsmBasketMilevsky1998, BsmBasket1Bm
 from .multiasset_mc import BsmNdMc, NormNdMc
 

pyfeng/bsm.py

@@ -308,7 +308,7 @@ def price_barrier(self, strike, barrier, spot, texp, cp=1, io=-1):
         return p
 
 
-class BsmDisp(smile.OptSmileABC):
+class BsmDisp(smile.OptSmileABC, Bsm):
     """
     Displaced Black-Scholes-Merton model for option pricing. Displace price,
 
@@ -321,17 +321,19 @@ class BsmDisp(smile.OptSmileABC):
         >>> import pyfeng as pf
         >>> m = pf.BsmDisp(sigma=0.2, beta=0.5, pivot=100, intr=0.05, divr=0.1)
         >>> m.price(np.arange(80, 121, 10), 100, 1.2)
-        >>> sigma = np.array([0.2, 0.3, 0.5])[:, None]
-        >>> m = pf.BsmDisp(sigma, beta=0.5, pivot=100, intr=0.05, divr=0.1) # sigma in axis=0
-        >>> m.price(np.array([90, 100, 110]), 100, 1.2, cp=np.array([-1,1,1]))
-        array([[ 5.75927238,  5.52948546,  2.94558338],
-               [ 9.4592961 ,  9.3881245 ,  6.45745004],
-               [16.812035  , 17.10541288, 14.10354768]])
+        array([15.9543935 ,  9.7886658 ,  5.4274197 ,  2.71430505,  1.22740381])
+        >>> beta = np.array([0.2, 0.6, 1])[:, None]
+        >>> m = pf.BsmDisp(0.2, beta=beta, pivot=100) # beta in axis=0
+        >>> m.vol_smile(np.arange(80, 121, 10), 100, 1.2)
+        array([[0.21915778, 0.20904587, 0.20038559, 0.19286293, 0.18625174],
+               [0.20977955, 0.20461468, 0.20025691, 0.19652101, 0.19327567],
+               [0.2       , 0.2       , 0.2       , 0.2       , 0.2       ]])
     """
 
-    pivot = None
     beta = 1  # equivalent to Black-Scholes
-    bsm_model = None
+    pivot = 0
+    sigma_disp = None
+    #_m_bsm = None
 
     def __init__(self, sigma, beta=1, pivot=0, *args, **kwargs):
         """
@@ -343,53 +345,80 @@ def __init__(self, sigma, beta=1, pivot=0, *args, **kwargs):
             divr: dividend/convenience yield (foreign interest rate)
             is_fwd: if True, treat `spot` as forward price. False by default.
         """
-        super().__init__(sigma, *args, **kwargs)
         self.pivot = pivot
         self.beta = beta
-        self.bsm_model = Bsm(sigma=beta * sigma)
+        #self._m_bsm = Bsm(sigma=beta*sigma, *args, **kwargs)
+        super().__init__(sigma, *args, **kwargs)
 
-    '''
     @property
     def sigma(self):
-        return self.sigma
+        return self.sigma_disp * self.beta
 
     @sigma.setter
     def sigma(self, sigma):
-        self.sigma = sigma
-        self.bsm_model.sigma = self.beta*sigma
-    '''
+        self.sigma_disp = sigma
 
-    def disp(self, x):
-        return self.beta*x + (1-self.beta)*self.pivot
+    def disp_spot(self, spot):
+        """
+        Displaced spot
 
-    def price(self, strike, spot, *args, **kwargs):
-        return (1/self.beta)*self.bsm_model.price(
-            self.disp(strike), self.disp(spot), *args, **kwargs)
+        Args:
+            spot: spot (or forward) price
 
-    def delta(self, strike, spot, *args, **kwargs):
-        return self.bsm_model.delta(
-            self.disp(strike), self.disp(spot), *args, **kwargs)
+        Returns:
+            Displaces spot
+        """
+        return self.beta*spot + (1-self.beta)*self.pivot
 
-    def cdf(self, strike, spot, *args, **kwargs):
-        return self.bsm_model.cdf(
-            self.disp(strike), self.disp(spot), *args, **kwargs)
+    def disp_strike(self, strike, texp):
+        """
+        Displaced strike
 
-    def vega(self, strike, spot, *args, **kwargs):
-        return self.bsm_model.vega(
-            self.disp(strike), self.disp(spot), *args, **kwargs)
+        Args:
+            strike: strike price
+            texp: time to expiry
 
-    def gamma(self, strike, spot, *args, **kwargs):
+        Returns:
+            Displace strike
+        """
+        return self.beta*strike + (1-self.beta)*self.forward(self.pivot, texp)
+
+    def price(self, strike, spot, texp, cp=1):
+        spot = self.disp_spot(spot)
+        strike = self.disp_strike(strike, texp)
+        return (1/self.beta) * super().price(strike, spot, texp, cp=cp)
+
+    def delta(self, strike, spot, texp, cp=1):
+        spot = self.disp_spot(spot)
+        strike = self.disp_strike(strike, texp)
+        return super().delta(strike, spot, texp, cp=cp)
+
+    def cdf(self, strike, spot, texp, cp=1):
+        spot = self.disp_spot(spot)
+        strike = self.disp_strike(strike, texp)
+        return super().cdf(strike, spot, texp, cp=cp)
+
+    def vega(self, strike, spot, texp, cp=1):
+        spot = self.disp_spot(spot)
+        strike = self.disp_strike(strike, texp)
+        return super().vega(strike, spot, texp, cp=cp)
+
+    def gamma(self, strike, spot, texp, cp=1):
         # need to mutiply beta because of (beta*sigma) appearing in the denominator of the bsm gamma
-        return self.beta * self.bsm_model.gamma(
-            self.disp(strike), self.disp(spot), *args, **kwargs)
+        spot = self.disp_spot(spot)
+        strike = self.disp_strike(strike, texp)
+        return self.beta * super().gamma(strike, spot, texp, cp=cp)
 
-    def theta(self, strike, spot, *args, **kwargs):
-        return (1/self.beta)*self.bsm_model.theta(
-            self.disp(strike), self.disp(spot), *args, **kwargs)
+    def theta(self, strike, spot, texp, cp=1):
+        spot = self.disp_spot(spot)
+        strike = self.disp_strike(strike, texp)
+        return (1/self.beta)*super().theta(
+            strike, spot, texp, cp=cp)
 
     def impvol(self, price_in, strike, spot, texp, cp=1, setval=False):
-        sigma = (1/self.beta)*self.bsm_model.impvol(
-            self.beta*price_in, self.disp(strike), self.disp(spot), texp, cp=cp, setval=setval)
+        spot = self.disp_spot(spot)
+        strike = self.disp_strike(strike, texp)
+        sigma = (1/self.beta)*super().impvol(self.beta*price_in, strike, spot, texp, cp=cp, setval=False)
         if setval:
             self.sigma = sigma
         return sigma
@@ -409,28 +438,29 @@ def vol_smile(self, strike, spot, texp, cp=1, model='bsm'):
             volatility smile under the specified model
         """
         if model.lower() == 'norm-approx':
-            fwd, _, _ = self._fwd_factor(spot, texp)
-            kkd = self.disp(strike) / self.disp(fwd)
+            fwdd = self.forward(self.disp_spot(spot), texp)
+            kkd = self.disp_strike(strike, texp) / fwdd
             lnkd = np.log(kkd)
-            sig_beta = self.beta * self.sigma
-            vol = self.sigma * self.disp(fwd) * np.sqrt(kkd) * (1 + lnkd ** 2 / 24) / \
-                (1 + sig_beta ** 2 * texp / 24)
+            # 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)
         elif model.lower() == 'bsm-approx':
-            fwd, _, _ = self._fwd_factor(spot, texp)
+            fwd = self.forward(spot, texp)
             kk = strike / fwd
             lnk = np.log(kk)
-            fwdd = self.disp(fwd)
-            kkd = self.disp(strike) / fwdd
+
+            fwdd = self.forward(self.disp_spot(spot), texp)
+            kkd = self.disp_strike(strike, texp) / fwdd
             lnkd = np.log(kkd)
-            sig_beta = self.beta * self.sigma
-            vol = self.sigma * (fwdd / fwd) * np.sqrt(kkd / kk)
-            vol *= (1 + lnkd ** 2 / 24) / (1 + lnk ** 2 / 24) * (1 + vol ** 2 * texp / 24) / \
-                (1 + sig_beta ** 2 * texp / 24)
+
+            # 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)
         else:
             vol = super().vol_smile(strike, spot, texp, model=model, cp=cp)
 
         return vol
 
-    def price_barrier(self, strike, barrier, spot, *args, **kwargs):
-        return (1/self.beta)*self.bsm_model.price_barrier(
-            self.disp(strike), self.disp(barrier), self.disp(spot), *args, **kwargs)
+    def price_barrier(self, strike, barrier, spot, texp, cp=1, io=-1):
+        fwd = self.forward(spot, texp)
+        return (1/self.beta)*super().price_barrier(
+            self.disp_spot(strike), self.disp_strike(barrier, strike), self.disp_spot(fwd), texp, cp=cp, io=io)

pyfeng/gamma.py

@@ -64,3 +64,42 @@ def cdf(self, strike, spot, texp, cp=1):
         x = strike/beta
         cdf = np.where(cp > 0, spst.gamma.cdf(1/x, a=alpha), spst.gamma.sf(1/x, a=alpha))
         return cdf
+
+
+class InvGauss(smile.OptSmileABC):
+    """
+    Option pricing model with the inverse Gaussian (IG) distribution.
+
+    The IG distribution with (gamma, delta) is modeled by scipy.stats.invgauss.invgauss(mu=1/(gamma*delta), scale=delta**2)
+    When, sig = gamma = delta, the IG variable has mean 1 and variance 1/sig^2. We match the momoent by
+
+        m2/m1^2 = 1/sig^2 = exp(sigma^2 T)
+
+    Examples:
+        >>> import numpy as np
+        >>> import pyfeng as pf
+        >>> m = pf.InvGauss(sigma=0.2, intr=0.05, divr=0.1)
+        >>> m.price(np.arange(80, 121, 10), 100, 1.2)
+        array([15.71924064,  9.70753358,  5.54459412,  2.95300168,  1.48019682])
+    """
+    sigma = None
+
+    def price(self, strike, spot, texp, cp=1):
+        fwd, df, _ = self._fwd_factor(spot, texp)
+        sig2_inv = np.exp(self.sigma**2 * texp) - 1
+        ig = spst.invgauss(mu=sig2_inv, scale=1/sig2_inv)
+        kk = strike / fwd
+        price = np.where(
+            cp > 0,
+            ig.cdf(1/kk) - kk * ig.sf(kk),
+            kk * ig.cdf(kk) - ig.sf(1/kk),
+        )
+        return df*fwd*price
+
+    def cdf(self, strike, spot, texp, cp=1):
+        fwd, df, _ = self._fwd_factor(spot, texp)
+        sig2_inv = np.exp(self.sigma**2 * texp) - 1
+        ig = spst.invgauss(mu=sig2_inv, scale=1/sig2_inv)
+        x = strike / fwd
+        cdf = np.where(cp > 0, ig.sf(x), ig.cdf(x))
+        return cdf

pyfeng/multiasset.py

@@ -1,4 +1,5 @@
-import scipy.stats as scst
+import scipy.stats as spst
+import warnings
 import numpy as np
 from . import bsm
 from . import norm
@@ -78,8 +79,8 @@ def price(self, strike, spot, texp, cp=1):
         d2 = (d3 - 0.5*std11 + self.rho*std12 + bb*(0.5*bb - 1)*std22) / std
         d3 = (d3 - 0.5*std11 + 0.5*bb**2*std22) / std
 
-        price = cp*(fwd1*scst.norm.cdf(cp*d1) - fwd2*scst.norm.cdf(cp*d2)
-                    - strike*scst.norm.cdf(cp*d3))
+        price = cp*(fwd1*spst.norm.cdf(cp*d1) - fwd2*spst.norm.cdf(cp*d2)
+                    - strike*spst.norm.cdf(cp*d3))
 
         return df * price
 
@@ -271,9 +272,9 @@ def price(self, strike, spot, texp, cp=1):
         price = np.zeros_like(strike, float)
         for k in range(n_strike):
             xx_ = xx + np.log(strike[k])/sig_std
-            term1 = fwd[0] * (scst.norm.cdf(yy[0]) - scst.multivariate_normal.cdf(np.array([xx_[0], yy[0]]), mu0, cor_m1))
-            term2 = fwd[1] * (scst.norm.cdf(yy[1]) - scst.multivariate_normal.cdf(np.array([xx_[1], yy[1]]), mu0, cor_m2))
-            term3 = strike[k] * np.array(scst.multivariate_normal.cdf(xx_ + sig_std, mu0, self.cor_m))
+            term1 = fwd[0] * (spst.norm.cdf(yy[0]) - spst.multivariate_normal.cdf(np.array([xx_[0], yy[0]]), mu0, cor_m1))
+            term2 = fwd[1] * (spst.norm.cdf(yy[1]) - spst.multivariate_normal.cdf(np.array([xx_[1], yy[1]]), mu0, cor_m2))
+            term3 = strike[k] * np.array(spst.multivariate_normal.cdf(xx_ + sig_std, mu0, self.cor_m))
 
             assert(term1 + term2 + term3 >= strike[k])
             price[k] = (term1 + term2 + term3 - strike[k])
@@ -283,3 +284,108 @@ def price(self, strike, spot, texp, cp=1):
         price *= df
 
         return price[0] if strike_isscalar else price
+
+
+class BsmBasket1Bm(opt.OptABC):
+    """
+    Multiasset BSM model for pricing basket/Spread options when all asset prices are driven by a single Brownian motion (BM).
+
+    """
+
+    def __init__(self, sigma, weight=None, intr=0.0, divr=0.0, is_fwd=False):
+        """
+        Args:
+            sigma: model volatilities of `n_asset` assets. (n_asset, )
+            weight: asset weights, If None, equally weighted as 1/n_asset
+                If scalar, equal weights of the value
+                If 1-D array, uses as it is. (n_asset, )
+            intr: interest rate (domestic interest rate)
+            divr: vector of dividend/convenience yield (foreign interest rate) 0-D or (n_asset, ) array
+            is_fwd: if True, treat `spot` as forward price. False by default.
+        """
+
+        sigma = np.atleast_1d(sigma)
+        self.n_asset = len(sigma)
+        if weight is None:
+            self.weight = np.ones(self.n_asset) / self.n_asset
+        elif np.isscalar(weight):
+            self.weight = np.ones(self.n_asset) * weight
+        else:
+            assert len(weight) == self.n_asset
+            self.weight = np.array(weight)
+        super().__init__(sigma, intr=intr, divr=divr, is_fwd=is_fwd)
+
+    @staticmethod
+    def root(fac, std, strike):
+        """
+        Calculate the root x of f(x) = sum(fac * exp(std*x)) - strike = 0 using Newton's method
+
+        Each fac and std should have the same signs so that f(x) is a monotonically increasing function.
+
+        fac: factor to the exponents. (n_asset, ) or (n_strike, n_asset). Asset takes the last dimension.
+        std: total standard variance. (n_asset, )
+        strike: strike prices. scalar or (n_asset, )
+        """
+
+        assert np.all(fac * std >= 0.0)
+        log = np.min(fac) > 0  # Basket if log=True, spread if otherwise.
+        scalar_output = np.isscalar(np.sum(fac * std, axis=-1) - strike)
+        strike = np.atleast_1d(strike)
+
+        with np.errstate(divide='ignore', invalid='ignore'):
+            log_k = np.where(strike > 0, np.log(strike), 1)
+
+            # Initial guess with linearlized assmption
+            x = (strike - np.sum(fac, axis=-1)) / np.sum(fac * std, axis=-1)
+            if log:
+                np.fmin(x, np.amin(np.log(strike[:, None] / fac) / std, axis=-1), out=x)
+            else:
+                np.clip(x, -3, 3, out=x)
+
+        # Test x=-9 and 9 for min/max values.
+        y_max = np.exp(9 * std)
+        y_min = np.sum(fac / y_max, axis=-1) - strike
+        y_max = np.sum(fac * y_max, axis=-1) - strike
+
+        x[y_min >= 0] = -np.inf
+        x[y_max <= 0] = np.inf
+        ind = ~((y_min >= 0) | (y_max <= 0))
+
+        if np.all(~ind):
+            return x[0] if scalar_output else x
+
+        for k in range(32):
+            y_vec = fac * np.exp(std * x[ind, None])
+            y = np.log(np.sum(y_vec, axis=-1)) - log_k[ind] if log else np.sum(y_vec, axis=-1) - strike[ind]
+            dy = np.sum(std * y_vec, axis=-1) / np.sum(y_vec, axis=-1) if log else np.sum(std * y_vec, axis=-1)
+            x[ind] -= y / dy
+            if len(y) == 0:
+                print(ind, y_vec, y)
+            y_err_max = np.amax(np.abs(y))
+            if y_err_max < BsmBasket1Bm.IMPVOL_TOL:
+                break
+
+        if y_err_max > BsmBasket1Bm.IMPVOL_TOL:
+            warn_msg = f'root did not converge within {k} iterations: max error = {y_err_max}'
+            warnings.warn(warn_msg, Warning)
+
+        return x[0] if scalar_output else x
+
+    def price(self, strike, spot, texp, cp=1):
+        fwd, df, _ = self._fwd_factor(spot, texp)
+        assert fwd.shape[-1] == self.n_asset
+
+        fwd_basket = fwd * self.weight
+        sigma_std = self.sigma * np.sqrt(texp)
+        cp = np.array(cp)
+        d2 = -cp * self.root(fwd_basket * np.exp(-sigma_std**2/2), sigma_std, strike)
+
+        if np.isscalar(d2):
+            d1 = d2 + cp*sigma_std
+        else:
+            d1 = d2[:, None] + np.atleast_1d(cp)[:, None]*sigma_std
+
+        price = np.sum(fwd_basket*spst.norm.cdf(d1), axis=-1)
+        price -= strike * spst.norm.cdf(d2)
+        price *= cp*df
+        return price