import pandas as pd
import numpy as np
bonds = pd.read_excel("Project_Data.xlsx")
liabilities = pd.read_excel("Project_Data.xlsx", sheet_name='Cash outflow (Liability)')
#We will only consider imp bond details and the trade details of bond is of no
#use to get the cashflows of the bond so we extract only curve based price, 
#coupon rate and time to maturity of the bond
cols = ['bond_id','current_coupon','time_to_maturity','curve_based_price','trade_size']
bonds = bonds[cols]
#the data has 50 bonds.
#the curve based price would be our price of the bond and will be used to 
#calculate YTM of the bond. 
#we assume that we get coupon anually and so for example if time to maturity 
#is 4.75 then coupon payments will be at 0.75,1.75,2.75,3.75 and 4.75.

#calculating certain values for easier calculations
#as we need to calculate cashflows for non integer values of time, so
#we seperate no of cashflows and ref_no so we can easily run for loop in 
#our calculations.
list = []
for i in range(len(bonds['time_to_maturity'])):
    list.append(int(bonds['time_to_maturity'].iloc[i])+1)
bonds['no_of_cashflows'] = list
list_2 = []
for i in range(len(bonds['time_to_maturity'])):
    list_2.append(float(bonds['time_to_maturity'].iloc[i])-int(bonds['time_to_maturity'].iloc[i]))
bonds['ref_no'] = list_2

#functions for calculating bond metrics
def bond_pv(coupon,int_rate,no_of_cashflows,ref_no):
    sum = 0
    for i in range(no_of_cashflows):
        sum = sum + (coupon*(1/((1+int_rate)**(i+ref_no))))
    sum = sum + (100*(1/((1+int_rate)**(no_of_cashflows-1+ref_no))))
    return sum
def find_ytm(coupon,no_of_cashflows,ref_no,price):
    from scipy import optimize 
    get_yeild = lambda int_rate : bond_pv(coupon,int_rate,no_of_cashflows,ref_no) - price 
    return optimize.newton(get_yeild, coupon/100)
def mod_dur(coupon,int_rate,no_of_cashflows,ref_no):
    sum = 0
    for i in range(no_of_cashflows):
        sum = sum + ((coupon*(1/((1+int_rate)**(i+ref_no))))*(i+ref_no))
    sum = sum + ((100*(1/((1+int_rate)**(no_of_cashflows-1+ref_no))))*(no_of_cashflows-1+ref_no))   
    macdur = sum / (bond_pv(coupon,int_rate,no_of_cashflows,ref_no))
    moddur = macdur/(1 + int_rate)      
    return moddur
def convexity(coupon,int_rate,no_of_cashflows,ref_no):
    num = bond_pv(coupon,int_rate-0.01,no_of_cashflows,ref_no)+bond_pv(coupon,int_rate+0.01,no_of_cashflows,ref_no)-(2*bond_pv(coupon,int_rate,no_of_cashflows,ref_no))
    den = 2*bond_pv(coupon,int_rate,no_of_cashflows,ref_no)*0.01*0.01
    return num/den

#updating bond metrics in table
#ytm
ytm_list = []
for i in range(len(bonds['time_to_maturity'])):
    ytm_list.append(find_ytm(bonds['current_coupon'].iloc[i],bonds['no_of_cashflows'].iloc[i],bonds['ref_no'].iloc[i],bonds['curve_based_price'].iloc[i]))
bonds['ytm'] = ytm_list

#modified duration
moddur_list = []
for i in range(len(bonds['time_to_maturity'])):
    moddur_list.append(-1*mod_dur(bonds['current_coupon'].iloc[i],bonds['ytm'].iloc[i],bonds['no_of_cashflows'].iloc[i],bonds['ref_no'].iloc[i]))
bonds['modified_duration'] = moddur_list

#convexity
convexity_list = []
for i in range(len(bonds['time_to_maturity'])):
    convexity_list.append(convexity(bonds['current_coupon'].iloc[i],bonds['ytm'].iloc[i],bonds['no_of_cashflows'].iloc[i],bonds['ref_no'].iloc[i]))
bonds['convexity'] = convexity_list

#next step is we need to find present value, modified duration and convexity of
#liabilities 
#to find interest rates to use for liabilities, we will find weighted average 
#of ytms based on liquidity of the bond (this is where we will use trade amount)
sum_tradesize = 0
for i in np.array(bonds['trade_size']):
    sum_tradesize = sum_tradesize + i
ytm_liabilities = 0
tradesize_list = np.array(bonds['trade_size'])
for i in range(len(bonds['time_to_maturity'])):
    w = tradesize_list[i]/sum_tradesize 
    ytm_liabilities += w*(bonds['ytm'].iloc[i])
    
#finding metrics of liabilities to be used (pv, modified duration and convexity)
#finding pv of liabilities
pv_liabilities = 0
for i in range(10):
    pv_liabilities += liabilities['Cashflow'].iloc[i]*((1 + ytm_liabilities)**(-liabilities['Due Time(In Years)'].iloc[i]))
#finding modified duration of liabilities
moddur_liabilities = 0
sum_modcal = 0
for i in range(10):
    sum_modcal += (liabilities['Cashflow'].iloc[i]*((1 + ytm_liabilities)**(-liabilities['Due Time(In Years)'].iloc[i])))*(liabilities['Due Time(In Years)'].iloc[i])
moddur_liabilities = (-sum_modcal/pv_liabilities)*(1/(1+ytm_liabilities))
#finding convexity of liabilities
pv_ = 0
pv1 = 0
for i in range(10):
    pv_ += liabilities['Cashflow'].iloc[i]*((1 + (ytm_liabilities-0.01))**(-liabilities['Due Time(In Years)'].iloc[i]))
for i in range(10):
    pv1 += liabilities['Cashflow'].iloc[i]*((1 + (ytm_liabilities+0.01))**(-liabilities['Due Time(In Years)'].iloc[i]))
convexity_liabilities = (pv_+pv1-(2*pv_liabilities))/(2*0.01*0.01*pv_liabilities)

#Now we have the required metrics of the assests and liabilities, we will 
#create reningtons equations and use it as constraints to perform linear 
#programming and find the units of each bond to buy at time 0 to immunise 
#our portfolio bearing minimum cost.
#we have reningtons principle to create three equations, the principles are:
#1) PV of Assests = PV of Liabilities
#2) Modified duration of assests = modified duration of liabilities 
#3) Convexity of Assests > Convexity of Liabilities
#no of units to be bough of each bond will be our variable to find. So we 
#assume we buy x1 units of bond 1,x2 units of bond 2,.....x50 units of bond 50.
#so pv of assests = units of bond*bond price 
#=x1*price of bond 1 + x2*price of bond 2.....
#this should be equal to pv of liablities making our first equation
#eq 1: 105.506x1 + 102.199x2 +.....+ 99.3448x50 = 341560957.014596
#second equation will be formed using second rule of renington immunisation
#which states that modified duration of assets and liabilities should be equal

#NOTE: modified duration of bond portfolio is weighted average of individual 
#bonds modified duration where weights are the price of the bonds

#105.506x1/(price of portfolio)*moddur of bond 1 + ...... = -8.210863321962403
#we know the denominator is total price of the portfolio which is equal to
#the pv of assets which is equal to pv of liablities so
#effectively eq 2: 105.506*moddur*x1 +.... = 341560957*-8.210863321962403
#similarly like modified duration, we will also calculate convexity 
#eq 3: 105.506*x1*convexity + ..... = 341560957*54.001684937246935

#once we have got our equations and variables, we will use linear programming 
#to get a portfolio with minimum cost and maximum profit with former being our
#priority cause our main objective is to immunise our portflio, our alternate 
#objective is to get max profit provided our portfolio is immunised.

from scipy.optimize import linprog

#linprog is used to minimize the objective function. our objective function is 
#price of the portfolio which is equal to the present value of the portfolio
#as this is fundamentals of bond pricing that its FAIR price is equal to PV of 
#expected cashflows from the bond.

obj = np.array(bonds['curve_based_price'])

#to make the equations, we will need to find out values of some variables for 
#each bond, namely price*modified duration and price*convexity 
list_Pmod = np.array(bonds['curve_based_price'])*np.array(bonds['modified_duration'])
list_Pcon = -np.array(bonds['curve_based_price'])*np.array(bonds['convexity'])

#setting up equations
lhs_eq = [obj,list_Pmod]
rhs_eq = [pv_liabilities,pv_liabilities*moddur_liabilities]
lhs_ineq = [list_Pcon]
rhs_ineq = [pv_liabilities*convexity_liabilities]
bnd = [(0,float("inf")) for q in range(len(obj))]

#linear programming using Scipy
#we use reverse simplex because this is the only method where the optimal 
#solution is found
opt = linprog(c=obj,A_ub = lhs_ineq, b_ub = rhs_ineq, A_eq = lhs_eq, b_eq = rhs_eq,bounds = bnd, method="revised simplex")
print(opt.success)  #wether or not the optimal solution is found
print(opt.fun) #value of ourobjective function
print(opt.x)  #gives the no of units of each bond to be bought
#1440407.58601837 units of bond 1 and 1443647.21515534 units of bond 7

#checking if the portfolio passes all constraints
#1) checking for equation 1
print(opt.fun == pv_liabilities) 
# output : True
#2) checking for equation 2
moddur_bond1 = bonds['modified_duration'].iloc[0]
moddur_bond7 = bonds['modified_duration'].iloc[6]
price_bond1 = bonds['curve_based_price'].iloc[0]
price_bond7 = bonds['curve_based_price'].iloc[6]
moddur_assets = (((opt.x[0]*price_bond1)/opt.fun)*moddur_bond1) + (((opt.x[6]*price_bond7)/opt.fun)*moddur_bond7)
print(moddur_assets == moddur_liabilities)
# output : True
#3) Checking for equation 3
convexity_bond1 = bonds['convexity'].iloc[0]
convexity_bond7 = bonds['convexity'].iloc[6]
convexity_assets = (((opt.x[0]*price_bond1)/opt.fun)*convexity_bond1) + (((opt.x[6]*price_bond7)/opt.fun)*convexity_bond7)
print(convexity_assets > convexity_liabilities)
# output : True

#so this portfolio satisfies all conditions of reningtons immunisation theory

#Our Portfolio will be:
#buy 1440407.5860183693 units of bond with bond id 1 and 1443647.2151553365 
#units of bond with bond id 7 at time 0

#bonds info:
    
#bond ID 1 
#curve_based_price : 105.506127
#current_coupon : 9.65% p.a.
#time_to_maturity : 0.652906 years
#coupon payment frequency : Annual

#Bond ID 7
#curve_based_price : 131.326496
#current_coupon : 6.125% p.a.
#time_to_maturity : 26.551478 years
#coupon payment frequency : Annual