-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathincome.py
198 lines (186 loc) · 6.74 KB
/
income.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
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
import numpy as np
import scipy.optimize as opt
import scipy.interpolate as si
from ogcore import parameter_plots as pp
from ogcore import utils
from ogcore.parameters import Specifications
import os
import json
import urllib.request
CUR_PATH = os.path.abspath(os.path.dirname(__file__))
OUTPUT_DIR = os.path.join(CUR_PATH, "OUTPUT", "ability")
def get_e_interp(E, S, J, lambdas, age_wgts, gini_to_match=41.1, plot=False):
"""
This function takes the calibrated lifetime earnings profiles
(abilities, e matrix) from OG-USA and then adjusts the shape of those
profiles to match the Gini coefficient for another economy. The
Gini coefficient to match is given in the argument gini_to_match.
Note that the calibrated OG-USA e matrix is of size (80, 10), where
80 is the number of ages and 10 is the number of ability types.
Users of this function specify their own number of age groups (S)
and ability types (J). The function will map the fitted functions
into these dimensions so long as the percentiles of the ability types
given in lambdas is not more refined at the top end than those in
OG-USA (which identifies up to the top 0.1%).
Args:
E (int): the age agents become economically active
S (int): number of ages to interpolate. This method assumes that
ages are evenly spaced between the beginning of age E
up to E+S, >= 3
J (int): number of ability types to interpolate
lambdas (Numpy array): distribution of population in each
ability group, length J
age_wgts (Numpy array): distribution of population in each age
group, length S
gini_to_match (float): Gini coefficient to match, default is
41.1, the Gini coefficient for MYS in 2019
plot (bool): if True, creates plots of emat_orig and the new
interpolated emat_new
Returns:
emat_new_scaled (Numpy array): interpolated ability matrix scaled
so that population-weighted average is 1, size SxJ
"""
assert lambdas.shape[0] == J
assert age_wgts.shape[0] == S
# Load USA e matrix as a baseline
usa_params = Specifications()
usa_params.update_specifications(
json.load(
urllib.request.urlopen(
"https://raw.githubusercontent.com/PSLmodels/OG-USA/master/ogusa/ogusa_default_parameters.json"
)
)
)
# Define a function that will find the "a" in the equation:
# e_Y = e_USA * exp(a * e_USA)
# such that the e_Y produces a gini coefficient in the model that
# gives the same ratio between the model implied Gini's in the USA
# and the target country and the empirical Gini's in the USA and given
# by gin_to_match for the target country
def f(
a,
emat_orig,
age_wgts,
abil_wgts,
gini_to_match,
gini_usa_data,
gini_usa_model,
):
gini_target_model = utils.Inequality(
emat_orig * np.exp(a * emat_orig),
age_wgts,
abil_wgts,
len(age_wgts),
len(abil_wgts),
).gini()
error = (gini_to_match / gini_usa_data) - (
gini_target_model / gini_usa_model
)
return error
# Note, USA gini in the World Bank data is 41.5
# See https://data.worldbank.org/indicator/SI.POV.GINI
gini_usa_data = 41.5
# Find the model implied Gini for the USA
gini_usa_model = utils.Inequality(
usa_params.e[0, :, :],
usa_params.omega_SS,
usa_params.lambdas,
usa_params.S,
usa_params.J,
).gini()
x = opt.root_scalar(
f,
args=(
usa_params.e[0, :, :],
usa_params.omega_SS,
usa_params.lambdas,
gini_to_match,
gini_usa_data,
gini_usa_model,
),
method="bisect",
bracket=[-1, 1],
xtol=1e-10,
)
a = x.root
e_new = usa_params.e[0, :, :] * np.exp(a * usa_params.e[0, :, :])
emat_new_scaled = (
e_new
/ (
e_new
* usa_params.omega_SS.reshape(usa_params.S, 1)
* usa_params.lambdas.reshape(1, usa_params.J)
).sum()
)
# Now interpolate for the cases where S and/or J not the same in the
# country parameterization as in the default USA parameterization
if (
S == usa_params.S
and np.array_equal(
usa_params.lambdas,
lambdas,
)
is True
):
pass # will return the e_new_scaled found above since dims the same
else:
# generate vector of mid points for the Malaysian ability groups
abil_midp = np.zeros(J)
pct_lb = 0.0
for j in range(J):
abil_midp[j] = pct_lb + 0.5 * lambdas[j]
pct_lb += lambdas[j]
# generate vector of mid points for the USA ability groups
M = usa_params.lambdas.shape[0]
emat_j_midp = np.zeros(M)
pct_lb = 0.0
for m in range(M):
emat_j_midp[m] = pct_lb + 0.5 * usa_params.lambdas[m]
pct_lb += usa_params.lambdas[m]
# Make sure that values in abil_midp are within interpolating
# bounds
if abil_midp.min() < emat_j_midp.min() or abil_midp.max() > (
1 - usa_params.lambdas[-1]
):
err = (
"One or more entries in abilities vector (lambdas) is outside the "
+ "allowable bounds for interpolation."
)
raise RuntimeError(err)
usa_step = 80 / usa_params.S
emat_s_midp = np.linspace(
usa_params.E + 0.5 * usa_step,
usa_params.E + usa_params.S - 0.5 * usa_step,
usa_params.S,
)
emat_j_mesh, emat_s_mesh = np.meshgrid(emat_j_midp, emat_s_midp)
newstep = 80 / S
new_s_midp = np.linspace(E + 0.5 * newstep, E + S - 0.5 * newstep, S)
new_j_mesh, new_s_mesh = np.meshgrid(abil_midp, new_s_midp)
newcoords = np.hstack(
(
emat_s_mesh.reshape((usa_params.S * usa_params.J, 1)),
emat_j_mesh.reshape((usa_params.S * usa_params.J, 1)),
)
)
emat_new = si.griddata(
newcoords,
emat_new_scaled.flatten(),
(new_s_mesh, new_j_mesh),
method="linear",
)
emat_new_scaled = (
emat_new
/ (emat_new * age_wgts.reshape(S, 1) * lambdas.reshape(1, J)).sum()
)
if plot:
kwargs = {"filesuffix": "_intrp_scaled"}
pp.plot_income_data(
new_s_midp,
abil_midp,
lambdas,
emat_new_scaled,
OUTPUT_DIR,
**kwargs,
)
return emat_new_scaled