1 files changed, 120 insertions, 0 deletions
diff --git a/python/analytics/option.py b/python/analytics/option.py
new file mode 100644
index 00000000..32f4f947
--- /dev/null
+++ b/python/analytics/option.py
@@ -0,0 +1,120 @@
+from .black import black
+from .utils import GHquad
+from yieldcurve import roll_yc
+from pandas.tseries.offsets import BDay
+
+class Option:
+    def __init__(self, index, exercise_date, strike, option_type="payer"):
+        self.index = index
+        self._exercise_date = exercise_date
+        self._forward_yc = roll_yc(self.index._yc, self.exercise_date)
+        self.exercise_date_settle = (pd.Timestamp(self.exercise_date) + 3* BDay()).date()
+        self._T = None
+        self.strike = strike
+        self.option_type = option_type.lower()
+        self._Z, self._w = GHquad(50)
+        self.notional = 1
+
+    @property
+    def exercise_date(self):
+        return self._exercise_date
+
+    @exercise_date.setter
+    def exercise_date(self, d : datetime.date):
+        self._exercise_date = d
+        self.exercise_date_settle = (pd.Timestamp(d) + 3* BDay()).date()
+        self._forward_yc = roll_yc(self.index._yc, self.exercise_date)
+
+    @property
+    def pv(self):
+        fp = self.index.forward_pv(self.exercise_date) / self.index.notional
+        T = self.T
+        tilt = np.exp(-self.sigma**2/2 * T + self.sigma * self._Z * math.sqrt(T))
+        rolled_curve = roll_yc(self.index._yc, self.exercise_date)
+        args = (fp, self.exercise_date, self.exercise_date_settle,
+                self.index, self._forward_yc, tilt, self._w)
+        eta = 1.1
+        a = self.index.spread
+        b = self.index.spread * eta
+        while True:
+            if calib(*((b,) + args)) > 0:
+                break
+            b *= eta
+
+        S0 = brentq(calib, a, b, args)
+
+        G = g(self.index, self.strike, self.exercise_date)
+        if T == 0:
+            pv = self.notional * (g(self.index, self.index.spread, self.exercise_date) - G)
+            if self.option_type == "payer":
+                return pv if self.index.spread > self.strike else 0
+            else:
+                return - pv if self.index.spread < self.strike else 0
+
+        Zstar = (math.log(self.strike/S0) + self.sigma**2/2 * T) / \
+                (self.sigma * math.sqrt(T))
+
+        if self.option_type == "payer":
+            Z = Zstar + np.logspace(0, 1.5, 300) - 1
+        elif self.option_type == "receiver":
+            Z = Zstar - np.logspace(0, 1.5, 300) + 1
+        else:
+            raise ValueError("option_type needs to be either 'payer' or 'receiver'")
+        S = S0 * np.exp(-self.sigma**2/2 * T + self.sigma * Z * math.sqrt(T))
+        a, b = strike_vec(S * 1e-4, rolled_curve, self.exercise_date,
+                          self.exercise_date_settle,
+                          self.index.start_date, self.index.end_date, self.index.recovery)
+        val = ((a - b * self.index.fixed_rate*1e-4) - G) * 1/math.sqrt(2*math.pi) * np.exp(-Z**2/2)
+        df_scale = self.index._yc.discount_factor(self.exercise_date_settle)
+        return self.notional * simps(val, Z) * df_scale
+
+    @property
+    def pv2(self):
+        G = g(self.index, self.strike, self.exercise_date)
+        fp = self.index.forward_pv(self.exercise_date) / self.index.notional
+        forward_annuity = self.index.forward_annuity(self.exercise_date)
+        DA_forward_spread = fp / forward_annuity + self.index.fixed_rate * 1e-4
+        strike_tilde = self.index.fixed_rate * 1e-4 + G  / forward_annuity
+        return forward_annuity * black(DA_forward_spread,
+                                       strike_tilde,
+                                       self.T,
+                                       self.sigma,
+                                       self.option_type) * self.notional
+
+    @property
+    def delta(self):
+        old_index_pv = self.index.pv
+        old_pv = self.pv
+        self.index.spread += 0.1
+        notional_ratio = self.index.notional/self.notional
+        delta = (self.pv - old_pv)/(self.index.pv - old_index_pv) * notional_ratio
+        self.index.spread -= 0.1
+        return delta
+
+
+    @property
+    def T(self):
+        if self._T:
+            return self._T
+        else:
+            return year_frac(self.index.trade_date, self.exercise_date) + 1/365
+
+    @property
+    def gamma(self):
+        pass
+
+    @property
+    def theta(self):
+        old_pv = self.pv
+        self._T = self.T - 1/365
+        theta = self.pv - old_pv
+        self._T = None
+        return theta
+
+    @property
+    def vega(self):
+        old_pv = self.pv
+        self.sigma += 0.01
+        vega = self.pv - old_pv
+        self.sigma -= 0.01
+        return vega