1 files changed, 161 insertions, 0 deletions
diff --git a/python/swaption.py b/python/swaption.py
new file mode 100644
index 00000000..4e81acb2
--- /dev/null
+++ b/python/swaption.py
@@ -0,0 +1,161 @@
+from pyisda.legs import ContingentLeg, FeeLeg
+from pyisda.flat_hazard import strike_vec
+from pyisda.curve import YieldCurve, BadDay, SpreadCurve
+import math
+from scipy.optimize import brentq
+from scipy.integrate import simps
+import datetime
+from tranche_functions import GHquad
+
+class Index():
+    """ minimal class to represent a credit index """
+    def __init__(self, start_date, end_date, recovery, fixed_rate):
+        self.fixed_rate = fixed_rate
+        self.notional = 1
+        self._start_date = start_date
+        self._end_date = end_date
+        self.recovery = recovery
+
+        self._fee_leg = FeeLeg(start_date, end_date, True, 1, 1)
+        self._default_leg = ContingentLeg(start_date, end_date, 1)
+
+    @property
+    def start_date(self):
+        return self._start_date
+
+    @property
+    def end_date(self):
+        return self._end_date
+
+    @start_date.setter
+    def start_date(self, d):
+        self._fee_leg = FeeLeg(d, self.end_date, True, 1, 1)
+        self._default_leg = ContingentLeg(d, self.end_date, 1)
+        self._start_date = d
+
+    @end_date.setter
+    def end_date(self, d):
+        self._fee_leg = FeeLeg(self.start_date, d, True, 1, 1)
+        self._default_leg = ContingentLeg(self.start_date, d, 1)
+        self._end_date = d
+
+    def pv(self, trade_date, exercise_date, yc, sc, h = None):
+        step_in_date = exercise_date + datetime.timedelta(days=1)
+        a = self._fee_leg.pv(trade_date, step_in_date, trade_date, yc, sc, False)
+        Delta = self._fee_leg.accrued(step_in_date)
+        if exercise_date > trade_date:
+            Delta *= math.exp(-h * year_frac(trade_date, exercise_date))
+        clean_forward_annuity = a - Delta
+
+        dl_pv = self._default_leg.pv(
+            trade_date, step_in_date, trade_date, yc, sc, self.recovery)
+        return dl_pv - clean_forward_annuity * self.fixed_rate, clean_forward_annuity
+
+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 DAforward_price(spread, trade_date, exercise_date, yield_curve, index):
+    cash_settle_date = (pd.Timestamp(trade_date) + 3* BDay()).date()
+    h, sc = flat_hazard(spread, yc, trade_date, cash_settle_date,
+                        start_date = trade_date,
+                        end_date = index.end_date,
+                        recovery_rate = index.recovery)
+    exercise_date_settle = (pd.Timestamp(exercise_date) + 3* BDay()).date()
+    forward_price = index.pv(trade_date, exercise_date, yield_curve, sc, h)[0]
+    fep = (1 - index.recovery) * (1-math.exp(-h*year_frac(trade_date, exercise_date)))
+    df = yc.discount_factor(exercise_date_settle)
+    price = fep + 1/df * forward_price
+    return price
+
+def flat_hazard(spread, yc, trade_date=datetime.date.today(),
+                cash_settle_date = None,
+                start_date = datetime.date.today(),
+                end_date = datetime.date(2021, 6, 20),
+                recovery_rate = 0.4):
+    step_in_date = start_date + datetime.timedelta(days=1)
+    if cash_settle_date is None:
+        cash_settle_date = (pd.Timestamp(trade_date) + 3* BDay()).date()
+    sc = SpreadCurve(trade_date, yc, trade_date, step_in_date,
+                     cash_settle_date,
+                     [end_date], array.array('d', [spread]), recovery_rate)
+    sc_data = sc.inspect()['data']
+    ## conversion to continuous compounding
+    hazard_rate = math.log(1 + sc_data[0][1])
+    return (hazard_rate, SpreadCurve.from_flat_hazard(trade_date, hazard_rate))
+
+def calib(S0, fp, forward_yield_curve, exercise_date_settle, index, tilt, w):
+    S = S0 * tilt
+    a, b = strike_vec(S, forward_yield_curve, exercise_date_settle,
+                      index.start_date, index.end_date, index.recovery)
+    vec = a - index.fixed_rate * b
+    df = forward_yield_curve.discount_factor(exercise_date_settle)
+    return 1/df*np.inner(vec - fp, w)
+
+def g(spread, forward_yield_curve, index):
+    """ computes the strike price using the expected forward yield curve """
+    exercise_date = forward_yield_curve.base_date
+    step_in_date = exercise_date + datetime.timedelta(days=1)
+    exercise_date_settle = (pd.Timestamp(exercise_date) + 3* BDay()).date()
+    sc = SpreadCurve(exercise_date, forward_yield_curve, exercise_date, step_in_date,
+                     exercise_date_settle,
+                     [index.end_date], array.array('d', [spread]), index.recovery)
+    a = index._fee_leg.pv(exercise_date, step_in_date, exercise_date,
+                          forward_yield_curve, sc, True)
+    dl_pv = index._default_leg.pv(
+        exercise_date, step_in_date, exercise_date, forward_yield_curve,
+        sc, index.recovery)
+    df = forward_yield_curve.discount_factor(exercise_date_settle)
+    return 1/df * (dl_pv - a * index.fixed_rate)
+
+def ATMstrike(spread, trade_date, exercise_date, yield_curve, index):
+    fp = DAforward_price(spread, trade_date, exercise_date, yc, index)
+    yc_forward = yc.expected_forward_curve(exercise_date)
+    closure = lambda S: g(S, yc_forward, index) - fp
+    eta = 1.1
+    a = spread
+    b = spread * eta
+    while True:
+        if closure(b) > 0:
+            break
+    return brentq(closure, a, b)
+
+def option(index, ref, trade_date, exercise_date, yield_curve, sigma, K, option_type="payer"):
+    """ computes the pv of an option using Pedersen's model """
+    fp = DAforward_price(ref, trade_date, exercise_date, yield_curve, index)
+    forward_yc = yield_curve.expected_forward_curve(exercise_date)
+    #expiry is end of day (not sure if this is right)
+    T = year_frac(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, forward_yc, exercise_date_settle, index, tilt, w)
+    ## atm forward is greater than spread
+    eta = 1.1
+    a = ref
+    b = ref * eta
+    while True:
+        if calib(*((b,) + args)) > 0:
+            break
+        b *= eta
+    S0 = brentq(calib, a, b, args)
+    S =  S0 * tilt
+    G = g(K, forward_yc, index)
+    handle = lambda Z: g(S0 * math.exp(-sigma**2/2 * T + sigma * Z * math.sqrt(T)),
+                         forward_yc, index) - 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))
+    df = forward_yc.discount_factor(exercise_date_settle)
+    a, b = strike_vec(S, forward_yc, 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)