# -*- coding: utf-8 -*-
"""
Created on Sat May 14 10:02:49 2022

@author: kuku
"""


import pandas as pd
import numpy as np

from scipy.optimize import linprog


""" Importing the dataset """
df_excel = pd.read_excel(r"C:/Users/kuku/Downloads/Project Data.xlsx",sheet_name="Bonds")

""" Deleting unnecessary columns """
df_bond = df_excel.drop(["weight","reporting_delay","trade_size","trade_type","trade_price_last1","trade_size_last1"],axis = 1)

""" Considering the time to maturity """

df_bonds_tm = df_bond[df_bond["time_to_maturity"] > 1]

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

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


""" Calculating the YTM for the bonds """

class Coupon_bond:

    # for 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
    
  # for pv of face value
    @staticmethod
    def get_face_value_pv(face_value,int_rate,years):
       fvpv = face_value/(1+int_rate)**years
       return fvpv
   
    # for 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
   
    # for 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)

bond_coupon_calc = Coupon_bond()

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

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


df_bonds_all_colmns = pd.merge(df_filtered_bonds,df_excel[["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_colmns["Investment_last_amt"] = df_bonds_all_colmns["trade_price_last1"]*df_bonds_all_colmns["trade_size_last1"]
last_invt = sum(df_bonds_all_colmns["Investment_last_amt"])
df_bonds_all_colmns["Weight_last_invst"] = df_bonds_all_colmns["Investment_last_amt"]/last_invt
df_bonds_all_colmns["Weight_YTM"] = df_bonds_all_colmns["Weight_last_invst"] * df_bonds_all_colmns["YTM"]
interest_rt_liab = sum(df_bonds_all_colmns["Weight_YTM"])

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

df_liabilities.sort_values("Due Time(In Years)",inplace=True)

df_liability = 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_liability.iloc[z,1])*(1+interest_rt_liab)**(-df_liability.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_liability.iloc[z,0]*(df_liability.iloc[z,1])*((1+interest_rt_liab)**(-(df_liability.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_liability.iloc[z,1])*(df_liability.iloc[z,0])*(df_liability.iloc[z,0]+1)*(1+interest_rt_liab)**(-(df_liability.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_volume = df_bonds_all_colmns
approx_modified_duration = []
for r in range(0,20):
    coupon_amt2 = df_volume.iloc[r,2]
    face_value2 = 100
    int_rate2 = df_volume.iloc[r,5]
    years2 = df_volume.iloc[r,3]
    trade_price2 = df_volume.iloc[r,1]
    pv_minus = bond_coupon_calc.get_price(coupon_amt2,face_value2,int_rate2 - 0.0001,years2)
    pv_plus = bond_coupon_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
    approx_modified_duration.append(app_mod_dur)
df_volume["Volatility_Bonds"] = approx_modified_duration

""" Computing approximate convexity of all bonds """

approx_convexity = []
for r in range(0,20):
    coupon_amt3 = df_volume.iloc[r,2]
    face_value3 = 100
    int_rate3 = df_volume.iloc[r,5]
    years3 = df_volume.iloc[r,3]
    trade_price3 = df_volume.iloc[r,1]
    pv_minus2 = bond_coupon_calc.get_price(coupon_amt3,face_value3,int_rate3 - 0.0001,years3)
    pv_plus2 = bond_coupon_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
    approx_convexity.append(app_conv)
df_volume["Convexity_Bonds"] = approx_convexity

""" Using linprog function to determine the number of units required to be 
purchased to satisfy Redington Immunization Theory """

# Setting the objective function to minimize the cost
obj = np.array(df_volume['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_volume['trade_price'])*np.array(df_volume['Volatility_Bonds'])
list_asset_conv = -np.array(df_volume['trade_price'])*np.array(df_volume['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

## 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.

"""