From 0f9eb185e23b79440ddb22d1a25103cf9da9b2cd Mon Sep 17 00:00:00 2001 From: Mridul Seth Date: Fri, 12 Jul 2019 12:40:19 +0200 Subject: [PATCH] comment out parallel SGU_solver, works for py3 --- .../TwoAssetCode/FluctuationsTwoAsset.py | 782 +++++++++--------- 1 file changed, 384 insertions(+), 398 deletions(-) diff --git a/HARK/BayerLuetticke/BayerLuetticke_code/TwoAssetCode/FluctuationsTwoAsset.py b/HARK/BayerLuetticke/BayerLuetticke_code/TwoAssetCode/FluctuationsTwoAsset.py index d494eecb9..ce446118f 100644 --- a/HARK/BayerLuetticke/BayerLuetticke_code/TwoAssetCode/FluctuationsTwoAsset.py +++ b/HARK/BayerLuetticke/BayerLuetticke_code/TwoAssetCode/FluctuationsTwoAsset.py @@ -10,7 +10,7 @@ import numpy as np from numpy.linalg import matrix_rank import scipy as sc -from scipy.stats import norm +from scipy.stats import norm from scipy.interpolate import interp1d, interp2d, griddata, RegularGridInterpolator, interpn import multiprocessing as mp from multiprocessing import Pool, cpu_count, Process @@ -1013,409 +1013,395 @@ def FF_2(range_, Xss,Yss,Gamma_state,indexMUdct,indexVKdct,par,mpar,grid,targets out_bl2.put(bl) -### This is the old parallel code that does not work in Python 3 -### We should fix this, potentially by using joblib Parallel instead of multiprocess -#def SGU_solver(Xss,Yss,Gamma_state,indexMUdct,indexVKdct,par,mpar,grid,targets,Copula,P_H,aggrshock): -# -# out_DF1 = mp.Queue() -# out_DF3 = mp.Queue() -# out_DF2 = mp.Queue() -# out_bl = mp.Queue() -# out_bl2 = mp.Queue() -# -# State = np.zeros((mpar['numstates'],1)) -# State_m = State.copy() -# Contr = np.zeros((mpar['numcontrols'],1)) -# Contr_m = Contr.copy() -# -# F = lambda S, S_m, C, C_m : Fsys(S, S_m, C, C_m, -# Xss,Yss,Gamma_state,indexMUdct,indexVKdct, -# par,mpar,grid,targets,Copula,P_H,aggrshock) -# -# -# start_time = time.clock() -# result_F = F(State,State_m,Contr.copy(),Contr_m.copy()) -# end_time = time.clock() -# print('Elapsed time is ', (end_time-start_time), ' seconds.') -# Fb=result_F['Difference'].copy() -# -# pool=cpu_count()/2 -# -# F1=np.zeros((mpar['numstates'] + mpar['numcontrols'], mpar['numstates'])) -# F2=np.zeros((mpar['numstates'] + mpar['numcontrols'], mpar['numcontrols'])) -# F3=np.zeros((mpar['numstates'] + mpar['numcontrols'], mpar['numstates'])) -# F4=np.asmatrix(np.vstack((np.zeros((mpar['numstates'], mpar['numcontrols'])), np.eye(mpar['numcontrols'],mpar['numcontrols']) ))) -# -# print('Use Schmitt Grohe Uribe Algorithm') -# print(' A *E[xprime uprime] =B*[x u]') -# print(' A = (dF/dxprimek dF/duprime), B =-(dF/dx dF/du)') -# -# #numscale=1 -# pnum=pool -# packagesize=int(ceil(mpar['numstates'] / float(3.0*pnum))) -# blocks=int(ceil(mpar['numstates'] / float(packagesize) )) -# -# par['scaleval1'] = 1e-5 -# par['scaleval2'] = 1e-5 -# -# start_time = time.clock() -# print('Computing Jacobian F1=DF/DXprime F3 =DF/DX') -# print('Total number of parallel blocks: ', str(blocks), '.') -# procs1=[] -# for bl in range(0,blocks): -# range_= range(bl*packagesize, min(packagesize*(bl+1),mpar['numstates'])) -# -# cc=np.zeros((mpar['numcontrols'],1)) -# ss=np.zeros((mpar['numstates'],1)) -# p1=Process(target=FF_1_3, args=(range_, Xss,Yss,Gamma_state,indexMUdct,indexVKdct,par,mpar, -# grid,targets,Copula,P_H,aggrshock,Fb,packagesize, bl,ss,cc, -# out_DF1,out_DF3, out_bl)) -# procs1.append(p1) -# p1.start() -# -# print('Block number: ', str(bl)) -# -# FF1 = [] -# FF3 = [] -# order_bl = [] -# -# for i in range(blocks): -# FF1.append(out_DF1.get()) -# FF3.append(out_DF3.get()) -# order_bl.append(out_bl.get()) -# -# -# print('bl order') -# print(order_bl) -# -# for p1 in procs1: -# p1.join() -# -# for i in range(0,int(ceil(mpar['numstates'] / float(packagesize)) )): -# range_= range(i*packagesize, min(packagesize*(i+1),mpar['numstates'])) -# F1[:,range_]=FF1[order_bl.index(i)].copy() -# F3[:,range_]=FF3[order_bl.index(i)].copy() -# -# end_time = time.clock() -# print('Elapsed time is ', (end_time-start_time), ' seconds.') -# -# # jacobian wrt Y' -# packagesize=int(ceil(mpar['numcontrols'] / (3.0*pnum))) -# blocks=int(ceil(mpar['numcontrols'] / float(packagesize))) -# print('Computing Jacobian F2 - DF/DYprime') -# print('Total number of parallel blocks: ', str(blocks),'.') -# -# start_time = time.clock() -# -# procs2=[] -# for bl in range(0,blocks): -# range_= range(bl*packagesize,min(packagesize*(bl+1),mpar['numcontrols'])) -# -# cc=np.zeros((mpar['numcontrols'],1)) -# ss=np.zeros((mpar['numstates'],1)) -# p2=Process(target=FF_2, args=(range_, Xss,Yss,Gamma_state,indexMUdct,indexVKdct,par,mpar, -# grid,targets,Copula,P_H,aggrshock,Fb,packagesize, bl,ss,cc, -# out_DF2, out_bl2)) -# procs2.append(p2) -# p2.start() -# print('Block number: ', str(bl)) -# -# FF=[] -# order_bl2 = [] -# -# for i in range(blocks): -# FF.append(out_DF2.get()) -# order_bl2.append(out_bl2.get()) -# -# print('bl2 order') -# print(order_bl2) -# -# -# for p2 in procs2: -# p2.join() -# -# for i in range(0,int( ceil(mpar['numcontrols'] / float(packagesize)) ) ): -# range_=range(i*packagesize, min(packagesize*(i+1),mpar['numcontrols'])) -# -# F2[:,range_]=FF[order_bl2.index(i)] -# -# end_time = time.clock() -# print('Elapsed time is ', (end_time-start_time), ' seconds.') -# -# FF=[] -# FF1=[] -# FF3=[] -# -# cc=np.zeros((mpar['numcontrols'],1)) -# ss=np.zeros((mpar['numstates'],1)) -# -# for Yct in range(0, mpar['oc']): -# Y=np.zeros((mpar['numcontrols'],1)) -# h=par['scaleval2'] -# Y[-1-Yct]=h -# Fx=F(ss.copy(),ss.copy(),cc.copy(),Y) -# F4[:,-1 - Yct]=(Fx['Difference'] - Fb) / h -# -# F2[mpar['nm']+mpar['nk']-3:mpar['numstates']-2,:] = 0 -# -# s,t,Q,Z=linalg.qz(np.hstack((F1,F2)), -np.hstack((F3,F4)), output='complex') -# abst = abs(np.diag(t))*(abs(np.diag(t))!=0.)+ (abs(np.diag(t))==0.)*10**(-11) -# #relev=np.divide(abs(np.diag(s)), abs(np.diag(t))) -# relev=np.divide(abs(np.diag(s)), abst) -# -# ll=sorted(relev) -# slt=relev >= 1 -# nk=sum(slt) -# slt=1*slt -# -# -# s_ord,t_ord,__,__,__,Z_ord=linalg.ordqz(np.hstack((F1,F2)), -np.hstack((F3,F4)), sort='ouc', output='complex') -# -# def sortOverridEigen(x, y): -# out = np.empty_like(x, dtype=bool) -# xzero = (x == 0) -# yzero = (y == 0) -# out[xzero & yzero] = False -# out[~xzero & yzero] = True -# out[~yzero] = (abs(x[~yzero]/y[~yzero]) > ll[-1 - mpar['numstates']]) -# return out -# -# if nk > mpar['numstates']: -# if mpar['overrideEigen']: -# print('Warning: The Equilibrium is Locally Indeterminate, critical eigenvalue shifted to: ', str(ll[-1 - mpar['numstates']])) -# slt=relev > ll[-1 - mpar['numstates']] -# nk=sum(slt) -# s_ord,t_ord,__,__,__,Z_ord=linalg.ordqz(np.hstack((F1,F2)), -np.hstack((F3,F4)), sort=sortOverridEigen, output='complex') -# -# else: -# print('No Local Equilibrium Exists, last eigenvalue: ', str(ll[-1 - mpar['numstates']])) -# -# elif nk < mpar['numstates']: -# if mpar['overrideEigen']: -# print('Warning: No Local Equilibrium Exists, critical eigenvalue shifted to: ', str(ll[-1 - mpar['numstates']])) -# slt=relev > ll[-1 - mpar['numstates']] -# nk=sum(slt) -# s_ord,t_ord,__,__,__,Z_ord=linalg.ordqz(np.hstack((F1,F2)), -np.hstack((F3,F4)), sort=sortOverridEigen, output='complex') -# -# else: -# print('No Local Equilibrium Exists, last eigenvalue: ', str(ll[-1 - mpar['numstates']])) -# -# -# -# -# z21=Z_ord[nk:,0:nk] -# z11=Z_ord[0:nk,0:nk] -# s11=s_ord[0:nk,0:nk] -# t11=t_ord[0:nk,0:nk] -# -# if matrix_rank(z11) < nk: -# print ('Warning: invertibility condition violated') -# -# z11i = np.dot(np.linalg.inv(z11), np.eye(nk)) # compute the solution -# -# gx = np.real(np.dot(z21,z11i)) -# hx = np.real(np.dot(z11,np.dot(np.dot(np.linalg.inv(s11),t11),z11i))) -# -# return{'hx': hx, 'gx': gx, 'F1': F1, 'F2': F2, 'F3': F3, 'F4': F4, 'par': par } -# return{'F1': F1, 'F2': F2, 'F3': F3, 'F4': F4, 'par': par } - -def SGU_solver(Xss,Yss,Gamma_state,indexMUdct,indexVKdct,par,mpar,grid,targets,Copula,P_H,aggrshock): # - - State = np.zeros((mpar['numstates'],1)) - State_m = State.copy() - Contr = np.zeros((mpar['numcontrols'],1)) - Contr_m = Contr.copy() - +def SGU_solver(Xss,Yss,Gamma_state,indexMUdct,indexVKdct,par,mpar,grid,targets,Copula,P_H,aggrshock): + + out_DF1 = mp.Queue() + out_DF3 = mp.Queue() + out_DF2 = mp.Queue() + out_bl = mp.Queue() + out_bl2 = mp.Queue() + + State = np.zeros((mpar['numstates'],1)) + State_m = State.copy() + Contr = np.zeros((mpar['numcontrols'],1)) + Contr_m = Contr.copy() + + F = lambda S, S_m, C, C_m : Fsys(S, S_m, C, C_m, + Xss,Yss,Gamma_state,indexMUdct,indexVKdct, + par,mpar,grid,targets,Copula,P_H,aggrshock) + + start_time = time.perf_counter() + result_F = F(State,State_m,Contr.copy(),Contr_m.copy()) + end_time = time.perf_counter() + print('Elapsed time is ', (end_time-start_time), ' seconds.') + Fb=result_F['Difference'].copy() + + pool=cpu_count()/2 + + F1=np.zeros((mpar['numstates'] + mpar['numcontrols'], mpar['numstates'])) + F2=np.zeros((mpar['numstates'] + mpar['numcontrols'], mpar['numcontrols'])) + F3=np.zeros((mpar['numstates'] + mpar['numcontrols'], mpar['numstates'])) + F4=np.asmatrix(np.vstack((np.zeros((mpar['numstates'], mpar['numcontrols'])), np.eye(mpar['numcontrols'],mpar['numcontrols']) ))) + + print('Use Schmitt Grohe Uribe Algorithm') + print(' A *E[xprime uprime] =B*[x u]') + print(' A = (dF/dxprimek dF/duprime), B =-(dF/dx dF/du)') + + #numscale=1 + pnum=pool + packagesize=int(ceil(mpar['numstates'] / float(3.0*pnum))) + blocks=int(ceil(mpar['numstates'] / float(packagesize) )) + + par['scaleval1'] = 1e-5 + par['scaleval2'] = 1e-5 + + start_time = time.perf_counter() + print('Computing Jacobian F1=DF/DXprime F3 =DF/DX') + print('Total number of parallel blocks: ', str(blocks), '.') + procs1=[] + for bl in range(0,blocks): + range_= range(bl*packagesize, min(packagesize*(bl+1),mpar['numstates'])) + cc=np.zeros((mpar['numcontrols'],1)) + ss=np.zeros((mpar['numstates'],1)) + p1=Process(target=FF_1_3, args=(range_, Xss,Yss,Gamma_state,indexMUdct,indexVKdct,par,mpar, + grid,targets,Copula,P_H,aggrshock,Fb,packagesize, bl,ss,cc, + out_DF1,out_DF3, out_bl)) + procs1.append(p1) + p1.start() + print('Block number: ', str(bl)) + + FF1 = [] + FF3 = [] + order_bl = [] + + for i in range(blocks): + FF1.append(out_DF1.get()) + FF3.append(out_DF3.get()) + order_bl.append(out_bl.get()) + + + print('bl order') + print(order_bl) + + for p1 in procs1: + p1.join() + + for i in range(0,int(ceil(mpar['numstates'] / float(packagesize)) )): + range_= range(i*packagesize, min(packagesize*(i+1),mpar['numstates'])) + F1[:,range_]=FF1[order_bl.index(i)].copy() + F3[:,range_]=FF3[order_bl.index(i)].copy() + + end_time = time.perf_counter() + print('Elapsed time is ', (end_time-start_time), ' seconds.') + + # jacobian wrt Y' + packagesize=int(ceil(mpar['numcontrols'] / (3.0*pnum))) + blocks=int(ceil(mpar['numcontrols'] / float(packagesize))) + print('Computing Jacobian F2 - DF/DYprime') + print('Total number of parallel blocks: ', str(blocks),'.') + + start_time = time.perf_counter() + + procs2=[] + for bl in range(0,blocks): + range_= range(bl*packagesize,min(packagesize*(bl+1),mpar['numcontrols'])) + cc=np.zeros((mpar['numcontrols'],1)) + ss=np.zeros((mpar['numstates'],1)) + p2=Process(target=FF_2, args=(range_, Xss,Yss,Gamma_state,indexMUdct,indexVKdct,par,mpar, + grid,targets,Copula,P_H,aggrshock,Fb,packagesize, bl,ss,cc, + out_DF2, out_bl2)) + procs2.append(p2) + p2.start() + print('Block number: ', str(bl)) + + FF=[] + order_bl2 = [] + + for i in range(blocks): + FF.append(out_DF2.get()) + order_bl2.append(out_bl2.get()) + + print('bl2 order') + print(order_bl2) + + + for p2 in procs2: + p2.join() + + for i in range(0,int( ceil(mpar['numcontrols'] / float(packagesize)) ) ): + range_=range(i*packagesize, min(packagesize*(i+1),mpar['numcontrols'])) + F2[:,range_]=FF[order_bl2.index(i)] + + end_time = time.perf_counter() + print('Elapsed time is ', (end_time-start_time), ' seconds.') + + FF=[] + FF1=[] + FF3=[] + cc=np.zeros((mpar['numcontrols'],1)) + ss=np.zeros((mpar['numstates'],1)) + + for Yct in range(0, mpar['oc']): + Y=np.zeros((mpar['numcontrols'],1)) + h=par['scaleval2'] + Y[-1-Yct]=h + Fx=F(ss.copy(),ss.copy(),cc.copy(),Y) + F4[:,-1 - Yct]=(Fx['Difference'] - Fb) / h + + F2[mpar['nm']+mpar['nk']-3:mpar['numstates']-2,:] = 0 + + s,t,Q,Z=linalg.qz(np.hstack((F1,F2)), -np.hstack((F3,F4)), output='complex') + abst = abs(np.diag(t))*(abs(np.diag(t))!=0.)+ (abs(np.diag(t))==0.)*10**(-11) + #relev=np.divide(abs(np.diag(s)), abs(np.diag(t))) + relev=np.divide(abs(np.diag(s)), abst) + + ll=sorted(relev) + slt=relev >= 1 + nk=sum(slt) + slt=1*slt + + s_ord,t_ord,__,__,__,Z_ord=linalg.ordqz(np.hstack((F1,F2)), -np.hstack((F3,F4)), sort='ouc', output='complex') + + def _sortOverridEigen(x, y): + out = np.empty_like(x, dtype=bool) + xzero = (x == 0) + yzero = (y == 0) + out[xzero & yzero] = False + out[~xzero & yzero] = True + out[~yzero] = (abs(x[~yzero]/y[~yzero]) > ll[-1 - mpar['numstates']]) + return out + + if nk > mpar['numstates']: + if mpar['overrideEigen']: + print('Warning: The Equilibrium is Locally Indeterminate, critical eigenvalue shifted to: ', str(ll[-1 - mpar['numstates']])) + slt=relev > ll[-1 - mpar['numstates']] + nk=sum(slt) + s_ord,t_ord,__,__,__,Z_ord=linalg.ordqz(np.hstack((F1,F2)), -np.hstack((F3,F4)), sort=_sortOverridEigen, output='complex') + else: + print('No Local Equilibrium Exists, last eigenvalue: ', str(ll[-1 - mpar['numstates']])) + elif nk < mpar['numstates']: + if mpar['overrideEigen']: + print('Warning: No Local Equilibrium Exists, critical eigenvalue shifted to: ', str(ll[-1 - mpar['numstates']])) + slt=relev > ll[-1 - mpar['numstates']] + nk=sum(slt) + s_ord,t_ord,__,__,__,Z_ord=linalg.ordqz(np.hstack((F1,F2)), -np.hstack((F3,F4)), sort=_sortOverridEigen, output='complex') + else: + print('No Local Equilibrium Exists, last eigenvalue: ', str(ll[-1 - mpar['numstates']])) + + z21=Z_ord[nk:,0:nk] + z11=Z_ord[0:nk,0:nk] + s11=s_ord[0:nk,0:nk] + t11=t_ord[0:nk,0:nk] + + if matrix_rank(z11) < nk: + print ('Warning: invertibility condition violated') + + z11i = np.dot(np.linalg.inv(z11), np.eye(nk)) # compute the solution + gx = np.real(np.dot(z21,z11i)) + hx = np.real(np.dot(z11,np.dot(np.dot(np.linalg.inv(s11),t11),z11i))) + return{'hx': hx, 'gx': gx, 'F1': F1, 'F2': F2, 'F3': F3, 'F4': F4, 'par': par } + + +# Sequential implementation of SGU_solver +# see https://github.com/econ-ark/HARK/issues/198 + +# def SGU_solver(Xss,Yss,Gamma_state,indexMUdct,indexVKdct,par,mpar,grid,targets,Copula,P_H,aggrshock): # + +# State = np.zeros((mpar['numstates'],1)) +# State_m = State.copy() +# Contr = np.zeros((mpar['numcontrols'],1)) +# Contr_m = Contr.copy() + + +# F = lambda S, S_m, C, C_m : Fsys(S, S_m, C, C_m, +# Xss,Yss,Gamma_state,indexMUdct,indexVKdct, +# par,mpar,grid,targets,Copula,P_H,aggrshock) + + +# start_time = time.clock() +# result_F = F(State,State_m,Contr.copy(),Contr_m.copy()) +# end_time = time.clock() +# print ('Elapsed time is ', (end_time-start_time), ' seconds.') +# Fb=result_F['Difference'].copy() + +# pool=cpu_count()/2 + +# F1=np.zeros((mpar['numstates'] + mpar['numcontrols'], mpar['numstates'])) +# F2=np.zeros((mpar['numstates'] + mpar['numcontrols'], mpar['numcontrols'])) +# F3=np.zeros((mpar['numstates'] + mpar['numcontrols'], mpar['numstates'])) +# F4=np.asmatrix(np.vstack((np.zeros((mpar['numstates'], mpar['numcontrols'])), np.eye(mpar['numcontrols'],mpar['numcontrols']) ))) + +# print ('Use Schmitt Grohe Uribe Algorithm') +# print (' A *E[xprime uprime] =B*[x u]') +# print (' A = (dF/dxprimek dF/duprime), B =-(dF/dx dF/du)') + +# #numscale=1 +# pnum=pool +# packagesize=int(ceil(mpar['numstates'] / float(3*pnum))) +# blocks=int(ceil(mpar['numstates'] / float(packagesize) )) + +# par['scaleval1'] = 1e-5 +# par['scaleval2'] = 1e-5 + +# start_time = time.clock() +# print ('Computing Jacobian F1=DF/DXprime F3 =DF/DX') +# print ('Total number of parallel blocks: ', str(blocks), '.') + +# FF1=[] +# FF3=[] + +# for bl in range(0,blocks): +# range_= range(bl*packagesize, min(packagesize*(bl+1),mpar['numstates'])) +# DF1=np.asmatrix( np.zeros((len(Fb),len(range_))) ) +# DF3=np.asmatrix( np.zeros((len(Fb),len(range_))) ) +# cc=np.zeros((mpar['numcontrols'],1)) +# ss=np.zeros((mpar['numstates'],1)) +# for Xct in range_: +# X=np.zeros((mpar['numstates'],1)) +# h=par['scaleval1'] +# X[Xct]=h +# Fx=F(ss.copy(),X,cc.copy(),cc.copy()) +# DF3[:, Xct - bl*packagesize]=(Fx['Difference'] - Fb) / h +# Fx=F(X,ss.copy(),cc.copy(),cc.copy()) +# DF1[:, Xct - bl*packagesize]=(Fx['Difference'] - Fb) / h +# if sum(range_ == mpar['numstates'] - 2) == 1: +# Xct=mpar['numstates'] - 2 +# X=np.zeros((mpar['numstates'],1)) +# h=par['scaleval2'] +# X[Xct]=h +# Fx=F(ss.copy(),X,cc.copy(),cc.copy()) +# DF3[:,Xct - bl*packagesize]=(Fx['Difference'] - Fb) / h +# Fx=F(X,ss.copy(),cc.copy(),cc.copy()) +# DF1[:,Xct - bl*packagesize]=(Fx['Difference'] - Fb) / h +# if sum(range_ == mpar['numstates'] - 1) == 1: +# Xct=mpar['numstates'] - 1 +# X=np.zeros((mpar['numstates'],1)) +# h=par['scaleval2'] +# X[Xct]=h +# Fx=F(ss.copy(),X,cc.copy(),cc.copy()) +# DF3[:,Xct - bl*packagesize]=(Fx['Difference'] - Fb) / h +# Fx=F(X,ss.copy(),cc.copy(),cc.copy()) +# DF1[:,Xct - bl*packagesize]=(Fx['Difference'] - Fb) / h +# FF1.append(DF1.copy()) +# FF3.append(DF3.copy()) +# print ('Block number: ', str(bl),' done.') + +# for i in range(0,int(ceil(mpar['numstates'] / float(packagesize)) )): +# range_= range(i*packagesize, min(packagesize*(i+1),mpar['numstates'])) +# F1[:,range_]=FF1[i] +# F3[:,range_]=FF3[i] + +# end_time = time.clock() +# print ('Elapsed time is ', (end_time-start_time), ' seconds.') + +# # jacobian wrt Y' +# packagesize=int(ceil(mpar['numcontrols'] / (3.0*pnum))) +# blocks=int(ceil(mpar['numcontrols'] / float(packagesize))) +# print ('Computing Jacobian F2 - DF/DYprime') +# print ('Total number of parallel blocks: ', str(blocks),'.') + +# FF=[] + +# start_time = time.clock() + +# for bl in range(0,blocks): +# range_= range(bl*packagesize,min(packagesize*(bl+1),mpar['numcontrols'])) +# DF2=np.asmatrix(np.zeros((len(Fb),len(range_)))) +# cc=np.zeros((mpar['numcontrols'],1)) +# ss=np.zeros((mpar['numstates'],1)) +# for Yct in range_: +# Y=np.zeros((mpar['numcontrols'],1)) +# h=par['scaleval2'] +# Y[Yct]=h +# Fx=F(ss.copy(),ss.copy(),Y,cc.copy()) +# DF2[:,Yct - bl*packagesize]=(Fx['Difference'] - Fb) / h +# FF.append(DF2.copy()) +# print ('Block number: ',str(bl),' done.') + + +# for i in range(0,int(ceil(mpar['numcontrols'] / float(packagesize) ))): +# range_=range(i*packagesize, min(packagesize*(i+1),mpar['numcontrols'])) +# F2[:,range_]=FF[i] + +# end_time = time.clock() +# print ('Elapsed time is ', (end_time-start_time), ' seconds.') + + +# FF=[] +# FF1=[] +# FF3=[] + +# cc=np.zeros((mpar['numcontrols'],1)) +# ss=np.zeros((mpar['numstates'],1)) + +# for Yct in range(0, mpar['oc']): +# Y=np.zeros((mpar['numcontrols'],1)) +# h=par['scaleval2'] +# Y[-1-Yct]=h +# Fx=F(ss.copy(),ss.copy(),cc.copy(),Y) +# F4[:,-1 - Yct]=(Fx['Difference'] - Fb) / h + +# F2[mpar['nm']+mpar['nk']-3:mpar['numstates']-2,:] = 0 + +# s,t,Q,Z=linalg.qz(np.hstack((F1,F2)), -np.hstack((F3,F4)), output='complex') +# abst = abs(np.diag(t))*(abs(np.diag(t))!=0.)+ (abs(np.diag(t))==0.)*10**(-11) +# #relev=np.divide(abs(np.diag(s)), abs(np.diag(t))) +# relev=np.divide(abs(np.diag(s)), abst) + +# ll=sorted(relev) +# slt=relev >= 1 +# nk=sum(slt) +# slt=1*slt + + +# s_ord,t_ord,__,__,__,Z_ord=linalg.ordqz(np.hstack((F1,F2)), -np.hstack((F3,F4)), sort='ouc', output='complex') + +# def sortOverridEigen(x, y): +# out = np.empty_like(x, dtype=bool) +# xzero = (x == 0) +# yzero = (y == 0) +# out[xzero & yzero] = False +# out[~xzero & yzero] = True +# out[~yzero] = (abs(x[~yzero]/y[~yzero]) > ll[-1 - mpar['numstates']]) +# return out + +# if nk > mpar['numstates']: +# if mpar['overrideEigen']: +# print ('Warning: The Equilibrium is Locally Indeterminate, critical eigenvalue shifted to: ', str(ll[-1 - mpar['numstates']])) +# slt=relev > ll[-1 - mpar['numstates']] +# nk=sum(slt) +# s_ord,t_ord,__,__,__,Z_ord=linalg.ordqz(np.hstack((F1,F2)), -np.hstack((F3,F4)), sort=sortOverridEigen, output='complex') + +# else: +# print ('No Local Equilibrium Exists, last eigenvalue: ', str(ll[-1 - mpar['numstates']])) - F = lambda S, S_m, C, C_m : Fsys(S, S_m, C, C_m, - Xss,Yss,Gamma_state,indexMUdct,indexVKdct, - par,mpar,grid,targets,Copula,P_H,aggrshock) - - - start_time = time.clock() - result_F = F(State,State_m,Contr.copy(),Contr_m.copy()) - end_time = time.clock() - print ('Elapsed time is ', (end_time-start_time), ' seconds.') - Fb=result_F['Difference'].copy() - - pool=cpu_count()/2 +# elif nk < mpar['numstates']: +# if mpar['overrideEigen']: +# print ('Warning: No Local Equilibrium Exists, critical eigenvalue shifted to: ', str(ll[-1 - mpar['numstates']])) +# slt=relev > ll[-1 - mpar['numstates']] +# nk=sum(slt) +# s_ord,t_ord,__,__,__,Z_ord=linalg.ordqz(np.hstack((F1,F2)), -np.hstack((F3,F4)), sort=sortOverridEigen, output='complex') - F1=np.zeros((mpar['numstates'] + mpar['numcontrols'], mpar['numstates'])) - F2=np.zeros((mpar['numstates'] + mpar['numcontrols'], mpar['numcontrols'])) - F3=np.zeros((mpar['numstates'] + mpar['numcontrols'], mpar['numstates'])) - F4=np.asmatrix(np.vstack((np.zeros((mpar['numstates'], mpar['numcontrols'])), np.eye(mpar['numcontrols'],mpar['numcontrols']) ))) - - print ('Use Schmitt Grohe Uribe Algorithm') - print (' A *E[xprime uprime] =B*[x u]') - print (' A = (dF/dxprimek dF/duprime), B =-(dF/dx dF/du)') - - #numscale=1 - pnum=pool - packagesize=int(ceil(mpar['numstates'] / float(3*pnum))) - blocks=int(ceil(mpar['numstates'] / float(packagesize) )) +# else: +# print ('No Local Equilibrium Exists, last eigenvalue: ', str(ll[-1 - mpar['numstates']])) - par['scaleval1'] = 1e-5 - par['scaleval2'] = 1e-5 - - start_time = time.clock() - print ('Computing Jacobian F1=DF/DXprime F3 =DF/DX') - print ('Total number of parallel blocks: ', str(blocks), '.') - - FF1=[] - FF3=[] - - for bl in range(0,blocks): - range_= range(bl*packagesize, min(packagesize*(bl+1),mpar['numstates'])) - DF1=np.asmatrix( np.zeros((len(Fb),len(range_))) ) - DF3=np.asmatrix( np.zeros((len(Fb),len(range_))) ) - cc=np.zeros((mpar['numcontrols'],1)) - ss=np.zeros((mpar['numstates'],1)) - for Xct in range_: - X=np.zeros((mpar['numstates'],1)) - h=par['scaleval1'] - X[Xct]=h - Fx=F(ss.copy(),X,cc.copy(),cc.copy()) - DF3[:, Xct - bl*packagesize]=(Fx['Difference'] - Fb) / h - Fx=F(X,ss.copy(),cc.copy(),cc.copy()) - DF1[:, Xct - bl*packagesize]=(Fx['Difference'] - Fb) / h - if sum(range_ == mpar['numstates'] - 2) == 1: - Xct=mpar['numstates'] - 2 - X=np.zeros((mpar['numstates'],1)) - h=par['scaleval2'] - X[Xct]=h - Fx=F(ss.copy(),X,cc.copy(),cc.copy()) - DF3[:,Xct - bl*packagesize]=(Fx['Difference'] - Fb) / h - Fx=F(X,ss.copy(),cc.copy(),cc.copy()) - DF1[:,Xct - bl*packagesize]=(Fx['Difference'] - Fb) / h - if sum(range_ == mpar['numstates'] - 1) == 1: - Xct=mpar['numstates'] - 1 - X=np.zeros((mpar['numstates'],1)) - h=par['scaleval2'] - X[Xct]=h - Fx=F(ss.copy(),X,cc.copy(),cc.copy()) - DF3[:,Xct - bl*packagesize]=(Fx['Difference'] - Fb) / h - Fx=F(X,ss.copy(),cc.copy(),cc.copy()) - DF1[:,Xct - bl*packagesize]=(Fx['Difference'] - Fb) / h - FF1.append(DF1.copy()) - FF3.append(DF3.copy()) - print ('Block number: ', str(bl),' done.') - - for i in range(0,int(ceil(mpar['numstates'] / float(packagesize)) )): - range_= range(i*packagesize, min(packagesize*(i+1),mpar['numstates'])) - F1[:,range_]=FF1[i] - F3[:,range_]=FF3[i] - - end_time = time.clock() - print ('Elapsed time is ', (end_time-start_time), ' seconds.') - - # jacobian wrt Y' - packagesize=int(ceil(mpar['numcontrols'] / (3.0*pnum))) - blocks=int(ceil(mpar['numcontrols'] / float(packagesize))) - print ('Computing Jacobian F2 - DF/DYprime') - print ('Total number of parallel blocks: ', str(blocks),'.') - - FF=[] - - start_time = time.clock() - - for bl in range(0,blocks): - range_= range(bl*packagesize,min(packagesize*(bl+1),mpar['numcontrols'])) - DF2=np.asmatrix(np.zeros((len(Fb),len(range_)))) - cc=np.zeros((mpar['numcontrols'],1)) - ss=np.zeros((mpar['numstates'],1)) - for Yct in range_: - Y=np.zeros((mpar['numcontrols'],1)) - h=par['scaleval2'] - Y[Yct]=h - Fx=F(ss.copy(),ss.copy(),Y,cc.copy()) - DF2[:,Yct - bl*packagesize]=(Fx['Difference'] - Fb) / h - FF.append(DF2.copy()) - print ('Block number: ',str(bl),' done.') - - for i in range(0,int(ceil(mpar['numcontrols'] / float(packagesize) ))): - range_=range(i*packagesize, min(packagesize*(i+1),mpar['numcontrols'])) - F2[:,range_]=FF[i] - - end_time = time.clock() - print ('Elapsed time is ', (end_time-start_time), ' seconds.') - - - FF=[] - FF1=[] - FF3=[] - - cc=np.zeros((mpar['numcontrols'],1)) - ss=np.zeros((mpar['numstates'],1)) - - for Yct in range(0, mpar['oc']): - Y=np.zeros((mpar['numcontrols'],1)) - h=par['scaleval2'] - Y[-1-Yct]=h - Fx=F(ss.copy(),ss.copy(),cc.copy(),Y) - F4[:,-1 - Yct]=(Fx['Difference'] - Fb) / h - - F2[mpar['nm']+mpar['nk']-3:mpar['numstates']-2,:] = 0 - - s,t,Q,Z=linalg.qz(np.hstack((F1,F2)), -np.hstack((F3,F4)), output='complex') - abst = abs(np.diag(t))*(abs(np.diag(t))!=0.)+ (abs(np.diag(t))==0.)*10**(-11) - #relev=np.divide(abs(np.diag(s)), abs(np.diag(t))) - relev=np.divide(abs(np.diag(s)), abst) - - ll=sorted(relev) - slt=relev >= 1 - nk=sum(slt) - slt=1*slt - - - s_ord,t_ord,__,__,__,Z_ord=linalg.ordqz(np.hstack((F1,F2)), -np.hstack((F3,F4)), sort='ouc', output='complex') - - def sortOverridEigen(x, y): - out = np.empty_like(x, dtype=bool) - xzero = (x == 0) - yzero = (y == 0) - out[xzero & yzero] = False - out[~xzero & yzero] = True - out[~yzero] = (abs(x[~yzero]/y[~yzero]) > ll[-1 - mpar['numstates']]) - return out - - if nk > mpar['numstates']: - if mpar['overrideEigen']: - print ('Warning: The Equilibrium is Locally Indeterminate, critical eigenvalue shifted to: ', str(ll[-1 - mpar['numstates']])) - slt=relev > ll[-1 - mpar['numstates']] - nk=sum(slt) - s_ord,t_ord,__,__,__,Z_ord=linalg.ordqz(np.hstack((F1,F2)), -np.hstack((F3,F4)), sort=sortOverridEigen, output='complex') - - else: - print ('No Local Equilibrium Exists, last eigenvalue: ', str(ll[-1 - mpar['numstates']])) - - elif nk < mpar['numstates']: - if mpar['overrideEigen']: - print ('Warning: No Local Equilibrium Exists, critical eigenvalue shifted to: ', str(ll[-1 - mpar['numstates']])) - slt=relev > ll[-1 - mpar['numstates']] - nk=sum(slt) - s_ord,t_ord,__,__,__,Z_ord=linalg.ordqz(np.hstack((F1,F2)), -np.hstack((F3,F4)), sort=sortOverridEigen, output='complex') - - else: - print ('No Local Equilibrium Exists, last eigenvalue: ', str(ll[-1 - mpar['numstates']])) - - - z21=Z_ord[nk:,0:nk] - z11=Z_ord[0:nk,0:nk] - s11=s_ord[0:nk,0:nk] - t11=t_ord[0:nk,0:nk] - - if matrix_rank(z11) < nk: - print ('Warning: invertibility condition violated') +# z21=Z_ord[nk:,0:nk] +# z11=Z_ord[0:nk,0:nk] +# s11=s_ord[0:nk,0:nk] +# t11=t_ord[0:nk,0:nk] - z11i = np.dot(np.linalg.inv(z11), np.eye(nk)) # compute the solution +# if matrix_rank(z11) < nk: +# print ('Warning: invertibility condition violated') - gx = np.real(np.dot(z21,z11i)) - hx = np.real(np.dot(z11,np.dot(np.dot(np.linalg.inv(s11),t11),z11i))) - - return{'hx': hx, 'gx': gx, 'F1': F1, 'F2': F2, 'F3': F3, 'F4': F4, 'par': par } +# z11i = np.dot(np.linalg.inv(z11), np.eye(nk)) # compute the solution + +# gx = np.real(np.dot(z21,z11i)) +# hx = np.real(np.dot(z11,np.dot(np.dot(np.linalg.inv(s11),t11),z11i))) + +# return{'hx': hx, 'gx': gx, 'F1': F1, 'F2': F2, 'F3': F3, 'F4': F4, 'par': par } ############################################################################### @@ -1427,7 +1413,7 @@ def sortOverridEigen(x, y): EX3SS=pickle.load(open("EX3SS_20.p", "rb")) - start_time0 = time.clock() + start_time0 = time.perf_counter() ## Aggregate shock to perturb(one of three shocks: MP, TFP, Uncertainty) @@ -1456,7 +1442,7 @@ def sortOverridEigen(x, y): plot_IRF(SR['mpar'],SR['par'],SGUresult['gx'],SGUresult['hx'],SR['joint_distr'], SR['Gamma_state'],SR['grid'],SR['targets'],SR['Output']) - end_time0 = time.clock() + end_time0 = time.perf_counter() print('Elapsed time is ', (end_time0-start_time0), ' seconds.') # \ No newline at end of file