import array import datetime import math import numpy as np import pandas as pd from .black import black from .utils import GHquad from .index import g, ForwardIndex from yieldcurve import roll_yc from pandas.tseries.offsets import BDay from pyisda.curve import SpreadCurve from pyisda.flat_hazard import strike_vec from scipy.optimize import brentq from scipy.integrate import simps def calib(S0, fp, exercise_date : datetime.date, exercise_date_settle :datetime.date, 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 ATMstrike(index, exercise_date : datetime.date): exercise_date_settle = (pd.Timestamp(exercise_date) + 3* BDay()).date() fp = index.forward_pv(exercise_date) / index.notional 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 Swaption(ForwardIndex): def __init__(self, index, exercise_date : datetime.date, strike : float, option_type="payer"): ForwardIndex.__init__(self, index, exercise_date) self._exercise_date = exercise_date self._forward_yc = roll_yc(index._yc, exercise_date) self._T = None self._strike = strike self.option_type = option_type.lower() self._Z, self._w = GHquad(50) self.notional = 1 self._G = g(index, strike, exercise_date, self._forward_yc) @property def exercise_date(self): return self._exercise_date @exercise_date.setter def exercise_date(self, d : datetime.date): self._exercise_date = d ForwardIndex.__init__(self, self.index, d) self._forward_yc = roll_yc(self.index._yc, d) self._G = g(self.index, self.strike, self.exercise_date, self._forward_yc) @property def strike(self): return self._strike @strike.setter def strike(self, K : float): self._strike = K self._G = g(self.index, K, self.exercise_date, self._forward_yc) @property def pv(self): T = self.T tilt = np.exp(-self.sigma**2/2 * T + self.sigma * self._Z * math.sqrt(T)) args = (self.forward_pv, self.exercise_date, self.exercise_date_settle, self.index, self._forward_yc, tilt, self._w) eta = 1.1 a = self.index.spread b = a * eta while True: if calib(*((b,) + args)) > 0: break b *= eta S0 = brentq(calib, a, b, args) if T == 0: pv = self.notional * ( g(self.index, self.index.spread, self.exercise_date, self._forward_yc) - self._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, math.log(4 / (self.sigma * math.sqrt(T)), 10), 300) - 1 elif self.option_type == "receiver": Z = Zstar - np.logspace(0, math.log(4 / (self.sigma * math.sqrt(T)), 10), 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, self._forward_yc, 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) - self._G) * 1/math.sqrt(2*math.pi) * np.exp(-Z**2/2) df = self.index._yc.discount_factor(self.exercise_date_settle) return self.notional * simps(val, Z) * df @property def pv_black(self): """compute pv using black-scholes formula""" df = self.index._yc.discount_factor(self.exercise_date_settle) strike_tilde = self.index.fixed_rate * 1e-4 + self._G / self.forward_annuity * df return self.forward_annuity * black(self.forward_spread * 1e-4, 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 ((self.exercise_date - self.index.trade_date).days + 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 @property def DV01(self): old_pv = self.pv self.index.spread += 1 dv01 = self.pv - old_pv self.index.spread -= 1 return dv01 @property def breakeven(self): pv = self.pv / self.notional G = g(self.index, self.strike, self.exercise_date) aux = lambda S: g(self.index, S, self.exercise_date) - G - pv eta = 1.1 a = self.strike b = a * eta while True: if aux(b) > 0: break b *= eta return brentq(aux, a, b)