"The aim of the project is to satisfy Redington's theory of immunization." 
"Here we have to determine the portfolio of bonds that immunizes the liabilities."

"Assumption  : Bonds pay interest annually"

"Step 1: determine the YTM of the bonds"
"Step 2: Calculate the effective duration and convexity of of the bonds using the calculated YTM."
"Step 3: Use a weighted average YTM to discount the liabilities." 
"Step 4: Calculate the modified duration and convexity of the liabilities."
"Step 5: Minimize portfolio value to satisfy Redington's Theory of Immunization."



import numpy as np
import pandas as pd
from scipy import optimize


#%%
newdata = pd.read_excel(r"C:/Users/Dell/Desktop/Bonds.xlsx" , sheet_name=[0, 1, 2])
bonds_data = newdata[0]
#%%

"As there are a list of 50 bonds , we have got it down to just 35"
"We have filtered using time to maturity with the max time to maturity as 21Yrs and min time to maturity 1 Yr"
bonds_data = bonds_data[bonds_data['time_to_maturity']<=21]
bonds_data = bonds_data[bonds_data['time_to_maturity']>=1]

#%%

liabilities_data = newdata[1]
liabilities_data = liabilities_data.rename({"Cashflow": "Amount", "Due Time(In Years)": "Time"}, axis="columns")

#%%

# STEP 1: 

def cashflows(time, coupon_payment):
    data = pd.DataFrame(columns=["Time", "Coupon", "Principal", "Total Cashflow"])
    coupon_payment = coupon_payment / 100
    if time < 0.5:
        data.loc[0, "Time"] = 1
        data.loc[0, "Coupon"] = coupon_payment * 100
        data.loc[0, "Principal"] = 100
        data["Total Cashflow"] = data["Coupon"] + data["Principal"]
        return data
    else:
        for i in range(int(round(time, 0))):
            data.loc[i, "Time"] = i + 1
            data.loc[i, "Coupon"] = coupon_payment * 100
            if data.loc[i, "Time"] == int(round(time, 0)):
                data.loc[i, "Principal"] = 100
            else:
                data.loc[i, "Principal"] = 0
        data["Total Cashflow"] = data["Coupon"] + data["Principal"]
        return data

def YTM(ytm, time, coupon_payment, original_price):
    cash_flows = cashflows(time, coupon_payment)
    price_of_bond = (cash_flows["Total Cashflow"] * ((1 + ytm) ** (-1 * cash_flows["Time"].astype(float)))).sum()
    return price_of_bond - original_price

#%%
#Step2
bonds_data["YTM"] = 0
for i in range(len(bonds_data)):
    try:
        time = bonds_data.loc[i, "time_to_maturity"]
        coupon_payment = bonds_data.loc[i, "current_coupon"]
        original_price = bonds_data.loc[i, "curve_based_price"]
        bonds_data.loc[i, "YTM"] = optimize.fsolve(YTM, 0.05, args=(time, coupon_payment, original_price))[0]
    except:
        continue
del original_price, i, time, coupon_payment

bonds_data["YTMWeight"] = (bonds_data["trade_price"] * bonds_data["trade_size"]) / (bonds_data["trade_price"] * bonds_data["trade_size"]).sum()


#%%

#Step 3

discounting = (bonds_data["YTM"] * bonds_data["YTMWeight"]).sum()
v = 1 / (1 + discounting)
bonds_data.drop(columns=["YTMWeight"], inplace=True)



#%%

#Step 4 and Step 5

liabilities_data["PV"] = liabilities_data["Amount"] * ((1 + discounting) ** (-liabilities_data["Time"]))
pv_liabilities_data = liabilities_data["PV"].sum()

def bond_price(ytm, time, coupon_payment):  # function to get bond price.
    cash_flows = cashflows(time, coupon_payment)
    bond_price = (cash_flows["Total Cashflow"]* ((1 + ytm) ** (-1 * cash_flows["Time"].astype(float)))).sum()
    return bond_price

bonds_data["ModifiedDuration"] = 0
bonds_data["Convexity"] = 0
for i in range(len(bonds_data)):
    try:            
        time = bonds_data.loc[i, "time_to_maturity"]
        coupon_payment = bonds_data.loc[i, "current_coupon"]
        ytm = bonds_data.loc[i, "YTM"]
        PV0 = bonds_data.loc[i, "curve_based_price"]
        PV1 = bond_price(ytm - 0.0001, time, coupon_payment)
        PV2 = bond_price(ytm + 0.0001, time, coupon_payment)
        bonds_data.loc[i, "ModifiedDuration"] = (PV1 - PV2) / (2 * 0.0001 * PV0)
        bonds_data.loc[i, "Convexity"] = (PV1 + PV2 - (2 * PV0)) / ((0.0001 ** 2) * PV0)
    except:
        continue

c = liabilities_data["Amount"]
t = liabilities_data["Time"]
ModDur_l = ((t * c * (v ** t)).sum() / (c * (v ** t)).sum()) / (1 + discounting)

liabilities_data["v"] = (1 + discounting) ** (-liabilities_data["Time"])
liabilities_data["v'"] = (-liabilities_data["Time"]) * ((1 + discounting) ** (-liabilities_data["Time"] - 1))
liabilities_data["v''"] = ((-liabilities_data["Time"]) * (-liabilities_data["Time"] - 1)) * ((1 + discounting) ** (-liabilities_data["Time"] - 2))
convexity_l = (liabilities_data["v''"] * liabilities_data["Amount"]).sum()
PVb = np.array(bonds_data["curve_based_price"])
lhs_eq1 = np.array([PVb, bonds_data["ModifiedDuration"] * PVb])
rhs_eq1 = np.array([pv_liabilities_data, ModDur_l * pv_liabilities_data])
opt = optimize.linprog(c=PVb, A_eq=lhs_eq1, b_eq=rhs_eq1)
X = opt.x
weights = (X * PVb) / (X * PVb).sum()
convexity_a = (opt.x * PVb * bonds_data["Convexity"]).sum()
bonds_data["Number of Bonds"] = opt.x
final = pd.DataFrame({"Bond ID": bonds_data["bond_id"],"Price": bonds_data["curve_based_price"],"Number Of Each Bonds to Buy": opt.x,})

#%%
print(final)
print('The Bond Portfolio for immunization has been printed above')
print('PV Of Assets is ' + str(int((opt.x * PVb).sum()))+ ' =' + ' PV Of Liabilities ' + str(int(pv_liabilities_data)))
print('Duration Of Assets is ' + str(int(((weights * bonds_data['ModifiedDuration']).sum()))) + ' =' + ' Duration Of Liabilities is ' + str(int(ModDur_l))) 
print('Convexity Of Assets > Convexity Of Liabilities as ' + str(int(convexity_a)) + ">" + str(int(convexity_l)))
#%%