diff options
Diffstat (limited to 'python')
| -rw-r--r-- | python/calibrate_tranches.py | 63 | ||||
| -rw-r--r-- | python/tranche_functions.py | 58 | ||||
| -rw-r--r-- | python/yieldcurve.py | 30 |
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() |
