# -*- coding: utf-8 -*-
"""
Created on Sat May 14 12:14:34 2022

@author: Sayali
"""

""" Libraries to import"""

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

"""Introducing the bonds dataset"""
bond= pd.read_csv(r"C:/Users/Sayali/Downloads/bond.csv")

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

"""  Picking a portfolio of bonds from a list of 50"""

#Bonds that mature over long periods of time have a longer duration than bonds 
#that pay off within a shorter time period. In a portfolio of 20 bonds,
# we consider only bonds that mature over one year.

bond_mm_1 = df_bonds[df_bonds["time_to_maturity"] > 1]

#Due to the fact that our last liability occurs at time 20,
# we are searching for a bond with a maturity of greater than 1 and less than 21.

df_bonds_new = bond_mm_1[(bond_mm_1["time_to_maturity"] > 1.00) & (bond_mm_1["time_to_maturity"] < 21.00)]
portfolio_20 = 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])]

""" The YTM is calculated for all 20 bonds in the portfolio"""

class portfolio:
    #  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
    
    #  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
   
    # 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
   
    #  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
bond_mat = portfolio()

""" By using a constructed for loop, we calculate the ytm for a portfolio of 20 Bonds"""
amount_invested = 0
coupon_amt = 0
time_coupons = 0
ytm = []
for i in range(0,20):
    amount_invested = portfolio_20.iloc[i,1]  
    coupon_amt = portfolio_20.iloc[i,2]  
    time_coupons = portfolio_20.iloc[i,3]
    redemption = 100
    ytm.append(bond_mat.get_ytm(bond_price = amount_invested, face_value = redemption, coupon = coupon_amt, years = time_coupons))
# adding column 
portfolio_20["YTM"] = ytm
# resetting index
portfolio_20 = portfolio_20.reset_index()
portfolio_20 = portfolio_20.drop(["index"],axis = 1)
# merging  data set 
final_portfolio = pd.merge(portfolio_20,bond[["bond_id","trade_price_last1","trade_size_last1"]],on = "bond_id",how = "left")



""" Utilizing the last investments as the basis for calculating the weighted average 
of YTM to Determine the interest rate that will be used to discount the liabilities
 """

final_portfolio["last_price"] = final_portfolio["trade_price_last1"]*final_portfolio["trade_size_last1"]
portfolio_last_inst = sum(final_portfolio["last_price"])
final_portfolio["Weight_last_invst"] = final_portfolio["last_price"]/portfolio_last_inst
final_portfolio["Weight_YTM"] = final_portfolio["Weight_last_invst"] * final_portfolio["YTM"]
interest_rt_liab = sum(final_portfolio["Weight_YTM"])

# The set of liabilities is loaded
df_liabilities = pd.read_csv(r"C:/Users/Sayali/Downloads/liabilities.csv")
# Sorting liabilities with respect to the time of cashflows
df_liabilities.sort_values("time",inplace=True)

# Liabilities occurring at the same time should be added together
liability_at_time = df_liabilities.groupby(['time'],as_index=False).sum()


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

""" volatility of liabilities (modified duration) """
volatility  = 0
volatility__ = 0
for n in range(0,8):
    volatility  = liability_at_time.iloc[z,0]*(liability_at_time.iloc[z,1])*((1+interest_rt_liab)**(-(liability_at_time.iloc[z,0]+1)))
    volatility__ = volatility__ + volatility 
    
vol_final = volatility__/total_pv

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


""" Approximate Modified Duration for all bonds """
df_bonds_vols = final_portfolio
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 = bond_mat.get_price(coupon_amt2,face_value2,int_rate2 - 0.0001,years2)
    pv_plus = bond_mat.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


""" 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 = bond_mat.get_price(coupon_amt3,face_value3,int_rate3 - 0.0001,years3)
    pv_plus2 = bond_mat.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

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


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 present value of assets is the same as present value of liabilities

opt.fun 

# hence The present value of assets is equal to present value of liabilities,
# hence it satisfies Redington Immunization conditions