# =============================================================================
#                   """ AOI - Python and FIP Project """
# =============================================================================

"""
Group:
431. Kshitij Srivastava
432. Priyanshu Mehta
433. Darshan Lodha
434. Sabyasachi Rathore
435. Ajay Dhandha
"""

# =============================================================================
# """
# Assumptions:
# 1. The periodicity of all bonds is assumed to be annual
# 2. There is no default risk
# 3. The rate at which the liabilities are discounted is assumed to be the weighted
#    average of the ytm of the selected bond portfolio.
# 4. Coupons are assumed to be paid in arrears
# """ 
# =============================================================================

""" Importing necessary libraries """
import pandas as pd
import numpy as np
# import scipy
# from scipy import optimize
from scipy.optimize import linprog

""" Importing the bonds dataset """
df_main = pd.read_excel(r"C:/Users/S145 DOIN/Downloads/Project Data.xlsx",sheet_name="Bonds")

""" Deleting redundant columns """
df_bonds = df_main.drop(["weight","reporting_delay","trade_size","trade_type","trade_price_last1","trade_size_last1"],axis = 1)

""" Selecting a portfolio of bonds out of given 50 bonds"""

#Selecting bonds in the portfolio that have a time to maturity of more than 1 """
df_bonds_high_tm = df_bonds[df_bonds["time_to_maturity"] > 1]

#Further, selecting bonds having a time to maturity between 1 and 21 """
df_bonds_new = df_bonds_high_tm[(df_bonds_high_tm["time_to_maturity"] > 1.00) & (df_bonds_high_tm["time_to_maturity"] < 21.00)]

#Further, removing the bonds having proximate time to  maturity """
df_bonds_filtered = df_bonds_new.loc[df_bonds_new["bond_id"].isin([50,19,36,21,17,37,34,35,8,6,40,24,46,27,44,45,5,48,15,25])]

""" Calculating YTM for all 20 bonds in the portfolio"""
# A class Bond_Coupon has been constructed consisting of user defined functions 
# to compute numerous bond parameters like price and ytm.

class Bond_Coupon:
    # Defining a function to compute bond price
    def get_price(self,coupon,face_value,int_rate,years,freq=1):
        total_coupons_pv = self.get_coupons_pv(coupon,int_rate,years,freq)
        face_value_pv = self.get_face_value_pv(face_value,int_rate,years)
        result = total_coupons_pv + face_value_pv
        return result
    
    # Defining a function to compute pv of face value
    @staticmethod
    def get_face_value_pv(face_value,int_rate,years):
       fvpv = face_value/(1+int_rate)**years
       return fvpv
   
    # Defining a function to compute pv of coupon payments
    def get_coupons_pv(self,coupon,int_rate,years,freq = 1):
       pv = 0
       for period in np.arange((years - int(years)),years+0.1,1):
           pv += self.get_coupon_pv(coupon, int_rate, period,freq)
       return pv
    
    @staticmethod
    def get_coupon_pv(coupon, int_rate,period,freq):
       pv = coupon/(1+int_rate/freq)**period
       return pv
   
    # Defining a function to compute ytm of a bond
    def get_ytm(self,bond_price,face_value,coupon,years,freq = 1,estimate = 0.05):
       import scipy
       from scipy import optimize
       get_yield = lambda int_rate: self.get_price(coupon,face_value,int_rate,years,freq) - bond_price
       return optimize.newton(get_yield,estimate)

# Storing the class as a variable for feasiblity to use
coupon_bond_calc = Bond_Coupon()

""" Computing the ytm of the selected portfolio of 20 Bonds using a constructed for loop """
invest_amt = coupon_amt = time_coupons = 0
irr_list = []
for i in range(0,20):
    invest_amt = df_bonds_filtered.iloc[i,1]  
    coupon_amt = df_bonds_filtered.iloc[i,2]  
    time_coupons = df_bonds_filtered.iloc[i,3]
    redemption = 100
    irr_list.append(coupon_bond_calc.get_ytm(bond_price = invest_amt, face_value = redemption, coupon = coupon_amt, years = time_coupons))

# Concatenating the list of ytm within the bonds dataset
df_bonds_filtered["YTM"] = irr_list
df_bonds_filtered = df_bonds_filtered.reset_index()
df_bonds_filtered = df_bonds_filtered.drop(["index"],axis = 1)

# Concatenating the columns of 'trade_price_last1' and 'trade_size_last1'
# within the bonds dataset for respective bonds
df_bonds_all_cols = pd.merge(df_bonds_filtered,df_main[["bond_id","trade_price_last1","trade_size_last1"]],on = "bond_id",how = "left")

""" Computing the weighted average of ytm on the basis of last investments to
 arrive at the interest rate to be used to discount the liabilities """
df_bonds_all_cols["Investment_last_amt"] = df_bonds_all_cols["trade_price_last1"]*df_bonds_all_cols["trade_size_last1"]
portfolio_last_inst = sum(df_bonds_all_cols["Investment_last_amt"])
df_bonds_all_cols["Weight_last_invst"] = df_bonds_all_cols["Investment_last_amt"]/portfolio_last_inst
df_bonds_all_cols["Weight_YTM"] = df_bonds_all_cols["Weight_last_invst"] * df_bonds_all_cols["YTM"]
interest_rt_liab = sum(df_bonds_all_cols["Weight_YTM"])

# Loading the set of liabilities
df_liabilities = pd.read_excel(r"C:/Users/S145 DOIN/Downloads/Project Data.xlsx",sheet_name="Cash outflow (Liability)")

# Sorting liabilities with respect to the time of cashflows
df_liabilities.sort_values("Due Time(In Years)",inplace=True)

# Adding the liabilties occuring at the same time
df_liab_unique = df_liabilities.groupby(['Due Time(In Years)'],as_index=False).sum()

""" Computing the present value of liabilities """
pv_liab = 0
total_pv = 0
for z in range(0,8):
    pv_liab = (df_liab_unique.iloc[z,1])*(1+interest_rt_liab)**(-df_liab_unique.iloc[z,0])
    total_pv = total_pv + pv_liab

""" Computing the volatility (Modified duration) of liabilities """
vol_num1 = 0
vol_num = 0
for n in range(0,8):
    vol_num1 = df_liab_unique.iloc[z,0]*(df_liab_unique.iloc[z,1])*((1+interest_rt_liab)**(-(df_liab_unique.iloc[z,0]+1)))
    vol_num = vol_num + vol_num1
    
vol_final = vol_num/total_pv

""" Computing the convexity of liabilities """
convex_num1 = 0
convex_num = 0
for c in range(0,8):
    convex_num1 = (df_liab_unique.iloc[z,1])*(df_liab_unique.iloc[z,0])*(df_liab_unique.iloc[z,0]+1)*(1+interest_rt_liab)**(-(df_liab_unique.iloc[z,0]+2))
    convex_num = convex_num + convex_num1
convexity_final = convex_num/total_pv 

""" Computing volatility (Approximate Modified Duration) for all bonds """
df_bonds_vols = df_bonds_all_cols
approximate_modified_duration = []
for r in range(0,20):
    coupon_amt2 = df_bonds_vols.iloc[r,2]
    face_value2 = 100
    int_rate2 = df_bonds_vols.iloc[r,5]
    years2 = df_bonds_vols.iloc[r,3]
    trade_price2 = df_bonds_vols.iloc[r,1]
    pv_minus = coupon_bond_calc.get_price(coupon_amt2,face_value2,int_rate2 - 0.0001,years2)
    pv_plus = coupon_bond_calc.get_price(coupon_amt2,face_value2,int_rate2 + 0.0001,years2)
    app_mod_dur_num = pv_minus - pv_plus
    app_mod_dur_den = 2*(0.0001)*trade_price2
    app_mod_dur = app_mod_dur_num/app_mod_dur_den
    approximate_modified_duration.append(app_mod_dur)
df_bonds_vols["Volatility_Bonds"] = approximate_modified_duration

""" Computing approximate convexity of all bonds """
approximate_convexity = []
for r in range(0,20):
    coupon_amt3 = df_bonds_vols.iloc[r,2]
    face_value3 = 100
    int_rate3 = df_bonds_vols.iloc[r,5]
    years3 = df_bonds_vols.iloc[r,3]
    trade_price3 = df_bonds_vols.iloc[r,1]
    pv_minus2 = coupon_bond_calc.get_price(coupon_amt3,face_value3,int_rate3 - 0.0001,years3)
    pv_plus2 = coupon_bond_calc.get_price(coupon_amt3,face_value3,int_rate3 + 0.0001,years3)
    app_conv_num = pv_minus2 + pv_plus2 - 2*trade_price3
    app_conv_den = trade_price3*((0.0001)**2)
    app_conv = app_conv_num/app_conv_den
    approximate_convexity.append(app_conv)
df_bonds_vols["Convexity_Bonds"] = approximate_convexity

""" Using linprog function to determine the number of units required to be 
purchased to satisfy Redington Immunization Theory """
from scipy.optimize import linprog

# Setting the objective function to minimize the cost
obj = np.array(df_bonds_vols['trade_price'])

# =============================================================================
# # Setting the first two conditions of Redington Immunization Theory as equalities
# # for the purpose of linear programming.
# # Additionally, setting the third condition of Redington Immunization Theory as inequality
# # for the purpose of linear programming.
# =============================================================================

list_asset_vol = np.array(df_bonds_vols['trade_price'])*np.array(df_bonds_vols['Volatility_Bonds'])
list_asset_conv = -np.array(df_bonds_vols['trade_price'])*np.array(df_bonds_vols['Convexity_Bonds'])
lhs_eq = [obj,list_asset_vol]
rhs_eq = [total_pv,total_pv*vol_final]
lhs_ineq = [list_asset_conv]
rhs_ineq = [total_pv*convexity_final]
bnd = [(0, float("inf")) for q in range(len(obj))]
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")

# Checking whether optimization has terminated successfully.
opt.success  # The output is True which indicates that optimization has terminated successfully.

# Checking whether present value of assets is the same as present value of 
# liabilities
opt.fun # The present value of assets is equal to present value of liabilities

# Determining the number of units to be invested in the portfolio of bonds
opt.x

""" 
The output shows that 551117.33921959 units of Bond ID 25 and 
2687361.44841479 units of Bond ID 50 will have to be purchased in order to satisfy
the Redington Immunization Theory.
"""