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): """Swaption class""" 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 intrinsic_value(self): V = self.df * (self.forward_pv - self._G) return max(V, 0) if self.option_type == "payer" else max(-V, 0) @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.05 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: return self.notional * self.intrinsic 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) return self.notional * simps(val, Z) * self.df @pv.setter def pv(self, val: float): if val < self.intrinsic_value: raise ValueError("{}: is less than intrinsic value: {}". format(val, self.intrinsic_value)) def handle(x): self.sigma = x return self.pv - val # use sigma_black as a starting point self.pv_black = val eta = 1.1 a = self.sigma while True: if handle(a) < 0: break a /= eta b = a * eta while True: if handle(b) > 0: break b *= eta self.sigma = brentq(handle, a, b) @property def pv_black(self): """compute pv using black-scholes formula""" strike_tilde = self.index.fixed_rate * 1e-4 + self._G / self.forward_annuity * self.df return self.forward_annuity * black(self.forward_spread * 1e-4, strike_tilde, self.T, self.sigma, self.option_type) * self.notional @pv_black.setter def pv_black(self, val: float): if val < self.intrinsic_value: raise ValueError("{}: is less than intrinsic value: {}". format(val, self.intrinsic_value)) def handle(x): self.sigma = x return self.pv_black - val eta = 1.01 a = 0.1 b = a * eta while True: if handle(b) > 0: break b *= eta self.sigma = brentq(handle, a, b) @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)