From e6538fee8b5f3c771f50ce16ea65bdf24ba9ce0e Mon Sep 17 00:00:00 2001 From: raja Date: Mon, 4 May 2020 21:17:16 +1000 Subject: [PATCH] Partial update to use the code in python 3 --- .gitignore | 6 + AAA_RUN_KANSM2_Est_LB.py | 327 ++++++++++++++++++----------------- AAC_KAGM_SingleLoop_file.py | 25 ++- AAE_KAGM_f_and_df_dx_file.py | 2 +- data_read.py | 46 ++--- functions.py | 6 +- 6 files changed, 213 insertions(+), 199 deletions(-) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..d397acc --- /dev/null +++ b/.gitignore @@ -0,0 +1,6 @@ +.idea/** +__pycache__ +venv/ +*.csv +*.pdf +*.jpg \ No newline at end of file diff --git a/AAA_RUN_KANSM2_Est_LB.py b/AAA_RUN_KANSM2_Est_LB.py index 1e58607..8fb33a0 100644 --- a/AAA_RUN_KANSM2_Est_LB.py +++ b/AAA_RUN_KANSM2_Est_LB.py @@ -14,82 +14,81 @@ import csv # Hyperparameters. -N = 2 # FIXED: Number of factors. -dTau = 0.01 # OPTION: Spacing for TauGrid, used to numerically obtain R. -KappaP_Constraint = 'Direct' # FIXED: KappaP matrix values are set directly, (but subject to an eigenvalue constraint in 'AAC_EKF_CAB_GATSM_SingleLoop'). -ZLB_Imposed = 1 # FIXED: 0=ANSM(2) or 1=K-ANSM(2). -DailyIterations = 200 # OPTION: Sets number of iterations between interim saves. -IEKF_Count = -1e-5 # OPTION: EKF if 0, IEKF steps if >0, tolerance if <0 (e.g. -1e-5). -FinalNaturalParametersGiven = 1 # OPTION: Full estimation if 0 (the Optimization toolbox is required), partial estimation with given parameters if 0. -HessianRequired = 0 # OPTION: Omits Hessian and standard errors if 0, calculates them if 1. +N = 2 # FIXED: Number of factors. +dTau = 0.01 # OPTION: Spacing for TauGrid, used to numerically obtain R. +KappaP_Constraint = 'Direct' # FIXED: KappaP matrix values are set directly, (but subject to an eigenvalue constraint in 'AAC_EKF_CAB_GATSM_SingleLoop'). +ZLB_Imposed = 1 # FIXED: 0=ANSM(2) or 1=K-ANSM(2). +DailyIterations = 200 # OPTION: Sets number of iterations between interim saves. +IEKF_Count = -1e-5 # OPTION: EKF if 0, IEKF steps if >0, tolerance if <0 (e.g. -1e-5). +FinalNaturalParametersGiven = 1 # OPTION: Full estimation if 0 (the Optimization toolbox is required), partial estimation with given parameters if 0. +HessianRequired = 0 # OPTION: Omits Hessian and standard errors if 0, calculates them if 1. # GSW US data file. Country = 'UK' DataFrequency = 'Monthly' -DataFileName=Country+'_GSW_Govt' +DataFileName = Country + '_GSW_Govt' -#load([DataFileName,'.mat']) +# load([DataFileName,'.mat']) # DailyDateIndex DailyYieldCurveData Maturities MonthlyDateIndex MonthlyYieldCurveData WeeklyDateIndex WeeklyYieldCurveData -Maturities=numpy.array([0.25,0.5,1,2,3,4,5,7,10,15,20,30]) -sub_datenum = datetime.date(1899,12,30).toordinal() + 366 -daily_file_name = Country+'_Daily.csv' +Maturities = numpy.array([0.25, 0.5, 1, 2, 3, 4, 5, 7, 10, 15, 20, 30]) +sub_datenum = datetime.date(1899, 12, 30).toordinal() + 366 +daily_file_name = Country + '_Daily.csv' daily_data = numpy.genfromtxt(daily_file_name, delimiter=',') DailyDateIndex = daily_data[:, 0].copy().astype('int') -DailyDateIndex=DailyDateIndex.__add__(sub_datenum) -DailyYieldCurveData = daily_data[:, 1:].copy() +DailyDateIndex = DailyDateIndex.__add__(sub_datenum) +DailyYieldCurveData = daily_data[:, 1:].copy() -weekly_file_name = Country+'_Weekly.csv' +weekly_file_name = Country + '_Weekly.csv' weekly_data = numpy.genfromtxt(weekly_file_name, delimiter=',') WeeklyDateIndex = weekly_data[:, 0].copy().astype('int') -WeeklyDateIndex=WeeklyDateIndex.__add__(sub_datenum) -WeeklyYieldCurveData = weekly_data[:, 1:].copy() +WeeklyDateIndex = WeeklyDateIndex.__add__(sub_datenum) +WeeklyYieldCurveData = weekly_data[:, 1:].copy() -month_file_name = Country+'_Monthly.csv' +month_file_name = Country + '_Monthly.csv' month_data = numpy.genfromtxt(month_file_name, delimiter=',') MonthlyDateIndex = month_data[:, 0].copy().astype('int') -MonthlyDateIndex= MonthlyDateIndex.__add__(sub_datenum) -MonthlyYieldCurveData = month_data[:, 1:].copy() - -#print numpy.shape(DailyDateIndex) -#print numpy.shape(DailyYieldCurveData) -#print numpy.shape(Maturities) -#print numpy.shape(MonthlyDateIndex) -#print numpy.shape(MonthlyYieldCurveData) -#print numpy.shape(WeeklyDateIndex) -#print numpy.shape(WeeklyYieldCurveData) +MonthlyDateIndex = MonthlyDateIndex.__add__(sub_datenum) +MonthlyYieldCurveData = month_data[:, 1:].copy() +# print(numpy.shape(DailyDateIndex)) +# print(numpy.shape(DailyYieldCurveData)) +# print(numpy.shape(Maturities)) +# print(numpy.shape(MonthlyDateIndex)) +# print(numpy.shape(MonthlyYieldCurveData)) +# print(numpy.shape(WeeklyDateIndex)) +# print(numpy.shape(WeeklyYieldCurveData)) FirstDay = DailyDateIndex[0] LastDay = DailyDateIndex[-1] -#FirstDay = datetime_to_matlab_datenum(datetime.datetime.strptime('30-Dec-1994', "%d-%b-%Y")) # Start of 30-year data. -#LastDay = datetime_to_matlab_datenum(datetime.datetime.strptime('31-Jul-2013', "%d-%b-%Y")) +# FirstDay = datetime_to_matlab_datenum(datetime.datetime.strptime('30-Dec-1994', "%d-%b-%Y")) # Start of 30-year data. +# LastDay = datetime_to_matlab_datenum(datetime.datetime.strptime('31-Jul-2013', "%d-%b-%Y")) SampleMaturities = numpy.array([0.25, 0.5, 1, 2, 3, 5, 7, 10, 30]) # Set starting parameters. # These lines set the estimated values from Krippner (2015). -#LoadName='BOOK_US_GSW_Govt_OIS_rL_30_K_AFNSM2_Monthly_IEKF_E-5_Final_2014_8_31_10_52_X'; -#load([LoadName,'.mat'],'FinalNaturalParameters'); +# LoadName='BOOK_US_GSW_Govt_OIS_rL_30_K_AFNSM2_Monthly_IEKF_E-5_Final_2014_8_31_10_52_X'; +# load([LoadName,'.mat'],'FinalNaturalParameters'); # FinalNaturalParameters -FinalNaturalParameters_country="FinalNaturalParameters_"+Country+".dat" +FinalNaturalParameters_country = "FinalNaturalParameters_" + Country + ".dat" FinalNaturalParameters = numpy.loadtxt(FinalNaturalParameters_country) InitialNaturalParameters = numpy.copy(FinalNaturalParameters) # Select the required data from the yield curve data file. IncludeMaturities = numpy.array([x in SampleMaturities for x in Maturities]) -if (DataFrequency == 'Daily'): - (StartT,) = numpy.where(DailyDateIndex==FirstDay) - (EndT,) = numpy.where(DailyDateIndex==LastDay) - YieldCurveDateIndex = DailyDateIndex[StartT:EndT+1] - YieldCurveData = DailyYieldCurveData[StartT:EndT+1,IncludeMaturities] +if DataFrequency == 'Daily': + (StartT,) = numpy.where(DailyDateIndex == FirstDay) + (EndT,) = numpy.where(DailyDateIndex == LastDay) + YieldCurveDateIndex = DailyDateIndex[StartT:EndT + 1] + YieldCurveData = DailyYieldCurveData[StartT:EndT + 1, IncludeMaturities] Dt = (YieldCurveDateIndex[-1] - YieldCurveDateIndex[0] + 1) / (len(YieldCurveDateIndex) * 365.25) Iterations = numpy.copy(DailyIterations) -elif (DataFrequency == 'Weekly'): +elif DataFrequency == 'Weekly': # Find earliest Friday consistent with FirstDay. ReferenceFriday = datetime_to_matlab_datenum(datetime.datetime.strptime('24-May-2013', "%d-%b-%Y")) WeeksToStepBackForStart = numpy.floor((ReferenceFriday - FirstDay) / 7) @@ -97,37 +96,38 @@ # Find latest Friday consistent with LastDay. WeeksToStepBackForEnd = 1 + numpy.floor((ReferenceFriday - LastDay) / 7) LastWeek = ReferenceFriday - 7 * WeeksToStepBackForEnd - (StartT,) = numpy.where(WeeklyDateIndex==FirstWeek) - (EndT,) = numpy.where(WeeklyDateIndex==LastWeek) - YieldCurveDateIndex = WeeklyDateIndex[StartT:EndT+1] - YieldCurveData = WeeklyYieldCurveData[StartT:EndT+1,IncludeMaturities] + (StartT,) = numpy.where(WeeklyDateIndex == FirstWeek) + (EndT,) = numpy.where(WeeklyDateIndex == LastWeek) + YieldCurveDateIndex = WeeklyDateIndex[StartT:EndT + 1] + YieldCurveData = WeeklyYieldCurveData[StartT:EndT + 1, IncludeMaturities] Dt = 7 / 365.25 Iterations = DailyIterations * 5 else: # Find earliest end-month consistent with FirstDay. - [FirstYear, FirstMonth] = [datetime.datetime.fromordinal(numpy.int(FirstDay-366)).year, datetime.datetime.fromordinal(numpy.int(FirstDay-366)).month] + [FirstYear, FirstMonth] = [datetime.datetime.fromordinal(numpy.int(FirstDay - 366)).year, + datetime.datetime.fromordinal(numpy.int(FirstDay - 366)).month] if FirstMonth == 12: - FirstYear = FirstYear+1 - FirstMonth = 0 - FirstMonth = numpy.int(datetime_to_matlab_datenum(datetime.datetime(FirstYear, FirstMonth+1, 1)) - 1) + FirstYear = FirstYear + 1 + FirstMonth = 0 + FirstMonth = numpy.int(datetime_to_matlab_datenum(datetime.datetime(FirstYear, FirstMonth + 1, 1)) - 1) # Find latest end-month consistent with LastDay. - [LastYear, LastMonth] = [datetime.datetime.fromordinal(numpy.int(LastDay-366)).year, datetime.datetime.fromordinal(numpy.int(LastDay-366)).month] + [LastYear, LastMonth] = [datetime.datetime.fromordinal(numpy.int(LastDay - 366)).year, + datetime.datetime.fromordinal(numpy.int(LastDay - 366)).month] if LastMonth == 12: - LastYear = LastYear+1 - LastMonth = 0 - LastMonth1 = numpy.int(datetime_to_matlab_datenum(datetime.datetime(LastYear, LastMonth+1, 1)) - 1) - if (LastMonth1 != LastDay): + LastYear = LastYear + 1 + LastMonth = 0 + LastMonth1 = numpy.int(datetime_to_matlab_datenum(datetime.datetime(LastYear, LastMonth + 1, 1)) - 1) + if LastMonth1 != LastDay: if LastMonth == 0: - LastYear = LastYear-1 - LastMonth = 12 + LastYear = LastYear - 1 + LastMonth = 12 LastMonth1 = datetime_to_matlab_datenum(datetime.datetime(LastYear, LastMonth, 1)) - 1 - - (StartT,) = numpy.where(MonthlyDateIndex==FirstMonth) - (EndT,) = numpy.where(MonthlyDateIndex==LastMonth1) - if ((numpy.size(StartT) > 0) and (numpy.size(EndT) > 0)): - YieldCurveDateIndex = MonthlyDateIndex[StartT:EndT+1] - YieldCurveData = MonthlyYieldCurveData[StartT:EndT+1,IncludeMaturities] + (StartT,) = numpy.where(MonthlyDateIndex == FirstMonth) + (EndT,) = numpy.where(MonthlyDateIndex == LastMonth1) + if (numpy.size(StartT) > 0) and (numpy.size(EndT) > 0): + YieldCurveDateIndex = MonthlyDateIndex[StartT[0]:EndT[0] + 1] + YieldCurveData = MonthlyYieldCurveData[StartT[0]:EndT[0] + 1, IncludeMaturities] else: YieldCurveDateIndex = numpy.array([]) YieldCurveData = numpy.array([]) @@ -137,44 +137,54 @@ Tau_K = numpy.copy(SampleMaturities) # Estimation. -if (FinalNaturalParametersGiven == 1): - print 'Finalizing model K-AFNSM(2) for ' + Country + " using %i" % Maturities[0] + "-%i" % SampleMaturities[-1] + ' year data at ' + DataFrequency + ' frequency for period ' + matlab_datenum_to_datetime(YieldCurveDateIndex[0]).strftime('%d-%b-%Y') + ' to ' + matlab_datenum_to_datetime(YieldCurveDateIndex[-1]).strftime('%d-%b-%Y') +if FinalNaturalParametersGiven == 1: + print('Finalizing model K-AFNSM(2) for ' + Country + " using %i" % Maturities[0] + "-%i" % SampleMaturities[ + -1] + ' year data at ' + DataFrequency + ' frequency for period ' + matlab_datenum_to_datetime( + YieldCurveDateIndex[0]).strftime('%d-%b-%Y') + ' to ' + matlab_datenum_to_datetime( + YieldCurveDateIndex[-1]).strftime('%d-%b-%Y')) FinalNaturalParameters = numpy.copy(InitialNaturalParameters) FINAL = 1 Exitflag = -1 - - [Fval, x_T] = AAC_KAGM_SingleLoop(numpy.matrix(YieldCurveData), Tau_K,N, FinalNaturalParameters, Dt, dTau, KappaP_Constraint, ZLB_Imposed, IEKF_Count, FINAL) + [Fval, x_T] = AAC_KAGM_SingleLoop(numpy.matrix(YieldCurveData), Tau_K, N, FinalNaturalParameters, Dt, dTau, + KappaP_Constraint, ZLB_Imposed, IEKF_Count, FINAL) Time0 = 0 Time1 = 0 Output = 'Final parameters given' rL = FinalNaturalParameters[0] KappaQ2 = FinalNaturalParameters[1] - KappaP = numpy.array([[FinalNaturalParameters[2], FinalNaturalParameters[3]], [FinalNaturalParameters[4], FinalNaturalParameters[5]]]) + KappaP = numpy.array([[FinalNaturalParameters[2], FinalNaturalParameters[3]], + [FinalNaturalParameters[4], FinalNaturalParameters[5]]]) ThetaP = numpy.array([[FinalNaturalParameters[6]], [FinalNaturalParameters[7]]]) Sigma1 = FinalNaturalParameters[8] Sigma2 = FinalNaturalParameters[9] Rho12 = FinalNaturalParameters[10] else: # Estimate final parameters. - print 'Estimating K-AFNSM(2) for ' + Country + " using %i" % SampleMaturities[0] + "-%i" % SampleMaturities[-1] + ' year data at ' + DataFrequency + ' frequency for period ' + matlab_datenum_to_datetime(YieldCurveDateIndex[0]).strftime('%d-%b-%Y') + ' to ' + matlab_datenum_to_datetime(YieldCurveDateIndex[-1]).strftime('%d-%b-%Y') + print('Estimating K-AFNSM(2) for ' + Country + " using %i" % SampleMaturities[0] + "-%i" % SampleMaturities[ + -1] + ' year data at ' + DataFrequency + ' frequency for period ' + matlab_datenum_to_datetime( + YieldCurveDateIndex[0]).strftime('%d-%b-%Y') + ' to ' + matlab_datenum_to_datetime( + YieldCurveDateIndex[-1]).strftime('%d-%b-%Y')) Exitflag = 0 while (Exitflag == 0): if (KappaP_Constraint == 'Direct'): InitialParameters = numpy.copy(InitialNaturalParameters) - tmp_func = lambda x: x/(1+numpy.abs(x))-InitialNaturalParameters[9] - InitialParameters[9] = scipy.optimize.fsolve(tmp_func,1) + tmp_func = lambda x: x / (1 + numpy.abs(x)) - InitialNaturalParameters[9] + InitialParameters[9] = scipy.optimize.fsolve(tmp_func, 1) elif (KappaP_Constraint == 'S/A'): - print 'Nothing here.' + print('Nothing here.') # Extended Kalman filter estimation. - #Time0 = matlab_now(datetime.datetime.now()) + # Time0 = matlab_now(datetime.datetime.now()) Time0 = datetime_to_matlab_datenum(datetime.datetime.now()) FINAL = 0 Max_IEKF_Count = 0 Max_IEKF_Point = 0 - [x_T, FinalParameters, Fval, Exitflag, Output] = AAB_KAGM_Estimation_NelderMead(YieldCurveData, Tau_K, N, InitialParameters, Dt, dTau, KappaP_Constraint, ZLB_Imposed, IEKF_Count, FINAL, Iterations) + [x_T, FinalParameters, Fval, Exitflag, Output] = AAB_KAGM_Estimation_NelderMead(YieldCurveData, Tau_K, N, + InitialParameters, Dt, dTau, + KappaP_Constraint, ZLB_Imposed, + IEKF_Count, FINAL, Iterations) Time1 = matlab_now(datetime.datetime.now()) if (KappaP_Constraint == 'Direct'): @@ -187,10 +197,11 @@ FinalNaturalParameters[11:] = abs(FinalParameters[11:]) # Calculate the state equation quantities based on parameter values. KappaQ = numpy.array([[0, 0], [FinalNaturalParameters[0], 0]]) - KappaP = numpy.array([[FinalNaturalParameters[1], FinalNaturalParameters[2]], [FinalNaturalParameters[3], FinalNaturalParameters[4]]]) + KappaP = numpy.array([[FinalNaturalParameters[1], FinalNaturalParameters[2]], + [FinalNaturalParameters[3], FinalNaturalParameters[4]]]) D, V = numpy.linalg.eig(KappaP) - d1 = D[0,0] - d2 = D[1,1] + d1 = D[0, 0] + d2 = D[1, 1] if ((d1.real < 0) or (d2.real < 0)): if (d1.real < 0): d1 = 1e-6 + d1.imag * 1j @@ -199,47 +210,49 @@ D = numpy.diag([d1, d2]) tmp1 = numpy.matrix(V) - tmp2 = tmp1 * numpy.matrix(numpy.reshape(D, (len(D),1))) + tmp2 = tmp1 * numpy.matrix(numpy.reshape(D, (len(D), 1))) KappaP = numpy.array(numpy.real(numpy.linalg.solve(tmp1.getH(), tmp2.getH()).getH())) - FinalNaturalParameters[2] = KappaP[0,0] - FinalNaturalParameters[3] = KappaP[0,1] - FinalNaturalParameters[4] = KappaP[1,0] - FinalNaturalParameters[5] = KappaP[1,1] + FinalNaturalParameters[2] = KappaP[0, 0] + FinalNaturalParameters[3] = KappaP[0, 1] + FinalNaturalParameters[4] = KappaP[1, 0] + FinalNaturalParameters[5] = KappaP[1, 1] elif (KappaP_Constraint == 'S/A'): - print 'Nothing here' + print('Nothing here') rL = FinalNaturalParameters[0] KappaQ2 = FinalNaturalParameters[1] - KappaP = numpy.array([[FinalNaturalParameters[2], FinalNaturalParameters[3]], [FinalNaturalParameters[4], FinalNaturalParameters[5]]]) + KappaP = numpy.array([[FinalNaturalParameters[2], FinalNaturalParameters[3]], + [FinalNaturalParameters[4], FinalNaturalParameters[5]]]) ThetaP = numpy.array([[FinalNaturalParameters[6]], [FinalNaturalParameters[7]]]) Sigma1 = FinalNaturalParameters[8] Sigma2 = FinalNaturalParameters[9] Rho12 = FinalNaturalParameters[10] # disp(Exitflag) - print [Max_IEKF_Point, Max_IEKF_Count] - print Fval - print FinalNaturalParameters[0:10] + print([Max_IEKF_Point, Max_IEKF_Count]) + print(Fval) + print(FinalNaturalParameters[0:10]) - #plotyy(1:length(x_T),x_T',1:length(x_T),sum(x_T)') - #pause 0.1 + # plotyy(1:length(x_T),x_T',1:length(x_T),sum(x_T)') + # pause 0.1 # Create figure figure = pyplot.figure(num=None, figsize=(8, 6), dpi=100, facecolor='w') - - subplot1 = figure.add_subplot(1,1,1, position=[0.15, 0.10, 0.75, 0.80], frame_on=True, zorder=0) - tmp1 = numpy.arange(0, max(numpy.shape(x_T))+1) + + subplot1 = figure.add_subplot(1, 1, 1, position=[0.15, 0.10, 0.75, 0.80], frame_on=True, zorder=0) + tmp1 = numpy.arange(0, max(numpy.shape(x_T)) + 1) subplot1.plot(tmp1, x_T.getH(), linewidth=2, marker='', markersize=3, zorder=1, label="") subplot2 = subplot1.twinx() subplot1.plot(tmp1, numpy.sum(x_T).getH(), linewidth=2, marker='', markersize=3, zorder=1, label="") - SaveName = AAL_CommonSaveName(DataFileName, ZLB_Imposed, IEKF_Count, SampleMaturities, N, DataFrequency, FINAL, -10); + SaveName = AAL_CommonSaveName(DataFileName, ZLB_Imposed, IEKF_Count, SampleMaturities, N, DataFrequency, FINAL, + -10); disp(SaveName) - + figure1 = pyplot.figure(num=None, figsize=(8, 6), dpi=100, facecolor='w') figure1.plot(tmp1, x_T.getH(), linewidth=2, marker='', markersize=3, zorder=1, label="") # Save final output in CSV file. NaturalParameterStandardErrors = -9.999 * numpy.ones(numpy.size(FinalNaturalParameters)) - numpy.savetxt(SaveName+".csv", FinalNaturalParameters, fmt='%25.15e', delimiter=',', newline='\n') + numpy.savetxt(SaveName + ".csv", FinalNaturalParameters, fmt='%25.15e', delimiter=',', newline='\n') # Reset InitialNaturalParameters for next iteration. InitialNaturalParameters = numpy.copy(FinalNaturalParameters) @@ -250,121 +263,117 @@ if (HessianRequired == 1): # Calculate Hessian and standard errors for natural model parameters. FINAL = 1 - NaturalHessian = AAF_FiniteDifferenceHessian(AAC_KAGM_SingleLoop, FinalNaturalParameters, 1e-10, YieldCurveData, Tau_K, N, Dt, dTau, KappaP_Constraint, ZLB_Imposed, IEKF_Count, FINAL) + NaturalHessian = AAF_FiniteDifferenceHessian(AAC_KAGM_SingleLoop, FinalNaturalParameters, 1e-10, YieldCurveData, + Tau_K, N, Dt, dTau, KappaP_Constraint, ZLB_Imposed, IEKF_Count, FINAL) NaturalParameterStandardErrors = numpy.sqrt(numpy.abs(numpy.diag(numpy.linalg.inv(NaturalHessian)))); else: NaturalParameterStandardErrors = -9.999 * numpy.ones(numpy.size(FinalNaturalParameters)) - # Display output. dTime = Time1 - Time0 -print dTime*24, 'hours (=', dTime*24*60, 'minutes)' -print Output -print Exitflag -print Fval -print InitialNaturalParameters[0:10] -print FinalNaturalParameters[0:10] -print NaturalParameterStandardErrors[0:10] -KappaQ = numpy.matrix([[0,0], [0,KappaQ2]]) -print KappaP, numpy.linalg.eig(KappaP), KappaQ-KappaP - -(T,K) = numpy.shape(YieldCurveData) -Residuals = numpy.ones((T,K)) * float('nan') +print(dTime * 24, 'hours (=', dTime * 24 * 60, 'minutes)') +print(Output) +print(Exitflag) +print(Fval) +print(InitialNaturalParameters[0:10]) +print(FinalNaturalParameters[0:10]) +print(NaturalParameterStandardErrors[0:10]) +KappaQ = numpy.matrix([[0, 0], [0, KappaQ2]]) +print(KappaP, numpy.linalg.eig(KappaP), KappaQ - KappaP) + +(T, K) = numpy.shape(YieldCurveData) +Residuals = numpy.ones((T, K)) * float('nan') PlotCurves = 0 - for t in range(0, T): - YieldCurveData_t = YieldCurveData[t,:] - (Fitted_R_t,tmp1) = AAD_KAGM_R_and_dR_dx(x_T[:,t], rL, KappaQ2, Sigma1, Sigma2, Rho12, Tau_K, dTau, ZLB_Imposed) + YieldCurveData_t = YieldCurveData[t, :] + (Fitted_R_t, tmp1) = AAD_KAGM_R_and_dR_dx(x_T[:, t], rL, KappaQ2, Sigma1, Sigma2, Rho12, Tau_K, dTau, ZLB_Imposed) Fitted_R_t = numpy.reshape(numpy.array(Fitted_R_t), (numpy.size(Fitted_R_t),)) if (PlotCurves == 1): figure = pyplot.figure(num=None, figsize=(8, 6), dpi=100, facecolor='w') - subplot = figure.add_subplot(1,1,1, position=[0.15, 0.10, 0.75, 0.80], frame_on=True, zorder=0) - subplot.plot(Tau_K, 100*Fitted_R_t, linewidth=2, marker='', markersize=3, zorder=1, label="") + subplot = figure.add_subplot(1, 1, 1, position=[0.15, 0.10, 0.75, 0.80], frame_on=True, zorder=0) + subplot.plot(Tau_K, 100 * Fitted_R_t, linewidth=2, marker='', markersize=3, zorder=1, label="") subplot.plot(Tau_K, YieldCurveData_t, linestyle='', linewidth=2, marker='o', markersize=3, zorder=1, label="") - #ylim([-2 10]); - #pause(0.3) + # ylim([-2 10]); + # pause(0.3) Residual_t = 0.01 * YieldCurveData_t - Fitted_R_t - print numpy.shape(Residuals) - print numpy.shape(Residual_t) - Residuals[t,:] = Residual_t + print(numpy.shape(Residuals)) + print(numpy.shape(Residual_t)) + Residuals[t, :] = Residual_t # -RMSE_Residuals = numpy.sqrt(numpy.sum(numpy.multiply(Residuals,Residuals)) / T) -print numpy.mean(Residuals) -print RMSE_Residuals +RMSE_Residuals = numpy.sqrt(numpy.sum(numpy.multiply(Residuals, Residuals)) / T) +print(numpy.mean(Residuals)) +print(RMSE_Residuals) Phi = KappaQ2 -(SSR,EMS_Q,ETZ_Q) = AAH_EMS_N23_function(Phi, x_T, dTau) +(SSR, EMS_Q, ETZ_Q) = AAH_EMS_N23_function(Phi, x_T, dTau) # Save final output in MatLab file. SaveName = AAL_CommonSaveName(DataFileName, ZLB_Imposed, IEKF_Count, SampleMaturities, N, DataFrequency, FINAL, -10); -print SaveName +print(SaveName) # Save final output in CSV file. -#numpy.savetxt(SaveName+".csv", FinalNaturalParameters, fmt='%25.15e', delimiter=',', newline='\n') -#save(SaveName) +# numpy.savetxt(SaveName+".csv", FinalNaturalParameters, fmt='%25.15e', delimiter=',', newline='\n') +# save(SaveName) # Save final output in Excel spreadsheet. -RangeName = 'A2:O%i' % (max(numpy.shape(YieldCurveDateIndex))+1) -#xlswrite([SaveName,'.xlsm'],[YieldCurveDateIndex-datenum('30-Dec-1899'),YieldCurveData,100*x_T',100*SSR,100*EMS_Q,ETZ_Q],'D. Monthly',RangeName) +RangeName = 'A2:O%i' % (max(numpy.shape(YieldCurveDateIndex)) + 1) +# xlswrite([SaveName,'.xlsm'],[YieldCurveDateIndex-datenum('30-Dec-1899'),YieldCurveData,100*x_T',100*SSR,100*EMS_Q,ETZ_Q],'D. Monthly',RangeName) SSR_pos = SSR.copy() SSR_pos[SSR_pos <= 0] = numpy.nan -fig,ax=pyplot.subplots() -pyplot.plot(SSR_pos*100,linewidth=2, marker='', markersize=3, zorder=1, label="slope", color='b') +fig, ax = pyplot.subplots() +pyplot.plot(SSR_pos * 100, linewidth=2, marker='', markersize=3, zorder=1, label="slope", color='b') SSR_neg = SSR.copy() SSR_neg[SSR_neg > 0] = numpy.nan -pyplot.plot(SSR_neg*100,linewidth=2, marker='', markersize=3, zorder=1, label="slope", color='r') - - +pyplot.plot(SSR_neg * 100, linewidth=2, marker='', markersize=3, zorder=1, label="slope", color='r') -a=ax.get_xticks() -y=[matlab_datenum_to_datetime(YieldCurveDateIndex[int(x)]).date().year for x in a[:-1]] -pyplot.xticks(a,y, rotation=90) +a = ax.get_xticks() +y = [matlab_datenum_to_datetime(YieldCurveDateIndex[int(x)]).date().year for x in a[:-1]] +pyplot.xticks(a, y, rotation=90) pyplot.ylabel('Percentage') fig.suptitle('SSR', fontsize=20) -fig.savefig('SSR_'+DataFrequency+'.jpg') +fig.savefig('SSR_' + DataFrequency + '.jpg') +fig, ax = pyplot.subplots() +pyplot.plot(EMS_Q * 100, linewidth=2, marker='', markersize=3, zorder=1, label="slope") -fig,ax=pyplot.subplots() -pyplot.plot(EMS_Q*100,linewidth=2, marker='', markersize=3, zorder=1, label="slope") - -a=ax.get_xticks() -y=[matlab_datenum_to_datetime(YieldCurveDateIndex[int(x)]).date().year for x in a[:-1]] -pyplot.xticks(a,y, rotation=90) +a = ax.get_xticks() +y = [matlab_datenum_to_datetime(YieldCurveDateIndex[int(x)]).date().year for x in a[:-1]] +pyplot.xticks(a, y, rotation=90) pyplot.ylabel('Percentage') fig.suptitle('EMS', fontsize=20) -fig.savefig('EMS_'+DataFrequency+'.jpg') +fig.savefig('EMS_' + DataFrequency + '.jpg') -fig,ax=pyplot.subplots() -pyplot.plot(ETZ_Q,linewidth=2, marker='', markersize=3, zorder=1, label="slope") +fig, ax = pyplot.subplots() +pyplot.plot(ETZ_Q, linewidth=2, marker='', markersize=3, zorder=1, label="slope") -a=ax.get_xticks() -y=[matlab_datenum_to_datetime(YieldCurveDateIndex[int(x)]).date().year for x in a[:-1]] -pyplot.xticks(a,y, rotation=90) +a = ax.get_xticks() +y = [matlab_datenum_to_datetime(YieldCurveDateIndex[int(x)]).date().year for x in a[:-1]] +pyplot.xticks(a, y, rotation=90) pyplot.ylabel('Years') fig.suptitle('ETZ', fontsize=20) -fig.savefig('ETZ_'+DataFrequency+'.jpg') +fig.savefig('ETZ_' + DataFrequency + '.jpg') pyplot.show() -output=numpy.concatenate((YieldCurveData, 100*x_T.T), axis=1) -output=numpy.concatenate((output, 100*SSR.reshape((SSR.shape[0],1))), axis=1) -output=numpy.concatenate((output, 100*EMS_Q.reshape((SSR.shape[0],1))), axis=1) -output=numpy.concatenate((output, 100*ETZ_Q.reshape((SSR.shape[0],1))), axis=1) +output = numpy.concatenate((YieldCurveData, 100 * x_T.T), axis=1) +output = numpy.concatenate((output, 100 * SSR.reshape((SSR.shape[0], 1))), axis=1) +output = numpy.concatenate((output, 100 * EMS_Q.reshape((SSR.shape[0], 1))), axis=1) +output = numpy.concatenate((output, 100 * ETZ_Q.reshape((SSR.shape[0], 1))), axis=1) date_str = [matlab_datenum_to_datetime(x).date().__str__() for x in YieldCurveDateIndex] -k=numpy.array(date_str) -k=k.reshape((len(k), 1)) -with open(SaveName+'_final.csv', 'wb') as csvfile: +k = numpy.array(date_str) +k = k.reshape((len(k), 1)) +with open(SaveName + '_final.csv', 'w') as csvfile: spamwriter = csv.writer(csvfile, delimiter=',') - header =['Date']+[str(x) for x in SampleMaturities]+['Level','Slope','SSR','EMS-Q','ETZ-Q'] + header = ['Date'] + [str(x) for x in SampleMaturities] + ['Level', 'Slope', 'SSR', 'EMS-Q', 'ETZ-Q'] spamwriter.writerow(header) - for i in range(len(date_str)): + for i in range(len(date_str)): data_to_write = list(k[i]) + list(output[i]) spamwriter.writerow(data_to_write) - -print 'Finished' + +print('Finished') diff --git a/AAC_KAGM_SingleLoop_file.py b/AAC_KAGM_SingleLoop_file.py index ee06711..cd107f2 100644 --- a/AAC_KAGM_SingleLoop_file.py +++ b/AAC_KAGM_SingleLoop_file.py @@ -73,7 +73,7 @@ def AAC_KAGM_SingleLoop(R_data, Tau_K, N, Parameters, Dt, dTau, KappaP_Constrain F = numpy.matrix(scipy.linalg.expm(-KappaP*Dt)) D, V = numpy.linalg.eig(F) if (numpy.any(numpy.abs(D)>1.0001)): - print 'SingleLoop:argChk abs(eig(F))>1' + print('SingleLoop:argChk abs(eig(F))>1') exit() Q = numpy.matrix([[G(2*d1,Dt), G(d1+d2,Dt)], [0, G(2*d2,Dt)]]) @@ -116,8 +116,8 @@ def AAC_KAGM_SingleLoop(R_data, Tau_K, N, Parameters, Dt, dTau, KappaP_Constrain # Observations for time t. y_Obs = numpy.matrix(0.01 * R_data[t,:]).getH() y_Missing = numpy.squeeze(numpy.array(numpy.isnan(y_Obs))) #numpy.array(numpy.isnan(y_Obs))#[:,0] - y_Obs = y_Obs[-y_Missing] - R = numpy.diag(numpy.power(Sigma_Nu[-y_Missing], 2)) + y_Obs = y_Obs[~y_Missing] + R = numpy.diag(numpy.power(Sigma_Nu[~y_Missing], 2)) x_Plus_i_Minus_1 = numpy.copy(x_Minus) x_Plus_i0 = numpy.copy(x_Minus) @@ -131,18 +131,18 @@ def AAC_KAGM_SingleLoop(R_data, Tau_K, N, Parameters, Dt, dTau, KappaP_Constrain (y_Hat, H_i) = AAD_KAGM_R_and_dR_dx(x_Plus_i0, rL, KappaQ2, Sigma1, Sigma2, Rho12, Tau_K, dTau, ZLB_Imposed) - y_Hat = y_Hat[-y_Missing] - H_i = H_i[-y_Missing,:] + y_Hat = y_Hat[~y_Missing] + H_i = H_i[~y_Missing,:] HPHR_i = (H_i * P_Minus * H_i.getH() + R) K_i = numpy.linalg.solve(HPHR_i.getH(), (P_Minus*H_i.getH()).getH()).getH() w_i = y_Obs - y_Hat - H_i * (x_Minus - x_Plus_i0) x_Plus_i1 = x_Minus + K_i * w_i if (IEKF_Count == 20): - # Using tolerance, so check for convergence. + # Using tolerance, so check for convergence. if (i > 15): # Large number of iterations, so print output to screen. - print t, i-1, numpy.matrix(x_Plus_i1).getH(), numpy.matrix(x_Plus_i0).getH(), numpy.matrix(x_Plus_i1).getH() - numpy.matrix(x_Plus_i0).getH() + print(t, i-1, numpy.matrix(x_Plus_i1).getH(), numpy.matrix(x_Plus_i0).getH(), numpy.matrix(x_Plus_i1).getH() - numpy.matrix(x_Plus_i0).getH()) if (numpy.all(numpy.abs(x_Plus_i1-x_Plus_i0)0) + numpy.linalg.solve(HPHR_i.getH(),w_i).getH() * w_i # Hold IEKF count. + if Max_IEKF_Count is None: + Max_IEKF_Count = IEKF_Count if (i-1 > Max_IEKF_Count): Max_IEKF_Count = i - 1 Max_IEKF_Point = t @@ -182,7 +184,7 @@ def AAC_KAGM_SingleLoop(R_data, Tau_K, N, Parameters, Dt, dTau, KappaP_Constrain # Negate the log likelihood value because fminunc minimizes. EKF_logL = -EKF_logL - print EKF_logL*1e-3, Parameters[0:10], Rho12 + print(EKF_logL*1e-3, Parameters[0:10], Rho12) tmp1 = numpy.sum(x_T, axis=0) figure = pyplot.figure(num=None, figsize=(8, 6), dpi=100, facecolor='w') @@ -193,7 +195,4 @@ def AAC_KAGM_SingleLoop(R_data, Tau_K, N, Parameters, Dt, dTau, KappaP_Constrain pyplot.savefig("plot.pdf") pyplot.show() - return (EKF_logL,x_T) - - - + return (EKF_logL,x_T) \ No newline at end of file diff --git a/AAE_KAGM_f_and_df_dx_file.py b/AAE_KAGM_f_and_df_dx_file.py index 2d94fc1..20e92c1 100644 --- a/AAE_KAGM_f_and_df_dx_file.py +++ b/AAE_KAGM_f_and_df_dx_file.py @@ -33,7 +33,7 @@ def AAE_KAGM_f_and_df_dx(x_t, rL, KappaQ2, Sigma1, Sigma2, Rho12, TauGrid, ZLB_I Omega = numpy.sqrt(Sigma1**2 * G_11 + Sigma2**2 * G_22 + 2 * Rho12 * Sigma1 * Sigma2 * G_12) # Calculate cumulative normal probabilities for N(0,1) distribution. - d = numpy.divide(GATSM_f-rL, Omega) + d = numpy.divide(GATSM_f-rL, Omega,where=Omega!=0) normsdist_erf_d = normsdist_erf(d) # Calculate gradiant and CAB_GATSM_f diff --git a/data_read.py b/data_read.py index 4112f0d..a792525 100755 --- a/data_read.py +++ b/data_read.py @@ -8,8 +8,8 @@ import os def write_to_csv(idx,idx1,val,val1,sub_datenum,name): - daily_values_data=map(list,zip(*val)) - daily_values_data1=map(list,zip(*val1)) + daily_values_data=list(map(list,list(zip(*val)))) + daily_values_data1=list(map(list,list(zip(*val1)))) daily_date_index=list(idx) daily_date_index[:] = [x - sub_datenum for x in idx] @@ -22,19 +22,19 @@ def write_to_csv(idx,idx1,val,val1,sub_datenum,name): for i in range(len(daily_values_data1)): daily_values_data1[i].insert(0,daily_date_index1[i]) - with open(name+'.csv', 'wb') as csvfile: + with open(name+'.csv', 'w') as csvfile: spamwriter = csv.writer(csvfile, delimiter=',',quotechar='|', quoting=csv.QUOTE_NONNUMERIC) spamwriter.writerows(daily_values_data1) spamwriter.writerows(daily_values_data) def filter_data_index(a,thr,method): - if method is 'ge': + if method =='ge': index=list([idx for idx,i in enumerate(a) if ithr]) - elif method is 'l': + elif method =='l': index=list([idx for idx,i in enumerate(a) if i>=thr]) - elif method is 'g': + elif method =='g': index=list([idx for idx,i in enumerate(a) if i<=thr]) return index @@ -62,16 +62,16 @@ def filter_data(a,index,dim): CurveType='OIS' PathName=os.getcwd() ExcelName=os.path.join(PathName, 'A_'+Country+'_All_Data_Bloomberg.xlsx') -if CurveType is 'OIS': +if CurveType =='OIS': ExcelSheetName='D. Live OIS data' -if CurveType is 'Govt': +if CurveType =='Govt': ExcelSheetName='D. Live Govt data'; Maturities=[0.25,0.5,1,2,3,4,5,7,10,15,20,30] wb=openpyxl.load_workbook(ExcelName) -sheet_names=wb.get_sheet_names() -sheet=wb.get_sheet_by_name('D. Live OIS data') +sheet_names=wb.get_sheet_names +sheet=wb['D. Live OIS data'] datenum = [[] for i in range(len(datelist))] values_data=[[] for i in range(len(datelist))] for col_num in range(len(datelist)): @@ -97,15 +97,15 @@ def filter_data(a,index,dim): diff_dates=list(set(common_datenum) ^ set(datenum[i])) for j in diff_dates: index=datenum[i].index(j) - print(str(j)+' date num at '+str(index)+'removed') + print((str(j)+' date num at '+str(index)+'removed')) value_to_remove = values_data[i][index] values_data[i].remove(value_to_remove) -ref_friday=datetime.date(2013,05,24).toordinal() +ref_friday=datetime.date(2013,0o5,24).toordinal() weeks_to_step_back = (ref_friday - common_datenum[0])/7 first_friday = ref_friday - (weeks_to_step_back*7) -week_date_index = range(first_friday,common_datenum[-1],7) +week_date_index = list(range(int(first_friday),common_datenum[-1],7)) week_values_data=[[] for i in range(len(datelist))] for i in range(len(datelist)): @@ -139,7 +139,7 @@ def filter_data(a,index,dim): #GOVT -if Country is 'EA': +if Country =='EA': wbG=openpyxl.load_workbook(os.path.join(PathName, 'A_GE_All_Data_Bloomberg.xlsx')) sheet=wbG.get_sheet_by_name('D. Live Govt data') datenumG = [[] for i in range(len(datelist))] @@ -168,7 +168,7 @@ def filter_data(a,index,dim): diff_dates=list(set(common_datenumG) ^ set(datenumG[i])) for j in diff_dates: index=datenumG[i].index(j) - print(str(j)+' date num at '+str(index)+'removed') + print((str(j)+' date num at '+str(index)+'removed')) value_to_remove = values_dataG[i][index] values_dataG[i].remove(value_to_remove) @@ -210,7 +210,7 @@ def filter_data(a,index,dim): diff_dates=list(set(common_datenumF) ^ set(datenumF[i])) for j in diff_dates: index=datenumG[i].index(j) - print(str(j)+' date num at '+str(index)+'removed') + print((str(j)+' date num at '+str(index)+'removed')) value_to_remove = values_dataF[i][index] values_dataF[i].remove(value_to_remove) @@ -227,14 +227,14 @@ def filter_data(a,index,dim): diff_dates=list(set(post_inter_index) ^ set(post_euro_fr_index)) for j in diff_dates: index=post_euro_fr_index.index(j) - print(str(j)+' date num at '+str(index)+'removed') + print((str(j)+' date num at '+str(index)+'removed')) value_to_remove = post_euro_fr_vlues[i][index] post_euro_fr_vlues[i].remove(value_to_remove) for i in range(len(post_euro_ge_vlues)): diff_dates=list(set(post_inter_index) ^ set(post_euro_ge_index)) for j in diff_dates: index=post_euro_ge_index.index(j) - print(str(j)+' date num at '+str(index)+'removed') + print((str(j)+' date num at '+str(index)+'removed')) value_to_remove = post_euro_ge_vlues[i][index] post_euro_ge_vlues[i].remove(value_to_remove) @@ -246,7 +246,7 @@ def filter_data(a,index,dim): values_data1[i] = pre_euro_ge_vlues[i] + post_euro_avg_values else: - sheet=wb.get_sheet_by_name('D. Live Govt data') + sheet=wb['D. Live Govt data'] datenum1 = [[] for i in range(len(datelist))] values_data1=[[] for i in range(len(datelist))] for col_num in range(len(datelist)): @@ -273,15 +273,15 @@ def filter_data(a,index,dim): diff_dates=list(set(common_datenum1) ^ set(datenum1[i])) for j in diff_dates: index=datenum1[i].index(j) - print(str(j)+' date num at '+str(index)+'removed') + print((str(j)+' date num at '+str(index)+'removed')) value_to_remove = values_data1[i][index] values_data1[i].remove(value_to_remove) -ref_friday=datetime.date(2013,05,24).toordinal() +ref_friday=datetime.date(2013,0o5,24).toordinal() weeks_to_step_back = (ref_friday - common_datenum1[0])/7 first_friday = ref_friday - (weeks_to_step_back*7) -week_date_index1 = range(first_friday,common_datenum1[-1],7) +week_date_index1 = list(range(int(first_friday),common_datenum1[-1],7)) week_values_data1=[[] for i in range(len(datelist))] for i in range(len(datelist)): f=interpolate.interp1d(common_datenum1,values_data1[i],'zero',bounds_error=False) diff --git a/functions.py b/functions.py index c759d81..b659783 100644 --- a/functions.py +++ b/functions.py @@ -1,15 +1,15 @@ import datetime import numpy + ####################################################################################################################### def datetime_to_matlab_datenum(dt): return 366 + 1 + (dt - datetime.datetime(1, 1, 1, 0, 0)).total_seconds() / datetime.timedelta(1).total_seconds() + ####################################################################################################################### def matlab_datenum_to_datetime(dn): # Convert a datetime object to its MATLAB datenum() equivalent - return datetime.datetime(1,1,1,0,0) + datetime.timedelta(dn-366-1) + return datetime.datetime(1, 1, 1, 0, 0) + datetime.timedelta(int(dn - 366 - 1)) ####################################################################################################################### - -