aboutsummaryrefslogtreecommitdiffstats
path: root/python
diff options
context:
space:
mode:
Diffstat (limited to 'python')
-rw-r--r--python/calibrate_tranches.py63
-rw-r--r--python/tranche_functions.py58
-rw-r--r--python/yieldcurve.py30
3 files changed, 128 insertions, 23 deletions
diff --git a/python/calibrate_tranches.py b/python/calibrate_tranches.py
new file mode 100644
index 00000000..743a0df7
--- /dev/null
+++ b/python/calibrate_tranches.py
@@ -0,0 +1,63 @@
+import numpy as np
+from tranche_functions import *
+from yieldcurve import YC
+import yaml
+import datetime
+import os
+import pandas as pd
+import pdb
+
+n_int = 500
+n_credit = 100
+Z, w = GHquad(n_int)
+
+with open("../R/index_definitions.yml") as fh:
+ indices = yaml.load(fh)
+indices['hy21']['maturity'] = datetime.date(1970, 1, 1) + datetime.timedelta(indices['hy21']['maturity'])
+hy21 = indices['hy21']
+hy21["startdate"] = datetime.date(2013, 9, 20)
+root = "/home/share/CorpCDOs"
+dates = [f[9:19] for f in os.listdir(os.path.join(root, "data", "Backtest")) if "survprob" in f]
+
+Rho = np.zeros((len(dates), 3))
+for i, d in enumerate(dates):
+ startdate = datetime.datetime.strptime(d, "%Y-%m-%d")
+ ts = YC(startdate)
+ with open(os.path.join(root, "data", "Backtest", "recov_{0}.csv".format(d))) as fh:
+ recov = np.array([float(e) for e in fh], dtype='double', order='F')
+
+ with open(os.path.join(root, "data", "Backtest", "survprob_{0}.csv".format(d))) as fh:
+ fh.readline() ##skip header
+ SurvProb = np.array([[float(e) for e in line.split(",")] for line in fh], dtype='double', order='F')
+
+ defaultprob = 1 - SurvProb
+ issuerweights = np.ones(100)/100
+
+ rho = 0.4
+ Ngrid = 101
+
+ K = np.array([0, 0.15, 0.25, 0.35, 1])
+
+ Kmod = adjust_attachments(K, hy21["loss"], hy21["factor"])
+ quotes = pd.read_csv(os.path.join(root, "Scenarios", "Calibration",
+ "hy21_tranches_{0}.csv".format(d)))
+ quotes = quotes["Mid"]/100
+ dK = np.diff(Kmod)
+ quotes = np.cumsum(dK * (1-quotes))
+ sched = creditSchedule(startdate, "5Yr", 0.05, ts, enddate=hy21["maturity"])
+ acc = cdsAccrued(startdate, 0.05)
+ for j, q in enumerate(quotes[:-1]):
+ def aux(rho):
+ L, R = BClossdist(defaultprob, issuerweights, recov, rho, Z, w, 101)
+ cl = tranche_cl(L, R, sched, 0, Kmod[j+1])
+ pl = tranche_pl(L, sched, 0, Kmod[j+1])
+ return cl+pl+q-acc
+ l, u = (0, 1)
+ for _ in range(10):
+ rho = (l+u)/2.
+ if aux(rho) > 0:
+ u = rho
+ else:
+ l = rho
+ Rho[i, j] = (l+u)/2.
+ print(Rho[i,:])
diff --git a/python/tranche_functions.py b/python/tranche_functions.py
index 6e095ed5..8cfbb116 100644
--- a/python/tranche_functions.py
+++ b/python/tranche_functions.py
@@ -1,6 +1,9 @@
import numpy as np
from ctypes import *
-import pdb
+from quantlib.time.schedule import Schedule, CDS
+from quantlib.time.api import Actual360, Period, UnitedStates, Following, today
+from quantlib.util.converter import qldate_to_pydate, pydate_to_qldate
+import pandas as pd
libloss = np.ctypeslib.load_library("lossdistrib", "/home/share/CorpCDOs/code/R")
libloss.fitprob.restype = None
@@ -66,20 +69,27 @@ def BClossdist(defaultprob, issuerweights, recov, rho, Z, w, Ngrid = 101, defaul
byref(c_int(Ngrid)), byref(c_int(defaultflag)), L, R)
return L, R
+def adjust_attachments(K, losstodate, factor):
+ """
+ computes the attachments adjusted for losses
+ on current notional
+ """
+ return np.minimum(np.maximum((K-losstodate)/factor, 0), 1)
+
def trancheloss(L, K1, K2):
- np.maximum(L - K1, 0) - np.maximum(L - K2, 0)
+ return np.maximum(L - K1, 0) - np.maximum(L - K2, 0)
def trancherecov(R, K1, K2):
- np.maximum(R - 1 + K2, 0) - np.maximum(R - 1 + K1, 0)
+ return np.maximum(R - 1 + K2, 0) - np.maximum(R - 1 + K1, 0)
def tranche_cl(L, R, cs, K1, K2, scaled = False):
if(K1 == K2):
return 0
else:
support = np.linspace(0, 1, L.shape[0])
- size = K2 - K1 - np.dot(trancheloss(support, K1, K2), L) -
- np.dot(trancherecov(support, K1, K2), R)
- sizeadj = 0.5 * (size + np.hstack([K2-K1], size[:-1]))
+ size = K2 - K1 - np.dot(trancheloss(support, K1, K2), L) - \
+ np.dot(trancherecov(support, K1, K2), R)
+ sizeadj = 0.5 * (size + np.hstack((K2-K1, size[:-1])))
if scaled:
return 1/(K2-K1) * np.dot(sizeadj * cs["coupons"], cs["df"])
else:
@@ -91,11 +101,41 @@ def tranche_pl(L, cs, K1, K2, scaled=False):
else:
support = np.linspace(0, 1, L.shape[0])
cf = K2 - K1 - np.dot(trancheloss(support, K1, K2), L)
- cf = np.hstack([K2-K1], cf)
+ cf = np.hstack((K2-K1, cf))
if scaled:
return 1/(K2-K1) * np.dot(np.diff(cf), cs["df"])
else:
return np.dot(np.diff(cf), cs["df"])
-def tranche_pv(L, R, cs, K2, K2):
- tranche_pl(L, cs, K1, K2) + tranche_cl(L, R, cs, K2, K2)
+def tranche_pv(L, R, cs, K1, K2):
+ return tranche_pl(L, cs, K1, K2) + tranche_cl(L, R, cs, K2, K2)
+
+def creditSchedule(tradedate, tenor, coupon, yc, enddate = None):
+ tradedate = pydate_to_qldate(tradedate)
+ start = tradedate - Period('4Mo')
+ enddate = pydate_to_qldate(enddate)
+ if enddate is None:
+ enddate = tradedate + Period(tenor)
+ cal = UnitedStates()
+ DC = Actual360()
+ sched = Schedule(start, enddate, Period('3Mo'), cal, date_generation_rule=CDS)
+ prevpaydate = sched.previous_date(tradedate)
+ sched = [d for d in sched if d>=prevpaydate]
+ df = [yc.discount(d) for d in sched if d > tradedate]
+
+ dates = pd.to_datetime([str(d) for d in sched if d > tradedate], "%m/%d/%Y")
+ coupons = np.diff([DC.year_fraction(prevpaydate, d) * coupon for d in sched])
+
+ return pd.DataFrame({"df":df, "coupons":coupons}, dates)
+
+def cdsAccrued(tradedate, coupon):
+ tradedate = pydate_to_qldate(tradedate)
+ start = tradedate - Period('3Mo')
+ end = tradedate + Period('3Mo')
+
+ start_protection = tradedate + 1
+ DC = Actual360()
+ cal = UnitedStates()
+ sched = Schedule(start, end, Period('3Mo'), cal, date_generation_rule=CDS)
+ prevpaydate = sched.previous_date(tradedate)
+ return DC.year_fraction(prevpaydate, start_protection) * coupon
diff --git a/python/yieldcurve.py b/python/yieldcurve.py
index b524a04c..50588aa6 100644
--- a/python/yieldcurve.py
+++ b/python/yieldcurve.py
@@ -4,10 +4,10 @@ import requests, zipfile
from io import BytesIO
import xml.etree.ElementTree as ET
import datetime
-from quantlib.time.api import Calendar, Date, Period, Days, Schedule, today
+from quantlib.time.api import Calendar, Period, Days, Schedule, today, Actual360
from quantlib.time import imm
from quantlib.util.converter import qldate_to_pydate, pydate_to_qldate
-from quantlib.market.market import libor_market
+from quantlib.market.market import libor_market, next_imm_date
import numpy as np
import matplotlib.pyplot as plt
@@ -45,30 +45,32 @@ def get_futures_data(date = datetime.date.today()):
def YC(date = datetime.date.today(), MarkitData=None, futures = None):
if not MarkitData:
- MarkitData = getMarkitIRData(date)
+ MarkitData = getMarkitIRData(date.date())
if not futures:
- calendar = Calendar.from_name('TARGET')## need to figure out which is the right calendar to use
- prev_day = calendar.advance(pydate_to_qldate(date), -1, Days)
- futures = get_futures_data(qldate_to_pydate(prev_day).date())
+ futures = get_futures_data(date.date())
m = libor_market('USD(NY)')
quotes = [('ED',i+1, v) for i, v in enumerate(futures)]
+ # if next_imm_date(date, 9) == pydate_to_qldate(date) + Period('2Yr'):
+ # quotes.pop(8)
quotes += [('SWAP', k, v) for k, v in MarkitData['swaps'].items()]
m.set_quotes(date, quotes)
- m.bootstrap_term_structure()
- return m
+ ts = m.bootstrap_term_structure(interpolator='linear')
+ return ts
if __name__=="__main__":
- m = YC()
+ date = datetime.date(2013, 12, 23)
+ ts = YC(date)
+ cal = Calendar.from_name('USA')
p1 = Period('1Mo')
p2 = Period('2Mo')
p3 = Period('3Mo')
p6 = Period('6Mo')
p12 = Period('12Mo')
- calendar = Calendar.from_name(m.calendar)
- sched = Schedule(m.settle_date, m.settle_date+Period('5Yr'), Period('3Mo'), calendar)
+ sched = Schedule(ts.reference_date, ts.reference_date+Period('5Yr'), Period('3Mo'), cal)
days = [qldate_to_pydate(day) for day in sched]
- f3 = [m.forward(d, p3, 1) for d in sched]
- f6 = [m.forward(d, p6, 1) for d in sched]
- f2 = [m.forward(d, p2, 1) for d in sched]
+ f3 = [ts.forward_rate(d, d+p3, Actual360()).rate for d in sched]
+ f6 = [ts.forward_rate(d, d+p6, Actual360()).rate for d in sched]
+ f2 = [ts.forward_rate(d, d+p2, Actual360()).rate for d in sched]
+
plt.plot(days, f2, days, f3, days, f6)
plt.show()