#import libraries
import math
import pandas as pd
import numpy as np
from scipy.optimize import fsolve
#define index change rate
u = 1.128208823
d = 1/u
#define current interest rate
r_0 = 0.0006
#define bond prices and volatilities
price_2 = 98.96589012
vol_2 = 0.013504395
price_3 = 97.4512835
vol_3 = 0.026782142
price_4 = 95.92012621
vol_4 = 0.038335293
price_5 = 94.06593247
vol_5 = 0.051799679
Period 1
def period_1 (r_1d):
#risk neutral probability
pr_0_1 = (math.exp(r_0) - d) / (u - d)
#interest-rates
r_1u = r_1d*math.exp(2*vol_2)
#From period 2 to 1 valuation
val_2_1 = pr_0_1 * (100 / (1+r_1u)) + (1 - pr_0_1) * (100 / (1 + r_1d))
#From period 0 to 1 valuation
val_0_1 = price_2*(1 + r_0)
return val_2_1 - val_0_1
#solve the equation
r_1d = fsolve(period_1, [1])[0]
r_1u = r_1d * math.exp(2*vol_2)
print(r_1d*100,r_1u*100)
#probability for the next period
pr_0_1 = (math.exp(r_0) - d) / (u - d)
Period 2
def period_2 (r_2d):
#risk neutral probabilities for two nodes
pr_1d_2 = (math.exp(r_1d) - d) / (u - d)
pr_1u_2 = (math.exp(r_1u) - d) / (u - d)
#interest-rates
r_1u1d = r_2d*math.exp(2*vol_3)
r_2u = r_2d*math.exp(4*vol_3)
#From period 3 to 2 valuation
val_r_1d = (pr_1d_2 * (100 / (1+r_1u1d)) + (1 - pr_1d_2) * (100 / (1+r_2d)))
val_r_1u = (pr_1u_2 * (100 / (1+r_2u)) + (1 - pr_1u_2) * (100 / (1+r_1u1d)))
#From period 2 to 1 valuation
val_2_1 = pr_0_1 * (val_r_1u / (1+r_1u)) + (1-pr_0_1) * (val_r_1d / (1+r_1d))
#From period 0 to 1 valuation
val_0_1 = price_3*(1 + r_0)
return val_2_1 - val_0_1
#solve the equation
r_2d = fsolve(period_2, [1])[0]
r_1u1d = r_2d * math.exp(2*vol_3)
r_2u = r_2d * math.exp(4*vol_3)
print(r_2d*100,r_1u1d*100,r_2u*100)
#probabilities for next period
pr_1d_2 = (math.exp(r_1d) - d) / (u - d)
pr_1u_2 = (math.exp(r_1u) - d) / (u - d)
Period 3
def period_3 (r_3d):
#risk neutral probabilities for three nodes
pr_2d_3 = (math.exp(r_2d) - d) / (u - d)
pr_1d1u_3 = (math.exp(r_1u1d) - d) / (u - d)
pr_2u_3 = (math.exp(r_2u) - d) / (u - d)
#interest-rates
r_2d1u = r_3d*math.exp(2*vol_4)
r_1d2u = r_3d*math.exp(4*vol_4)
r_3u = r_3d*math.exp(6*vol_4)
#From period 4 to 3 valuation
val_r_2d = (pr_2d_3 * (100 / (1+r_2d1u)) + (1 - pr_2d_3) * (100 / (1+r_3d)))
val_r_1d1u = (pr_1d1u_3 * (100 / (1+r_1d2u)) + (1 - pr_1d1u_3) * (100 / (1+r_2d1u)))
val_r_2u = (pr_2u_3 * (100 / (1+r_3u)) + (1 - pr_2u_3) * (100 / (1+r_1d2u)))
#From period 3 to 2 valuation
val_r_1d = pr_1d_2 * (val_r_1d1u / (1 + r_1u1d)) + (1 - pr_1d_2) * (val_r_2d / (1 + r_2d))
val_r_1u = pr_1u_2 * (val_r_2u / (1 + r_2u)) + (1 - pr_1u_2) * (val_r_1d1u / (1 + r_1u1d))
#From period 2 to 1 valuation
val_2_1 = pr_0_1 * (val_r_1u / (1+r_1u)) + (1-pr_0_1) * (val_r_1d / (1+r_1d))
#From period 0 to 1 valuation
val_0_1 = price_4 * (1 + r_0)
return val_2_1 - val_0_1
#solve the equation
r_3d = fsolve(period_3, [1])[0]
r_2d1u = r_3d * math.exp(2*vol_4)
r_1d2u = r_3d * math.exp(4*vol_4)
r_3u = r_3d * math.exp(6*vol_4)
print(r_3d*100, r_2d1u*100, r_1d2u*100, r_3u*100)
#probabilities for next period
pr_2d_3 = (math.exp(r_2d) - d) / (u - d)
pr_1d1u_3 = (math.exp(r_1u1d) - d) / (u - d)
pr_2u_3 = (math.exp(r_2u) - d) / (u - d)
Period 4
def period_4 (r_4d):
#risk neutral probabilities for four nodes
pr_3d_4 = (math.exp(r_3d) - d) / (u - d)
pr_2d1u_4 = (math.exp(r_2d1u) - d) / (u - d)
pr_1d2u_4 = (math.exp(r_1d2u) - d) / (u - d)
pr_3u_4 = (math.exp(r_3u) - d) / (u - d)
#interest-rates
r_3d1u = r_4d*math.exp(2*vol_5)
r_2d2u = r_4d*math.exp(4*vol_5)
r_1d3u = r_4d*math.exp(6*vol_5)
r_4u = r_4d*math.exp(8*vol_5)
#From period 5 to 4 valuation
val_r_3d = (pr_3d_4 * (100 / (1+r_3d1u)) + (1 - pr_3d_4) * (100 / (1+r_4d)))
val_r_2d1u = (pr_2d1u_4 * (100 / (1+r_2d2u)) + (1 - pr_2d1u_4) * (100 / (1+r_3d1u)))
val_r_1d2u = (pr_1d2u_4 * (100 / (1+r_1d3u)) + (1 - pr_1d2u_4) * (100 / (1+r_2d2u)))
val_r_3u = (pr_3u_4 * (100 / (1+r_4u)) + (1 - pr_3u_4) * (100 / (1+r_1d3u)))
#From period 4 to 3 valuation
val_r_2d = (pr_2d_3 * (val_r_2d1u / (1+r_2d1u)) + (1 - pr_2d_3) * (val_r_3d / (1+r_3d)))
val_r_1d1u = (pr_1d1u_3 * (val_r_1d2u / (1+r_1d2u)) + (1 - pr_1d1u_3) * (val_r_2d1u / (1+r_2d1u)))
val_r_2u = (pr_2u_3 * (val_r_3u / (1+r_3u)) + (1 - pr_2u_3) * (val_r_1d2u / (1+r_1d2u)))
#From period 3 to 2 valuation
val_r_1d = pr_1d_2 * (val_r_1d1u / (1 + r_1u1d)) + (1 - pr_1d_2) * (val_r_2d / (1 + r_2d))
val_r_1u = pr_1u_2 * (val_r_2u / (1 + r_2u)) + (1 - pr_1u_2) * (val_r_1d1u / (1 + r_1u1d))
#From period 2 to 1 valuation
val_2_1 = pr_0_1 * (val_r_1u / (1+r_1u)) + (1-pr_0_1) * (val_r_1d / (1+r_1d))
#From period 0 to 1 valuation
val_0_1 = price_5 * (1 + r_0)
return val_2_1 - val_0_1
#solve the equation
r_4d = fsolve(period_4, [1])[0]
r_3d1u = r_4d * math.exp(2*vol_5)
r_2d2u = r_4d * math.exp(4*vol_5)
r_1d3u = r_4d * math.exp(6*vol_5)
r_4u = r_4d * math.exp(8*vol_5)
print(r_4d*100, r_3d1u*100, r_2d2u*100, r_1d3u*100, r_4u*100)
#probabilities for next period
pr_3d_4 = (math.exp(r_3d) - d) / (u - d)
pr_2d1u_4 = (math.exp(r_2d1u) - d) / (u - d)
pr_1d2u_4 = (math.exp(r_1d2u) - d) / (u - d)
pr_3u_4 = (math.exp(r_3u) - d) / (u - d)
#create dictionary for interest rates
result = [{'Period' : 0, '0':r_0},
{'Period' : 1, '0':r_1d, '1':r_1u},
{'Period' : 2, '0':r_2d, '1':r_1u1d, '2':r_2u},
{'Period' : 3, '0':r_3d, '1':r_2d1u, '2':r_1d2u, '3':r_3u},
{'Period' : 4, '0':r_4d, '1':r_3d1u, '2':r_2d2u, '3':r_1d3u, '4': r_4u}]
#create dataframe for interest rates
interest_rates = pd.DataFrame(result)
interest_rates.set_index('Period', inplace= True)
interest_rates.fillna('-', inplace=True)
#create dictionary for risk neutral probabilities
probs = [ {'Period' : 0, '0':pr_0_1},
{'Period' : 1, '0':pr_1d_2, '1':pr_1u_2},
{'Period' : 2, '0':pr_2d_3, '1':pr_1d1u_3, '2':pr_2u_3},
{'Period' : 3, '0':pr_3d_4, '1':pr_2d1u_4, '2':pr_1d2u_4, '3':pr_3u_4}]
#create dataframe for risk neutral probabilities
risk_neutral = pd.DataFrame(probs)
risk_neutral.set_index('Period', inplace= True)
risk_neutral.fillna('-', inplace=True)
#export dataframes to excel
risk_neutral.to_excel('risk_neutral_probabilities.xlsx')
interest_rates.to_excel('interest_rates.xlsx')