from paper_functions import * ##import paper_functions as paper from numpy import * from scipy import stats as st from matplotlib.pyplot import * ioff() # to turn off interactive mode in spyder so that it can accurtely picture this import xlrd import pdb # to debug # put pdb.set_trace() in a place where you want to debug #seterr(divide='ignore', invalid='ignore') # this is for ignoring warning given by dividing by zero from datetime import datetime start_time = datetime.now() # here are the parameters workbook=xlrd.open_workbook('input.xlsm') sheet=workbook.sheet_by_name('Inputs') numScenarios = int( sheet.cell_value(0,1)) Months = int(sheet.cell_value(1,1)) dt = sheet.cell_value(2,1) iniRate = sheet.cell_value(3,1) iniTarget = sheet.cell_value(4,1) meanTarget = sheet.cell_value(5,1) targetVol = sheet.cell_value(6,1) adjVol = sheet.cell_value(7,1) Reversion = sheet.cell_value(8,1) alpha = sheet.cell_value(9,1) beta = sheet.cell_value(10,1) DiGamma_alpha = sheet.cell_value(21,1) DiGamma_beta = sheet.cell_value(22,1) tau = sheet.cell_value(23,1) # getting necessary matrix in order to make random elements non random # so that we can test if our code is correct # you don't really need this in order to run the program... You need thiis block only if you want to make random values non random ## another_sheet=workbook.sheet_by_name('another test') #declaring necessary numpy arrays # first random index starts from 0 to 9 # second random index starts from row of 0 to 9, colum of 0 to 720 ## first_random=zeros(10) ## second_random=zeros((10,731)) ## third_random=zeros((10,731)) ## fourth_random=zeros((10,731)) ## for i in range(10): ## first_random[i]=another_sheet.cell_value(i+1,0) ## ## for i in range(10): ## for j in range(731): ## second_random[i,j]=another_sheet.cell_value(i+1,j+2) ## third_random[i,j]=another_sheet.cell_value(i+13,j+2) ## fourth_random[i,j]=another_sheet.cell_value(i+25,j+2) targetRates=zeros(Months+1) #Calculating Drift Terms driftTerms=[0] for i in range(1,Months+1): driftTerms.append(Drift(Reversion, alpha, beta, dt, targetVol, adjVol ** 2, i * dt)) driftTerms=array(driftTerms) #plot(driftTerms) #show() #savetxt('drift_terms.txt', driftTerms, delimiter=',') # saving Drift Terms in a text file # Declaring numpy matrix for interest AllRate=zeros((numScenarios+1,Months+1)) intRates=zeros(Months+1) # this will change in each scenario #Declearing numpy matrix for volatility ( required for calculating spread. ) volatility_random_values=zeros(Months+1) volatility_matrix=zeros((numScenarios+1,Months+1)) long_rate_targets=zeros(Months+1) long_rate_targets_matrix=zeros((numScenarios+1,Months+1)) # Declaring numpy matrix and variables for calculating statistics ScenarioStat=zeros((10,numScenarios+1)) LnDeltaIntRates=zeros(Months+1) # ===========================BEGIN OF BIG LOOP================================ a = math.log(meanTarget) - ( targetVol ** 2 ) / 2 currentTarget = iniTarget # Holds the current target rate for the scenario previousTimeChange = 0 # Holds the previous time change for the weighted average calculation for the equation after the do loop nextTimeChange = 0 nextTarget = iniTarget ##"""I might not need the following three lines either ##Initializations for worst scenario calculation ##""" for i in range(1,numScenarios+1): #I might not need the following three lines either #Initializations for worst scenario calculation CumulativeBelow = 0 CumulativeAll = 0 WorstGroup = 0 #Set the current target back to the initial target at the beginning of each scenario currentTarget = iniTarget nextTarget = iniTarget ##' Calculate first target change for the scenario # the first random variable nextTimeChange = RandomdBar(alpha, beta) / dt # nextTimeChange = first_random[i-1] previousTimeChange = nextTimeChange #===========================BEGIN OF SMALL LOOP================================ for j in range(1,Months+1): ## 'If the target changed last period, begin calculating the weighted average ## 'Otherwise, skip this. ## 'Note that j=1 corresponds to time = 0 ## ' that means If t > nextTimeChange ## ' that mean t > t_j+1 ## ' that means If theCurrentTime > nextTimeChange ## ' Then, it is time to change the target ## ' The line below is the first element from the equation that calculates T_t ( Target at time t) ## ' So, it is obvious that the first element in calculating T_t is sure to exist ## ' and currentTime - nextTimeChange is less than dt ## ' while currentTime-previousTimeChange is greater than dt ## ' So, the second, third and all the other terms that came later are supposed to be 1 most of the time ## ' The only time when the second and third terms exist is when we enter this While Loop ## ## if j > 1 + nextTimeChange: currentTarget = currentTarget ** ((((1 - Reversion) **((j - 1 - nextTimeChange) * dt) - ((1 - Reversion) ** dt)) / (1 - (1 - Reversion) ** dt))) previousTimeChange = nextTimeChange # second random variable second_random_variable = st.gamma.ppf(random.uniform(0,1), alpha,0, beta) / dt # second_random_variable = second_random[i-1,j-1] nextTimeChange = nextTimeChange + second_random_variable ## 'Enter loop if the target changed more than once last period while j > 1 + nextTimeChange: J=j print 'go into evil loop' pdb.set_trace() nextRate = exp(st.norm.ppf(random.uniform(0,1), a, targetVol)) currentTarget = currentTarget * nextRate ** ((((1 - Reversion) ** ((j - 1 - nextTimeChange) * dt) - (1 - Reversion) ** ((j - 1 - previousTimeChange) * dt)) / (1 - (1 - Reversion) ** dt))) previousTimeChange = nextTimeChange nextTimeChange = nextTimeChange + st.gamma.ppf(random.uniform(0,1),alpha,0,beta) / dt # third random variable nextTarget = exp(st.norm.ppf(random.uniform(0,1), a, targetVol)) # nextTarget = third_random[i-1,j-1] currentTarget = currentTarget * nextTarget ** (((1 - (1 - Reversion)** ((j - 1 - previousTimeChange) * dt)) / (1 - (1 - Reversion) ** dt))) ## 'Set the current period's target to the weighted average targetRates[j - 1] = currentTarget ## ' The following two lines prepare for the next loop ( or the next month in loop ) currentTarget = nextTarget targetRates[j] = currentTarget ## 'Calculate interest rates ## 'note that j=1 corresponds to time 0 if j == 1: # Here is the initial month. We only need to initialize runningTotal and intRates() intRates[j] = iniRate currentRate = iniRate else: # Here is the following interest rates calculation ## 'Calculate current month's interest rate Target = targetRates[j - 1] previousRate = intRates[j - 1] # fourth random variable fourth_random_variable=DiGammaInv(random.uniform(0,1), DiGamma_alpha, DiGamma_beta, tau) # fourth_random_variable = fourth_random[i-1,j-1] Dln_rt = (1 - (1 - Reversion) ** dt) * (math.log(Target) - math.log(previousRate)) + (1 - Reversion) ** dt * driftTerms[j - 1] * dt + (1 - Reversion) ** dt * fourth_random_variable currentRate = previousRate * exp(Dln_rt) intRates[j] = currentRate LnDeltaIntRates[j] = log(intRates[j] / intRates[j - 1]) volatility_random_values[j]=fourth_random_variable long_rate_targets[j]=Target #===========================END OF SMALL LOOP================================ AllRate[i,j] = intRates[j] # j is months and i is scenarios volatility_matrix[i,j] = volatility_random_values[j] long_rate_targets_matrix[i,j] = long_rate_targets[j] #========================== Calculating Statistics ========================= ## if j != 1: ## Appaverage = average(intRates[1:]) ## AppStDev_P = std(intRates[1:]) ## DLnAve = average(LnDeltaIntRates[1:]) ## DLnStDev_P = std(LnDeltaIntRates[1:]) ## ## third = 0 ## fourth = 0 ## sixth = 0 ## dlnthird = 0 ## DLnFourth = 0 ## DlnSixth = 0 ## ## for s in range(1, Months+1): ## third = third + ( intRates[s] - Appaverage ) ** 3 ## fourth = fourth + ( intRates[s] - Appaverage ) ** 4 ## sixth = sixth + ( intRates[s] - Appaverage ) ** 6 ## dlnthird = dlnthird + ( LnDeltaIntRates[s] - DLnAve ) ** 3 ## DLnFourth = DLnFourth + ( LnDeltaIntRates[s] - DLnAve ) ** 4 ## DlnSixth = DlnSixth + ( LnDeltaIntRates[s] - DLnAve ) ** 6 ## ## #print 'it is good '+str(i) ## #Store ScenarioStats ## ScenarioStat[1][i] = Appaverage ## ScenarioStat[2][i] = AppStDev_P ## ScenarioStat[3][i] = third / Months / AppStDev_P ** 3 ## ScenarioStat[4][i] = fourth / Months / AppStDev_P ** 4 ## ScenarioStat[5][i] = sixth / Months / AppStDev_P ** 6 ## ScenarioStat[6][i] = DLnStDev_P ## ScenarioStat[7][i] = dlnthird / Months / DLnStDev_P ** 3 ## ScenarioStat[8][i] = DLnFourth / Months / DLnStDev_P ** 4 ## ScenarioStat[9][i] = DlnSixth / Months / DLnStDev_P ** 6 ## #===========================END OF BIG LOOP================================ # Reformmating the AllRate (interest rate array) # Since we calculated using array index of 1, index 0 has no value in both rows and colums # so, I am now deleting the first row and the first column AllRate=AllRate[1:,1:] volatility_matrix=volatility_matrix[1:,1:] ScenarioStat=ScenarioStat[1:,1:] long_rate_targets_matrix=long_rate_targets_matrix[1:,1:] #############final_statistics=ScenarioStats(numScenarios,ScenarioStat) #==================================== Generating the Yield Curve ===========================================# # AAA values are ##a 0.02685 ##b 0.0002 ##c 0.0414 a=sheet.cell_value(32,1) b=sheet.cell_value(33,1) c=sheet.cell_value(34,1) theta=sheet.cell_value(35,1) rho=-0.19197 # rho is the correlation coefficient between 1Zt and 2Zt, which is used in calculating spread initial_spread = sheet.cell_value(37,1) # The difference between long rate and short rate at April 1953 target_spread_at_t= sheet.cell_value(38,1) # The average of historical spreads until Feb 2014 ( this should probably not be constant ) Long_rates_matrix=AllRate Lower_bound_for_short_rate=sheet.cell_value(39,1) # this lower bound and kappa are here to make sure that short rate doesn't get ridiculous values kappa =sheet.cell_value(40,1) # since spread could be negative spread_matrix=generating_spread(Long_rates_matrix,a,b,c,target_spread_at_t, initial_spread,long_rate_targets_matrix,theta,volatility_matrix ,DiGamma_alpha,DiGamma_beta,tau,rho) short_rates_matrix=short_rate_generator(Long_rates_matrix,spread_matrix,kappa,Lower_bound_for_short_rate) #Plotting ##random_scenario=int(round(random.uniform(0,numScenarios-1))) ##m=range(1,Months+1) ##plot(m,Long_rates_matrix[random_scenario,:],label='Long Rate') ##plot(m,short_rates_matrix[random_scenario,:],label='Short Rate') ##plot(m,spread_matrix[random_scenario,:],label='Spread') ##xlabel('Months') ##ylabel('Interest Rates') ##title(' Scenario '+ str(random_scenario)) ##legend() ##show() #==================================== Saving Excel Files ===========================================# # Saving all the interest rates savetxt("long_rate.csv", AllRate, delimiter=",") savetxt("short_rate.csv", short_rates_matrix, delimiter=",") savetxt("spread_rate.csv",spread_matrix, delimiter=",") ##savetxt("scenario_stat.csv", ScenarioStat, delimiter=",") #==================================== Plotting ===========================================# if numScenarios<=50: photos_to_print=numScenarios else: photos_to_print=50 for i in range(0,photos_to_print): m=range(1,Months+1) plot(m,Long_rates_matrix[i,:],label='Long Rate') plot(m,short_rates_matrix[i,:],label='Short Rate') plot(m,spread_matrix[i,:],label='Spread') xlabel('Months') ylabel('Interest Rates') title(' Scenario '+ str(i+1)) legend() savefig( ' Scenario '+ str(i+1)+'.png') clf() # calculting the program run time end_time = datetime.now() print('Duration: {}'.format(end_time - start_time))