-
Notifications
You must be signed in to change notification settings - Fork 0
/
plotPensions.py
165 lines (133 loc) · 4.79 KB
/
plotPensions.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
"""
Computes the Euler errors and K regions for the benchmark pension savings
model by Druedhal and Jorgensen (2017) using RFC and compares to G2EGM.
Author: Akshay Shanker, University of New South Wales, Sydney
Email: akshay.shanker@me.com
Acknowledgements:
-----------------
The pensions module (G2EGM) is a modification
of the original G2EGM mod written by Druedhal and Jorgensen. The code original
was modified to:
i) invert the Euler equation according to the necessary conditions
in Dobrescu and Shanker (2024) and
ii) implement roof-top cut algorithm.
"""
import numpy as np
import dill as pickle
from G2EGM.G2EGMModel import G2EGMModelClass
from G2EGM.figs import euler_hist, Kregions, decision_functions, segments
from timingPensions import timing
# Suppress all numpy errors (overflows, divisions by zero, etc.)
np.seterr(all='ignore')
if __name__ == '__main__':
# Simulation parameters
Nm = 800 # Grid points for each axis
T = 20 # Periods
k = 65 # Nearest neighbors points for RFC search
s = 0.08 # Proportion of the grid to be used in the RFC at each iteration of RFC
rho_r = 0.33 # Max radius for RFC to eliminate points
rho_rI = 0.475 # Radius for RFC to search for intersections
J_bar = 1 + 1e-5 # Jump detection threshold
k1 = 30 # Neighbors for intersection point search
k2 = 1 # Neighbors of uniform grid to construct triangulation for interpolation
segplot_t = 3 # t for plotting constrained regions
do_print = True
p_L = 2 # Upper bound for pension contributions
# File paths for saving plots and data
scrpath = '/scratch/tp66/dcdp/data/' # drive where raw egmgrids are saved
plotpath = 'plots/pensions/'
# Initialize and configure the RFC model
model_RFC = G2EGMModelClass(
name='RFC',
par={
'solmethod': 'RFC',
'T': T,
'do_print': do_print,
'k': k,
'Nm': Nm,
'rad': rho_r,
'rad_I': rho_rI,
'J_bar': J_bar,
'k1': k1,
'k2': k2,
'intersection': False,
'interp_intersect': False,
's': s,
'correct_jumps': True,
't_save': segplot_t,
'do_print': do_print,
'p_L': p_L,
'save_data': False
}
)
#
model_RFC.precompile_numba()
model_RFC = timing(model_RFC, rep=1)
#segments(model_RFC, 2,'test')
# Initialize and configure the G2EGM model
model_G2EGM = G2EGMModelClass(
name='G2EGM',
par={
'solmethod': 'G2EGM',
'T': T,
'do_print': do_print,
'Nm': Nm,
'p_L': p_L,
}
)
model_G2EGM.precompile_numba()
model_G2EGM = timing(model_G2EGM, rep=1)
#segments(model_G2EGM, 3,'testG2E')
# Load endogenous grid data
egrids_intersect = pickle.load(open(f"{scrpath}e_grids_intersect.pkl", "rb"))
egrids_clean = pickle.load(open(f"{scrpath}/e_grids_clean.pkl", "rb"))
egrids_raw = pickle.load(open(f"{scrpath}/e_grid_raw.pkl", "rb"))
# Create a histogram of the Euler errors
#smodelDict = {'RFC': model_RFC, 'G2EGM': model_G2EGM}
#euler_hist(modelDict, Nm, plotpath)
# Plot K regions
Kregions(model_RFC, egrids_raw, egrids_clean, egrids_intersect, plotpath, segplot_t)
decision_functions(model_RFC,3,'RFC_t3')
models = [model_RFC, model_G2EGM]
# tab
postfix = '_G2EGM_vs_RFC'
# b. euler errors
lines = []
txt = '| All (average)'
for i,model in enumerate(models):
txt += f' | {np.nanmean(model.sim.euler):.3f}'
txt += ' |\n'
lines.append(txt)
txt = '| 5th percentile'
for i,model in enumerate(models):
txt += f' | {np.nanpercentile(model.sim.euler,5):.3f}'
txt += ' |\n'
lines.append(txt)
txt = '| 95th percentile'
for i,model in enumerate(models):
txt += f' | {np.nanpercentile(model.sim.euler,95):.3f}'
txt += ' |\n'
lines.append(txt)
txt = '| Median'
for i,model in enumerate(models):
txt += f' | {np.nanpercentile(model.sim.euler,50):.3f}'
txt += ' |\n'
lines.append(txt)
# c. timings
txt = '| Total time (min)'
for model in models:
txt += f' | {np.sum(model.par.time_work)/60:.2f}'
txt += ' |\n'
lines.append(txt)
txt = '| Inversion time (min)'
for model in models:
txt += f' | {np.sum(model.par.time_invert)/60:.2f}'
txt += ' |\n'
lines.append(txt)
txt = '| RFC time (min)'
for model in models:
txt += f' | {np.sum(model.par.time_rfc)/60:.2f}'
txt += ' |\n'
lines.append(txt)
with open(f'tabs_euler_errors{postfix}.md', 'w') as txtfile:
txtfile.writelines(lines)