from pyisda.flat_hazard import strike_vec from quantlib.time.api import Date import array from scipy.optimize import brentq from scipy.integrate import simps import datetime import numpy as np import pandas as pd from pandas.tseries.offsets import BDay from tranche_functions import GHquad from yieldcurve import roll_yc from scipy.stats import norm def year_frac(d1, d2, day_count_conv = "Actual/365"): """ compute the year fraction between two dates """ if day_count_conv.lower() in ["actual/365", "act/365"]: return (d2-d1).days/365 elif day_count_conv.lower() in ["actual/360", "act/360"]: return (d2-d1).days/360 def calib(S0, fp, exercise_date, exercise_date_settle, index, rolled_curve, tilt, w): S = S0 * tilt * 1e-4 a, b = strike_vec(S, rolled_curve, exercise_date, exercise_date_settle, index.start_date, index.end_date, index.recovery) vec = a - index.fixed_rate * b * 1e-4 return np.inner(vec, w) - fp def g(index, spread, exercise_date, use_rolled_curve = True): """ computes the strike clean price using the expected forward yield curve """ if use_rolled_curve: rolled_curve = roll_yc(index._yc, exercise_date) else: rolled_curve = index._yc step_in_date = exercise_date + datetime.timedelta(days=1) exercise_date_settle = (pd.Timestamp(exercise_date) + 3* BDay()).date() sc = SpreadCurve(exercise_date, rolled_curve, index.start_date, step_in_date, exercise_date_settle, [index.end_date], array.array('d', [spread * 1e-4]), index.recovery) a = index._fee_leg.pv(exercise_date, step_in_date, exercise_date_settle, rolled_curve, sc, True) return (spread - index.fixed_rate) * a *1e-4 def ATMstrike(index, exercise_date): exercise_date_settle = (pd.Timestamp(exercise_date) + 3* BDay()).date() fp = index.forward_pv(exercise_date) closure = lambda S: g(index, S, exercise_date) - fp eta = 1.1 a = index.spread b = index.spread * eta while True: if closure(b) > 0: break b *= eta return brentq(closure, a, b) 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 def d1(F, K, sigma, T): return (np.log(F / K) + sigma**2 * T / 2) / (sigma * np.sqrt(T)) def d2(F, K, sigma, T): return d1(F, K, sigma, T) - sigma * np.sqrt(T) def black(F, K, T, sigma, option_type = "payer"): chi = 1 if option_type == "payer" else -1 if option_type == "payer": return F * norm.cdf(d1(F, K, sigma, T)) - K * norm.cdf(d2(F, K, sigma, T)) else: return K * norm.cdf(- d2(F, K, sigma, T)) - F * norm.cdf(- d1(F, K, sigma, T)) def option(index, exercise_date, sigma, K, option_type="payer"): """ computes the pv of an option using Pedersen's model """ fp = index.forward_pv(exercise_date)/index.notional #forward_yc = yield_curve.expected_forward_curve(exercise_date) #expiry is end of day (not sure if this is right) T = year_frac(index.trade_date, exercise_date) Z, w = GHquad(50) tilt = np.exp(-sigma**2/2 * T + sigma * Z * math.sqrt(T)) exercise_date_settle = (pd.Timestamp(exercise_date) + 3* BDay()).date() args = (fp, exercise_date, exercise_date_settle, index, tilt, w) ## atm forward is greater than spread eta = 1.1 a = index.spread b = index.spread * eta while True: if calib(*((b,) + args)) > 0: break b *= eta S0 = brentq(calib, a, b, args) S = S0 * tilt G = g(index, K, exercise_date) handle = lambda Z: g(index, S0 * math.exp(-sigma**2/2 * T + sigma * Z * math.sqrt(T)), exercise_date) - G Zstar = brentq(handle, -3, 3) if option_type.lower() == "payer": Z = Zstar + np.logspace(0, 1.1, 300) - 1 elif option_type.lower() == "receiver": Z = Zstar - np.logspace(0, 1.1, 300) + 1 else: raise ValueError("option_type needs to be either 'payer' or 'receiver'") S = S0 * np.exp(-sigma**2/2 * T + sigma * Z * math.sqrt(T)) a, b = strike_vec(S, index._yc, exercise_date, exercise_date_settle, index.start_date, index.end_date, index.recovery) val = ((a - b * index.fixed_rate)/df - G) * 1/math.sqrt(2*math.pi) * np.exp(-Z**2/2) return simps(val, Z) * yield_curve.discount_factor(exercise_date_settle) if __name__ == "__main__": import datetime from analytics import Index from swaption import Option ig27_5yr = Index.from_name('ig', 27, '5yr', datetime.date(2016, 10, 24)) ig27_5yr.spread = 74 payer = Option(ig27_5yr, datetime.date(2016, 12, 21), 65) payer.sigma = 0.428 payer.notional = 10e6