-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathps1_felipe
464 lines (464 loc) · 99.2 KB
/
ps1_felipe
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
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Assignment 1 - Econ262A"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Author: David Henning"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import statsmodels.api as sm\n",
"import linearmodels as ln\n",
"import pandas as pd\n",
"from statsmodels.sandbox.regression.gmm import IV2SLS\n",
"from pystout import pystout\n",
"from tabulate import tabulate\n",
"import matplotlib.pyplot as plt\n",
"from stargazer.stargazer import Stargazer\n",
"from sklearn.neighbors import KernelDensity\n",
"from sklearn.model_selection import GridSearchCV"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 1"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"#Simulate data\n",
"np.random.seed(1)\n",
"X = np.random.normal(loc=0, scale=1, size=10000)\n",
"Z = np.random.randint(2, size=10000)\n",
"u, ε = np.random.multivariate_normal([0, 0], [[1, 0.25], [0.25, 1]], (10000)).T\n",
"\n",
"T = np.empty(10000)\n",
"for i in range(len(T)):\n",
" T[i] = 1 if -0.5 + 0.25 * X[i] + Z[i] + u[i] > 0 else 0\n",
"\n",
"Y = 0.5 * X + T + ε\n",
"\n",
"data = np.vstack((Y, X, T, Z, u, ε)).T"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Question 1.1"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"data = np.vstack((Y, X, T, Z, u, ε)).T\n",
"df = pd.DataFrame(data=data, columns=[\"Y\", \"X\", \"T\", \"Z\", \"u\", \"ε\"])"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"#1) OLS regression of Y on X and T\n",
"x = df[['X','T']]\n",
"y = df[['Y']]\n",
"m1 = sm.OLS(y, x).fit()\n",
"\n",
"#2) OLS regression of T on Z and X\n",
"x = df[['X','Z']]\n",
"y = df[['T']]\n",
"m2 = sm.OLS(y, x).fit()\n",
"\n",
"#2) OLS regression of T on Z and X\n",
"x = df[['X','T']]\n",
"z = df[['X','Z']]\n",
"y = df[['Y']]\n",
"m3 = IV2SLS(Y, x, instrument=z).fit()\n",
"\n",
"#stargazer = Stargazer([m1, m2, m3])"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"#Produce table with regression results\n",
"pystout(models=[m1, m2, m3],\n",
" file='simulated_iv1.tex',\n",
" digits=3,\n",
" endog_names=['Y','T','Y'],\n",
" exogvars=['T', 'Z', 'X'],\n",
" #varlabels={'price':'Price','prom':'Promtion'},\n",
" stars=False,\n",
" mgroups={'OLS':[1,2], 'IV':3},\n",
" modstat={'nobs':'Obs','rsquared_adj':'Adj. R\\sym{2}'}\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Question 1.2"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"α, u, ε = np.random.multivariate_normal([1, 0, 0], [[1, 0.25, 0.25], [0.25, 1, 0.25], [0.25, 0.25, 1]], (10000)).T\n",
"\n",
"T = np.empty(10000)\n",
"for i in range(len(T)):\n",
" T[i] = 1 if -0.5 + 0.25 * X[i] + Z[i] + u[i] > 0 else 0\n",
"\n",
"Y = 0.5 * X + α * T + ε"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"data = np.vstack((Y, X, T, Z, u, ε, α)).T\n",
"df = pd.DataFrame(data=data, columns=[\"Y\", \"X\", \"T\", \"Z\", \"u\", \"ε\", \"α\"])"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"#1) OLS regression of Y on X and T\n",
"x = df[['X','T']]\n",
"y = df[['Y']]\n",
"m1 = sm.OLS(y, x).fit()\n",
"\n",
"#2) OLS regression of T on Z and X\n",
"x = df[['X','Z']]\n",
"y = df[['T']]\n",
"m2 = sm.OLS(y, x).fit()\n",
"\n",
"#2) OLS regression of T on Z and X\n",
"x = df[['X','T']]\n",
"z = df[['X','Z']]\n",
"y = df[['Y']]\n",
"m3 = IV2SLS(Y, x, instrument=z).fit()\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"#Produce table with regression results\n",
"pystout(models=[m1, m2, m3],\n",
" file='simulated_iv2.tex',\n",
" digits=3,\n",
" endog_names=['Y','T','Y'],\n",
" exogvars=['T', 'Z', 'X'],\n",
" #varlabels={'price':'Price','prom':'Promtion'},\n",
" stars=False,\n",
" mgroups={'OLS':[1,2], 'IV':3},\n",
" modstat={'nobs':'Obs','rsquared_adj':'Adj. R\\sym{2}'}\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot the dsitribution of $Y$ for always takers and never takers."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"a = df[(df['Z']==0) & (df['T']==1)] #Define always takers\n",
"n = df[(df['Z']==1) & (df['T']==0)] #Define never takers\n",
"c = df[(df['Z'] == df['T'])] #Define compliers\n",
"N = df.shape[0]\n",
"\n",
"φ_n = n.shape[0]/N\n",
"φ_a = a.shape[0]/N\n",
"φ_c = 1 - φ_n - φ_a"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"c0 = df[(df['Z']==0) & (df['T']==0)]\n",
"y_c0 = c0['Y']\n",
"c1 = df[(df['Z']==1) & (df['T']==1)]\n",
"y_c1 = c1['Y']\n",
"\n",
"y_n = n['Y']\n",
"y_a = a['Y']"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"distributions = [y_a, y_n, y_c0, y_c1]\n",
"y_min = min(df['Y'])\n",
"y_max = max(df['Y'])\n",
"mylist = []\n",
"\n",
"for i in distributions:\n",
" y_d = np.linspace(y_min, y_max, 2000)\n",
" kde = KernelDensity(bandwidth=0.1, kernel='gaussian')\n",
" kde.fit(i[:, None])\n",
"\n",
" # score_samples returns the log of the probability density\n",
" y = np.exp(kde.score_samples(y_d[:, None]))\n",
" mylist.append(y)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"g_a = mylist[0]\n",
"g_n = mylist[1]\n",
"\n",
"g_c0 = ((φ_n + φ_c) / φ_c) * mylist[2] - (φ_n / φ_c) * g_n\n",
"g_c1 = ((φ_a + φ_c) / φ_c) * mylist[3] - (φ_a / φ_c) * g_a"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD7CAYAAABkO19ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydd3Sc1bW3nzNdM+rVkiVbsi03ybJs3MCFYhA1EAMXTAkthBBCbgqkQMplcROSexcXCPeGsCBAki90OwZCCWBCcQE3bGy5y7Ks3ttoRqNp5/vjHY0kaySNZFX7PGt5SXPec+bdI0u/2bPPPnsLKSUKhUKhOH3RjbUBCoVCoRhZlNArFArFaY4SeoVCoTjNUUKvUCgUpzlK6BUKheI0Rwm9QqFQnOaEJfRCiEuEEIeFEEVCiJ/1M2+xEMInhLh2sGsVCoVCMTIMKPRCCD3wB+BSYC5wgxBibh/z/gt4f7BrFQqFQjFyGMKYswQoklIWAwghXgGuAg6cNO97wHpg8RDW9iAxMVFmZmaGY79CoVAogF27dtVLKZNCXQtH6CcDZd0elwNLu08QQkwG1gAX0FPoB1wbiszMTHbu3BmGaQqFQqEAEEKc6OtaODF6EWLs5LoJTwA/lVL6hrBWmyjEXUKInUKInXV1dWGYpVAoFIpwCMejLwcyuj1OBypPmrMIeEUIAZAIXCaE8Ia5FgAp5TPAMwCLFi1SBXgUCoVimAhH6HcA2UKILKACWAvc2H2ClDKr83shxJ+Bt6WUbwghDAOtVSgUCsXIMqDQSym9Qoh70bJp9MDzUsr9Qoi7A9efHuza4TFdoVCMBB6Ph/Lyclwu11ibogiBxWIhPT0do9EY9hoxHssUL1q0SKrNWIVibDh+/DhRUVEkJCQQCMcqxglSShoaGrDb7WRlZfW4JoTYJaVcFGqdOhmrUCh64HK5lMiPU4QQJCQkDPrTlhJ6hULRCyXy45eh/N8ooVf0SYfXx56yZlqcnrE2RaFQnAJK6BV98snhOj4+VMu6L8vx+PxjbY7iDEIIwX333Rd8/Oijj/LQQw+N2v1LSkp46aWXBpz3ySefcMUVV4yCRaeGEnpFSDq8Po5U24m1Gmlt93C42j7WJinOIMxmM3//+9+pr68flft5vd4ej8MV+uG850iihF4RkuoWF16/5ILZycTbTByobB2W5/X4/IzHTC/F+MJgMHDXXXfx+OOP97pWV1fHNddcw+LFi1m8eDFbtmzB7/eTmZlJc3NzcN6MGTOoqakJOR/goYce4q677qKgoIBbbrmlxz1+9rOfsWnTJvLz83n88ccpKSlh5cqVLFy4kIULF7J169Zedu3YsYMFCxZQXFzMrl27OPfccznrrLO4+OKLqaqqAuC8887jwQcf5Nxzz+X3v/89r7/+Orm5ucyfP59Vq1YN54+wB+EcmFKcgVQ2uxACJsVYmJMazZaielraPcREhJ+7ezKbj9az80Qj6XFWrpyfhsmg/IzxzieHa6mzdwzrcyZFmTlvVvKA87773e+Sl5fHT37ykx7j3//+9/nhD3/IihUrKC0t5eKLL+bgwYNcddVVbNiwgdtvv51t27aRmZlJSkoKN954Y8j5ALt27WLz5s1ERET0uMfvfvc7Hn30Ud5++20AnE4nH374IRaLhaNHj3LDDTf0qMe1detWvve97/Hmm2+SmprKzTffzJtvvklSUhKvvvoqP//5z3n++ecBaG5u5tNPPwVg3rx5vP/++0yePLnHm9Rwo4ReEZJGh5uYCCNmg55ZKVFsKarnSI2dxZnxPeZJKfn8WAMHqlrJiLdywexkjPreAl7T6mJHSSOT4yIob3Lyr0O1XJI7abRejmICEh0dzS233MKTTz7ZQ4g3btzIgQNdBXBbW1ux2+1cf/31PPzww9x+++288sorXH/99f3OB7jyyit7iXwoPB4P9957L3v27EGv13PkyJHgtYMHD3LXXXfxwQcfkJaWRmFhIYWFhVx00UUA+Hw+UlNTg/M77QJYvnw5t912G9dddx1XX331YH9EYaOEXhGSRqebeJsJgBirkbRYC4eqewv9tuONbDuuCfjBqlb0QnDh3JRez7evvAWjXnDl/DR2lzbzRXEDmYlWZk+KHpXXoxga4XjeI8kPfvADFi5cyO233x4c8/v9fP75570E+uyzz6aoqIi6ujreeOMNfvGLX/Q7H8Bms4Vlx+OPP05KSgpfffUVfr8fi8USvJaamorL5WL37t2kpaUhpSQnJ4fPP/885HN1v+fTTz/Ntm3beOedd8jPz2fPnj0kJCSEZdNgUJ+dFb3w+yXNDjdxVlNwbGZKFPX2Durbuj7G19pdbCtuZPakKP7trHTOmhrHvooWalt7HuZwe/0crrEzIzkKi1HP0qx4JsdG8NHBWiqb20ftdSkmHvHx8Vx33XU899xzwbGCggL+7//+L/h4z549gJaps2bNGn70ox8xZ86coGD2Nb8/oqKigl4/QEtLC6mpqeh0Ov7f//t/+HxdhXpjY2N55513ePDBB/nkk0+YNWsWdXV1QaH3eDzs3x+68suxY8dYunQpDz/8MImJiZSVlYWcd6oooVf0wt7hxeuXvYReCDgSyL6RUvLxoVosRh3nzUpGCMHizHjMRh1fHG/s8XxHauy4vX7mpccAoNMJLpk3iQijntd2lvGPryppdroHNsznhYpdUL4LfCq3/0zhvvvu65F98+STT7Jz507y8vKYO3cuTz/dVW7r+uuv529/+1uP8Eh/8/siLy8Pg8HA/Pnzefzxx7nnnnv4y1/+wrJlyzhy5EivTwIpKSn84x//4Lvf/S67d+9m3bp1/PSnP2X+/Pnk5+eH3LwF+PGPf8y8efPIzc1l1apVzJ8/f7A/nrBQtW4UvahqaeeV7WVclZ/GtKTI4Pjfvyyn0eHm1nMyOVxt58MDNRTkpJCTFhOcs/VYPduKG7l52VSSoswAvLK9FLfPzzeWTe1xqs/l8fHliSZ2lzVj0AluXjYVm7mPaKLfB3tegpZy7XFsBuStBb2KPg43Bw8eZM6cOWNthqIfQv0fqVo3ikHhdGsfS62mniK6ODMeu8vLW3sq+eRwLelxEcxN7RljX5ARh8mgY2eJ5tXX2TuoanGRkxbT6+i2xajnnBmJXL84gw6vn23HG/o26tjHmsjPuQLmfA2ay6Bs2zC8WoXi9EcJvaIX7Z1Cb9b3GM+It3L29ATKm9pJjDRz2bzUXuIdYdKTlx7D4Ro7jQ43X5Y2YdQLctL63nRNjDQza1IUB6vsdHhPblIG1B6E8h2QvhgmzYNJuZA0C0q3gqvl1F+wQnGao4Re0YtOjz7CqO91bdm0BL53wQzWLpnSZ5hl4RTNq1+/q5yDVa3MS4/FEuK5upOTFo3b66e0wdnzQnMZHHwbYibD9PO7xqdfoDWlLP50UK9NoTgTUUKv6IXD7cVk0IXMhwdtM7U/bGYDl+WmIgRkJtg4e9rA6WJpMRGYjTqK6x3dDKmHwnVgiYbca0DX7c0iIhYyFkPNfmitCut1KRRnKkroFb1od/uwmvr3wAciM9HGnSun8fUFkzF1NEH9UW1DtQ90OkFGnJWKpkC6ZYcd9r4KQg9514EpRL7zlLPBZIVjH8E4TCpQKMYLYQm9EOISIcRhIUSREOJnIa5fJYTYK4TYI4TYKYRY0e1aiRBiX+e14TReMTI4h0HogzSXws7nYN86Tbj9fVfBTIu10NLuwenqgP0bwNOuiXxEXOgFBjNkrtTCO/VHQs9RKBQDC70QQg/8AbgUmAvcIISYe9K0j4D5Usp84A7gTyddP19Kmd9X6o9ifNHu9hJhGoa0RSmhaCOYo2DaedB0Aip39zl9Uox2crGxaAe0VMCsSyFqgDIJqflgS4TiT/r9xKBQnMmE49EvAYqklMVSSjfwCnBV9wlSyjbZlZBvQ9smU0xQnG4ftuHw6O1VYK+BjCUwZRnEToGyL/r06pOjzOjx4ynerOXJJ5/sT4RAp9PeRJyNUL331G1WKE5DwhH6yUD3c7nlgbEeCCHWCCEOAe+gefWdSOADIcQuIcRdp2KsYuTx+yXtHh8RwyH0DUUghCbYQkD6InC1QmNxyOlGvY6popoORwtkLNXWhEPCDC0r58TWfkNDiolFdXU1a9euZfr06cydO5fLLrusRzGxkSAyUjsgeM4554zofQ4ePEhWVhb+wO+r3++noKCAv/71ryNyv3CEPtRfWy+PXUq5QUo5G/g68J/dLi2XUi5EC/18VwgRsuiyEOKuQHx/Z11dXRhmKUaCdo8PKXsflhoSTSe00IsxUEwqYQYYLVB3sM8lUzzHafGZIH5a+PcRQsuxd7VC0/FTNFoxHpBSsmbNGs477zyOHTvGgQMHeOSRR6ipqRmV+/dVsiAUUsqgYIfLnDlzmD17drAM8oMPPsisWbN61cUfLsIR+nIgo9vjdKCyr8lSys+A6UKIxMDjysDXWmADWigo1LpnpJSLpJSLkpKSwjRfMdx0nYo9RY/e64bWSoid2jWm00P8dGg4Ftrz9vtJ9FRSbUjHNdhwe0K29iZSe2DguYpxz8cff4zRaOTuu+8OjuXn57Ny5Uoee+wxcnNzyc3N5YknngheLykpYfbs2dx5553k5uZy0003sXHjRpYvX052djbbt28Pzrn11lvJy8vj2muvxel09rp/p2cP8Le//Y0lS5aQn5/Pt7/9bXw+HyUlJcyZM4d77rmHhQsXUlZWhsPh4PLLL2f+/Pnk5uby6quv9vsaf/jDH/LHP/6R9evXs2XLFh577LFh+MmFJhy3bQeQLYTIAiqAtcCN3ScIIWYAx6SUUgixEDABDUIIG6CTUtoD3xcADw/rK1AMK+39HJYaFI5akH6ISe85npgdyH0v12L23bFXEan30mKZTIPDzeTYgeuEB9EbtE8BjcXaJnC4YR9F/xzdCG3D7EVHpkD2hf1OKSws5Kyzzuo1vmvXLl544QW2bduGlJKlS5dy7rnnsmDBAgCKiop4/fXXeeaZZ1i8eDEvvfQSmzdv5q233uKRRx7hiSee4PDhwzz33HMsX76cO+64g6eeeor7778/pB0HDx7k1VdfZcuWLRiNRu655x5efPFFVq1axeHDh3nhhRd46qmnAFi/fj1paWm88847gFbxsj8KCgq47777eOCBB/j0008xGofe1GcgBvTopZRe4F7gfeAg8JqUcr8Q4m4hROfb7TVAoRBiD1qGzvWBzdkUYLMQ4itgO/COlPKfI/FCFMOD06P1sTxlj75THCJPqmcelwVCFzpO33Qcq9lAiyWd+qF0NYqfDm4n2KsHv1YxIdi8eTNr1qzBZrMRGRnJ1VdfzaZNm4LXs7KymDdvHjqdjpycHFavXo0Qgnnz5lFSUgJARkYGy5cvB+Dmm29m8+bNfd7vo48+YteuXSxevJj8/Hw++ugjiou1392pU6eybNmy4Nx58+axceNGfvrTn7Jp0yZiYmL6etog55xzDj/60Y+CjUkcDge33nor3/rWt3jxxRcH/fPpi7ACsVLKd4F3Txp7utv3/wX8V4h1xcDI1N1UjAidoZs+q0iGS1utFkoxn1TjxmjRNk4bi7Vsme40lWCKTcPQZutR9z5s4jK1r82lEJ3a71RFmAzgeY8UOTk5rFu3rtf4QNV2zWZz8HudThd8rNPpgs24T67PdPLjk+9366238tvf/rbHeElJSa9SxTNnzmTXrl28++67PPDAAxQUFPCrX/2qX3sPHDjQo6nK3//+d6699lq+9rWvcf3113PTTTf1uz5c1MlYRQ/a3T50QmA+1X6ubTXaR/RQf0Tx07S0y462rjFvB7RUIOIySYw0Da1PqTkSLDHQWjF0uxXjggsuuICOjg6effbZ4Fhn8+033ngDp9OJw+Fgw4YNrFy5clDPXVpaGmwK8vLLL7NixYo+565evZp169ZRW1sLQGNjIydOnAg5t7KyEqvVys0338z999/Pl19+GXyOiorQv5P79+8nNzc3+Li8vJyMDG1LVK8fpkOLKKFXnISjw4vVpO/XyxkQKbW8dmti6OudGTXdM2SaS7WYfnwWiVFmGhzuAb23kESnafn7igmNEIINGzbw4YcfMn36dHJycnjooYdIS0vjtttuY8mSJSxdupQ777wzGJ8Plzlz5vCXv/yFvLw8Ghsb+c53vtPn3Llz5/LrX/+agoIC8vLyuOiii6iqCv37tW/fvuCm7W9+8xt+8Ytf4Pf7KSoqIj4+vtf8srIyYmNje2z8pqenU16u9VwYbCZPf6jGI4oevLmnArvLy83Lpg48uS/cTtjye5ixWjssdTJSwtb/hbipMDdw9u7IB1D9FSz/IYXVDj48UMPtyzOJ7dblKizKtkPRR3DO9zQPXzFoTufGIyUlJVxxxRUUFhaOyv0KCwt5/vnnw86ocTgc3HvvvVgsFlasWNFn6GawjUdUex5FD4alzk17k/a1rxo1QmhefUORlmapC2zOxmaC3hDsTFVn7xi80Hdu/jpqldArxpzc3NxBpU3abDZeeOGFYbdDhW4UPRgWoXc1a1/7EnrQhN7TroVZnI3am0N8lnbJZkIIqBvKhqwtcAbDUd//PMUZSWZm5qh58+MJ5dErgkgpaXd7T/1UbKdHb+knvSwuU/PsG4u1dEuAhOmAVgoh3jbEDVmTTfvnUKerFYpOlNArgnh8Eo9PDk/oxhwF+n4OgJis2oGpil2A1OL13T4BJEaaqWpxDe3+tiQl9ApFN1ToRhHE6dbyjE+5oFl7c/9hm06yVoHPrZVLyOpZAikx0kxruweXZwilhzuFfhwmGigUY4Hy6BVBuurcDEPoJhCG6ZeYdFh8p/a9tWf6WfcN2Yx46+Dub40Hnxc6WvsPHykUZwjKo1cEGZaCZl43uB3hefSgibK1d47xpGgLAJXN7YO3ofP5nI2DX6tQnIYooVcECRY0OxWhDyfjJgwiTHoSI01UtgxB6Dvv3bkprFCc4SihVwTpjNFbT6VyZTDjJvaU7ZkcF0Flswu/f5CxdnM06AzQrjz6icyGDRsQQnDo0KHgWElJSY+SAaPNE088EbKs8clkZmZSXz9+UnyV0CuCOD0+zEYdBv0p/FoMdFhqEEyOteL2+qluHWT2jRAQEattCismLJ11aF555ZWxNiVIuEI/VDoLrw03SugVQZwdvlPz5kETV2OEVqXyFJmaYEUnBMV1jsEvtsarGP0Epq2tjS1btvDcc8/1KfSXXXYZe/dqfYIXLFjAww9rrS5++ctf8qc//Ym2tjZWr17NwoULmTdvHm+++Wbw+u9///vg8/z85z/nySefpKqqilWrVpGfn09ubm6P8scATz75JJWVlZx//vmcf/75AHznO99h0aJF5OTk8B//8R+9bGxvb+eSSy7h2WefxeFwcMcdd7B48WIWLFgQtOfPf/4z//Zv/8bXvvY1CgoKBrRjKKisG0UQ53Adloo49bANgMWoJyM+gqJaO8tnJAyu0FpEXFcnK53yZ4bK5orN1LcPbwgiMSKRFZP7rhgJ8MYbb3DJJZcwc+ZM4uPj+fLLL1m4cGGPOatWrWLTpk1kZmZiMBjYsmWLZvPmzdx8881YLBY2bNhAdHQ09fX1LFu2jCuvvJJvfvObXH311Xz/+9/H7/fzyiuvsH37dv785z9z8cUX8/Of/xyfz9fLc//3f/93HnvsMT7++GMSE7WCfb/5zW+Ij4/H5/OxevVq9u7dS15eHqC9Wa1du5ZbbrmFW265hQcffJALLriA559/nubmZpYsWcKFF2ploD///HP27t1LfHw8//M//9OvHUNB/QUoggxLU/D2pmEJ23QyIzmSJqeH+jb34BZGxIPfBx39d/lRjE9efvll1q5dC8DatWt5+eWXe81ZuXIln332GZs3b+byyy+nra0Np9NJSUkJs2bNQkrJgw8+SF5eHhdeeCEVFRXU1NSQmZlJQkICu3fv5oMPPmDBggUkJCSwePFiXnjhBR566CH27dtHVFTUgHa+9tprLFy4kAULFrB//34OHOhqZXnVVVdx++23B/vAfvDBB/zud78jPz+f8847D5fLRWlpKQAXXXRRsMLlUOwYCOXRK4I43T4mx56C0Pt9Wu56xPBtls1IjuRfh2o5WmMP5taHRffMm2F84znTGMjzHgkaGhr417/+RWFhIUIIfD4fQgj++7//u8e8xYsXs3PnTqZNm8ZFF11EfX09zz77bLAF4YsvvkhdXR27du3CaDSSmZmJy6Xt99x55538+c9/prq6mjvuuAPQPiF89tlnvPPOO3zjG9/gxz/+cb/Nuo8fP86jjz7Kjh07iIuL47bbbgs+P8Dy5ct57733uPHGGxFCIKVk/fr1zJo1q8fzbNu2rUcTk8HaEQ5hefRCiEuEEIeFEEVCiJ+FuH6VEGKvEGKPEGKnEGJFuGsV4wO/X9Lu9p1a6Ka9WTuNOozCajUZmBJv5XCNfXD16YO59CrFcqKxbt06brnlFk6cOEFJSQllZWVkZWX1avlnMpnIyMjgtddeY9myZaxcuZJHH3002IikpaWF5ORkjEYjH3/8cY+GIWvWrOGf//wnO3bs4OKLLwbgxIkTJCcn861vfYtvfvObwcYh3YmKisJutwPQ2tqKzWYjJiaGmpoa3nvvvR5zH374YRISErjnnnsAuPjii/nf//3f4O/x7t27Q77+cOwYLAMKvRBCj9YH9lJgLnCDEGLuSdM+AuZLKfOBO4A/DWKtYhzQ7hmGw1LDmHHTnZkpUTQ7PdQOpsiZKVJrGK5y6SccL7/8MmvWrOkxds011/DSSy/1mrty5UpSUlKwWq2sXLmS8vLyoNDfdNNN7Ny5k0WLFvHiiy8ye/bs4DqTycT555/PddddF+zk9Mknn5Cfn8+CBQtYv3493//+93vd76677uLSSy/l/PPPZ/78+SxYsICcnBzuuOOOYB/a7jzxxBO4XC5+8pOf8Mtf/hKPx0NeXh65ubn88pe/DPn6w7FjsAzYeEQIcTbwkJTy4sDjBwCklL/tZ/7zUso5g13biWo8MvrU2l28+EUpV+Slkp0yxJhg2Q4o2gjL/12rIDlMuDw+nvmsmPyMWFbNTAp/4Y4/gTkG8v5t2Gw5EzidG4904vf7WbhwIa+//jrZ2dljbc6gGWzjkXBCN5OBsm6PywNjJ99kjRDiEPAOmlcf9trA+rsCYZ+ddXWq8uBoMyynYtubwGAG4yBr0wyAxahnaoKVI4MN30TEKY9e0YsDBw4wY8YMVq9ePSFFfiiEE5ANldPW669NSrkB2CCEWAX8J3BhuGsD658BngHNow/DLsUwMiwFzTo3Pk+l32wfTE+KpLjOQYPDTWJkmJuyEfEqxVLRi7lz51JcXDzWZowq4fz2lwMZ3R6nA5V9TZZSfgZMF0IkDnatYuwYloJmI5jhkhGnfUqoaBpE7ZuIOJViqVAQntDvALKFEFlCCBOwFnir+wQhxAwROM0ihFgImICGcNYqxgdOtxe9TmA2DNHz9fvA1TJiQh8dYSDKYqBiMNUsVXEzhQIII3QjpfQKIe4F3gf0aBut+4UQdweuPw1cA9wihPAA7cD1Ugumhlw7Qq9FcQo4OrResYM6fdodZwNIP9gSh9ewAEII0mIjBle2WAm9QgGEeWBKSvku8O5JY093+/6/gP8Kd61i/OF0e7GZTyE+39m6z5Y8PAaFIDnKzOFqOy6PD0s4NXnMUSrFUqFAlUBQBHC4facWn2+rBZ0+ZBOR4aJ716mwECKQeaOqWCrObJTQKwBod3uxnUrGjaMOrAma2I8QnUI/qINTEXGqiqXijEcJvQK/X+J0+7CahyjSUkJrJUSmDK9hJ2E1GYg0G8L36EETelezlmKpmFBUV1ezdu1apk+fzty5c7nssss4cuTIiN4zMjISgHPOOWdE73Pw4EGysrLwB34v/X4/BQUF/PWvfx2R+ymhV9Du8SElQ/fonY3gaYfYjIHnniKJUSbq2wYp9J3F1hQTBikla9as4bzzzuPYsWMcOHCARx55hJqamlG5/9atW8OeK6UMCna4zJkzh9mzZ/P2228D8OCDDzJr1qxTLl7WF0roFTgCLQRtQ/Xom0u0rzEjL/RxVhPNTnf4J2QjAnsGakN2QvHxxx9jNBq5++67g2P5+fmsXLmSxx57jNzcXHJzc3niiSeC10tKSpg9ezZ33nknubm53HTTTWzcuJHly5eTnZ3N9u3bg3NuvfVW8vLyuPbaa0PWe+/07AH+9re/sWTJEvLz8/n2t7+Nz+ejpKSEOXPmcM8997Bw4ULKyspwOBxcfvnlzJ8/n9zcXF599dV+X+MPf/hD/vjHP7J+/Xq2bNnCY489Ngw/udCoMsUKnB2neCq27rAWnx+FcsDxNhMen8Te4SXaYhx4QY8Uy6wRte10pG3TJrx1w9t4xJCUSGSg8FhfFBYWBssNd2fXrl288MILbNu2DSklS5cu5dxzz2XBggUAFBUV8frrr/PMM8+wePFiXnrpJTZv3sxbb73FI488whNPPMHhw4d57rnnWL58OXfccQdPPfUU999/f0g7Dh48yKuvvsqWLVswGo3cc889vPjii6xatYrDhw/zwgsv8NRTTwGwfv160tLSeOeddwCtemZ/FBQUcN999/HAAw/w6aefYjSG8fs8RJRHr6CtI9AUfChZN456aC6F5DkjUvrgZOKsJgCaHGE2IjFHqUbhpxGbN29mzZo12Gw2IiMjufrqq3u02svKymLevHnodDpycnJYvXo1QgjmzZtHSUkJABkZGcFKkzfffHOv8sfd+eijj9i1axeLFy8mPz+fjz76KFg+YerUqSxbtiw4d968eWzcuJGf/vSnbNq0iZiYmAFfzznnnMOPfvQjUlNTASguLuab3/wm11577aB/Nv2hPHpFUOgjw82jlxLqj2oCX38E9CaY3Nv7GgnibZrQNzrcTE0Io0KmahR+SgzkeY8UOTk5rFu3rtf4QCE7s7mrDpJOpws+1ul0wcbbJx8K7O+QoJSSW2+9ld/+tmfB3ZKSkh7NQgBmzpzJrl27ePfdd3nggQcoKCjgV7/6Vb/2HjhwgNtvvz34eNq0aTz33HPDLvTKo1dgd3mxmfUY9GH+OpzYAoXroWoPmKyQdzI3elAAACAASURBVJ32dRSwmvSYjTqanINoLTgGVSz90s+hxkO8dvg11h1ZN+x9V093LrjgAjo6Onj22WeDYzt27GDBggW88cYbOJ1OHA4HGzZsCNafD5fS0lI+//xzQKt9v2JF3120Vq9ezbp166itrQWgsbGxRwOT7lRWVmK1Wrn55pu5//77gw1DVq9eTUVFRcg1+/fvJzd3+Dqy9YXy6BXYXR6iwol3A3TY4cRWLVQz58pRrwophCDeaqLJ4Ql/kTUeGotHvIplp7iX28upclTh8DhIsCTg8Dp4t/hdbph9A0b9yMVhTyeEEGzYsIEf/OAH/O53v8NisZCZmckTTzzBbbfdxpIlSwCtJWBnfD5c5syZw1/+8he+/e1vk52dzXe+850+586dO5df//rXFBQU4Pf7MRqN/OEPf2DSpEm95u7bt48f//jH6HQ6jEYjf/zjH/H7/RQVFQX7wXanrKyM2NjYHhu/I8WAjUfGAtV4ZHT5y9YSEiJNXJGXNvDkks3av6XfHrNerP8srKa8ycmdK6eFt6B6Hxx8G5Z8a8Rq8Ugp2Vi6kaNNR7EZbaTYUpgZN5Os6CwqHZW8WfQm56afS05izojcfzg5nRuPlJSUcMUVV1BYWDgq9yssLOT5558PO6OmoaGBn//853z44YfceeedPPDAAyHnDbbxiPLoz3CklNhdHjITw+wIVX8EotPGtOF2vM3EwapW3F4/pnCqbXYe5LJXj5jQl9vLOdp0lEUpi1g8aXGPuG+aLY0ESwIHGg9MCKFXDB+5ubmDSptMSEjg6aefHnjiIFEx+jMcl8ePxyeJsoTxnu92gr0GEmaMvGH9EGfVwh/N4cbprQla5k3byB22+ar+K2xGGwtTFobc7MuOy6bOWYfD4xgxGxQDk5mZOWre/HhCCf0ZTqtLi3VHhyP0nUIZlTqCFg1MXGfmTbhCr9Nrnnxb7YjY4/F7qLBXMC1mGgZd6J/jlOgpAJTZy0JeVyhGEiX0Zzid2Suxgfz0fukUyhGuaTMQsRFGhGBwG7JRk8BeNSI1b6rbqvFJX1DMQ5FgScBmtCmhV4wJSujPcJocHoTQxHNA2mq0A0ijlErZFwa9jmiLcXAplrFTwNsBbdXDbk+pvRS90JMW2fdmthCCSbZJVDuG//4jwXhM0lBoDOX/JiyhF0JcIoQ4LIQoEkL8LMT1m4QQewP/tgoh5ne7ViKE2CeE2COEUKk044xmp5toizG8HPq2mjH35juJsw1S6OMyta+Nx4fdllJ7KamRqRh1/b9ZplhTsLvt4z5Ob7FYaGhoUGI/DpFS0tDQgMViGdS6AQOzQgg98AfgIrRm3zuEEG9JKQ90m3YcOFdK2SSEuBR4Blja7fr5Ukp1YmQc0uT0EGcLw5v3ebQqlUmzRt6oMIizmqhsbkVKGV77Q5MNolOh9gBMPUcrq1xTCInZEB9mmmYI7G47Ta4m5sTPwVNTg2PzZnzNzVjmzsW6bFkP2ybZtNzrakc102OnD/meI016ejrl5eXU1dWNtSmKEFgsFtLT0we1Jpz0yiVAkZSyGEAI8QpwFRAUeill95qeXwCDs0IxJkgpaXK6SY2NHniyo17rCTtePHqrCbfXT1uHN/zDXpPP0vLp972uHaCSEqq+gkV3DDntsqzqENaKRpKqT9BctAWd1YohOQXnzl0IoxHroq605qSIJPRCT42jZlwLvdFoJCtLFYA7nQgndDMZ6L6DVB4Y64tvAu91eyyBD4QQu4QQdw3eRMVI0ez04Pb6SYo0Dzy5M7YdOXI9YQdDZ82bZucgNmSTc7RPJI3FkDwXlt2t1cIp3zHo+/vsdpo3vEHrX18kfVsJpqJyIublEnfTjURfcTnm7Bk4tm/H29RVekGv05NkTaLaOTHi9IrTh3A8+lCfi0MG74QQ56MJfffiEcullJVCiGTgQyHEISnlZyHW3gXcBTBlSt/ZC4rho7MlX3J0OEJfCwYTWGJH2KrwiA3k0jc63GTEh7k5rNNB7tVaI5LOlocJ2VqBtpmXhF1909/RQcsbb+JzOiidFUfS9BUk5l2F6FZm1rZiJe6SEpw7dhBdUBAcT7GmUFhfiM/vQ6/T49yxg46SEmxnn4MpvT//SaEYOuF49OVA944S6UDlyZOEEHnAn4CrpJQNneNSysrA11pgA1ooqBdSymeklIuklIuSkpLCfwWKIVNrd6HXCRJs4Qh9YCN2FEoRh0Ok2YDJMMjiZp1072ubMB3cjrAPU0kpsW/ciK+1BfcFS6mfmUT6tPk9RB5AH2nDkjuPjiNHe3j1ydZkfNJHg6sB1+EjOL7Yhre6BvsHHyA9oT+dePweml2q+qZi6IQj9DuAbCFElhDCBKwF3uo+QQgxBfg78A0p5ZFu4zYhRFTn90ABcOYdSxun1LR2kBBpQq8bQLyl1Dz6cRKfBy1dMdY6yMybUMROweX34gqRjeNrbcVTWYnf3XWP9p07cRcfJ3L5ckpt7QghyIgK3VnLuiAfYdDj7Fa3KcWm/Qxr26pxbtuGISmJmK9fhd/hoKO4tw1SSt4+9jYvHXqJwnr1p6MYGgOGbqSUXiHEvcD7gB54Xkq5Xwhxd+D608CvgATgqUCWgTdQXCcF2BAYMwAvSSn/OSKvRDEofH5JdUs7uZMHbo5Ae5OWdTNO4vOdxFlNVLe4Tuk5an0u3nIWI4vWc2XiNFJsKUgpcWzZSvvu3dokvQ5TRgbCZKbjyBHMs2ZimT+fksOvkGZLw2IIneqms9mw5OTS/tVXeM86C0N8PFHGKCIMETQU7ialxU705ZdhTE9HFxVJx5EjWGbN7PEcnZUwAb6s+ZKchJzwsowUim6ElUcvpXxXSjlTSjldSvmbwNjTAZFHSnmnlDJOSpkf+LcoMF4spZwf+JfTuVYx9tS0uvD4JJNjIwae3BnWGEcePWhC3+ry4PUN/bTr1qrPEeZojG4Hn5Z/ipSS9l27aN+9G0tODtFXXI5/9jQKj25h75f/xDk7najVq6lx1tDkahowe8Z61kKE0Yjziy8A7ZNISkQyHbv3YEhMwJSVhRACc1YWnooKpM/XY31RcxEmvYnzMs6jzdNGrXNkyjgoTm/UydgzlIrmdgAmx4Up9EIH1pGp/DhU4mxGpITm9kFk3nSjzd1GZVslCxJzWaaPod5ZR0nZPhzbtmHOnkHk+efBlDT+mVLDwYtmcOzK+byTVMmmqi18Xvk5Jr2JmXEz+72HzmolYkE+HceKcZdrzSdS6r14GxvRz88NeufGyZORHg/e2p5CXt5WTnpUOlOitASFTu9eoRgMSujPUCqa2om3mcJrCN5WC7YE0I+vqtbxg+0fexKddWemxs8i2xCFTeg4/uEbCL2ByJUrEUKwrWobTo+Tr037GjfOuZHcxFwK6wupclSxPG05Jv3ANYKs+fnoY2Kwf/gh7hMniN1zHE90BM1pXecXjGla+QRPt05ETo8Tu9vOJOskIk2RRJuiJ0wJBcX4Qgn9GYjPL6lobu8K21R9BZv+Bw69o228nsw4Kn3QnTibCSGgvm3oQm8z2oiPzUIvdMx3ROIqPoY7Zxo6m42qtioK6wuZlzSPFFsKRp2RVemruG7WdaydvZY5CeE15xAmE9GXXoL0emh56x/Y/EbqF0+nxtXlveusVvRxcXiqu7J/OsM0ydbk4Ne6dnVaVTF4xpeLphgVqlracXv9ZCZatabZRz7QcuSr9kLsVJjUrYdlhx062iCyd+u0scao1xFnNVHX1jHotX7pp8xeRmZMJsKagJSSyfsqqIywcCDNT5LfwyflnxBpjGTppKU91iZGDD6EZUhKIv7GG3GXV2BMS8Va8XaveLshOQlPeZdHX9deh0CQFKGlG8db4ilqLsLtc4f1SUKh6ER59GcgJxqc6IQgPc6qtdmTPlh4q1YGoPTznl69PRAqiBp/Qg+QFGWmzj54oa9z1tHh69BSI01WPE0eqKoi4exVHLUX80bRGzS5mjgv47xh6/Oqs9mwzJqJPiqKZGsyNc6aHoXDjCkp+B0OfG1tANQ4a4i3xAfvHx+h9R1tdDUOiz2KMwcl9GcgJxqcpMZYsBj1UHcIYjIgIhYylmo1bZpLuybbq7VDUuMwdAOa0Le2e3B5fANP7kaZvQyBID0yHenz0Xa0Eb0ZFqy6hvSodJpdzZyTdk6/NeZPhUnWSbi8LlrdrcExQ7IWovHW1iKlpNZZS5K16/BggiUBUEKvGDxK6M8wnG4vNa0upiZYoa1OE/ak2drF5DlgtEDll10LWiu1VnyG8Rkq6KzTM1ivvsxeRqI1EavRinPbNnxOH7a5KZiMFr42/Wt8K+9b5Cfnj4TJQFfcvcbZFZM3JCaCTuCtraXV3YrL6yLF2vUGG22KxqAzKKFXDBol9GcYJxqcAFoz8LpDmrfeWXpYb4RJ86DuiBaX93mhpVSL249TkqICQt9PnN7vcuGpqcVdXo6npgZXaxPVjioyrJNx7tyJc9eXWGZnY47Vj0gHqlAkRCRg0Bl6xOmF0YghIQFvTU2vjVgInAY2x9LcocohKAaH2ow9wzhe78Bq0pMcZYZDgbCNObJrQtpCKNuhZeLETNbEvrNpxzjEZjZgM+upbe19QrajqAjnzl14T6qr3uRqYkrrMZIT23DorJizs4nMSYSiD8BtB0sYp4VPEZ3QkRSRRI2jZ40dQ3IKHceKqHUkoBd64i3xPa7HWeJUiqVi0CihP4Pw+vwcr3cwKyUK4ajXwjbZBT0nWeM1Ya/aA446LWQzjoUeIDUmgsrmnkLv2L4d57bt6BPisS1bij4hAWEyIz1uDhZ9THuth8SMZVimTsU4ZQqiKVBnxtUyKkIPWt2bvXV78fg8wQ1XQ0oyrv37qa85QVJ0EvruBdiAWHMsRU1FePyeATtaKRSdqNDNGcSJRidur58ZyZG9wzbdmXwWuFqh9iCk5o/b+HwnabERtLR7sLu0E7LukhKc27Zjnj2LuLVrsS5ejHnaNEzpkzFMnUJRKsSuWEX0qlWYpk7VTqd2ll9uH72wSEZkBn7pp9LRVQzWOGkSfumnraKkR9imkzhzHBJJS0fLqNmpmPgooT+DOFJtx2zUkREXoYn4yWGbThKzYdYlWsu9rFWjb+ggSQ+UcahsdmkZNJs2o4+NJer88xE6HW6fG49fexM41nIMl9fFjNgZPZ/EHK298blGT0AnRU5CL/SU28uDY/q4ONqFB31DC5OsvVNaYwNvSE2upl7XFIq+UKGbMwSn28vR2jbmpcegd1SDswEyQrYG0AQvbcHoGngKJEWaMRl0VDa3M7WlEl9zM9GXX44wGNhfv5/PKj5Dh46smCyqHFXEWeKYGn3SBrPeAKbIURV6o85IamRqsBQDgNDpsMeYMNe2BXvMdifWHItAqA1ZxaBQHv0Zwq4TTfj8kvnpsVCzH3SGrrTKCY5OJ5gcG8Hxegft+/ahj4nBlJVJQ3sDmyo2kWpLZXb8bErtpQgEF065MHSpX0vMqAo9QHpkOo2uRtrcbcGxxiiIsnuxid4NYQw6A1GmKOXRKwaF8ujPAGpaXewubWZuWjTxEXqoPaB1VjKGrqM+EZmeFEnF8cO0nCgn+fxV+KWff5X9C5PeRMHUAqxGK6vStTBUn/XcLdHauYFRJCsmiy+qvqCouYj85HwtZm/tYJrRhreuLljsrDtxljgl9IpBoTz605xWl4c391RgMxtYmZ0ITSXgdkJK7oBrJxLTkmzElB2j0enFMns2u2t3U+es49z0c7EatZ6yQoj+m3aYo7XaPqEKu40QcZY4EiMSKWouArRCZq0xBmLMsXiqQ6dRxpnjaO5oxi9HJ+dfMfFRQn8a4/dL3tlbhccn+Xp+mlaSuHofGCM0j/40wqqTTGmq4HhkMuXeJnbU7GBG7AwmRUwNvzyCJVprHN4tjDIazIybSa2zlhpHDUeajkBEBPGTMnGXloacH2uJxSd92N32UbVTMXEJS+iFEJcIIQ4LIYqEED8Lcf0mIcTewL+tQoj54a5VjByHqu1Ut7hYPSeZhEgzeFxQfxSS5/ZskD2BkH4/ntpapNfbY7xkzyY8vuNsjTnB/+14GZ/HgqM5m2c3FfPc5uPBRiv90pliOcpx+rkJczHrzWws3ciBhgNkx2VjzZyOt6oqZMPwOHMcgNqQVYTNgEIvhNADfwAuBeYCNwgh5p407ThwrpQyD/hP4JlBrFWMEIUVLSREmpiVEqUN1B8GvxdScsbWsCEi/X5a3niT5ldfo+mll/HZNY+2tLWUvZ+swx1rIG1mNv72qTgbzqKq2Ud+RiwRRj0f7K/G7x8gJGMONAJxtfY/b5gx6U1cMOUCnB4nUaYolqUuwzQlA+n19WhE0olKsVQMlnA2Y5cARVLKYgAhxCvAVcCBzglSyq3d5n8BpIe7VjEytLo8VDS3s3xGYldcuv6oFp6I7r3BNxFwHTiIp6KCiPz5uA4cpOWtt4i5+mp27niL2DYvi6/5NpfMy+N4vQOvX5KZYMNk0JEeZ+cfX1VxvMHB9KQQ5wY6sQSEvmN0hR60Tdk7cu9ACIFO6JCTJyPMZlxHjmDKzOwxN8IQgcVgUUKvCJtwQjeTgbJuj8sDY33xTeC9Ia5VDBPljVqoIivRpg34fdpGbPx0LU9+AuLatxdDUhK2FSuIvvxyfC0tlL36V4xb95A6eTaROVoP1mlJkcxMicJk0H69sxIjsZn1HKoaIKZtMGv/Rtmj70Sv06MTms3CYMCcPQN3cTF+V+86Pp0bsgpFOIQj9KFUIeRnYCHE+WhC/9MhrL1LCLFTCLGzrk61SztVKpvbMRt1JEYGyhe0lIHPA/HTxtawIeJracFb34BlzmyEEJjSJxN98cXUNpxAmExkfv1GhD70voNeJ8hMsHGi0TFw+MYSPSYefSgi5s1Den04t2/vdS3OEkdTh/LoFeERTuimHMjo9jgd6JVsLITIA/4EXCqlbBjMWgAp5TMEYvuLFi0avfy205SqlnZSYyxdYZvmMs2Tjxu/JYf7w12ulQkwZnT9OpmnT+foZTnEW+KJSO6/A1Zmoo39la3U2F2kxkT0PdESC67x4SkbEhOJmJdL+1d78dbWghDoIqOwLV9OrDkWl9dFu7edCEM/r0ehIDyPfgeQLYTIEkKYgLXAW90nCCGmAH8HviGlPDKYtYrhx+Pz0+BwMym6mwC0VmitAg29T1tOBDzlFehsNvRxccExp8dJS0cLqbbUAdenBRqhV7X0DoP0wBw9ZqGbUNhWrMC6eHHgkcB9vJjWt98m1qjtJzSPkzclxfhmQI9eSukVQtwLvA/ogeellPuFEHcHrj8N/ApIAJ4KeJBeKeWivtaO0GtRBGh2epAS4m2BsI2U2onP5Dlja9gQkVLiKS/HmJHe48BTlaMKIGRNmJOJNBuIshioGUjoLdHg7dD+jYM3RaHXY1u2FNAalHcUFdH63j+xndDCm00dTaRGDvxGpzizCasEgpTyXeDdk8ae7vb9ncCd4a5VjCxNTjcAcbZAvXJngyZc0RNzH9zX2Ijf6cSUnt5jvLKtEoPOQFJEUh8rezIpxhKeRw+aVx8Z3vOOJqbp0zEkJeLfX4Q+R6cybxRhoU7GnoY0OgJCbw149HbN8yVqfHp+npoa3CUlyD5KD3g64/MnCX2Vo4oUa0qv5hx9kRJtoWWgRuKdTUdG+dBUuAghsOTk4m9sIsllVP1jFWGhhP40pMnhJjrCiFEf+O911GsnYa3x/S8cA9r3FdL82uu0/ONt2v71r5Bz3GXl6GNi0EdHd4353DS0N4QVn+8kIRDK6nwjDEkwl358Cj2Aefo00Akm1XiocdaomjeKAVFCfxrS5PQQZ+3WZs5RDxFx467sga+1FceWzZimZBAxPw/XgYO4T5zoMUf6/XgqKnp58zWOGiRycEIfqcXc6/tpJI4pUvs5jaMN2ZPRWa2YMjKIr2zD7e2gvr1+rE1SjHOU0J9mSClpcrqJs3Vr/+esB9v4ize3796N9PuJvOACbOecgy4qEueuL3vM8dbVId1ujOk99xcqHZUIIUixpYR9v2iLAZNBR0NbPx69EGCOGje59H1hzs4myq3D3Oigoq13mQSFojtK6E8z2jq8uL1+4jvj81631gfVlji2hp2EdLtxHTqMecYM9FFRCIOBiLw8PBUVeLsdmHOfOAFCYMrI6LG+ylFFYkQiJn34/WyFECTYTP179DDuUixDYZo2DZPRQmqth2PNx8baHMU4Rwn9aUaTQ6t2GEytdAbOro0zj76juBjpdhOR21UX3zJ3LsJooP2rr4JjnrIyDMlJ6CK6zgT4/D5qHDWDCtt0khhppsHh7nPjVzNk9DtNDRad2Yxp6lTS6yS1bTWqWbiiX5TQn2Y0BlMrA0LvCHjH1vHl0buLi9FFRmJI7RJrncWCefZsXEeO4Hc68TsceKqrMU3peZq3rr0On/QNSegTIk20u3043f1l3kRrNen9YdaxHyPM2dnE+yOIqLdTWF841uYoxjFK6E8zmpxuTAYdNlNg49UZyLiJiOt/4SgiPR7cpaWYsjJ7dXyKmD8ffH7a9xXiOnwY/BLzzOweczoPSg3Vowf6j9Obo7VDZh3ju7GHOSuLiKg4ZpX6OFC/n3ZvGDX3FWckSuhPM5ocbuKspi4BdTRoaZW68fNf7S4vR3q8mKf1LrBmiIvDNC0L566dOLdvx5iejiG+Z1poVVsVMeaYYIvAwZAQKPJW11+cfgzLFQ8GYTRiXXQWk+1GbAdOsKNy21ibpBinjJ+/fsWw0OhwE2/rnlpZN+7CNp7yCoRBH7LxNUDUeedhSExEZ4sk8rxze1yTUlLlqBqSNw9gNRmwmfU09Cf05s5DU+Nb6AEseXnEzZnPrGI3zS++THXhjrE2STEOCasEgmJi4Pb6sbu8XSdivW5tUzF1fv8LRxlPZSWG5BSEIfSvn85mI+6660Jea3A10OHrYHLk0Ms5JNjM1PcXuun06Mf5hixomURRFxeQPW0K1W8/zaH1LxAblYRlauZYm6YYRyiP/jSiObAR25VxEzhIM45SK6Xbjbe+DuPkoXW56swZT4scepeshEgTjY6OvjNv9EYwWcd96KYTIQTRM+cw8xvfocXs49B7r/SfVaQ441BCfxrR4DhJ6B2dQj9+Uis91dXglxhThxZ6qWyrJNoUTZQpasg2JEaa8fgkLe29G28HmQC59CeTnTSb6CXLqCo/RPWxvWNtjmIcoYT+NKLR4UYnBLHWbqmVOoPWTGOc4KmsBCF6pFWGi8/vo6Kt4pS8eejKvOn34NQ46jQ1GJaefQ0Gk4Vtn72qSiMogiihP41odLiJtRrR6wIZN87xl3HjqazSNlpN4Z9o7aTKUYXb5yYzOvOUbIi3mRCC/uP05hit09QEC4FYI6LJPetiIioaWX9oHZ+UfaLVBZpgr0MxvKjN2NMILeOmm4A66iAmve8Fo4z0evHWVGPpdhp2MBxvOY5e6MmIyhh4cj+YDDpiIowDePQx4POCp12L108g4mbkMK+oBLM/iSNNRzjQcIC0yDQumnoRNqNtrM1TjAHjx9VTnBJen59mp6dL6L0dWox5HMXnvXV1SK+vz7TK/vBLP8UtxWREZWDUGwdeMACToi1UNbv69nQnSC59KExTMjDqTSz2ZnBrzq2smLyCWmct7xS/g8ffz76E4rQlLKEXQlwihDgshCgSQvwsxPXZQojPhRAdQoj7T7pWIoTYJ4TYI4TYOVyGK3pS19aBX0pSogPt7zpLH9iSx86ok/BUan3hh7IRe6L1BA6Pg9kJs4fFlrTYCNo6vLS6vKEndO80NcHQRURgSEnGXXoCs95MXlIeBVMLqG+v58uaLwd+AsVpx4BCL4TQA38ALgXmAjcIIeaeNK0R+Hfg0T6e5nwpZb6UctGpGKvom+pAi7yUaIs20FarfR1HqZWeykr0cXHorIMLhUgp2V27G5vRdsrx+U46m4VXNPVRNmACe/QApilT8dbU4m/XXl9mTCbZcdnsqd2Dw+MYY+sUo004Hv0SoEhKWSyldAOvAFd1nyClrJVS7gDU58IxoqrFFWiAHQhrOOrBYOpqjTfGSL8fT2XVkMI2Rc1FVDuqWTRpEToxPNHGBJsJq0lPcX1b6AlGK+gN2obsBMSUmQlS9mjksmTSEvzSz776fWNnmGJMCOevZjJQ1u1xeWAsXCTwgRBilxDirr4mCSHuEkLsFELsrOtWj1wxMD6/pKTBwdSEbp6yo06Lz59UNGys8NbXaw1EBnlQqtpRzafln5JiTWF23PCEbQB0OsHMlCiO1zloD1XJUohA5s3E9OgNyUnorFbcJSXBsRhzDFkxWeyv34/Hp3yyM4lwhD6UUgwmV2u5lHIhWujnu0KIVaEmSSmfkVIuklIuSkoaPxuIE4HSRicdHj/TkgIZFVJ2Cf04wdsZn58cvo9Q0lLCP479A4veQkFmQdhNwMMlLz0Gv4SNB2sormujsrm95+bsBM2lB+20rCkrE/eJUqSv641sftJ8OnwdHG46PHbGKUadcIS+HOiez5YOVIZ7AyllZeBrLbABLRSkGCY8Pj9fFDcQaTaQlRipDbrbtLTAcST07ooK9DHR6CMjw5p/vOU47x1/j1hLLGuy15zSSdi+SIg0c/b0BIpq23hzTyWv7ijj719W4PEFmm1PwNOx3TFlZiLdbjzl5cGxSbZJJFmT2Fu3V+XWn0GEI/Q7gGwhRJYQwgSsBd4K58mFEDYhRFTn90ABoDokDBNNDjcvby+lptXFubOSug5K2au1r1GTxs64bkivF09ZOcaM8PLfnR4nH5V+RKI1ka9P//qI5n4vyYrnjuVZ3LBkCufOSqK00cnWY4GuXJZocDu0fPoJiGnKFHTWCNr3df3JCSHIS8yjuaOZMntZP6sVpxMDHpiSUnqFEPcC7wN64Hkp5X4hxN2B608LISYBO4FowC+E+AFahk4isCFQG90AvCSl/OfIvJQzbBKAnwAAIABJREFUC6fby+u7yvBLWLNgMlMTuolhayUIHUSG3zh7JPGUlyM9HsxZWWHN31mzE4/fw4VTLhyWnPmBiLEaicHIpBgLDW1u9pQ2s3BKLFGdG9kdrdoJ4wmGMBiw5OTi3LkTd1kZ+thYPJWVZBj1WPUR7K3fy5ToKWNtpmIUCOtkrJTyXeDdk8ae7vZ9NVpI52RagfFVI/c04fNjDbg8fm5YMoWkKHPPi/YqLa1yFEQyHFxHjiBMJozpA5/SdfvcHGo8xKy4WcRZRr8r1pLMePZXtvBVWQsrEruVK56AQg9gXbiAjqNHaXnjzR7jeZNNfJF1giZX05j8nBWjizoZOwFxeXwcrGplTmp0b5GXUhP66FMr/DVc+Ox2OoqKsMye1Wf9+e4UNRfh9XuZm3DyUY3RIcZqZFpSJAeqWvB3CqCzcUxsGQ6EyUTsNVdjO+dsIs9dRez11xGRn8+kMge26la+rFUHqM4EVK2bCcjxegcenyQnLbr3RUcdeFwQPfTGHKeK9PvxOxz4HQ4cWz9H6HRELFgQ1trilmJizDGkWMcu7DRnUhTHatsod+iZYrSAo3bMbBkOdFYr1rPOCj42JCTgLikh50Qp2ycdIichh0m2vvdzfH4fQohhO8OgGH2U0E9ASuodWE16UmMsvS82BQ7IxE0dXaMCdBw9StvmLfjbAgeR9DqiLrgAfXSIN6WT8Pg9VNgryEnM6dU0fDTJSrRhMug4XNvGFFtSVzmJ0wSh1xORn0/6xw0U2eH9kve5cvqVPUI4Hr+HstYyDjUe4oRdK6Vw4ZQLVUx/gqKEfoIhpaSkwcm0JFtoMWwq0eLJY3Aitn3/ftr+9TGGlGSsixahs1kxJCeHnVJZYa/AJ31MjRqbN6lODHodmQk2jte3IVMSEbX7tZDYODl8NhyYZ83EsGULK9uSeD+++f+3d57BcV1pen5O54ycc2ZOYqYSJY0kKq5XE6SxPbteb01N2bu2y167Nnir/MehynG2vOXd2fU67XjkGY0taWYkKpKURFKMIAmCJEAARA6N2Dnf4x8XlEiiG2gSjQYavE8ViyDuube/CzTfPvc733k/ftr1UyodlQgh8Ef9zEZmkVJiM9jYWryVEd8IHw58yHc3fPeBmrJrrC6a0OcYM4Eo4ViCqnmvlrtIxGFuAMoezAZ4OUQHBvAfO46prhbXiy8i9Pe/uWnAO4BRZ6TC8WDdpzJJQ7Gd7gkfMyKfotu9d61rp4HLctGZTJhbmqGnl289+W3aZ64wEZwAwGVy0ZDXQIW9gmpnNTqhYzY8y5s33qRzupM95XtWOXqN+0UT+hxjbN68rDKZ0M/egkQMiluyGlPC68X74YcYigpxPf/8A4k8wJBviCpHFQbd6r8t64ttCAEDEQdFoO5NWEdCD2BuaiJ87Tqm8Vkeb0y6Yf0rCiwFVDuruTZ9jd1lu1c1taZx/2irKznGmCeMxainwJakdHLyBhgtUFCftXhkPI736FFQJK4jRxAP0DkKwBPx4I16l91UJFPYTAbKXBZuBu1qO0bv8NIn5RjGmhqExUyktyet8c35zQRiAabD0yscmUam0YQ+xxj3hKjIsyycUSkJmLoJRS2QYU+YxQieP098wo3zG8+gz3/wGe+wXxXSaufa6YhVX2RnzBcjZi8Dz/oTeqHXY25sJNp3CxlfevdvnUtdOxnwDiwxUmOtoQl9DhGJJ5gORClPVm0zc0vtKlWSOYfHpYjPzBC8cAHzhjbMjY3Lutawbxi70U6+ee2kR+qKbEgJ46IUfBMQDa52SBnH3NS0wA8nFTajjWJrMcO+9feht97RhD6HmPRFkPKO5iJ34r6mpm0K07MZyATBs2cRBiOOQ4cWHRdX4ihSSXlcSsmIf4RqZ/Wayv2WuyyYjTpuUQlSgembqx1SxjFWVyNMJiJ9fWmNr7BXMBGcWPT3qbH20IQ+h3D71GbWpffuhk3EYKobituylraJz84S6enFum1ryo5RilQ4PnScv7jyF/zV1b/ikvtSUsfEqdAU4XiYasfaSduA6llfW2ijO2BHWvJgonO1Q8o4wmDAVF9PtK8PqSwt3uX2cuJKnOmQlqfPJTShzyHc3ggOswG7+Z6qlOleVexLN2Ytlsj16yDAum1byjEnR05ybfoaG4s2UmGv4NToKb4Y+WKB2N/y3EIg1sxC7J3UFdrxRRJ4C7eom9Fm119+2tzUiBIKExsdW3Ls7R2044HxlQ5LI4NoQp9DTPrClLrMSQ5cB5Md8rOz0UhKSbi7W7XBtSe3EHYH3VydusqW4i08WfMkLzS8wPaS7XRMdXBt+tpdY/s8fVQ4KtbkRpza+a5dfaY2tbzy+rswt77sfU21tQiDnmhf75JjHUYHdqOdscDSHwoaawdN6HOEWEJhOhBdaGKmKOpCbFEz6LLz64yPjqL4/Jhb21KOOT9+HrPBzL6KfYDqg36w8iA1zho+H/mciYC6OWciMMFMeIbm/OasxH6/5FmNFNpN9M/FYMs3Vfvn9r+Gy2+qC7TrAGEyYaytJdLbt2QzEiEE5fbyrzZXaeQGmtDnCFN+dSG21HnPQqxvVK22KVxe1cv9EO7qRhiNmBvqkx6fC88x4B1gS9EWzPqvP5iEEDxT9wx2o52j/UfxRX2cnziPSW+itaA1O8E/ALVFNkZmQ8StRbDnt6HpsLqB6uL/gMnu1Q4vI5ibmlD8fuLupQ3cymxl+KI+ArFAFiLTyASa0OcIbu/8Quy9qZuZPtWDJUsmZjIeJ9LTg6mxIeXmqI6pDoQQbCleaMVgNVg50nCEmBLjf137Xwx4B9hdthuT/sE2WmWDukIbsYRkdC4MBjPU7oe93wdHKVx/BwJTqx3isjHV14NOEO1dOn1zO09/+6lMY+2TltALIZ4XQnQJIXqEEL+f5PgGIcRpIURECPF793OuRnq4fRGsJj3OexdiZ/pU73ljEkuEFSA6OIiMRLC0Jp+BJ5QEN+du0pDXkDLnXmwt5rWW19havJUna55ke8na7k1TXWBDrxMMzNwxgzXZ1FSOzgDdH6imZzmMzmLBVF1NpKd3yfRNsbUYndBpC7I5xJJCL4TQA38KHEFtD/iGEOLerhAzwD8A/t0DnKuRBm5fmFKn+e4682hATSFkMW0TuXEDnc2KsTa5Xe2gb5BwPExbQer8PajeKY9VP8amok1rqnY+GSaDjsp8K/3T92yYMjug4QmYG4Tp9GwE1jKmxiYSHg+JmcUbrRh0BkptpVqePodIZ0a/F+iRUvZJKaPAm8Crdw6QUrqllOeA2P2eq7E08YTCtD+6MD8/26/OJLMk9EooRKS/H3NrGyLFwm/3bDdWg3VNlkouh7oiG1O+CP7IPVYBFTtUS+iBUzk/qzc3NoAQRNJI35TZynAH3SSURBYi01gu6Qh9FXBnPdnw/PfSYTnnaswz4YuQUORC64OZPjVl40jdHSiTRG7ehISCZWNym4VIIkK/p5/m/Gb0WfTbyQZ1hWoaavDeWb1OB7X71Ibsc4OrEFnm0NntGCvKiaaxS7bcXk5CJpgK5f76xMNAOkKf7Lk63alL2ucKIb4vhDgvhDg/Obm+Ovosl9G5EACV+XcIvZRqWWVBfVbKKqWiELrSgaG0FENxcdIxA94BEjJBS0F2bZKzQYnTjM2kZ2A6SaVJ+TY1Zz90JvuBZRhTYyPxySkSc3OLjrvd6nE8qOXpc4F0FGIYuPM5vBoYTfP6aZ8rpfyRlHK3lHJ3SUlJmpd/OBidC1FgM2Iz3bEQ63erOfospW2ifX0kZmex7Urd+/WW5xY2g21V+72uFEII6opsDM4EFy5W6o1QtVvdoezP7UmKuakJgEjfrUXHOUwOXCYXI76RbISlsUzSEfpzQIsQokEIYQJeB95N8/rLOVcDUBS1rG9Bo5HZ+f+IWTAxk9EogZMn0RcWYJoXgnuJK3EGvYPU59Wv+cXVB6W20E4wmmBy3nPoLip3gt4Aw2ezH1gG0btcGEpK0tolW+uqZdg/TEy5d2lOY62xpNBLKePA7wAfANeBn0opO4UQPxBC/ABACFEuhBgG/jHwz4UQw0IIV6pzV+pm1iOjnhDhWIKG4nusBqZ7wFECZueKxxA4d46E14fzySdTLsKO+EeIKTEa8rLnnplt6m7bIUwlSd+YbFC+XTU+i/iyHFlmMTc1EhsbJ+FffENUvaueuBJn1J/uA77GapFWcldK+Z6UslVK2SSl/Jfz3/szKeWfzX89LqWsllK6pJT58197U52rkZxwLIEndPfsqG8ygF4nvvJcUQd61UYYxYuXMGaCmNtN6NIlLJs3YaxKvY7e7+nHqDNS5Vi/a+12s4GqfCs3J1IIefVu1c545EJ2A8swt5/alprVVzoqMelN9MzmfmnpekfbGbtGGJwO8pef9/FXX9zig85x4gmFeELh+piXuiIbZsMdVSzu6+pibOnKbkmQioL/2HF0Fiv2gwdTj5OSW55b1Lpq10S/15WkpczBlD/KTCC68KCtEIpbYeQixJMczxEMhYXoCwqI9C5efWPQGWgtaKVnrodwPJyl6DQeBE3o1wCxhMLRzjFcViO76gq4Nurl/14c4VjXJMFogl21BV8PTsRh5Dzk14C9aEXjivb2Ene7sT/6KDpLkmYn80wEJwjGg9S76lc0nrVAc6kDgO5Us/qafar30NjlLEaVeczNTcRGR1CCi3fV2lS0iYRM0DHVkaXINB4ETejXAFdHPAQiCZ7aUMoTrSW8uK0Cty/M1REP22vyqJmv4UZK6P1UTd3ULd7VablIKQlebEefn4+5ZXFnyVueW2pViis7fjuridNiVNM3bn/yAXlVkFcNw+dUZ9EcxdzSAookfO3aouOKrcU05jVyyX2JYGz9tVpcL2hCvwa4OuqlPM9CdYEq6K1lTv7uo41870Adh9tK1cbfIxfh7I/U/G/1nhWvtomNjBB3u7Hu3JlyAfY2tzy3qHJUYTGknvWvJ1rKHEz5IsnTN6CanoU9ap+AHMVQVISptobQ5StLNg7fV7GPhExwavRUlqLTuF80oV9lZgNRpnwR2srvrp6xmvQUOcwIJQEdb6nGWQYzbHoVmp9e8bhCFy+is1mxbFh8wXc2PMtcZI4G1/qttrmXJdM3Rc1gK4LBL3PaFsG6YwdKMEi4q2vRcQWWAnaV7qJ7tpsh7/pqyrJe0IR+lbmdAmiZF48FDHyhWh20Pge7fgPKNqm2xCtIfGqK6MAg1m3bEIbFF1f7vf0A1OfVr2hMa4kl0zdCqLYIfrfqR5SjGGtrMZSXETxzFhlbvFZ+V9ku8sx5nBg+odXVr0E0oV9luid8VOZbcFqMCw9GAzB0Fsq3QNWuFRf42wTb2xFGI5atW5cc2zvXS7G1GKdp5ev51xJLpm9KN6vuljlsiyCEwHHoEEogQOjSpUXHGnQGnqx5Em/Uy/nx8wuOeyIeftH7C358/cd0TmtbabKNJvSryFwwyqQvQnNpCpGc6FTz8zX7sxZTYm6OSHc3ls2bFq20AZgJz+AOutd0d6iV4nb6JmVNvd6g7paduQWhxX1j1jLGykrMTY0EL1wg4fEsOrbKUcWGwg1cmrx0l9lZJBHhV32/YiI4gVlv5sTQCbpmFk8HaWQWTehXkduP/s3J0jZSqiV6rkp1B2yWCJw6hdAbsO3ateTY69PXEUI8lEK/ZPoGoGyL+hQ2ntulh/bHHgME/hMnlmxKcrDyIBa9heNDx1GkgiIVPh74GE/UwwsNL/DrLb9OpaOSE8MnmAvn7gdgrqEJ/Spyc8JPeZ6FPGuStI1vTG1RV7Eta/HERkeJ9PZh3bUTnd2+6NhgLMi16Ws05TWl7CS13mkqdTDpi+AJpshJW/NVd9HxjpxelNU7ndgP7Cc6MEik++aiYy0GC4eqDuEOujk1eopjQ8cY8A7wWNVjVDoq0Qkdz9Q+g17o+WTwExSZuyWouYQm9KuEJxRjwhtOvQg7dkV9/C/ZmJV4pJT4v/gCnd2ObWdqh8rbY0+NniKuxNlTvicr8a1FmkvU313P5CLeNuVb1VLLHF6UBbBs3YqhrJTAF5+jhBffBduS38KGwg1cmbxC10wXu8t239U/2GFy8Hj140wEJ7gwkdt2EbmCJvSrRI9bFYeWZPn5RAzcnVCyAYzZqU2P3LxJfMKN/cB+hDHJE8Y83qiXt3vepnu2m93luymwFKQcu97JsxkpcZrpdS9i/lXcCgaTut6SwwidDufhwyjhMIFTi9fLCyE4XHOYb7V+izc2vMHeir0LxrQUtNBS0ML5ifO4g+6VCltjHk3oVwEpJZ3zm6TybElEdbJL9Uopz07aRsbjBE+fxlBSjLktdd38zdmb/LTrp0yHp3mq9il2l+3OSnxrmeZSB6OeEIF7WwzeRm9Un8omb+S0/w2AoaQE644dhDuvERtZ3IdeCEGJrWTRicBjVY9hM9j4eOBjrSRzhdGEfhUYmQsx7Y+ytSov+YDxK2p+Nz95A+5ME7pyhYTXh/3QoZS7YC9MXOCjgY8oshTx7bZvs6Fww7r1nb8fmkocSKm6jKakbLP6lDbVnb3AVgj7nj3oXU58x44vuWN2KSwGC0/VPsVcZI6TIyczFKFGMjShXwXO9c9gMeoX7IYFIDQLswPqbD4LQqqEwwTPX8BUX4epJnlD7xszNzgzdobWglZeaXoFl8m14nHlCsUOE/k2Y+pdsqB+YFvycj59AyBMJhxPPEFidpbgxYvLvl6Ns4adpTu5Nn2NGzM3MhChRjI0oc8yA9MB+qeC7G0owKhP8uMf71AFvnzLwmMrQOjSJWQkgv3AgaTHR/2jHB86TrWzmsM1h9dd0+/lIoRgY4WLwZkgc8EUqRkh1Fn97K2cb0oCYKqvx9zSQvD8eeKzs8u+3r6KfWrJ5dAJZsPLv57GQjShzyLBaJwPOycosBnZVp2/cICiqNU2BfXqDHCFUUIhQpcuY25pTtrw2xPxcLT/KC6Ti+fqn9NEPgVbqvLQCcGloUXqwsu2qCWWE4u7QeYKjsceRej0BM8sv3WiTuh4tu5Z9Do9nw1/tmStvsb9k5bQCyGeF0J0CSF6hBC/n+S4EEL8yfzxK0KIXXcc6xdCdAghLgkhFu6NfkiQUvJB5zjhWIIXtlUkn83P9KozvsrFyxvvRQmF8J88iffoUSK9S/f6vE2ovR0Zj2Pbu7AqYiY8wzs97yCl5MXGFzHrzfcV08OEw2xgQ4WTjmFP6pp6exG4KmAitzdP3UZnt2PdtpVIT09GZvU2o439FfsZ8Y9wc27xWn2N+2dJoRdC6IE/BY4Am4A3hBD3tjY6ArTM//k+8F/uOX5YSrlDSvnQlmlcHJylfyrIE20llDpTlEyOtqv+KEUtaV9XiUTwvP02oUuXiI2O4X3vfQJfLu2vogSDhK5cwdzSgqGwEFA/jCYCE5wYOsHPun5GQiZ4pekV8swr/3SR6xxsKkKnExztHCOWSLEJqGwr+CfBO5bd4FYI644dCIOeUAZy9aA2MSm2FnN69DSxhFaFk0nSmdHvBXqklH1SyijwJvDqPWNeBf6nVPkSyBdCVGQ41pzF7Qtzsmea5lJH6kqb4IzqUlmxA5bwf7+TwOnTxKdnyHv5ZQp/8zewbN5E8Ny5lNayilTom+vjwkc/ZtQzzERrEX2ePr4c+5Kf3PgJP7/5c7pmu2guaOY7bd+hxJY9+4Vcxmkx8uymMsY8Yd5uHyESTywcVLZZ3RfR/3n2A1wBdDYb5rY2It3dS26iSut6QsdjVY8RiAVod7dnIEKN26TT4LMKuNNkehjYl8aYKmAMkMCHQggJ/LmU8kcPHm7uIaXk2A03FqOOZzaWpS5JHDwNQn9faZv4zAzhq51Yt23FVKuWYjqefJLE7Cz+Y8cxlJZiKPi6jjmmxDh66yij7l5qLlzGX1PIlOcMeNRFxSpHFdtLt9Oc36ylah6AljInz0vJB1cneOfSKH9jZ9XdKTqjRTWo6zsOUz1QvHjnrlzAunUr4audhK9fX3JHdTpUOCpoKWih3d3OhqINWoVXhkhH6JMp072rJYuNOSSlHBVClAIfCSFuSCk/W/AiQnwfNe1DbW126sezQe9kgNG5MM9sLMNqSrGYGZqF8atQ9YiaukmT0JUrCL0O256vbQiETofz2WeZffNNfB98SP43X/vKU/7z4c8Z9g3z2HgepaWP4HjtOwStAkUq5JvzMelNy7pXDdhQ7kIgeP/qGJ/fnOSpDWV3D6jeo+567voVOH8LzLlt72woLsZYWUG446qayslASfD+iv3c8tzi9Ohpnqt/LgNRaqSTIxgG7iywrgZG0x0jpbz9txv4f6ipoAVIKX8kpdwtpdxdUrJ+0gUXBmbItxnZXLnIzGTgNAid2qwiTZRolMiNLswtLeis1ruO6Z1OnE8/Q3xyksBJdSNKn6ePGzM3eESpoWw4gG3nTmyFJRRbiym1lWoin0Hayp3srC3g8pCHMU/o7oN6A2x8FRJR6HxbtaHOcSxbtpLweIgNZaa7lNPkZGfpTnrnehnxL74DVyM90hH6c0CLEKJBCGECXgfevWfMu8D35qtv9gMeKeWYEMIuhHACCCHswLPA1QzGv6aZ9kcYnQuzrToPnS7FTCfshYmrqkvlfczuIl1dyFgsZXMQc2MD1h07CF3pwHfyJKeHTlIaNNJwYQx9fj72PQ+vGVk2ONBYhNWk53Tv9MKDjhJoewE8w9B3LPvBZRhzUyM6m23J5iT3w47SHbhMLj4e+Bhv1Jux6z6sLCn0Uso48DvAB8B14KdSyk4hxA+EED+YH/Ye0Af0AH8B/L3575cBXwghLgNngV9JKY9m+B7WLNfGvOiEYEP5IrP54bNqfXVN+rN5KSWhjg41B19amnKc/dBBLJs2MvDF++T/7FO2n5pAZzTieulFhEmbwa8kJoOO3XUFDEwHF87qQV2Yrd4NQ+dUb6McRhgMWLdtJTowSHxqaukT0sCoM/J8w/PElBhvdb/F9enrxJXlWS48zKSTo0dK+R6qmN/5vT+742sJ/P0k5/UB25cZY06iKJLrY17qi23YzSl+zLEQjF6C0o2qt02axEdHSUzP4Hz6qUVzokKnw/jEIS4n2imbcVHZdBjrli0LUj0aK8O26nzO9c9yrn+WV7Yn+Zk3PQVzg9DzCRQ2qWmdHMWydSvBCxcJnP5SnUhkIFdfbC3mtZbX+HjwY44NHePk6Ema85tpLWilwl6heS3dB7n7zlrj9E8HCEQSi+fmRy6oZle199cqMNRxFWE2Y25Zut7+6tRVPKU2nn30W9i1UsmsYjLo2F6Tx5m+Gab9EYoc91Qy6fSq2F9+U30v3McazVpDZ7Fg27uXwMmThC5dwrpjB8TjJLxeZCSCsFrROxxgMKAEAsRGRomNjqjpx02bMVVXJb1ugaWAb7Z8k2H/MF0zXXTPdnNt+houk4snqp+gxpXcn0njbjShXyGujXmxmvQ0FKeooolHYfg8FDWDI3X65V4Sfj+R3h6s27Yv6hsPEEvEuDJ1hTpXnVYPv0rsqMnn4sAsFwZmeXZz+cIBhQ1Q2AiDp6Bie9b6D6wE1h3biY2NEvjiJMEvv0Qm2UsgjAZkTE3BCJMJodcR6erG8eQTWFOsNwkhqHHWUOOsIZaI0efpo93dzi/7fslzDc/RmNe4ove1HtCEfgUIRRP0TQbYVp2HPtUi7PgVNXVzn7P5cEcHSLBuS/6f4k46pzsJx8M8UvbIfb2GRuawmQxsrsqjY9jD/qYiXJYkH86NT8D5/6au1zQ8nv0gM4TQ6XAdOaI2sXFPorOY0bny0FktKMEgis+HEo6gdzkxlFdgKCmGRALvBx/iP/EZ+vz8lA6qtzHqjbQVttGY18g7ve/wycAnFLUVabu3l0AzNVsBro97SSiSzZUp3nxKAobOQF415Kf/6KmEw4SuXsXUUI8+b/E3dlyJc3nyMpWOSsrtSWaSGlljV20BUkL7YArTM2c5lG6AobMQXcTXPgcQOh2WtjYcjz2Kbc8eLG2tmGprsWzYgG3PHhyPPYp1+3aMZaUInQ5hNOJ69hvoC/LxffgRSjCY1usY9Uaer38eIYRmhJYGmtBnGCkl10a9lLrMlDhT7C51X1PLKmuTWwOnInjuHDISxb5v6Vxu12wXgVhAm82vAfKsRtrKnVwd8RCKpqibr38clDj0f5Hd4NYAwmTC9dxzKJEwvk+PpS3aDpODfeX7GPIN0TuXvpnfw4gm9Blmwhth0hdJPZuXEga/BHsxFDWlvM5seJb3b73PW91vcXXqKuGeHkKXLmPZsjmppfCdKFKhfaKdUlsp1Y7q5dyORobYXV9APCE51uVOLmT2IqjaDSMXYfCM2lDcNw5TN1UPpHWwsWoxDMXF2A8cIHrrFuHO9K2cNxdvpthazKnRU1o7wkXQcvQZ5srwHCaDjo0VKTY/TfdAYAo2vpyyg9RUaIq3e95GICgI6bj22X8nPm2mrmE7jkOHloyhZ64Hb9TLwaqDWgnaGqHYYeZAUxEne6awmfQ80Vqy8HfTdBjCc9D7qfrnTuzFsOnV+1q4zzWsO3YQHRgg8MXnGKsq7/JpSoVO6DhUdYh3et7hkvsSe8q1jYDJ0Gb0GSQUTdA94aOtzInZkMTXRkrVvMySB6X3Oj2rRBIR3ut7D6svyku9eTx6xk/dlKCrNE7wmb1LVtpIKbk4cZFCSyENroZM3JZGhthTX8CO2nzaB+d49/LoQodLnR62vAbbX4e2I7Dl1+GR34DNvwbxMLT/NcxlxmZgLSKEwPnMM6DX4/vwI2QivaeYKkcVTflNtLvb8UVzv4PXSqAJfQY52z9DXJHsrE2x+WluEDwj6i7YFFbEZ0a/RHetl4Nn/ejds9j37WPH7/wx0X1bODl5FkWm8Dqfp8/Tx0x4hl1lu7TZ/BpDCMGTrSUc3lBK/1SQN88OMROI3jtILbms3AElbeCqVDfU7foemBxw5f+o76N1it55yRnHAAAQjklEQVThwHn4MHG3G+/Ro2k3ID9QeQApJSdHT36VGpNS4g66uThxkc6pTgKx3F7oXg5a6iZDDM0EaR+cZXNl3sKNMaDO5vs/V90pK7YlvcaEd5SJ93/BhklB4fa9OJ9+6qtdrAcqDvDhwIdcn7nO5qLNSc+XUnJh4gJ55jya83PfAnc9IoRgR00+RXYT73WM8ea5Qd7YU0uBfQlLCkse7PguXP6JKvZbvw0FddkJOsuYm5txPPE4/hOf4X3/fVxHjnzlwJoKl8nF7vLdnBk7w1Gptr/sn7tFrL8fQyBCuMTFZwUOGvIa2Fy8mWpH9UM1EdKEfplE4wqXh+c40zdNgc3E460pFkpn+oj1dxH0FBPv/jF6Vx6mulpMjY0YCgqIzkxz5X//kAL3HM1HfhPX/kN3vRGb8puomKrgzNiZlH7xvXO9TIWmeKr2KXRCe1hby9QU2nh9Ty0/OTfILzvG+O7e2tR7Lm5jdsD2N1Sx7/gpbP2W2l94HWLdtg10evzHj+P55S/Je/HFJdOWu0p3IaVUm5ZEorSeGaUkaCDfXEh8JM5EtaSjdZg+Tx92ox2T3kQ0EUUgaC1s5ZGyRzDqFn+NXEWsxfrT3bt3y/Pn13Z72YQiOdc/Q/vgHOFYgsYSO99oMGHzDaplcq4KcFWrKZpokMi7/x7fxVuIpkMYa2pIzMwSn5wE1E49Y1O36A+P0PbS36TlkaeTvuZkcJK3ut9ia8lWHq169K5jMSXGT67/BLPezLfavqUJfY7QN+nnnUujHGouZm9DYXonRQNw6X+rC7dbvqmmetYp4evX8X3yKaa6urQ9dJRYDM+775Jwu3E89TTGqirCVy4TvHARfXUl7oOtDIRGkEjMOjPhRJhbnltUOip5sfHFnBV7IcSFVO1atRn9AzAXjPL+1XHGPWGai80ccM5R7DsHF+/JnRotyIIGgucvEOzowbj9aVyvvYHOZgMg4fMR7evDNzpAp2MA++YXaN78VMrXLbGVsLFoIx1THbQVtN1la3By5CT+mJ9n6p7RRD6HaCxx0FLm4EzfNG1lTvJsaYiMyf51GqfjLdj6mmqjsA6xbNyIjMXwn/iM4Llz2JM0sr8TKSX+T48RHx3D+dyzWFpbAbAfPIi+oADfJ59S/oWk9aWX0Jm/firumuni08FP+WTwE56rey6tD5TZ8CwD3gHyzfnUuerWdCrooRT6208x9/uLkVLSOerlRPckjugUr7sGqPANwVwMrAXIhsehdBPCYIa5QZSRTvwffkJkwo/50Es4X/7OXblGvdOJZds2jjlH8PireL7lG0vGtK9iHwPeAd6/9T6vNr+K0+TkwsQFrk1fY2fpTiodlff/A9FYVZ5oLWFgOsjxbjev7khu7rUAkx223xb7n6sVOovsy8hlLFu3EhsfJ3j2HMaKikVtEoJnzhDp7sZ+YP9XIv/VdTZuRBgMeD/6iLmf/5y8l19G71TLoNsK2wjFQ5waPcVZ81n2VaTelJhQErS72zk/cf6r4oh6Vz3P1j+LQbc2JfWhSt1E4wonuie5MeZFrxdsq8pnf2MhBv3SM+ChmSCn+6YZmQnwSOIS+ww9mC02KN1INF5I8MYwsYkJUCQ6ixlhsaL4vEhFwX7gINadydustbvbOT16moOVB9lRuiOt+3AH3fyi9xfElBgWvYVgPEhrQauWm89hLgzM8Fn3FK/uqKSxJP12kkSDcOVN8LtVn5ya/ffVXD5XkNEosz/7GTIcoeD176Cz2xeMCV+7hu+TT7Fs3oTj8OGUk6bo0BDe994HIbDv3YO5tRWdzYaUkuNDx7k+c52na5+mrbBtwbmTwUmODR1jKjRFU34TByoP0DfXx6nRU7QVtvF0bfK0azZYLHXz0Ah9KJrg7UsjTHjDbKnMIxJX6J7wUZ5n4cVtFcnNpoChkSnaL3Yz5o1iKbDzpPkytcKNqNpFzLGJ4Pl2ooND6F1OTE3NCIMeJRxGhkLoHE4smzdhKFyYe40rcdrd7ZwbP0dTfhPP1j17X08YvqiPjskOQvEQda46mvKb1vSjo8biJBTJj88MEE9IvnegLq3Jx1fEo9D1Hrivqxurqh5RSzNNC8Uwl4lPTzP3s59hKCsn79VXEHd8oEX6+vC+/z7GqiryXn4ZoU/Rn/n2tWZn8R87TmxEbVUojEZ0VgvSZKI91sdAnYXnH3mdKof6hBWKhzg/fp6r01ex6q08Xv04jflfp8vOjp3l/MR5nqp9ig2FG1bg7pfmoRd6TyjG2+0jeEMxjmytoLlUnTH1uH180DmBXid4bnM5DcV2pJTEJiYY6eii79INgqPjGPWCaqegJNiLXhfH0PQIiqWUyPQkMaPAtXc/eTseWVACFklEGPYNM+QbYi4yR0yJoSgKCZkgEAsQU2K0FrRyuOYwet3ib0yN9c/QTJC3Lgyzv7GIA01F93eylGqnqv7P1Z3XQkB+rSr4pZtz2v74Tr6atW/ZjOPxx0GnI9x5Df9nJzAUF5P3a7+GLs3uaVJK4pOTxEZGUPx+dYIWDhMaHqRz4grupiJse/cgDHpG/aNIKdlcvJm95XuxGO7+ecYjYY5e+AnToWle2PFtSopqV+L2F2XZQi+EeB74IaAH/lJK+W/uOS7mj78ABIHflFJeTOfcZDyI0CuK5K0LQ/hnvJiCXvIMkGcxYLRZ6ZyLE7PYeHlXDdUFtrvOmwlE+dWFfkL9g1T5pyiYHSfk9ROIKsiSEmq3ttBWHiM+fBLvXICArRFfMMxM3MtkvsBbX4I06im0FFLpqCTfnE8gFmA8MM54cBwpJSa9iSJLEUa9EYMwoBM6LAYLjXmNVDs1LxqNrzl6dYwb4z5+bUcV9cUPMCOXUk3jTN5QhT84rXauKt0EFTvUDVg5/OQnpSR4+jTBCxfVoga9DsXnx1Rbg/O559BZlv+BpkQizH12jIHzx5m2K/gObKaksolNRZsotHz9dC6jUSL9/UR7eogODBCNhrk23YlO6NhUtwdbZS3GinLVkrm4aMmnjOWyLKEXQuiBbuAbwDBqs/A3pJTX7hjzAvC7qEK/D/ihlHJfOucm40GEXiYSfPGv/gRjLIKUknBMIRRLkFAkeVYDDQVGbDYdepsJvdWIzqQnHgoyOznO3PQ07kCQWSmZKXAQKcvHXOOgyJ5ACYzjiwaIWZxQ1AJGCzqho8xWRrWzmkJLIXOROUb9o0wEJ9S6XCEosZZ81SyhzFamzdg10iIaV/g/54fwBKN8Y1M5rWWO5aXkfOMw2g4TnWo3M0ep2q/WUQqOspxN70T7+wl3dYNUMDU0YG5puSuVkwkifX34P/0UJRrF3NiEsbIC9HoUn5+4e4Lo8DAkFHR2O+bmJky1tbjDU3zW8UvyfQk2xkqwRwWKVAgqEfxmBY85QcRqIO6wgNMBTjvku5A6QUyJYdQZOdJw5IHiXa7QHwD+hZTyufl//wGAlPJf3zHmz4HjUsqfzP+7C3gSqF/q3GQ8aOqm+6Ofo1hNSJcDZBglMEliehjd7DjS40WGYhCMEg1GCIcjeA0KYZuecL4BY4kTU74d3Z25UZ0Rg6MUZ+lmXIUtuMx5OE1OXGZX0lpbRSpEEhGMOuOaXX3XWPsEInHevTzKuCdMsdNMbaENl8WAQadDrxMPtNYqEhHM0zewTF7BEJz86vuK0UbCUqD+MbmQejNSZwBxx8Qkh58AlosMhVGudyMHh5GRCDBfredwICrKEFXliOKiuz5k3JFpTs6eJ5QIYwonMM8GsHrCGINRnBGBLSzRxeJIKZFAQg+RfDux4jx0ZaV865V/9kAf7suto68C7nRSGkadtS81pirNc28H+X3g+wC1tQ+W3zoe+5B4OArTUVDmPWH0JijPg/pq1StEb8RgtJJnLqDCWcUGRzUV1mKsEtU4Sirq46/ZoY6/jx+4TuiwGrTG2xrLw2428J3dNXSOerk+5qVjeI5YIhNraWXANzAYw9hi09ij09gCM1jiHiyxQYxKKAOvsU4pkxBLIBSJNOlBF4PILPTdgL6Fw3ejMI4fn4hikDry7Gby7Q7MtyU3lkBEYohQDJ0/hG7Mg+hxg3UYXsl8+OkIfTKlu/ddl2pMOueq35TyR8CPQJ3RpxHXAr5dfRgQCJMVYS9FuCrRWQsRQocQAoFACIFJZ0rxiblII28NjSyi0wm2VuextTrvq1RkXFFIKJIVq59Q4pCIIhIx5ALzvHtedA0Wcaw1lm72eTcyEgG/f0Wq59IR+mHgzh0K1cBommNMaZybMfK3vb5Sl9bQWDWEEFhNetR6hpXEBNiWHKWRe6ST7TsHtAghGoQQJuB14N17xrwLfE+o7Ac8UsqxNM/V0NDQ0FhBlpzRSynjQojfAT5AnVL8lZSyUwjxg/njfwa8h1px04NaXvl3Fjt3Re5EQ0NDQyMpD8WGKQ0NDY31zmJVN+vPFENDQ0ND4y40odfQ0NBY52hCr6GhobHO0YReQ0NDY52jCb2GhobGOmdNVt0IISaBgWVcohiYylA4axHt/nIb7f5ym7V6f3VSypJkB9ak0C8XIcT5VGVG6wHt/nIb7f5ym1y8Py11o6GhobHO0YReQ0NDY52zXoX+R6sdwAqj3V9uo91fbpNz97cuc/QaGhoaGl+zXmf0GhoaGhrzrGuhF0L8nhBCCiGKVzuWTCKE+LdCiBtCiCtCiP8nhMhf7ZgygRDieSFElxCiRwjx+6sdTyYRQtQIIY4JIa4LITqFEP9wtWNaCYQQeiFEuxDil6sdS6YRQuQLId6a/793fb7Nak6wboVeCFGD2pR8cLVjWQE+ArZIKbehNl//g1WOZ9nMN5L/U+AIsAl4QwixaXWjyihx4J9IKTcC+4G/v87u7zb/ELi+2kGsED8EjkopNwDbyaH7XLdCD/xH4J+RonVhLiOl/FBKGZ//55eonbtynb1Aj5SyT0oZBd4EXl3lmDKGlHJMSnlx/msfqkhUrW5UmUUIUQ28CPzlaseSaYQQLuBx4L8CSCmjUsq51Y0qfdal0AshXgFGpJSXVzuWLPBbwPurHUQGSNVgft0hhKgHdgJnVjeSjPOfUCdX9zacXQ80ApPAf5tPTf2lEMK+2kGlSzo9Y9ckQoiPgfIkh/4I+EPg2exGlFkWuz8p5TvzY/4INSXw42zGtkKk3Ug+lxFCOICfA/9ISuld7XgyhRDiJcAtpbwghHhyteNZAQzALuB3pZRnhBA/BH4f+OPVDSs9clbopZTPJPu+EGIr0ABcnu+mXg1cFELslVKOZzHEZZHq/m4jhPgN4CXgabk+amTTaUKf0wghjKgi/2Mp5f9d7XgyzCHgFSHEC4AFcAkh/lpK+bdWOa5MMQwMSylvP4W9hSr0OcG6r6MXQvQDu6WUa9GE6IEQQjwP/AfgCSnl5GrHkwmEEAbUheWngRHUxvLfXS89hoU66/gfwIyU8h+tdjwryfyM/veklC+tdiyZRAjxOfDbUsouIcS/AOxSyn+6ymGlRc7O6B9y/jNgBj6af2r5Ukr5g9UNaXk8BI3kDwF/G+gQQlya/94fSinfW8WYNO6P3wV+LIQwAX3A31nleNJm3c/oNTQ0NB521mXVjYaGhobG12hCr6GhobHO0YReQ0NDY52jCb2GhobGOkcTeg0NDY11jib0GhoaGuscTeg1NDQ01jma0GtoaGisc/4/Da39CKhG40EAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"f = plt.figure()\n",
"plt.plot(y_d, g_n, alpha=0.5, label='Never takers')\n",
"plt.plot(y_d, g_c0, alpha=0.5, label='Compliers, $Y_0$')\n",
"plt.plot(y_d, g_a, alpha=0.5, label='Always takers')\n",
"plt.plot(y_d, g_c1, alpha=0.5, label='Compliers, $Y_1$')\n",
"plt.legend()\n",
"f.savefig(\"Y.pdf\", bbox_inches='tight')"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"x = df['X']\n",
"x_n = n['X']\n",
"x_a = a['X']\n",
"\n",
"distributions = [x_a, x_n, x]\n",
"x_min = min(df['X'])\n",
"x_max = max(df['X'])\n",
"\n",
"x_distributions = []\n",
"\n",
"for i in distributions:\n",
" x_d = np.linspace(x_min, x_max, 2000)\n",
" kde = KernelDensity(bandwidth=0.1, kernel='gaussian')\n",
" kde.fit(i[:, None])\n",
" \n",
" # score_samples returns the log of the probability density\n",
" x = np.exp(kde.score_samples(x_d[:, None]))\n",
" \n",
" x_distributions.append(x)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"g_xa = x_distributions[0]\n",
"g_xn = x_distributions[1]\n",
"g_x = (1 / φ_c) * x_distributions[2] - (φ_n/ φ_c) * g_xn - (φ_a / φ_c) * g_xa"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOy9d3TcV533/7rTNCozI2nUm1Vs2ZZs2ZYtp7mkEDuBkJCEDSGweULZnCwty/LAcmD5bU6ewz48nOwCgd2TXcpmC6RCSCAJGEh1Yty7LcmyutVG0mhmVKbf3x9fSZaskWbUrHZf5+Q4+n7v3O+VLb3nzud+Pu+PkFKiUCgUiqWPbqEXoFAoFIq5QQm6QqFQLBOUoCsUCsUyQQm6QqFQLBOUoCsUCsUywbBQD05LS5OFhYUL9XiFQqFYkhw9erRbSpke6d6CCXphYSFHjhxZqMcrFArFkkQI0TTZPRVyUSgUimWCEnSFQqFYJihBVygUimXCgsXQFQrFwhIIBGhtbcXr9S70UhQRMJvN5OXlYTQaY36NEnSFYoXS2tqKxWKhsLAQIcRCL0cxBiklPT09tLa2UlRUFPPrVMhFoViheL1e7Ha7EvNFiBACu90+7U9PStAVihWMEvPFy0z+bVTIRaGIhn8AHDUQ9EJCGqSsAkPcQq9KoZiA2qErFFPRWw8H/w1qfw/1b8OZX8L7T8LZX4Pr0kKvbskjhOArX/nK6NdPPPEEjz322FV7fmNjI7/4xS+ijnvrrbe44447rsKKZocSdIViMtztcPqXYLbCtk/Drv8NWz4BWRXgbIBj/wXnXoGgf6FXumSJi4vjV7/6Fd3d3VflecFgcNzXsQr6XD5zPlGCrlBEIhyG2tfBaIZNHwdLJuiNkFwApXvh2s9D4Q3QdQ7OvayNV0wbg8HAww8/zPe+970J9xwOB/feey9VVVVUVVXx3nvvEQ6HKSwspK+vb3Tc6tWr6ezsjDge4LHHHuPhhx9mz549PPjgg+Oe8fWvf513332XzZs3873vfY/GxkZ27txJZWUllZWVvP/++xPWdfjwYbZs2UJ9fT1Hjx5l9+7dbN26lb1799Le3g7AjTfeyDe+8Q12797ND37wA1544QU2bNjApk2b2LVr11z+FY5DxdAVikhcOgqeTii/G0yJE+8bTFC0C4yJcGEftJ+A3Mqrv8454q2aLhwe35zOmW6J48a1GVHHff7zn6eiooKvfe1r464/+uijfPnLX2bHjh00Nzezd+9ezp8/z1133cVLL73Epz71KQ4ePEhhYSGZmZk88MADEccDHD16lP379xMfHz/uGd/5znd44okn+O1vfwvA4OAgf/jDHzCbzVy4cIGPf/zj4zyn3n//fb74xS/y8ssvk52dzSc/+Ulefvll0tPTee655/jmN7/Jz372MwD6+vp4++23Adi4cSO///3vyc3NHfdmNNcoQVcorsTrhoa3wV4C6WunHptbCY5qaHgHsjZqu3jFtLBarTz44IM8+eST4wT3j3/8I+fOnRv92u124/F4+NjHPsbjjz/Opz71KZ599lk+9rGPTTke4M4775wg5pEIBAJ84Qtf4MSJE+j1empra0fvnT9/nocffph9+/aRk5PDmTNnOHPmDLfeeisAoVCI7Ozs0fEj6wK44YYbeOihh7jvvvu45557pvtXFDNK0BWKK6n7I0gJa26FaKljQkDhDjjxC+g8Azlbrs4a55hYdtLzyd/8zd9QWVnJpz71qdFr4XCYAwcOTBDi6667jrq6OhwOB7/+9a/5+7//+ynHAyQmRviUFYHvfe97ZGZmcvLkScLhMGazefRednY2Xq+X48ePk5OTg5SS8vJyDhw4EHGusc986qmnOHjwIK+++iqbN2/mxIkT2O32mNY0HVQMXaEYS3edlqJYeAPEp8T2muQCSEqH9lPzu7ZlTGpqKvfddx8//elPR6/t2bOHH/3oR6NfnzhxAtAyY+6++27+9m//lvXr148K42Tjp8JisYzu4gFcLhfZ2dnodDr++7//m1AoNHovOTmZV199lW984xu89dZbrF27FofDMSrogUCAs2fPRnzOxYsXueaaa3j88cdJS0ujpaUllr+WaaMEXaEYIeCF2t9BYhrkXxP764SAjDJwt4HXNX/rW+Z85StfGZft8uSTT3LkyBEqKiooKyvjqaeeGr33sY99jP/5n/8ZF9aYavxkVFRUYDAY2LRpE9/73vf43Oc+x3/+539y7bXXUltbO2Fnn5mZyW9+8xs+//nPc/z4cV588UX+7u/+jk2bNrF58+aIh6gAX/3qV9m4cSMbNmxg165dbNq0abp/PTEhpJTzMnE0tm3bJlWDC8WiQUo4+xJ010Llg2DNmd7rB3u1fPXVt0D+9vlZ4xxz/vx51q9fv9DLUExBpH8jIcRRKeW2SOPVDl2hAGh6Twu1FN80fTEHSEjVdvY9F+d+bQpFjChBVygcNdDwLmRtmN3uOqUIXK0QunqFJArFWJSgKxYFUkrerO7ip/sbqOnwRH/BXDHQDedf0XblpbdHz2qZipRVEA6Cu3Xu1qdQTAMl6IpFQXWHhxMtfQRCYfad7aBvcG7L6T1+D06vk5Ezo94BP3WdbkLnfws6I2y4B/SzzOJNLgChA+ekPXwVinlF5aErFgVHm5ykJZm4uzKP/9jfwNEmJ7esz5z1vP6Qnzdb3uRinxbbtpqs5MRt5nRDHKn9tXj7q1l308cxxFlm/SwMcZpFgGt+UtIUimioHbpiwXEO+HF4fJTn2kiKM7A2y8L5djf+4Oz8UaSU7GvaR72rnq2ZW7kx/0aCIcF/nXoJr76aG80X6AxbOebLm6PvBLDmgqddebsoFgQl6IoFp6FnAICStCQA1mdbCYQk9d39s5q3ureaZnczO3N3ck32NZTZyygw3UKqoRiT4STnQzXIgiqOtfQRDM2RAFtztEPRAcfczLcCeOmllxBCUF1dPXqtsbGRDRs2LNiavv/97zM4OBh1XGFh4VVzioyFmARdCHGbEKJGCFEnhPj6FOOqhBAhIcRH526JiuVOc88g9iQTtgTNByUvJR6L2TCrw9FAOMDB9oNkJWZRbi8HwBsIcb6tn1sLb+bG+DTqgm6arA48vkHqHLN78xhlJOXRrbzSY+WZZ55hx44dPPvsswu9lFFiFfSZMl+WulEFXQihB/4FuB0oAz4uhCibZNz/A34/14tULG+6PF4yrZc9M4QQlGZaaOweZMgfmuKVk1PrrGUwOMg12deMtvI62dKHPxhm26pktvj83Jq7k4DeTZP/LU5dap+T7wVzMpgStKpRRVT6+/t57733+OlPfzqpoH/wgx/k1CnNVmHLli08/vjjAHzrW9/iJz/5Cf39/dxyyy1UVlayceNGXn755dH7P/jBD0bn+eY3v8mTTz5Je3s7u3btYvPmzWzYsIF333133POefPJJ2trauOmmm7jpppsA+Ou//mu2bdtGeXk5//AP/zBhjUNDQ9x22238+Mc/ZmBggE9/+tNUVVWxZcuW0fU8/fTT/MVf/AUf/vCH2bNnT9R1zIRYDkW3A3VSynoAIcSzwF3AuSvGfRH4JVA161UpVgwDviADvhDplvEt3dZlWTja5KSuq5+NebZpzSml5LTjNGnxaeQkajvmQCjMiZY+itISSQ87wD/AmjW3kphgo7nnBfa1vMINpZ8k15IdZfYoCKHF0ZeaoF/4I/R3zu2cSZmw5gNTDvn1r3/NbbfdRmlpKampqRw7dozKyvE2xLt27eLdd9+lsLAQg8Ew6nO+f/9+PvnJT2I2m3nppZewWq10d3dz7bXXcuedd/KZz3yGe+65h0cffZRwOMyzzz7LoUOHePrpp9m7dy/f/OY3CYVCE3biX/rSl/jnf/5n3nzzTdLS0gD49re/TWpqKqFQiFtuuYVTp05RUVEBaG9K999/Pw8++CAPPvgg3/jGN7j55pv52c9+Rl9fH9u3b+cDH9D+Hg4cOMCpU6dITU3ln/7pn6Zcx0yIJeSSC4w9tm8dvjaKECIXuBuIbp6gUIxhxIM7PWm8oKdb4khNNFHd4Y74Oikl59vd/Ol8J8eaneMOUDsHO+n19rIhbcPo7vx8u5tBf4itq1K0QiKdAVJLyEnK4WPrPgpSz8/P/orq3mpmbYdhzYHBHggMzW6eFcAzzzzD/fffD8D999/PM888M2HMzp07eeedd9i/fz8f+tCH6O/vZ3BwkMbGRtauXYuUkm984xtUVFTwgQ98gEuXLtHZ2UlhYSF2u53jx4+zb98+tmzZgt1up6qqiv/4j//gscce4/Tp01gs0TOcnn/+eSorK9myZQtnz54dZ9N711138alPfWq0eca+ffv4zne+w+bNm7nxxhvxer00NzcDcOutt5Kamgowo3VEI5YdeqRKiyt/4r8P/J2UMjRVp2ohxMPAwwAFBQWxrlGxjOnuHxb0K3boQgjWZlk4cLEHtzeA1XzZZzwUlvz2VBv1jgFMBh3+YJiTLX3csyUPW4KRC84L6IWekuQSAMJhydEmJ1k2M3nJZqiuBnux1qQCWJ+RzbrEDyBDp3ij+Q2a3E3sztuN2WBmRozs8j3tkFo8szmuNlF20vNBT08Pb7zxBmfOnEEIQSgUQgjBd7/73XHjqqqqOHLkCMXFxdx66610d3fz4x//mK1btwLw85//HIfDwdGjRzEajRQWFuL1egH47Gc/y9NPP01HRwef/vSnAW3H/8477/Dqq6/yl3/5l3z1q1+d0MloLA0NDTzxxBMcPnyYlJQUHnroodH5QfM6f/3113nggQcQQiCl5Je//CVr14730j948OA4s6/priMWYtmhtwL5Y77OA678PLkNeFYI0Qh8FPhXIcRHrpxISvnvUsptUspt6enpM1yyYjnh8PiwmA2YjfoJ99ZlaTuW822Xd+lSSv50vpN6xwC7StP53I0lfHRrHt5AmBeOtuAc8HGx7yKrrKuI02tvEjWdHvoGA1QVpiA8beDrh/R1o3OaDDoKUlLINezg2uxrqXfV83zN81zqn+HB5oigu+coLr9MefHFF3nwwQdpamqisbGRlpYWioqK2L9//7hxJpOJ/Px8nn/+ea699lp27tzJE088wc6dOwHN8jYjIwOj0cibb75JU9Plwq67776b3/3udxw+fJi9e/cC0NTUREZGBn/1V3/FZz7zGY4dOzZhbWNtdd1uN4mJidhsNjo7O3n99dfHjX388cex2+187nOfA2Dv3r388Ic/HP2kd/z48YjffyzrmC6xCPphYI0QokgIYQLuB14ZO0BKWSSlLJRSFgIvAp+TUv561qtTLHsc/b4Ju/MRkhNMFKYlcGL4MBPgQH0PZ9vcXFOcytZVKQghyE9N4KNb8wiEJE8fOobT62F18moAgqEwB+t7SLPEUZKepHUX0hnAvnrcs/JS4+ny+ClL3cQ9q+9Br9PzSt0r1PTWTP+bMpohwa7t0BWT8swzz3D33XePu3bvvfdGbNq8c+dOMjMzSUhIYOfOnbS2to4K+ic+8QmOHDnCtm3b+PnPf866dWPerE0mbrrpJu677z70em3T8NZbb7F582a2bNnCL3/5Sx599NEJz3v44Ye5/fbbuemmm9i0aRNbtmyhvLycT3/609xwww0Txn//+9/H6/Xyta99jW9961sEAgEqKirYsGED3/rWtyJ+/7GsY7rEZJ8rhPggWlhFD/xMSvltIcQjAFLKp64Y+zTwWynli1PNqexzFYFQmH95s47thalcvzot4ph21xDPHmqhNNOCNd7AkUYn5TlWbi3L5MrwXofLyz+/9xL94Rb+7y1fxBYfz7sXHBxpdPKRLbkU2RPgz/+qHdZtHJ9Z29I7yItHW7l7Sy6FaYkEQgFea3iNtoE27ii+g3xLPtPi/G+gtwGu/+Ls/GHmkZVgnxsOh6msrOSFF15gzZo1C72caTMv9rlSyteklKVSyhIp5beHrz11pZgPX38ompgrFKD5qUg5MX4+lmxbPDvWpFHb6RkV8w+snyjmABlWE1lpbuJkJs8dbuOXR1s50uikIs9GUVricAMKd8Q+oRnWOISAdpcWGzXqjXyw6IOkxKXwp6Y/MRiYZgaCJQf8A+CLfKirmH/OnTvH6tWrueWWW5akmM8E5eWiWDBGM1ymEHSAqsJUVqcnEZYSe9LkY9sH2jGbQty/ZRttjjjcQ0GuLbZzTVHq8AOrQacH+8Rf7jiDHntSHO2uy5kpRr2RW1fdygu1L/B+2/t8YNU0Dg6tY+Lo5umlXSrmhrKyMurr6xd6GVcVJeiKBcPh8WEy6LDFG6OOTUk0RR3T4GpAL/RUZpdyTd4Vc0qppSumFGox7gjk2MzUdHqQUo5+ArDH29mSsYWjnUcps5eRkxRj84vEDO3Nw9MGGeuij1co5gDl5aJYMBweH2lJpojhk+kipaS+r548Sx5GfYQ3CPclrd9nxuQx4yybGV8gTO/AeOveyoxKEo2JHGg7EHuOut4ASRkq00VxVVGCrlgQpJRTZrhMF8eQg/5AP8W2SfK+O89p2S1ppZPOkW2LBy7H0Ucw6o1sy9pG52AnzZ7m2BdlyVHOi4qrihJ0xYLgHgriD4ZJT5ph8c4VNLgaEAgKbYUTb4YC0HUO0lZrnuWTkJJgxGzUTxB0gHUp67CarBzqOBT7Lt2arT17sCfG70KhmB1K0BULgqNfE8252qHXu+rJTsom3hA/8WbnWa0MP6dy4r0xCCHItpnpcE0s2dfr9GzN3Ipj0EGjuzG2RVmG4+2eJebrchXp6Ojg/vvvp6SkhLKyMj74wQ9SW1s7L8966623uOOOOwB45ZVX+M53vjMvz1lIlKArFoQujw8hwJ4U/bAzGn3ePpxeZ+RwSzgELYe0eHZydLuJLJuZngE/3sBEl8fSlFJscTYOtce4S09I1ewFVBw9IlJK7r77bm688UYuXrzIuXPn+Md//Ec6O+fYJCwCd955J1//+qRO4BOYL7vbuUYJumJB6O73k5Jgwqif/Y/ghb4LCERkQW87roU8inbFVOCTY4tHSq1I6Ur0Oj3bMrfR4+0ZbWk3JUJcjqMrJvDmm29iNBp55JFHRq9t3ryZHTt28NWvfpUNGzawceNGnnvuOUDbYe/evZv77ruP0tJSvv71r/Pzn/+c7du3s3HjRi5e1P5NHnroIR555BF27txJaWkpv/3tbyc8++mnn+YLX/gCAA6Hg3vvvZeqqiqqqqpG3Rwfe+wxHn74Yfbs2cODDz7I2bNn2b59O5s3b6aiooILFy7M91/RtFFpi4oFweHxkW2bffxcSkmts5acpBySTEnjbzqb4OIbmkHWFaX+k5Fpu1xgVJiWOOH+mpQ1HOs6xqGOQxQnF6MTUd6QrNnQfFDrYjTbJtTzyP5L++kemtvOO2nxaezI3THp/TNnzowabI3lV7/6FSdOnODkyZN0d3dTVVXFrl27ADh58iTnz58nNTWV4uJiPvvZz3Lo0CF+8IMf8MMf/pDvf//7gNbx6O233+bixYvcdNNN1NXVTbqORx99lC9/+cvs2LGD5uZm9u7dy/nz5wE4evQo+/fvJz4+ni9+8Ys8+uijfOITn8Dv9xMKzcyrfz5ZvD9himWLNxDCPRSgYpo+55HoHOzE5XNRmXFFfLzfAWd+CfEpUHZnzOX3kQqMxqITOq7JuobfNf6OWmct61Kj5JhbckCGNa9xW+7UYxWA5nP+8Y9/HL1eT2ZmJrt37+bw4cNYrVaqqqrIztaKtkpKStizZw8AGzdu5M033xyd47777kOn07FmzRqKi4vHtbe7kj/+8Y/j7HDdbveoMdedd95JfLx2LnPdddfx7W9/m9bWVu65555FWX2qBF1x1ZnMA30m1Dpr0Qs9xcljwi0+D5x+XktT3PgXYIxwUDoFOTYz1R3jC4zGUmQrIj0hnYPtBym2FWPST3EOYB1jpbuIBX2qnfR8UV5ezosvTnQJmep8Ii7u8s+MTqcb/Vqn042Lc1/57zZVrUM4HObAgQOjwj2WsXa3DzzwANdccw2vvvoqe/fu5Sc/+Qk333zzpPMuBCqGrrjqjHigp80yw8Uf8lPrrKUkuWTUKpdQAE6/qGW1VNwH8cnTnjfbFo8/GKbnigKjEYQQ7MzdyWBgkMMdh6eeLM6i/bfUOhhdBW6++WZ8Ph8//vGPR6+NeI4/99xzhEIhHA4H77zzDtu3b5/W3C+88ALhcJiLFy9SX18/wZt8LHv27OFHP/rR6NcnTpyIOK6+vp7i4mK+9KUvceedd462xVtMKEFXXHUcHh8JJj2Jpoke6NOhprcGf8jPxrSNly/W/Qk8HVB2F1iyZjTvSGy/vW/iwegIWYlZrLev51T3KRyDjqkntGSpg9EICCF46aWX+MMf/kBJSQnl5eU89thjPPDAA1RUVLBp0yZuvvlmvvvd75KVNb1/y7Vr17J7925uv/12nnrqKczmyc9rnnzySY4cOUJFRQVlZWU89VTkxmvPPfccGzZsYPPmzVRXV8+6GcV8EJN97nyg7HNXLj8/2ES8Uc89lXkzniMswzxT/Qxx+jg+Wjpsheu6BMf+C/KrYPXMO/BIKfm3d+opTktkT/nkQuINenm2+lkSjAncu+Ze9LpJ3qCa3of6t2HH30w7/DOfLFf73Iceeog77riDj370o9EHL3LmxT5XoZgrQmFJT79/1gVFtc7a8YehUsLFP0FcEhTunNXcIwVGkSpGx2I2mNmVt4vuoW5OOCJ/TAfGdDBSYRfF/KIORRVXld4BP6GwJG0WB6LBcJAjHUdIi0+jyFakXXS1aDv00j1TlvfHSrYtnnrHAEP+EPFThIaKk4spTi7mcMdhim3FpJhTJg6y5mrOi31NYC+Z9doUU/P0008v9BIWDLVDV1xVujzarjdjFjv0413HcfvdXJd93eXsheaDWjgjq2IulkleihYaaXVGb2yxK3cXRp2RN1veJCwjGHEZTGDN0fLiFxkLFXJVRGcm/zZK0BVXFYfHh1EvSEmYWcm/0+vkaOdRVievJt863BZuqA966iB3K0Syzp0BmVYzJoOOlhgEPcGYwI7cHXQMdHCm+0zkQcmrtFz0QOT89oXAbDbT09OjRH0RIqWkp6dnysPcSKiQi+Kq0uXRLHN1uul7oAfCAfY17sOoM47Pm+48q/2ZtTHyC2eAXifIS4mnuSe21nOlKaVc6LvAn9v/TKGtEKvJOn5Ayipo3A99LZA+uYXv1SQvL4/W1lYcjihZOooFwWw2k5c3vcQBJeiKq4aUEofHR1m2NfrgCK99u+Vter29fKj4QyQYE0ZuQMdpTTBnkHM+FXkpCdQ7BvB4A1jMU+/8hRDsztvNL87/gqMdR7mp4KbxA6y5Wum/s3HRCLrRaKSoqGihl6GYQ1TIRXHV6BsMaB7o04yfSyl5v+19ap21VGVVUWAd45roaoUh55zuzkfIT9Xi6M29se3SLSYL6+3rqXZW4/F7xt/U6SG5UAsNqRCHYp5Qgq64anQNl/xnWGMX9GA4yBstb3DScZKK9Aq2Zl5h5tRxWoubp01eCThT0pPiiDfpaemNPe69JWMLAsGJrghpjGlrtDZ4/V1zuEqF4jJK0BVXjQ63F4NOYE+MTdAHA4O8cvEVanprqMqq4oacG8Z7coSC4KiG9LVaJskcI4SgIDWBpp6BmA8OLSYLJckl1DhrCIQC42+mrdFMwrrnp4GDQqEEXXHVaO8bItNqRh/DgWj3UDcv1r5I91A3ewv3UpVVNdFgqbcegj7IKJunFUNRWiKD/hCdbl/Mrym3l+MP+anru8Ky1ZQItjztTUiFXRTzgBJ0xVUhGArT5fGRnRw9DavF08Kv636NRHL36rspSZ6kGKfrLJgSIKVwbhc7hkJ7IkJAfXd/zK/JTswm1ZzK2Z6zE29mlMFAt/J2UcwLStAVV4VOj49QWJJtm9rLpGOgg9cbXsditHDvmntJT0iPPDDog+46SF+vHTjOE/EmPdk2M43dsR2MghaqWW9fT9dgF06vc/zNzHIt26VtCqsAhWKGKEFXXBXa+7SDxZwpduiDgUFeb3idRGMiHy758MQORGPproVwEDLnL9wyQlFaEp1uL/2+2PtKrk5ejUBwwXlFmzJDHGSUa58ufLHv+hWKWFCCrrgqXOobIjnBSIIpcumDlJI3Wt7AH/Jze9Htl/PMJ6PzHJhtWn73PFM03IqusXsg5tckGhPJScrhQt+FiQeqBddCOAzNB+ZymQqFEnTF/BMKS1qdQxSkTi7SdX11NLubuS7nOlLNqVNP6OvXCnQy1sfcWm42pCWZsJgN1E9D0EHrP+ryuXAMXVGJmZAK2Zvg0jHNUEyhmCOUoCvmnQ63F38wzCp7ZEEPhoP8uf3PpMWnsSFtQwwTntb6dM6REVc0hBAUpSXS3DNAIBTBfGsSim1aE+kJYReA4hs1q99zvwave87WqljZKEFXzDtNPQMIoZXSR+Jczzk8fg/X5VyHTkT5kZQS2k9Ccj4k2udhtZFZk2EhEJI09cS+SzcbzORb8ql31U8MuxjNsOFezazr1HMQmNp7XaGIBSXoinmnsXuQbJsZs3FiNkpYhjnpOEl2Yjb5lvzok/U1aaX+2ZvnYaWTk5cST7xJz4XO6R1kFtmK8Pg99Hh7Jt60ZGmiPuSEMy9qhVIKxSxQgq6YV1yDATrdXkrSI2esXOy7iMfvYVPGptgmvHRU292mz32p/1TodIKS9CTquwcITiPsUmgtRCBocDVEHpBaBOvu0FwY6/44R6tVrFSUoCvmlQtdmknVmgxLxPtnus9gi7NRZI3B9W+gGxy1c+p7Ph1KM5PwB8M0xWjWBZpXelZiFvV99ZMPyiyD/O3QdnxRNsFQLB2UoCvmDSkl1R0eMq1mbAkTBdjlc9E+0M661HUTy/oj0XxAE/LciP1x5528lATMRj0XOj3RB4+hyFZEj7cHl881xaBdEGeBhneULYBixihBV8wbzb2DODw+KvJsEe/XOmsRCEpTYvAHH3Jquec5m7Vy/wVArxOszUriQmc/3kAo5teN9D2dNOwC2htVwXWaHXBf82yXqlihxCToQojbhBA1Qog6IcTXI9y/SwhxSghxQghxRAixI9I8ipXF4UYnSXEG1mVNDLdIKSKMon4AACAASURBVKl11pJrycViihyOGUfzQS3nPP+aeVhp7GzItREMS862xZ5qaIuzYTfbpxZ00HLTjWYt9KJQzICogi6E0AP/AtwOlAEfF0JcWW/9J2CTlHIz8GngJ3O9UMXSot01REvvIJWrUjDoJ/6YOX1OXD4XJbZJjLfG4vNAxykt7zwuBvGfRzIsZnKSzZxo6SMUHh8a6en38fyRFv7z/UbOXSH4xcnFdAx0MBiYIv6uN0DmBs3WwB97nF6hGCGWHfp2oE5KWS+l9APPAneNHSCl7JeXE20TARUEXOEcaujFbNSzMTdyuGVkt1poK4w+WcshrZCoYGF35yNcU2THPRTgePNl462zbS6eOdRM74Afg17w+7Md43LWi2xFSCQN7ii79KwKCIegu2a+lq9YxsQi6LlAy5ivW4evjUMIcbcQohp4FW2XPgEhxMPDIZkjqjHt8qW730e9Y4DN+cmYDJF/xBpdjWQkZJBoTJx6soAX2k9A+jqIT5mH1U6fVfYESjKSeK+uh0MNvbx6qp19ZzvJtJr55LWruG9bPikJRt6udYwWFNnNdmxxtqmzXQCSMrTvsztCdalCEYVYBD1S+sGEHbiU8iUp5TrgI8D/iTSRlPLfpZTbpJTb0tMnsUVVLHmONPZiMujYnB+5afNgYJDOwU4KrYXRJ2s/AUG/Zmi1SBBCsKcsk5xkM+/VdVPv6Of6Ejv3VuaRFGfAqNdxbYmdnn4/jT2Do68pshXR2t+KNzhFVagQWhNpZ6OqHlVMm1gEvRUYW8KXB7RNNlhK+Q5QIoRIm+XaFEsQ12CAmo5+NuTaiDdF9ilv7W8FGN/sORLhMLQe0RpYWLLmeKWzw2zU89GteXx6RxEP7y7mmmI7ujGdmNZkWDAb9ZxvvxxLL7YVI6WkyR0l1zytVAu79NRNPU6huIJYBP0wsEYIUSSEMAH3A6+MHSCEWC2GE4mFEJWACYhQ66xY7hxrcSIEVBZE3p0DtPe3Y9KbSIuP8p7f16gdiOZc3TL/WBFCYIs3EmeY+Mal1wlKM5Ood/SPGnplJmSSaEyk3hUl7GLN1VIze6OMUyiuIKqgSymDwBeA3wPngeellGeFEI8IIR4ZHnYvcEYIcQItI+ZjMtauuoplgzcQ4lybm9JMCxbz5JWcbQNtZCZkRjfi6jitpfHZ18zxSq8OxelJBEKSS06tuYcQgmJbMc3u5okNpMcihPapxNmoiowU0yKmPHQp5WtSylIpZYmU8tvD156SUj41/P//T0pZLqXcLKW8Tkq5fz4XrVicnG1z4w+G2TLF7nwwMIjT6yQnKWfqyUIBLX0vo0xL51uC5KXEY9CJcVYBxcnFhGQoetglpQj8A9DfNc+rVCwnVKWoYk6QUnLmkoucZDOZ1snbzHUMdACQkxhF0J1Nmvtg2tLcnQMY9TpykuPHpS9mJ2YTb4jnouvi1C9OHfa2cUZJc1QoxqAEXTEnOPp99A74WZ9tnXJc20AbeqGfvPnzCD11Wjm8LcrB6SInPzWBnn4/Q37NKkAndJQkl9DkbsIX8k3+wjgLJKZBrxJ0RewoQVfMCTUdHnRCTOqqOEJbfxuZiZkYdFOEUaTUBD21aMmGW0bItmmfVtpdQ6PXSlNKCYaDXOyLYZfualU+6YqYUYKumDVSSmo7+1llT5g0VRHAH/LTM9QTPdwy5NSyW1JisNRd5GRazeiEoMN1Oac8MyGT5LhkanqjVIMmr4JwENyq76giNpSgK2ZN74Af91Bg0iYWI3QMdCCRZCdlTz2ha7gwOXlph1sATAYdaRYTbWMEXQjButR1tA+0T22pa8vTMl5cLZOPUSjGoARdMWtGqiFXpU1ta9s20IYQgqyEKEVCfS1gjIeEq9czdD7JscXT6fYSHmPmVZpSikBwrufc5C80xkNiuvb3oVDEgBJ0xaxp6hnAnmTCOkXuOWgFRRnxGRijdRtytVzenS4Dsmxm/MEwPQP+0WtJpiSKbEWc6zlHIDxFTnpyAbhbtcpRhSIKStAVsyIQCnPJOcQq+9QmW4FwgM7BzujhFp8HhvqWRbhlhAxLHABdnvHeLBXpFfhCPmp7ayd/sS1fOxT1tM/nEhXLBCXoillxyTlEMCxZlTp1uKVrsIuwDEc/EHUPC5d1gqHnkiUlwYRRL3B4xqcpZidmkxafxknHSSYtrE4etlFSYRdFDChBV8yKNtcQQkBOcvzU4/rbEIjoO/T+Di3UkpQxh6tcWHQ6Qboljq4rBF0Iwab0TfT5+ib3STclavno6mBUEQNK0BWzosPlJS0pblLf8xHaB9qxx9uJ08dNPaGnUzsMjRZnX2JkWMw4PL4JO/E1KWuwxdk40nFk8l26LV8T9HD4KqxUsZRRgq6YMVJK2l3e0eKZyQiGg3QMdJCdGGV3DtoOfZFZ5c4F6ZY4/MEwfYPjD0B1QkdlRiXdQ900uhsjvzg5X/OEH1C+LoqpUYKumDE9A378wTBZUQTdMeggGA6SmxQlLu7zgK8fkpafoI8cjDr6J5b7l6aUYjVZOdI5yS7dOnzu4J60DYFCAShBV8yC9j4tayPHNnX8vLW/FYGI7rA44ixoyZyL5S0q7Elx6HWCLvdEQdfr9GzN3Ipj0BHZhdGcrPmjK0FXREEJumLGtLuGiDfpSU6YOt7d1t+GPd6O2TD1Th6P5sRI0vITdL1OYE8yTUhdHGHKXboQWtaPEnRFFJSgK2ZMh1uLn4spCoBG4udRd+egxc8TUsEQ5eB0iZJhMdMV4WAUtF36lowtdA12jbboG4c1FwZ7IDA08Z5CMYwSdMWM8AZC9PT7yZrC+xw0/5aQDMUm6J7OZbk7HyHdEseQP4THF9k9cW3qWhKNiRztPDrxpoqjK2JACbpiRoy4B2ZHiZ83uZvQCR15SXlTTxgYAq9rWQv66MGoJ7IPukFnYHPGZtr622jvv6Iy1JKthV6U86JiCpSgK2bESEFRpm3q8Eiju5HcpFxMetPUE/Z3an8uo4KiK0lLikMIIh6MjlCWWka8IZ4jnUfG3zCYNKMutUNXTIESdMWM6HB5sSfFRex4P4LT68Tlc1FoLYw+Yb9D+3MZ79BNBh2piZMfjAIY9UYq0ito8bTQM9Qz/ubIwahqHK2YBCXoimkjpaTD7SUnSv55XV8dAkGhrTD6pP2dWmpe3NSe6kudDEvcpCGXEcrsZeiFnjPdZ8bfsOZA0KcdjioUEVCCrpg2PQN+fIGpC4qklNT01pCTlIPFNHVbOkCrglzGu/MR0i1xeLxBBv2Tt5WLN8SzJmUNNc6a8X1HRw5GlfOiYhKUoCumTSwHoq39rbj9btalros+YTgEA93LOn4+QoZFexOMtkuvSK8gGA5yvuf85YsJdi2WruLoiklQgq6YNm19Q5iNelKmKCg60XWCBEMCJckl0Scc7NFEPXH5C3r6qDf61IKeFp9GdmI2Z3vOXs5bF0LLdlGCrpgEJeiKaROtoKitv40WTwsV6RUYdIboE46U/K+AkIvZqMcWb5wy02WE8rRyXD4XbQNjBNySDQMOremFQnEFStAV02KkoGgyh8VQOMTbrW9jMVnYmL4xtkn7O0FnWDY9RKORYY2bMtNlhCJbESa9iere6ssXrTnap5mRNE+FYgxK0BXTIlr8/ITjBE6vk115uzDqYvQ07++CRDvoVsaPY4bFTN9gAG9g6j6hRp2R1cmrqe+rJxAatt21DFsQq4NRRQRWxm+QYs5od3knLShy+Vwc6ThCcXIxq6yrYptQyhWT4TJCepSK0bGsS11HIBzgouuidsFs1VI7VRxdEQEl6Ipp0e4amrSg6L1L76ETOnbk7Ih9Qn8/+AdXlKBP5Y1+JZkJmdjibNQ6xzSStmSrHboiIkrQFTETDmsdinKTJ8bPm9xNNLob2Za1jSTTNIqDlrFl7mQkxhlIijPEdDAqhGB18mou9V9iMDCoXbTmwGCvcl5UTEAJuiJmugd8+IPhCfFzKSXvt71PclwyFWkV05vU3QZCt6IEHbSDUUcMB6MAJcklSCkvt6hTcXTFJChBV8TMZB2K6l31OL1OqrKq0Osm93aJiKdd62pviGLetcxIt8SNtvCLht1sxxZn42LfcBx9RNDdStAV41GCroiZdtcQSXEGrPGXc8ullBztPIotzhZbEdFYpNQE3RJD8+hlRpbVjJTElL4ohKDYVkxrfyveoBeMZi3FU+3QFVegBF0RM5f6vGQnjy8o6hzspHuom80Zm9GJaf44DTkh4AXryhP0kbDVSBpoNCaEXazZynlRMQEl6IqY6PcFcQ8FJsTPq3urMegMrEleM/1JR1LvLDF0M1pmjPRibY9R0NPj00k0Jl5uIm3JAf8A+NzzuErFUkMJuiImLjm1jIrc5MuCHggHqOuro8RWEr2BRSRcLZcbN6xAsm1m2l1DEXuMXokQggJLAS2eFkLh0OVPNSqOrhhDTIIuhLhNCFEjhKgTQnw9wv1PCCFODf/3vhBi09wvVbGQNPcOEmfUjeZQA9T31eMP+Vlnj8FRMRJ9zWArWDEVoleSZYtnwDd5j9ErWWVdhT/kp2OwQzMy0+nBowqMFJeJ+pskhNAD/wLcDpQBHxdClF0xrAHYLaWsAP4P8O9zvVDFwiGlpKlngPyUBHS6y/HzGmcNVpOVnMQZhEx8Hi2XOrlgDle6tBjxw4k1jp5nyUMndFrYRW/Q7IbVDl0xhli2RtuBOillvZTSDzwL3DV2gJTyfSmlc/jLPwNROgIrlhKuoQAeb5CC1ITRa26/m1ZPK2tT107qujglfc3anytY0NOS4jDoRMxxdJPeRE5Szvg4uqcdwtFTHxUrg1gEPRdoGfN16/C1yfgM8HqkG0KIh4UQR4QQRxwOR+yrVCwoTT1aheJYQa/prUEgYmtgEYmei1r63QorKBqLXifItJrpcMVe8VlgKcDpdeL2u7WK0VBAtaRTjBKLoEfafkU8xRFC3IQm6H8X6b6U8t+llNuklNvS01fmQdhS5KKjn+QEI8nDDS2klFT3VpNryY2tvdyVhEPQexHsq1ds/HyELJuZLrePYCi2XfZIw+0mV9OYlnQqjq7QiOW3qRXIH/N1HjDhJ0gIUQH8BLhLSqm2DMsEbyBES+8QqzOSRkMrl/ov4fF7Zr47d7Vq+edppXO40qVJTrKZYFhG7WA0gi3Ohi3ORrOnGeJTwBCn4uiKUWIR9MPAGiFEkRDCBNwPvDJ2gBCiAPgV8JdSytoIcyiWKPWOAcJSsjrjsuFWdW81Jr2JIlvRzCZ1VGsNLVJm+PplxEhef1tfbGGXkfTFS/2XCMigtktXO3TFMFEFXUoZBL4A/B44DzwvpTwrhHhECPHI8LD/D7AD/yqEOCGEODJvK1ZcVc63u7HGG8myahkZvpCPi30XKU0pjb2BxVhCQeg8C+mlK86/JRKJcQZSE01cilHQQUtfDIaDtPcP2yb0O7RYumLFE0PDR5BSvga8dsW1p8b8/2eBz87t0hQLTd+gn+beQa4vsY+GWy44LxCSIdanrp/ZpI5qCPoga5qujMuYnOR4LnR5kFLGlDGUnZSNQWegyd1EgS0PZFjLGrJP00tHsexY2SdSiik53tyHTgjKcqyj1871nCMtPo20+LTpTygltPxZM5ZKKZy7hS5xcpPj8QXCdPf7Yxpv1BnJTcqlyd2EtBWA3gjdF+Z5lYqlgBJ0RUQ83gCnL7koz7FiMWuhlc4BzYhrXeq6meWed1/QwgOrroOZvH6ZkpuixdGnG3Zx+924ggOQWgw9F5RRl0IJuiIy+y90A1BVlDp67VjXMUx608yyW0IBuPgnbXeeUT5Xy1wWWM0GLGbDqF9OLORbtMSzJk8TpK0BX7+WPaRY0ShBV0ygrquf6g4P24tSscVru/OeoR4aXA1UpFXMzIir4W0Y6oPSPSs+9/xKhBDkJsfT1hebURdo6Ysp5hSa3c2QtlZLX2w7Ns8rVSx21G+WYhzeQIg3qjtJt8RRVXh5d3686zhGnZGN6RunP2nXeWg5DLmVKnY+Cbkp8fT7griGYs9WWWVZpaUvCqEdMjtqtDdNxYpFCbpiHG/XOhjyh9lTlol+2IjL5XNxwXmBcns58Yb4KDNcQb8Dql8FWy6s/sA8rHh5MGJL3DqNsEuBtYCwDNPa3wr5Vdq5xIV9E71d+h1Q/zYcfRoO/huceh46z6mY+zIkprRFxcqgqWeAc21urilKJWM47xy03blO6NiUMU1X5IAXzv4K9CYov1uze1VEJDXRRLxJz6W+ITbk2mJ6TXZiNkadkWZ3M0X5RVByM9TugzMvQvparSNU9wUY6NbE3pYHSTbwdMC5l6HjFJR9RPPUUSwLlKArAAiHJW/VOEhJMLJ9zEFov7+f6t5qyuxlJBoTY59QSqj+rRYC2PwAxM3A82UFMRJHb+kdjDkfXa/Tk2/J19IXpUTkbtX+3hve1szPhE4T8dI9Wpw9brjaV0poOw51f4QTP4fNn1CivkxQgq4AoLrDQ++Anw9vysagvxyJO+E4gUSyOWPz9CZsO67tDld/AJLzo49XsMqeQF1XP32DAVISYzt4LrQVUu+qp3Owk6zELMjbBjlbtNZ0xsTI1bhCaOcZZhucfhHOvgQV96lPUMsAFUNXIKXkYEMPGdY4StIve7YMBgY513OO0pRSrCbrFDNcgdcN9W9qB6B52+Z+wcuUEXvipt7BmF9TZCtCJ3Rc7Lt4+aJOP2zcFeVNwV4Ca28DZyM0vTeDFSsWG0rQFVzqG6JvMEBlQcq4j/qnuk8RCoeozKic3oQX39DK0Uv3qgKiaWCLN2KNN9I8DUGP08dRYCmgrq8u5pTHcWRvgqwN0HTgctNuxZJFCbqC8+0eTAbduN25L+TjTPcZipKLSDGnxD6Zp0NLU8zbDgmp0ccrRhFCsCo1gZbeQcLh2MW5JLmEgcAAHQMdM3vw6lvBlAg1r6vuR0scJegrnEAoTG2nh9UZSZgMl38cznSfwR/yszVj6/QmrH9bO2DLv2aOV7oyKLAn4A+G6XDH1pYOtLCLXui50DdDPxejGUpugv4u6DwzszkUiwIl6CucescA/mCYsuzLMfJAKMBJx0kKrAWkJ0yjs5SzCXrroeB6lTUxQ/JTEhACGnsGYn7NiDf9BecFguHgzB6cUQbWbC1DJjTDORQLjhL0Fc65dhcWs4G8lMsFQ+d6z+ENeqe3Ox9Jl4uzaBkUihkRb9KTY4vnoiN2QQdYb1+PL+SjwdUwswcLAUW7NU+YztMzm0Ox4ChBX8H0+4I09QxSlm0dPQwNhoOc6DpBTlIO2UnZsU/WWw+uS7Dqes3OVTFjVmcm0e3x4RyIzU4XIC8pD4vJwvne8zN/cEohWLKg+aCKpS9RlKCvYGo63EgJ68aEW2qcNQwEBma2O49P1rImFLNipN1fnaM/5tcIIVifup5WTysun2tmDxYCCq7TKkx7lL/6UkQJ+gpFSsm5dg/ZNjOpw0UsYRnmRNcJMhIyyLPkxT6ZowY8nVC4QxWnzAFWs5Esm5nqDs+0UhHXpa5DIDjbc3bmD08rBbNVKwxTLDmUoK9QHP0+uj0+1o/Zndf31ePyuajMqIy9gUUoCPVvQWKa8jmfQzbk2Oj2+Gh3xZ7tkmRKoii5iPM95wnMtMeoTqc5NzobtZ26YkmhBH2Fcr7dg14nWJuleaxIKTnhOIEtzkahrTD2iVoOar/4q29RPudzyNosCyaDjhMt07PD3ZS2CV/IR42zZuYPzx7u99p+auZzKBYE9Ru4AgmFJefb3RSnJ2I2aiGS9oF2uga72Jy+GZ2I8cdiyAnN72vOfqnF87jilYfJoGNzfjI1HR46prFLz0rMIiMhg5OOkzOrHAXN4yW1WHNjVIejSwol6CuQekc/Q/4QG3Iu27SecpzCbDBTmloa2yThMJx7BYRe250r5pxthSkkmPTsO9eBLxiK6TVCCDalb8Llc9Hobpz5w7M2Dre1a575HIqrjhL0FcipVi33fMQMajAwSIO7gXUp6zDqYkw5bNqveX+U3qbt6BRzTpxBz+0bsnEOBHj2UAs1HR6Coeg75mJbMUnGJI53HZ/5Lt2+Wks/7Tw3s9crFgQl6MsQKeWkv8gdLi/NvYNsyk9GN9yRqMZZg5SS9fb1sT2grxma3td2cZllc7VsRQQK7AncvSWXsJS8drqdf3unnt+d6aDLM3kYRq/TU5lZScdAB62eGTaO1hu1UJqjWlWOLiGUH/oywh8M8+4FB+fb3VpecraF60vSRuPkUkr213VjNuqpyLONXjvXc47sxOzYTLgCXjj/GzAnw5pb5/PbUQxTYE/gf11XSItzkJoOD3WOfuq6PNy5KZcCe0LE16xLXcexzmMc7jxMniUv9qylsWSUQccZrWgsPcZQnGJBUTv0ZUIoLHnlZBunL7kozbRQkp7E6VY3T7/fyJlLLkJhyYH6Hlp6B7lhtZ04gybynYOduHyu2HbnUkLt61pstexOrdO84qqg0wlW2RPZU57FQ9cXYos38urpdjzeyOmJBp2BrZlb6RjooMXTMrOHphSBKQG6ZpHXrriqKEFfJhxu7KWld5A9ZVnsKc/itg1ZPHBNASkJRv5wrpMn/3SBg/W9rM+2sHFMz8oGVwNCCAqthdEf0nEauqqhaCdYc+bvm1FMSYLJwB0VOQRDYfZf6J503NrUtVhNVt5re49QOLZD1XHodJC+HrrrIOibxYoVVwsl6MsA54CfQw29rM2yUJZzuVAo3RLHfdvy+fCmbLYXpXJHRTZ7y7NGP35LKal31ZOXlIfZEMUdsb8LLvwekgsg/9r5/HYUMZCSaGJrYQrVHR66JrHaNegM3JB7A06vk9PdMzTcyiyDcBC6a2exWsXVQgn6MuBQYy8C2F060epWCMHqDAs3rE5jTaZlXCy119uLy+eiyFY09QMCQ3DmV2AwQ9ldqoBokVBZkILJoONI0+QVnYXWQgqsBRzpPDIzjxdrrpbFpLJdlgTqN3OJ4xoKUN3uYUOejcS46Z1xj1itTinoQR+cel5rOlz2kcud4xULzsjhdm2nB9dg5Fi6EIJdebsAeKP5DcJymoVCQmi7dGcj+Kdn6au4+ihBX+IcbepFCNi2ahpt4oapd9WTlZhFojEx8oDBXjj+31pbubKPQHL+LFermGs25ycjEBxvmXyXbjVZ2ZW3i/aBdg51HJr+QzLKtB6xjupZrFRxNVCCvoTp9wU5e8lNWbYVi3l6HuRuv5vuoe7Jd+c9F+HYf4LPAxs/qtLWFikWs5HSzCTOtrmnrCYtTSllfep6jnUe44Jzmta4SRma+ZoKuyx6lKAvYY42OQlLqCqcfjPmkXBLsS2CB8ulY3D6BYizQuX/AnvJbJeqmEe2FKTgD4Y52+aectyuvF1kJ2bzZsubOAYd03tIZjm4WsE7Q691xVVBCfoSZcgf4nRrH2uzkrAlTL9DUIOrgVRzKra4K8r2ey7ChX2QWgKVD0LC9N8sFFeXLJuZ3OR4TjT3EQ5PXuqv1+nZW7gXs97Maw2vMRgYjP0hGcN1Cl2z6IikmHeUoC9RjjU7CYbljHbng4FB2vvbJ+7Ogz6oeR0S7FD+EdVKbgmxpSAZ11CA+u6pDy4TjAncXnQ73qCXPzX/KXavl/gUrfagUxUZLWZiEnQhxG1CiBohRJ0Q4usR7q8TQhwQQviEEP977pepGIs3EOJESx9rMizYk6ZfrdnobkQiJ8bPLx3VYuZrP6jEfIlRkp6ENd7I8eboTSnSE9LZkbeDFk8LRzuPxv6QzHKtHmFg8mImxcISVdCFEHrgX4DbgTLg40KIKx2ZeoEvAU/M+QoVEzjW5MQfDHNN8czCIQ2uBiwmC2nxaZcvhgLQckiLl9ty52iliquFTifYnG+j1Tk0aaHRWMpSyyhNKeVwx2E6Bzpje0j6Oi2NUe3SFy2x7NC3A3VSynoppR94Frhr7AApZZeU8jAww75XilgZ8oc43tLHmswk0mawO/eH/LR6Wim2FY83bOqu1QqI8rfP4WoVV5PyHBsmg45jzdG7HAkh2Jm3kwRjAm+2vBmbNUBcEiSv0uLoM7XlVcwrsQh6LjDW3ad1+Nq0EUI8LIQ4IoQ44nBM85RdAcCB+m6CIcl1xfYZvb7R3UhIhibGzzvOaM2Bk1fNwSoVC4HZqKc8x0p1h5t211DU8XH6OHbn7abX28uxrmOxPSSzTOtU5Wmf5WoV80Esgh7Jd3NGb89Syn+XUm6TUm5LT59Ypq6Ymi6Pl1OtLirybTOKnYPWCDrRmEhWYtbli75+cDZA5gbtI7ViyXJtsZ2kOAOvne7APYkT41gKbYWsSVnDsc5jsVkDpK0FnV7lpC9SYhH0VmBsiWAe0DY/y1FMRiAU5vdnOog36me8Ow+EAjR7mimyFY0Pt/TUaR+hM2JscKFYtJiNeu6oyMEbCPHMwWaONjmjCvt12dchhOBA24HoDzCatX6jjvOq3+giJBZBPwysEUIUCSFMwP3AK/O7LMWVvF3joLvfz97yrNGGFdOlyd1EMBykJPmKQqHeixBngUT1qWk5kGUz87GqfFISTLxT6+Cn7zbw339u4lRr5Dz1JFMSlRmV1LvqudR/KfoDMstVv9FFSlRBl1IGgS8AvwfOA89LKc8KIR4RQjwCIITIEkK0An8L/L0QolUIYZ18VsV0qO30cPqSi6rCVArTJvFdiYEaZw2JxkSyE7MvXwyHNOOl1GIVbllGpCXFcV9VPg9dX8iu0jQMOsGfznfxm1NtEfuSbsrYhMVkYX/r/ugGXqrf6KIlpjx0KeVrUspSKWWJlPLbw9eeklI+Nfz/HVLKPCmlVUqZPPz/U9chK2LCNRjgD+c6yUk2c13JzEItAAOBAZrdzZSmlKITY/7Z3Zcg6Ffl/cuUlEQTW1elcn9VPjeuTafeMcBbNRMTEow6I9dlX0ePt4fzvVGqQfVGSCtV/UYXIapSdBETCktePd2OEHDbhmz0upnv9+S1sAAAF2RJREFUoGt6a5BI1qdeESfvuQhCp7JbljlCCLYUpLC9KJXTl1zUO/onjClJLiE7MZtD7YfwhaJ0KMos1yqLnQ3ztGLFTFCCvoh5r66bTreXPWWZ2OJnXrkZCoc43X2anKQcks3J42/21muFRMYoHYsUy4Jri+2kJpp4u9YxIfQihOCG3BvwBr3RK0hTCsEYr4qMFhlK0BcpLb2DHGt2UpFnY3WGZVZz1TprGQgMsCVjy/gbPo9Wyp2qwi0rBb1OsLs0nb7BQER3xoyEDEpTSznlODV1GqNOr1WO9lzQQnaKRYES9EWILxji92c7SEkwsStCW7npEAwHOdp5lLT4NAosBeNv9tZrf6ZGsNBVLFtW2RPItpk50uQkFCHr5drsa9ELPQfao6QxZpZpMfSeafqrK+YNJeiLkEMNvXi8QfaUZ2LUz+6f6HjXcdx+N9fnXD8+9xw0QY9L0hoYKFYMQgiqilJxDwWo6fBMuJ9oTGRLxhbq+6KkMdrytXRXle2yaFCCvshwDvg53txHeY6VbFv8rOZqdjdzpPMIq5NXk2fJG38zHIbeBpWuuEIpTkskLcnE0WZnRAvdTRmbSDIm8d6l9yZPYxRCK0brrQf/NLzVFfOGEvRFxoH6HvQ6wQ2r06IPnoLuoW72Ne3DbrZzY/6NEwe4W7UsBRU/X5EIIahclUK3x0dz70QxNuqMXJ9zPd1D3RzvOj75RJnlWr/R7pp5XK0iVpSgLyKcA35qOz1U5NlIjDPMeJ5+fz+v1r+KUWfkg0UfxKQ3TRw0kq6YUjjzBSuWNOuyrCTFGTjaFNlDvSS5hNXJqznUcYiOgY7IkyRlag1RVNhlUaAEfRFxuLEXvRBUFqTMeA5/yM9rDa8RCAf4UPGHSDIlRR7YUwfJ+SpdcQWj1wk2FyTT1DOIwzMx71wIwe783ViMFl5reA2nN4Lwj4RdXC3gVbWEC40S9EXCgC9IdYeHDbkz352HZZh9Tfvo8fawZ9We8Q0sxuJ1aV1nVLhlxbMxV/NQn2yXHqeP48MlH0YgeKnuJdr6I/jyZZZr5m6q3+iCowR9kXCu3U0oLNmUnxx9cASklLzb+i7N7mZ25e2iwFow+eDu4TQz++oZPUuxfDAb9ZTlWKnp8OCZxJXRFmfj7tV3E6eP4+W6l3n/0vsEwmPGJqSCJQs6z6jGFwuMEvRFgJSS060u8lLiSU2MEO+OgZOOk5ztOcuWjC2U28unHtx1DpLSIfH/b+9OY+O6rgOO/8+bfcgZ7psoaqckM7I2y7IdJbEdxY4j70mQOkWcNkURBEiDBGhapAkK9GOBAm0DtEBgpIuDpA0ML4ntOlbsxHLiRdZmy7asjRLFRSQlrjMkZ593+uFStiRzGYqkRqTuTxgMxXkzPI/kHN53l3OvvDaMtXhsbapAUQ53Tr6QqDxYzpfXfpmWqhbe6XuHJ44/cenWdQ0bzSK1kUn62q2rwib0a0D7QIJYMsvGpVfWOm8dauWN7jdYXb6aWxtunfrg5DDEzkLt5dvCWtersrCP5toI754dJp2bfCu6gCfA7U238+CaB8m7eZ5pfYbjg+OzW2o/AR4vdE8xI8aadzahXwPePRsj7PewpnaSAcwpdMQ7eLnjZRpKGti5bOfHFw9d7kLtjZr1VxCptVhtXV5OOutOWA7gco2ljXxl3VdYUrqE33f8nmODx8zgem2LufrLTVPYy5o3NqEX2Wg6R1vfGC1LojOuptgeb+fFMy9SGaxk16pdeJ1pBlNd17SgKleafk/LGtdQFqKxPMTbHRNvgnG5oDfIrpW7aIw0sqdzj5nW2LAZ8lnTl24VhU3oRfZBdxxXlRsbywp+jqpyoPcAL5x+gfJAOfetuo+Ap4A9RvuOmoJcjTfNImJrsdq6vIJ4MsuxCcoBTMTreLl7+d1E/BF2n9lNuqQKog3Quc9uT1ckNqEXkary/tkYTZVhysOFDYZm8hl2n9nNvt59NFc083Dzw4R94emf6Oah7Q9mMNTObrEmsLqmhLpokDdO9U+4q9FEgt4gdy2/i0QuwWvdr8Oy28w4TZ+dwlgMNqEXUddQklgyy4bGwnbri6VjPHXyKdribexo3MHOZTvxOQXWSW9/w7zRVt1pa7dYExIRPt1czUgqx+Gu4YKfVxuuZUvtFo4PHqfdHzArRzvetFMYi8Am9CJ6/2yMoM/DmprpB0O7R7t58sSTJHNJ7l91P5tqNk0/AHrBcIdJ6HWfsFvNWVNqqgyzojrMvrYhUtnJZ7xcblvdNiqDlezpepVM080w2mc3vygCm9CLJJnJc/L8KOsbIninKZF7dvQsz516jpA3xJeav/TxyolTGT0P7z0JoXJovnuWUVvXg0+tqSGdy7P/zGDBz/E6Xu5ouoNENsHefBwiddD2qhkkta4am9CL5MLK0E8smbq7pT/Zz2/afkM0EOXh5ocpCxQ+eMpwB7z9c7Op78Y/sXVbrILURALc0BDlnY5hYonCE3J9ST0bazby/sARuhs2mNounW/NY6TW5WxCLwLXVd7pHKaxIkRtZPIkm86nebHtRXyOj/tX3U/IO4P66OePwuFfgr8UtjxqWuiWVaBPrq5CBF5r7Z/R87bXbyfqj/JK/CTZ6rWmq2/0/DxFaV3OJvQiONU3SjyZnbKqoqrySucrjGRH+PyKz09eNXEinfvhg1+b+hpbvmaTuTVjkaCPm5ZXcuLcCGeHkwU/z+fxcXvT7cTSMQ6UVYE3CEefnXzfUVWz0UrbH+DMaxDvmaMzuD7ZhH6VqSoH2ocoC/lYVV0y6XFHBo5wevg0tzbcSn1JfaEvDq2/g9aXoboZNn0V/AVMabSsCdy0vIJI0Murx/sm3NVoMk2RJm6ovIF3Bo/SuWybqez5wa8u7U9XNTsdHXrcXEm2v2ES+sH/hqPPm71KrRm78l0UrCtyqm+U3liKu1rqcCZZGdqX6OO1s6+xLLqMzTWbC3vhXBqOPmcqKTbeBGs+B479e21dOb/XYceaal58v5ejPSO0TDPec7EdjTvoTfTyUuwYX175KaJtr5nkvXS7OaD3PTPGEyyD9btM2QA3D517of1NSMfhxq+Y+jBWwew7/irKu8obpwaoLPHT0jDxmyOTz/Db9t8S8oYKq80CMHIODv3M7ELUfDc032WTuTUn1tdHqC8L8nprP5lc4as//R4/X1jxBVx1eS7RTmL9vaaFfuz/zC05aH5Pt38TGjaZgXtfEFbdAevvhaF201Vj57LPiP3zdxUdbB9iYDTDA5uXTNg6V1X2dO4hnonz0OqHph8EzWXgzB+g66B5M2x6BCqWz1P01vVIRPjM2hqe2N/JwfYhbltdeMnl8mA59668l+dOP8ezw0d4YMujhHMp82C4avIFbg0bIZuEU7+H03tg9Z2zP5HrhG3GXSU9sSR7Tw+wti7C6kkWEh04d4DW4VZuqb+FhtKGqV+w7wTse8wMgDZsMi0dm8ytedBYHmJdfYSD7YOTboIxmYbSBnat3EU8E+epU88w6HGgpHr61cpN22HJFujYCz2HZxH99cUm9KtgaCzD84d7KA14+ez62o89rqocOneI/b37WVe5ji21WyZ/sVTMLBR6/ynTKt/6KKy7B3wzmNJoWTO0Y001qvD6DKcxAiyNLOWhNQ+Rd/M8ffJpOuOd0z9JxHTJVK6EE7tNF4w1LZvQ59mZ/jGeONCJq8oDm5cQ8nsueTyWjrH7zG729uyluaKZO5beMXG/eT5nBov2PQZDbbD6s3DTN6BsBqtGLesKlYV8bF1ewdGeEXpjqRk/vzZcyxfXfpFSXynPn36eg+cOTj9zxvFAy0MQqoAjT0Oi8JWr1yuZyXSkubRt2zY9cOBAUb72fFFVMnmXsXSenliSYz0jnB4YJhxMcVtzKUG/SyafIZlLkswlGUgN0J/oxxGH7Q3b2Vyz+ePJ3HWh/7jpS0wOm+mIaz5n55ZbV106l+fxN84Q8nt55OYmfNOUrJhINp9lT9ceTg6dZEV0BTuX75y+9HNyCA4+bq5Ct379ur8aFZGDqrptwsdsQr8yqkrHYIK2/jHOx9OMpnMks3nimRjxXA+j+T7yMkRlRKmPBi8ZBBURwt4wUX+UpZGltFS1UOIrufjFITEA/Seg513zC11SbRJ55coinK1lGW39Y/zq7bO0LIlyd0vdhFeTfSNp3u4YonMoSTqXJxL0sbq6hE1N5ZQEvGYP3f73eL37daL+KDuX7Zx+rcVwJxz+X7NYbsOXwD/5Go7Fzib0OaSqnOob5Y8n+xlOZPF5hIoSIea2MZBpJ60xfB6HmpIymiubqAnXUBYoI+KPEPAE8Hv8+B3/BC3xvJl2OHgKhs6Y1jhAWaOZu1u91k5FtK4Jb54aYO/pAVqWRLlzXS1+r4Oq0jmY5FDHEG39Y/i9DiuqSgj7PQyOZegcSuDzOHxydRWbm8oREXpGe3ip/SXGsmNsqtnEzfU34/NMUQ76/DGz1sIbMDNf6jZcl6WgbUKfI+dHUrx6vI+uoSTVpX5uaPQyoq2cjJ0g5+aoL6lndflqVpatJOovcBFGPgdd+6DrAGTGwOuH8uWmJV7VDMHCF3NY1tWgqrx5eoC3Tg8S9HmojQQYTmaJJ7OE/R42N5WzqamcoO+j8aKhsQx7TpznTH+CpRUhPr+hnmjQRyaf4c3uNzkycISgN8jG6o00VzRPXoRu9LyZxz7Sa7odl2wxZaEDkat09sVnE/oEsnmXobEMmbxL0OchEvQS8HomPDaWyLLvzCBHuk398rVL8qQ8p2iPn0FEWFuxlo01G6kOVRcegKrpUmn9nZm5UrUalmw1idyZOA7LupZ0Dyd5tytGLJkh7PeyqqaEdXWTl4NWVY50x3n1RB8i8Nn1tayvNw2Wc2Pn2Ne7j84RMwOm1FdKebAcv+PHM/5+kPF/oHjiPVTGe6lNxKn1hnHKlpqNz6tWm0HUCy13VXBz5goYTOt+gbfqbUIfl3dNd8l7XTG6hpK4l517ScBDRdhPRdhPJOglnXM5F09xdjhJThPUVMbxh3sZSvcT8ATYUL2BG6s2EFaF1LBZDCGO+YW6+JfqciO9ZtHEULvpG2++CypWzP83wLKuAbFElheP9NA9nGJ9fYQ719d+2JqPpWO0xdroS/YRT8fJulnyapLxhVylKFk3SyqXgkyCQCrGsmyGZTllqbeUEl/YDJyqmveke1FdGI/PvDdLa6G07qP7BTTQOuuELiL3AD8GPMBPVfUfL3tcxh/fBSSAP1fVQ1O95pUm9AvxFrxbDxBLZjlyNsb73THG0nmiIR/r6iLURQP4vQ6prMtwMsX5kVH6xkY4PzZMLB0nRwK/P0UgMEYkmCGgaaodP+v8lbR4SvClYmYq1URF/D0+iC6BaKO59wbMsf0nTFEibxBWfhoatti+ceu647rK/jOD7D09SEnAw03LK1hXHyHsL2zxuqoylh2jZ6yHjngHHSMdJBMDkI4RcZVabwnl3jDRQDmRYCVhX5iQ4yOYTSPJQRg7D+nRj14wGDWJvaTmo0Q/VaPs4wFBNmFqKoljbo5n/GPPRf+f/dXBrBK6iHiAE8BdQBewH/iqqn5w0TG7gO9gEvotwI9V9ZapXvdKE3rr+RFeeK+XgNch6PNQEvBSeuEW9FIa8BDyeYilE/TE45zuH6QrFiOnKaoi0FDhEA25ZNwM6Xz6w1vOzX10eZbP4GYS+HIpom6eynyeBldo9JRQ5QmaH0qwDEKVEK78qEXuC5vnJwdNKzx+1mzFpRfVwAiUmpWdS7fbDSes6965eIpXjp2nZ3xueyToJejzkHeVTM4l67rk8yZHhQPm/R0N+igL+YiGxu+DPkoDHgbS/XSPdnMucY6+RB8j2ZGPzXUXEUKeED6PD8llcLJjOOlRJD2KpEdwMmMI4CD4vH78wUr84Sp8gSh+bwi/ePCriy+XIZDLEMil8WdTBHNpfDp5Q1NVyaPkxUM+FMW7ZCv+Zbde0fdsqoReyJ/D7UCrqp4ef7FfAg8CH1x0zIPAz9R89/aKSLmINKjqnBc3Hhs8AIP/w5i6xF0l6+bJ5l1yrouLi2K+cRd+kH6PEPZ7qAp48Q0L6ZiXUcdLEC/ljge/OATwEBAPASCgQtTxE3V8hJwypLTKdIuU1Jj6E+Fqk7ynqgJX3mSSNpjW++g5cx+Imj8AC7wPz7LmSl00yCPbl3EunqJjMMHAaJp0zsXncfB5HLwewec4uKokMnlG0zm6YylOnBu9pMvUESHkd3AkikgZjqwliEtGE6TzY2Q1RdZNkdU0WTdFXnOAD6UM1Qg6/g9vFm9uFF9uFCc1jAz1QP4kec1c8vVcx4srfvKOD/fDmx/HCSAK4CKqKC6qeRQX1EXUxeum2JLx8cgVJvSpFJLQG4GL1+p2YVrh0x3TCFyS0EXkm8A3AZYtWzbTWAGoKytnXf0yHBEcHDNMImawJO9CLg95V4j4AlQEQpT5g4Q8AcKOn6DHjyOeC8EActHHmEsif8S0ooPlJonPtnynx2dXc1rWNOqiQeqihV+xuq4yksoRT2WJjc+wSWTyuKq4alrErgKEgRlMVpiIKqJZNJckJ0oGyJEn52bIaZasmybnZshqhpxrul8/bLKJ4ODBEQ+OOOYeD0vr5ycnFJKtJmpOXt5PU8gxqOpjwGNgulwK+NofU9+whfqGKWqdWJa16DmOUBb2URb20VTsYK4hhYzGdcEl37OlQPcVHGNZlmXNo0IS+n6gWURWiogfeAR49rJjngW+LsatQGw++s8ty7KsyU3b5aKqORH5K2A3Ztrif6rqERH51vjjPwFewMxwacVMW/zG/IVsWZZlTaSgET9VfQGTtC/+3E8u+liBb89taJZlWdZM2BUtlmVZi4RN6JZlWYuETeiWZVmLhE3olmVZi0TRqi2KSB8wm51fq4GZ71i7cNjzW9js+S1s1/L5LVfVmokeKFpCny0ROTBZgZrFwJ7fwmbPb2FbqOdnu1wsy7IWCZvQLcuyFomFnNAfK3YA88ye38Jmz29hW5Dnt2D70C3LsqxLLeQWumVZlnURm9Aty7IWiUWR0EXk+yKiIjLLrUmuLSLyTyJyTETeFZFnRKS82DHNBRG5R0SOi0iriPyg2PHMJRFpEpFXROSoiBwRke8WO6b5ICIeEXlbRJ4vdixzbXwLzSfH33tHReS2YsdUqAWf0EWkCbOBdUexY5kHLwEbVHUjZqPuvytyPLM2vun4vwNfAFqAr4pIS3GjmlM54K9V9QbgVuDbi+z8LvgucLTYQcyTHwMvqup6YBML6DwXfEIH/gX4WybY8m6hU9Xfqmpu/L97MTtBLXQfbjquqhngwqbji4Kq9qjqofGPRzDJoLG4Uc0tEVkK3Av8tNixzDURiQKfAf4DQFUzqjpc3KgKt6ATuog8AJxV1cPFjuUq+AvgN8UOYg5MtqH4oiMiK4AtwFvFjWTO/SumEeUWO5B5sAroA/5rvEvppyJSUuygCjXLLe3nn4i8DNRP8NCPgB8Cd1/diObWVOenqr8eP+ZHmEv5X1zN2OZJQRuKL3QiUgo8BXxPVePFjmeuiMh9wHlVPSgidxQ7nnngBbYC31HVt0Tkx8APgL8vbliFueYTuqp+bqLPi8iNwErgsIiA6Y44JCLbVbX3KoY4K5Od3wUi8mfAfcBOXRyLBhb9huIi4sMk81+o6tPFjmeO7QAeEJFdQBCIisjPVfVrRY5rrnQBXap64arqSUxCXxAWzcIiETkDbFPVa7VC2oyJyD3APwO3q2pfseOZCyLixQzw7gTOYjYh/1NVPVLUwOaImNbF48Cgqn6v2PHMp/EW+vdV9b5ixzKXROSPwF+q6nER+QegRFX/pshhFeSab6Ff5/4NCAAvjV+F7FXVbxU3pNmZbNPxIoc1l3YAjwLvicg745/74fi+vNbC8B3gFyLiB06zgDa9XzQtdMuyrOvdgp7lYlmWZX3EJnTLsqxFwiZ0y7KsRcImdMuyrEXCJnTLsqxFwiZ0y7KsRcImdMuyrEXi/wFd75PAdZobCQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"f = plt.figure()\n",
"plt.plot(y_d, g_xn, alpha=0.5, label='Never takers')\n",
"plt.plot(y_d, g_xa, alpha=0.5, label='Always takers')\n",
"plt.plot(y_d, g_x, alpha=0.5, label='Compliers')\n",
"plt.legend()\n",
"f.savefig(\"X.pdf\", bbox_inches='tight')"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"α = x = df['α']\n",
"α_n = n['α']\n",
"α_a = a['α']\n",
"\n",
"distributions = [α_a, α_n, α]\n",
"α_min = min(df['α'])\n",
"α_max = max(df['α'])\n",
"\n",
"α_distributions = []\n",
"\n",
"for i in distributions:\n",
" α_d = np.linspace(α_min, α_max, 2000)\n",
" kde = KernelDensity(bandwidth=0.1, kernel='gaussian')\n",
" kde.fit(i[:, None])\n",
" \n",
" # score_samples returns the log of the probability density\n",
" α = np.exp(kde.score_samples(α_d[:, None]))\n",
" \n",
" α_distributions.append(α)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"g_αa = α_distributions[0]\n",
"g_αn = α_distributions[1]\n",
"g_α = (1 / φ_c) * α_distributions[2] - (φ_n/ φ_c) * g_αn - (φ_a / φ_c) * g_αa"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"ename": "NameError",
"evalue": "name 'g_αc1' is not defined",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m<ipython-input-20-716b964b4151>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m 2\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0my_d\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mg_αn\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0malpha\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0.5\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlabel\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'Never takers'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 3\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0my_d\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mg_αa\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0malpha\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0.5\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlabel\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'Always takers'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 4\u001b[1;33m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0my_d\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mg_αc1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0malpha\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0.5\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlabel\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'Compliers'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 5\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlegend\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 6\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msavefig\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"alpha.pdf\"\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbbox_inches\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'tight'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mNameError\u001b[0m: name 'g_αc1' is not defined"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3Dk6Xng9+/buRsdkHOcnHZnw2wguUwiKS55FFfR4krWHU9SsdY+yqeyZR99Z7tcdeUqq5xOVcczRfHk4lmyKXopihS1FPOJS3J3ObNxZifsYDAzyBmd0fn1Hy8wg5kBMI3u368bDTyfWlQPOry/txfAgxdveB6ltUYIIUTjc9S7A0IIIawhAV0IIfYICehCCLFHSEAXQog9QgK6EELsEa56Xbi9vV0PDw/X6/JCCNGQXn311UWtdcdmj9UtoA8PD3Pu3Ll6XV4IIRqSUurmVo/JlIsQQuwREtCFEGKPkIAuhBB7hAR0IYTYIySgCyHEHiEBXQgh9ggJ6EIIsUdIQBdiGzOxVa7OJagozXQ2CeOvwPxlkDTVogbqdrBIiN3u2kKSv31zGq3hzHAL7z286eG8zaWW4PX/G/Kr5vP+M3D4I/Z0VIg1MkIXYhOlkubH7yzQ1uThaHeI125Gia3my3xxEd7+a1AKzvyuCeaT52Dlhq19FkICuhCbGF9OE03nefJAG08dbkejuTgdL+/Fk2chtQhH/xGEuuDAB8EXhrF/sLfTYt8rK6ArpZ5WSl1RSo0qpT63zfMeU0oVlVK/bl0Xhai9y7MJvG4HI+1NhH1u+lsCXJ1P3P+FxTxMvAKtB6D9kLnP6YKBJyE+DfEZezsu9rX7BnSllBP4PPAx4ATwrFLqxBbP+2PgO1Z3UohaKpY01xaSHOoI4nKaH5GR9gBLyRzxzH2mXWbPQy4Ng0/eeX/XSRPYZ96wqddClDdCfxwY1VqPaa1zwFeAZzZ53h8AXwPmLeyfEDU3G8+QK5Q40NF0676hNvPv8aX01i8slWDi5xDugebBOx9z+6DjOMxfNKN4IWxQTkDvAyY2fD65dt8tSqk+4FeAL2zXkFLqM0qpc0qpcwsLCzvtqxA1MbGcRinobwmYBU6taWvy0OR1MrmyTUBffAdWV8z0ilL3Pt59Cgo5WLxqX+fFvlZOQN/kO5O7N9X+G+BfaK2L2zWktf6i1vqM1vpMR8cOtoAJUUMTy2k6gh58kz+DF/83eOVPUYkZuiN+ZmKZzV+kNUy8DP4WaD+y+XOah8zi6NyFrS++cAXe/CuYfLX6NyL2nXIC+iQwsOHzfmD6ruecAb6ilLoB/Drw75RSv2xJD4WoQCyd59vnZ/jR5Xky+W3HGXfIF0vMxDIcdUzBjZ+YxU1dgrf+igF/jmg6TzpX2OSCE2bBc+AxcGzxY6WUmUtfHjOHju62cgPe/jrEJ+Hqd2F2m8AvxCbKCehngcNKqRGllAf4FPDNjU/QWo9orYe11sPA88B/rrX+G8t7K0QZ8sUSX3ttktH5JG9Nxnjh/EzZJz1nohlKxQIjsVfMlsOTvwqnPwWlAkPLPwGtNx+l33wJPAHofnD7C3Q9YEbz8xfvvL9UgqvfA18zvOuzZh7++j+YKR8hynTfgK61LgCfxexeuQR8VWv9tlLqOaXUc3Z3UIidemPCHAL65Yf7+OCxDm4upXlnbpMR8SYmVtK0ZcaJOFZh+L1mtB1ohZEPEEmP05EZYyZ6V0BPzJpRd/9j4HRvf4GmNhOsZ8/fef/cebN3/cAHwOWFoacgE4el0bLftxBl7UPXWr+gtT6itT6otf6f1u77gtb6nkVQrfWntdbPW91RIcqhteatyRiDrQEGWgOc6o3QHvLy8thSWaP0ieU0R4pXcQVaoPXg7Qf6HsUZ7uFk+ixzyyt3vujGT8Dlgd5Hyutk1wOQnDcfYHa9XH8Rwr3QcdTc13oAvMF7A78Q25CTomJPmVxZJb6a52RfGACHQ/HoYAvLqRyTK6vbvjZbKLK0EqNHz5u57o1z4Q4HHP0Yza48vvEXKZXWfjksjppdK0PvMVsTy9F5HBxOmHlrrdPnIJswo/P13TEOh3ne8pjZGSNEGSSgiz3lnbkEHpeDgx3BW/cd6Qri9zh5czK67WsnV1aJpG8Q8bk236kS6sY5+Bit8YtEr/4UohNw+VvQ1A59Z8rvpCdggvX0azBx1ozw2w9Dy9Cdz2s7ZObQV66X37bY1ySgiz1Da82NpTQDrQHcztvf2i6ngxM9Ya7Np0hlN9mhsmZiOU175gbBSCuEujd9TvjER1jxD5K//D14/S/A4YJTv2ZOge7EwV8AbwhGv29ujzx973MiA2Y+XebRRZkkfa7YM1bSeeKreR4bbrnnsQf6Irx6c4WLM3EeG27d9PWTSwkeZg5nx/s2PxgERJp8jPd+nLBvnq6hJjOSd/t33llPE5z5PUjMQKjHzMHfzeE0c+lLo2ZnzBZ9EmKdjNDFnnFjKQXAUGvTPY+1NHnob/FzfjK26eJoOlcguzROs89hgugWlFL0NPsZ073Qc7qyYL7O5THTLJsF83WtB0xumKRk1BD3JwFd7BmTK6s0B9xEAptvHTzVFyG2mmdi+d7F0bGFFJHMNJGAF5oHNnn1bV1hH0vJ3I4OLFWsZdjcyjy6KIMEdLEnaK2Zia7S27z1iPlwZxCf28mF6dg9j11bSNJVmqWpvd/MW2+jJ2J2s8zFt0gDYCVf2Cy6SnEMUQYJ6GJPMEfyi/RGtg7oLqeD4z0hRueTdxzfzxaKTC5E6XdFUesj4m10hU1An90qr4vVWkbMjpri1gu6QoAEdLFHTMfMNEpv8/Z7wU/1RSiWNBemblcfujSTILA6TUeT596tg5vwuZ20BT3M1mKEDmbapVQw+WKE2IYEdLEnTEczeN0OWpu2WWAE2oNeDnQ0cfbGMqlsgUKxxOvjKww7lwj6fRDuL+t6PRE/U9HV2weM7NQ8aHa8yDy6uA8J6GJPmImt0hvxo8rY2vfewx0US5pvvTXN9y7OEU3neSCUMHvPy9xPPtQWIJsvMVOLUbrLA+E+WJaALrYnAV00vEy+yFIyt+2C6EatTR4+erKbhUSWy7MJHh8K065XTC6VMg22BnAoxY3FVKXd3pmWYbN1MVej64mGJAeLRMObjpr58/XdJ+U42h1iqC1AvlgilJ2DqSJEyptuATOP3tvsY3Q+ybsPtpX1l0FVWkfg+o/Nbpeuk/ZeSzQsGaGLhjcTy+BQiu4dBHQwQTnkc0N8rV5LuG/7F9zleE+Y5VSO6Vrsdgl2m+Rfsn1RbEMCumh409FVOsPeO/K37EhsEvzNJl3tDhzpCuF1Ozh3Y7my6+6Ew2FK2C1fN2kAhNiEBHTR0IolzVw8s6PpljtoDfGpHY/OATwuB0+MtDK2kOLidPz+L6hW64hJs5uuwS8Q0ZAkoIuGtpjMki9qerY5ULStTMzU96wgoAM8NNDCQGuA716c5cfvLJAvlirrRzkkDYC4DwnooqGt1/fsuc+Boi0lZs1tuKeilzsdik+e7uVUr8nm+P+8Mk48k6+sL/fjbzEfsn1RbEECumhoM9FVQj4XIW+FG7aSs6Ac0NRZcR88LgcfPtHFrz7SRypX4BuvT1Gwa6TeMgyxcVNUWoi7SEAXDW06lqE74qt822BizhRu3mmBik0MtTXxsVM9LCZzvD6xfXWkijUPmpJ0yTl72hcNTQK6aFjJbIH4ar7y+XOtzQg9uHl1okqMtDcx1Bbg9fEVinakBVhP7Rsdt75t0fAkoIuGNVtmQq4tZROmeMQW5eYq9chgC6lskdH5pKXtAqZcXaBVEnWJTUlAFw1rOprB6VB0BLfPX76l9WmLYJd1ncKkBQh6Xbwzl7C03VsiA2aELvPo4i4S0EXDmomt0hX24qr0QFFi1tTptDigOxyKQ51BbiymyBVsCLrNg1DIQmrB+rZFQ5OALhpSvlhiLp4tOyHXppJz4G/dvqZnhQ52BCmUNFPRe8vdVU3m0cUWJKCLhjQdXaVY0gy0BCpvJDkHIWtH5+t6mn24HIrx5bT1jfsipjRdfMr6tkVDk4AuGtLkyioOpSofoReykIlDU4e1HVvjdjroafbbE9DBpPpNzNjTtmhYEtBFQ5pYTtMd8eJxVfgtnFo0t4F26zp1l8HWAIuJ7B31Sy0T6oXVqORHF3eQgC4aTiZfZC6epb+a6Zb0WkBvsi+g97eYvx6mozak111PVRCXUbq4TQK6aDjjy2lKWjPc3lR5I6lFcLjA12xdx+7SGfLidChmYjYsjAa7zQ6dxLT1bYuGJQFdNJxr80n8Hic94QoPFIEJ6IFWk2fcJi6ng86Q91YCMWsb95i/LmSELjaQgC4aSqFY4vpSipH2JhyOKsq+pRdtnW5Z1x3xMR/P2JMGINxnRuhS8EKskYAuGsq1hRTZfIlj3aHKG7F5h8tGPRE/+aJmMZm1vvFQN+QzsLpifduiIUlAFw1Da82bE1HCfjeDrVUsiNZgh8u69Tzt03YcMAr1mtv1nO5i35OALhrG6HySqegqZ4ZaKk+XCzXZ4bIu5HUR9LqYtWMevandLOzKwqhYU30SaCFsUCiWOHtjhanoKl6XA5/byZXZON0RH6f6ItU1XoMdLuuUUvQ0++xZGHU4IdgpI3Rxi4zQxa6jtebvzs/w8tgSuUKJ5VSOq/MJBtua+OTpXpzVLIZCTXa4bNQT8RFbzZPK2nHAqMcEdMm8KJARutiFLs7EGVtI8YGjHTw82GL9BdKLEOm3vt0tdK8V4JiJZTjUGbS28XAPTL0Kq8s1mUISu5uM0MWuUixpXrq2RE/Ex0MDNkyJ1HCHy7qutQNGtsyjh9ZOjEpeF0GZAV0p9bRS6opSalQp9blNHn9GKfWWUuoNpdQ5pdRT1ndV7AfXFpIkMgUeH2mtbuFzKzXc4bLO5XTQEfIybceJUX8rON1ywEgAZQR0pZQT+DzwMeAE8KxS6sRdT/sBcFpr/RDwu8CXrO6o2B/WtyUOt1VxrH87NdzhspFtB4wcDrMfXUbogvJG6I8Do1rrMa11DvgK8MzGJ2itk1rfOq7WBMjRNbFji8kskyurnO6PVHcKdDs13OGyUe/aAaMluw4YJeehVLS+bdFQygnofcDGirSTa/fdQSn1K0qpy8DfYUbp91BKfWZtSubcwoKUzxJ3Oj8Zw+lQnOytclvidmq8w2Vdd2TtgJEt8+i9UCrcnk4S+1Y539WbDZXuGYFrrb+utT4G/DLwrzdrSGv9Ra31Ga31mY6O2i1Kid0vXyxxaTbO4c4gfo/TvgulF2u6ILou7HPR5HUya8c8eqjb3Mq0y75XTkCfBAY2fN4PbHk0TWv9Y+CgUkr2UImyvTOXIJsv8UC/jaPzWztcav+tqZSiJ+K354CRvwVcXgnooqyAfhY4rJQaUUp5gE8B39z4BKXUIbW2JUEp9QjgAZas7qzYm7TWvDYepS3ooa+aos/3sz4lUYcROpgDRtF03voKRkqtHTCSgL7f3Tega60LwGeB7wCXgK9qrd9WSj2nlHpu7Wm/BlxQSr2B2RHzmxsWSYXY1jtzSRYTWR4btmmr4rrU2rpNoM2+a2yjp3m9gpFN0y6pRSjacBpVNIyyTopqrV8AXrjrvi9s+PcfA39sbdfEXqO15vpiiomVVSJ+NwMtZufHj67M0xn2crSripS45UgvgtNlpijqoDvsw+NyML6c5lCnxe813Gt2uaTmzb/FviRH/0XNnL2xwk9HF3EoRWnDH3BNXicfP9Vj31bFdaklMzq386+AbTgdiv4WPzeX0tY3vnFhVAL6viUBXdTEYjLLz64tcrQ7xEdPdpPMFBhfTlPUmqNdIXt3tqxLLUDLkP3X2cZQWxNjCyli6TyRgNu6hr1h8AQk8+I+JwFd1MQrY8u4nQ4+eLQTp0MRCbh5IGDjjpa75TOQTdT0yP9m1gtz3FhKcTpg4eGm9YXRuORG388kOZewXSZf5NpCkpO94dqMxDeTru8Ol3UtATctATej80nrG48MmIXRbML6tkVDkIAubHdlNkGxpDnRG65fJ5Jz5jZY34CulOJIV4iJlbT12xdbR8ztyg1r2xUNQwK6sN3ofJK2oIeOoLd+nUjMgdtv5prr7HBXCK3h6pzFo/Rgl5lHX75ubbuiYUhAF7bKF0tMR1cZamuyd4/5/SRnzU6QevZhTXvQQ1vQw6WZuLUNKwUtw7ByXSoY7VMS0IWtpqOrFEr61mJgXZSKZm452FW/PmyglOJUX4SZWIa5uMWpADqOQS5tgrrYdySgC1vdXErjdCh7j/TfT2rBBPX1vdq7wImeMB6XgzcmotY23HrQ5HWZe9vadkVDkIAubDUdXaU7Yk5I1s36Vr5dMkIH8LmdnOgJc2U2QSydt65hpwu6TsLCZZOITOwrEtCFbYolzUIiS3fYV9+OxCbB01S3I/9bOTPcggJeGrM4j93AE6A1XP8Ha9sVu54EdGGbxWSWQknfKu5QN7EJaB7YFQuiG4V8bk4PNHN5Ns5CwsJKRv5mGHwSZi/A+MuyQLqPSEAXtllf8Ouq5wg9EzNTD5GB+z+3Dh4bbsXndvLDy3NYmqB0+CnoOArXfgQ//1O48VOZgtkHJKAL28zGMvg9TsK+OmaYWB4zt831zeGyFb/HyXsPtzMdzXBhysKA63DCyV8xH94wXP8xvPKnMHfRumuIXUdyuQjbzMUzdId99d1/vjhqpiDqUKWoXCd6wlycjvPi6AKDrQHrknYpBZ3HzEd6Ga68AJf+FrwhMwUl9hwZoQtb5AolllK5+k635FfNMfi2w7tu/nwjpRS/eMJsqfy78zMUijbMeQda4dSvgy8M7/y9zKvvURLQhS3mExm0pr4LorPnoVSA7gfq14cyRQJufvFEF3PxDN++MEupZEPBL7cPDn7IHLKau2B9+6LuJKALW9xeEK1R/pZS0ZyQXF9YzK/C+EvQPAih3bP/fDuHOkO8/2gHo/NJnn9tkpVUzvqLtB8200+TZ2//vxJ7hsyhC1vMxrKE/W4Cnhp8i829De98BwpZMz/cMmQKPRSycOhD9l/fQo8MtuB3O/nh5Xm+/NINDnUGOd3fTH+L35q1CKWg/wxc+XuIT0Gkv/o2xa4hAV3YYnZtQdR20Qmz0Bfug/YjJkgtj4HLD6d+bVcd9y/X8Z4wg60BXh+Pcn4qxtW5JB0hLx850WXNmkTnCbj6fZi/JAF9j5GALiyXzhWIr+Z5aMDmikRaw9XvmlH5g78JLo+916uhJq+Lpw6388SBVq7MJnh5bImvnp3gkw/1MtTWVF3jLi+0HTAB/eCHwCEzr3uFfCWF5ebi5tRjZ8jmEfrKdUjOw/B791Qw38jtdHCqL8JvPzFES5OHb701Y83ceucJyKXMKVqxZ0hAF5abjWVQqgYnRKdfNwUdOk/Ye51dwO9x8sxDvTiU4jtvz1Z/qrT1gDl8tHzNmg6KXUECurDcXDxDW5PH3gyLhRwsjUHHcZNhcB8I+dy870g7M7EM71Rb7cjlNfPn6ydpxZ4gAV1YSmvNbDxDp92j85XrZo95xxF7r7PLHO8O0xHy8tK1RWtG6ckFyfGyh0hAF5aKrxZYzRXpsftA0cIVc1AmMmjvdXYZh0NxZriFlXSeG0vp6hprPWhuZZS+Z0hAF5aaXTtQZOuWRa3Nkf7WA/tyh8bhzhBBr4s3Jlaqa6ipHbxBiN60pmOi7vbfT4Ow1UxsFZdD0Ra08YRoetns0GjeX6PzdU6H4mRfmJtLaZLZQuUNKWXSCkcn5NToHiEBXVhqLp6hK+zD6bAxGdb6iHKXpsSthaNdIbSGd+YS1TXUPADZBKxWOdoXu4IEdGGZYkkzH8/SZff8eXTcTBXsspJytdQW9NIZ9nJlttqAvvZLUfaj7wkS0IVl1kvO2bogqrUJPpHdV1Ku1g53hpiNZaqbdgm0mb380XHrOibqRgK6sMxsrAYl51ZXIJvct/PnG420mxQANxZTlTeycR5dNDwJ6MIys/EMAbtLziVmzG24z75rNIj2oIeQz8X1agI6mF+OmRisRq3pmKgbCejCMtPRVXqaLUrzupXEDDhcu7qkXK0opRhpb2J8OV1dlaP1AtqxSWs6JupGArqwRDJbIJrO09fst/dC8RlTsMLhtPc6DWKorYlcoXRr/39FmjpMKgBZGG14EtCFJaZWVgHob7ExoJdKkJyFUI9912gwpvAFTCyvVt6Iw2HyusgIveFJQBeWmIqm8bgcdNh6oGgJioWGLFphF5/bSUfIy8RKlWkAIgOm1miuynZEXUlAF5aYWlmlt9mHw84DRYlpcxvqte8aDWigJcBsLEO+qnn0tcpFMkpvaBLQRdVWc0UWkzn6mgP2XigxawpZBFrtvU6D6W/xUyxpZqJVzKOHesxic0z2ozeysgK6UupppdQVpdSoUupzmzz+20qpt9Y+fqaUOm19V8VuNRWtwfw5QHzaBJ59fqDobn0tfhxKMVnNtIvTBeEeGaE3uPsGdKWUE/g88DHgBPCsUuruEjHXgfdrrR8E/jXwRas7KnavyZU0bqey90BRsQCpBZk/34TX5aQzbMU8ej8k5kzxENGQyhmhPw6Maq3HtNY54CvAMxufoLX+mdZ6PbvPy4CUEt9Hbi6l6Wvx25uQKzUPpaLMn2/BzKNnyRWq3I+uSxCfsq5joqbKCeh9wMYNqpNr923l94Bvb/aAUuozSqlzSqlzCwsL5fdS7FrxTJ7lVI7B1ior0d/P+glRGaFvqr/FT0lrZmJVbF+M9JvpLNmP3rDKCeibDbs2TZ6slPogJqD/i80e11p/UWt9Rmt9pqOjo/xeil1rfK1qzlBbDRZE3X7wRey9ToPqbV6fR68ioLu8EOyUefQGVk5AnwQGNnzeD0zf/SSl1IPAl4BntNZL1nRP7HY3l9KEfC7amjz2XigxIwui2/C4HHSFvdUtjIKZdolPmekt0XDKCehngcNKqRGllAf4FPDNjU9QSg0Cfw38jtb6Heu7KXajUkkzvpxmsDVgb/6WYh5SSzLdch8DrRbNoxcL5i8i0XDuG9C11gXgs8B3gEvAV7XWbyulnlNKPbf2tP8BaAP+nVLqDaXUOdt6LHaNuUSGTL7IUJvN8+fJObNYJ0f+t7U+jz4drXIeHWQevUGVledUa/0C8MJd931hw79/H/h9a7smdrubS2mUgsFWu+fP58ytjNC31RO5PY8+3F7hL1lv0Bzcknn0hiQnRUXFbiym6Az58HtsznyYmAFPE3hD9l6nwXlcDrojVsyj95sRuhSObjgS0EVFUtkCs/EMBzpsnm4BWRDdgf6WAHPxLNlCFYuakQHIZ0yyLtFQJKCLiowtpNAaDnYE7b1QIWeyLMp0S1kGWgJmP3o1eV3W59GlzmjDkYAuKjK2mCTsd9MetHm7YnLO/OkvAb0sPc0+nI4q96P7W8zH0lXrOiZqQgK62LFcocT4UpqDHU32bleE29vnJKCXxe100B32VZfXRSnoOAIrN83Ui2gYEtDFjo0vpyiUtP3TLWAOuXhDsiC6A4NtAebiGVLZQuWNtB8xW0WXRq3rmLCdBHSxY5dnEwQ8Tvvrh4IJ6JHtUgeJux3oaEJruL6YqryRcB/4wjB73rqOCdtJQBc7ki0Uub6Q4kh3yN7qRACZGGTit6vSi7J0BL2EfC7GqgnoSkHvw7Byw5zSFQ1BArrYkdH5JIWS5lh3DaZAYmtpXMMyQt8JpRQHOpoYX0pRqKYsXfeDporRzZ9a1zlhKwnoYkcuzSSI+N1021nMYl18ylTSCXbaf609ZqQ9SL6omahmt4s3CAOPwdzbsHTNus4J20hAF2VbTuWYWE5zqi9i/+4WMPugw33gsPkk6h400OLH63ZwZTZRXUOD74ZgB1z8G1iUbYy7nQR0Uba3JqM4HYqTvWH7L5ZNQnIeWobtv9Ye5HI6ONIZ4tpCsrrsiy4PPPCfgL8Vzj8Pl1+QrYy7mAR0UZZcocTFmTiHO4M0ecvK6VadlevmtvWA/dfao471hMgVSowtJqtryBeGh38HBp+E2bfg3L+XtAC7lAR0UZZ35hJk8yUe6K9RxaDl6+AJQLCrNtfbg/qa/YR8Li7NxKtvzOmCgx80gb1UgAt/bfKmi11FArq4L601r09E6Qh5a7P3vFSE5WtmdC4JuSqmlOJ4T5ibS2nimbw1jUb64NgnTH6dKSl7sNtIQBf3NbmyymIiy0MDzbVZDF25YeZpO47bf6097lSf+YvqwmTMukbbDppfthOvSKm6XUYCurivNyai+D1OjtZi7zmYbXIuL7SO1OZ6e1jE72akvYnzUzGKJQvzm/efgVwaFq5Y16aomgR0sa14Js+1hSSneiO4nTX4dsnEYOEydJ2S7YoWebC/mXSuyOh8lYujG7UeAF/E/PIVu4YEdLGtS9NxtKZ2i6Hjr5jbwSdqc719YLgtQMTv5s3JqHWN3srIeB0KWevaFVWRgC62pLXm8myCvhY/Eb/b/gtmkzDzphmd+2r0C2QfUErxYH+EqZVVFpMWBt/2I2sL2GPWtSmqIgFdbGkunmU5leNETw0OEoFZZNNFs99ZWOpEbxinQ3HeysXRcD+4fRLQdxEJ6GJLV+cTOB2KQ501yHueS8P069B5wlSdF5YKeFwc6QpycSZe3cnRjRwOaB40hTDEriABXWzp+mKKvmY/PncNFicnz0IxD0Pvtv9a+9QD/c3kCqXq87ts1DxsFrJXV6xrU1RMArrYVCydZymZY6Sjyf6L5TPmkErHUWhqt/96+1RvxEd7yMtbU1G0tmgLY/OguZWC0ruCBHSxqfX8HwfaaxDQZ96AQk5G5zZTSvFgX4T5eJbZuEUJtpraTYqG6IQ17YmqSEAXm7q+mKK1yUNzwGPvhUpFmDwHLUNSCLoGjvWE8LgcvGXV4qhSJsVxfNqa9kRVJKCLexRLmunoKoNtAfsvtnAZsgnof9z+awm8LifHukO8M5sgW7Do2H6ox+R2yVdRTENYQgK6uMdcPEO+qBlosTkRl9Zmq2KgzeQHETVxojdMoaS5OmfRydFwr7mVUXrdSUAX95hcK1vWa3dmxdgEJOZMmTPJqlgz3WEfLQE3F61IqwsmoCslAX0XkIAu7jEVTdMW9BDw2F7ugCkAABeSSURBVFzIYvIsuP3mZKioGaUUJ3rNydFY2oK0ui6v+SsrMVN9W6IqEtDFHUolzXQ0Q7/d0y2rUVOjsvchcNYgrYC4w7GeEEph4Si9zxT1tmo7pKiIBHRxh/lEllyhRF+zzQuiU+cABb2P2Hsdsamwz01/S4BLM3Fr9qSHe815AjlgVFcS0MUdpqJpgOpH6PnVrUuUZZNrx/yPmXqVoi5O9ISJreaZilqwO+XWwuhU9W2JitWg2q9oJJMrq7QE3JUXgs6l4dLfmoRNThf0PATD7zVJnNbd+AmUSuZ+UTeHOoP86IqDSzMJ+luq/Iss0G6mzuIz0P2ANR0UOyYjdHFLqaSZiq5W/sNdKsGFr5lj4MPvMYm2pl6Fs392u7LNwhUzOu97VJJw1ZnH5eBgR5B35hLki1Um7HI4zH50GaHXlYzQxS2LqSzZfIm+SqdbZl6H2CQc/yXoXtu50vcoXP47UyXe32wWQ8O9cOD91nVcVOxET5hLM3HGFlLVlxgM95hTv8WC+etM1JyM0MUtU2v7zysK6KUijL9sqsJ3nbx9f6gbHv00HPowBDth+Ck4/azsbNkl+lv8hHwuLlmx2yXcZ74PknPVtyUqIr9GxS1T0VXCfjdhXwXBduEKZOJw+KP3HhJyOM3hoYHHrOmosIzDoTjWHebVmyuksoXK107ATLmAOWAU6bOmg2JHZIQuAFNubmpllb5KT4fOvQ3ekBzhb0DHe0KU1soNVsUXNt8DCTkxWi9lBXSl1NNKqStKqVGl1Oc2efyYUuolpVRWKfVH1ndT2G0lnSedK1a2XTGXNrtaOo/LEf4G1Bb00h3xWTTt0iMpAOrovgFdKeUEPg98DDgBPKuUOnHX05aB/wL4Xy3voaiJW/PnlYzQl6+BLpmALhrS8Z4wC4ks84kq86SH+8zCdy5lTcfEjpQzQn8cGNVaj2mtc8BXgGc2PkFrPa+1PgtYkBhC1MNUNE2T10lzoIL58+XrpsjB+hyqaDhHu0I4HYpLM1VOu9yaR5e8LvVQTkDvAzaWI5lcu2/HlFKfUUqdU0qdW1hYqKQJYZPJlVX6mgOonU6ZaA0r16FlWKZbGpjf42S4vYkrs3FKpSpSAYR6zPeBzKPXRTkBfbOf0oq+4lrrL2qtz2itz3R0dFTShLBBLJ0nkSlUtl0xOWfm0FsPWN8xUVMnesKkssVb5Qcr4vKYsnQyj14X5QT0SWBgw+f9gHy19pCby2a+c7C1ghOiKzfMbcuIdR0SdXGgvYmQz8UbE1WWp1svSVeq8vSp2LFyAvpZ4LBSakQp5QE+BXzT3m6JWrq5lCbsd9NSyfx5bNIc4fcGre+YqCmHQ/FgfzMTy2mWktnKG2oehEJWDhjVwX0Duta6AHwW+A5wCfiq1vptpdRzSqnnAJRS3UqpSeC/BP47pdSkUkrS6DWAYkkzvpxmqLXC+fPYpBmRiT3hVF8Yl0Px5mS08kZahs3tynVL+iTKV9axMK31C8ALd933hQ3/nsVMxYgGMxvPkCuUGKqkIHR62aTJjciXfq8IeFwc6Q5xaSbBuw+243M7d96Ip8mkeVi+DkPvtr6TYktyUnSfu7mUQikYqGT+PD5pbiWg7ykPDzaTK5R4c6KKUXrriMm8WMhZ1zFxXxLQ97mxhRS9EX9lI7HYlMlzHmizvmOibjpDPkbam3h9IkquUOHCZusBk6hreczazoltSUDfx2LpPAuJLAc7K1zQjE9BuF/2n+9Bj4+0sporcn6qwh0vkUEz9TJ/0dqOiW1JQN/HRhfMqcBDHRUE9PwqpBYlq94e1dvsZ6A1wGs3VyhUUvzC4TCpIJaume8VURMS0Pex0fkkHSEvkYq2K65VppEdLnvWEyOtJLMF3qh0Lr3nNJQKMPWatR0TW5KAvk8lswVmYhkOVTPdohy3iwOLPWegNcBIexM/v7HMaq648waCnSad8uRZGaXXiAT0ferKbBytTVKmisSnIdghlYf2uKcOt5MrlHjl+lJlDYy8zxwyuvKCnBytAalYtE9dnEnQE/HR0uTZ+YtLJTNCl+rue1570MsDfRHemIhyrDtMd8S3swZC3XDwgzD6A3jty2sFUJSZilEOs+W19YAsrFtEAvo+NJ/IsJjI8gvHOitrIL0IxbzMn+8T7znUzthCiu9dnOXZxwdxOXf4h/3A42bHy82fwY2fmvscTnPSWJfMX3rHP2mmaERVJKDvQ5dmEjgdiiMVT7esL4jK/Pl+4HM7+dDxTr7xxjQ/urLAh4937jxNRNdJ86HXErUqBcUCLF65PXo/+atSwrBKMoe+z+SLJS5OxznYEcTvqeAwEZj5c7cf/C3Wdk7sWgc6gjw+0sqFqRhvTlaRjVGp29MrTpcJ8md+F/ytcOFrpti4qJgE9H3mymyCTL7I6YFI5Y3Epsx0i8x77ivvOtDGgY4m/uOVea4vWlhizhuEh34Lgl3w9tdh5k3r2t5nJKDvI1prXp+I0h7yVlY7FMz2s/SSHCjahxwOxcdO9dAR8vLC+Rnm41XWH93I7YfTz5pMjZdfgKvfM4VTxI5IQN9Hri2kWExkeWSweedzoOvWK9HI/Pm+5HE5eOahPrwuB994Y5p4xsIywi4PPPAb0H8Gpl6Fl/4tvPpluPJtmH4D8hb+AtmjJKA3sInlNH/z+hR/8fJNvvv27LYV20slzUtjS7QE3BzvriJVfWzSbDeTgtD7VtDr4pmH+sgVS3zjjWmyhQoOHW3F4YTDH4Ezvwd9j4LDBQuXTVB/5QsmJa/YkgT0BnV5Ns7XXptkMZkl6HVxdT7JX748zg8uzZHJ3/sD9ur4CouJLO851I7DUcXcd3Tc7C12eavovWh0HSEvn3iwh6Vklp+NVnjoaDvBDjj0IXj4t+E9fwiP/hOz9fHC81KvdBsS0BvQcirH9y/O0dvs5x+/a5hffriP33tqhEeGWjg/FePLP7vB+cnYrertV+cS/Gx0icNdwcqP+oPJbZ2YMSXGxL431NbE6f5m3pyMbvvXYdWUMlN8D/0WuANw6Vtmy6O4hwT0BvTi1QUcDsXHH+jB4zJfQp/byfuPdPBbjw/SEvDw/Utz/NmLY/yHl27wrbdm6Ap7+ciJrsrnzsHsPy8VJaCLW951sA2vy8lL12wYpd/N0wRHnjaL8lPn7L9eA5KDRQ1mOrrK2EKKpw63E/Te++XrDPv4jTP9jC2muDqXJFsocrI3wun+yM5P+N0tOn77uLYQmIHEo0Mt/HR0kbl4hq7wDlMD7FTbQVMNaeIV6H3ELKSKW2SE3mBevbmC3+PkdH/zls9RSnGwI8jTp7p55qE+Hh1qqT6YA0RvQqhL5s/FHU4PRPC5nbxyfbk2Fxx+ymxplP3q95CA3kBS2QJjCylO9IRvTbXUTC5tFqNa5Wi2uJPX5eTB/ghjC0liaQu3MW4l0m/OQUy9ejuVgAAkoDeUizNxSlpzqq+KU56VWr5mfnjaDtX+2mLXe7A/gkLxxmQVhaV3ou9RWF2RmqV3kYDeILTWXJiK0dfip7WSlLfVWho1R7RD3bW/ttj1Qj43h7uCvD0dq7yw9E50HDOLpFIN6Q4S0BvExPIq0XSeB+oxOi8WzEio9aDkbxFbeniwmWy+xMWZuP0XczhNibvla2akLgAJ6A3jwnQMn9vJ4Wr2kVdqadTsQe88Vvtri4bRE/HTE/HxxvgKuhZz270PA0pG6RtIQG8A6VyB0fkkx3pC1uxW2am5C2a6pXm49tcWDeXhwRZW0nlrszFuxReGjiNmt0shZ//1GoDsQ28Al2biFEu68umWfAamXzfTJqUCtAxB/+PgCdz/tdmkeV3fI+CQ3/9ie4c6g4R8Ll4bj3KgowZ/TfadgfnLMP/22oh9f5Of0F1Oa835yRh9zX7agzvc/10qwsTP4ZX/E8b+IxRzJtnR+Mtw7s8hMXf/NqZeNWXCeh+pqP9if3E6FKcHmplYTrOQyNp/wUi/ORsxeU62MCIBfde7uZRmJZ3f+VbF5TETtEd/YDIjPvppOPNPTbKjRz9tFjff/H8htc2R7ULOjOzbDkGgtZq3IfaRB/oiuJ2K18drsFiplBmlpxZh6Zr919vlJKDvcj+/sUzI5+Jod5n1PxNz8OZfmY9S0eSXfvA3Ibwh3W2o2xQTUAre+gpkE5u3NfGyKWgx+GT1b0TsGz63kxO9YS7PJkhla5BEq+ukKYd4/R/2/ShdAvoudnMpxdTKKmeGW3HeL+VtLm2y0J37c0hMw8FfgMd+H9oPbb7VMNBqAn1+Fd76q3uLB6QWTb6MzuOSu0Xs2MMDLWhtBiS2czhNOoDkPMxfsv96u5gE9F0qVyjxw8vztATcnOq9T0GK5TE4+yWYe9uMpp/4z2DwCVOEdzuhbjj1q2ba5cLzt0t+ZZOmtqPTY3JSC7FDLU0eTvSGOT8Zs7aq0VY6T5i59NHvm0HKPiUBfRfKF0t8+8IMsdU8Hzretf1WxalX4a3/z9RkfPTTcPCD4N5BxrvWA3D8l0yelrN/Bhe/YX45ZKJw8lfAW+ZUjxB3eeKAWXepSWpdhwOOftwE86vf3bdTL7JtcRcplTRX55O8cn2J5VSOXzjWyUDrFlsLSyW49kOYPGsWLU98svIsiF0nzBzkjZ9AbMpMsYy8D4Kdlb8Zse+FfW4eHmzm3I0VTvSEt/5etkqoG4bfA9dfNBsBBh6393q7kAT0XaBU0lyZS/DK2BIr6TxtQQ/PPNTHSHvT5i8oZOHiN80Jzv7HzHx5tXvEwz3w4G9U14YQd3nyQBtX55J8/9Iczz4+iM/ttPeCQ+8xc+mjPzBz632P2nu9XUYCep3FVvP8/YUZpqMZOkJeful0Dwc7gltXFkotwcWvm9sjv7jvvmFFY3E7HXz0VDfPn5vkhfMz/NLpXtzbTCFqrUlkCxSKGgUEvE68rh38ElAKjn8S9N/AO981eV5GPrD5epLWkJiF5Kw5a9HUAeH+hj5AJwG9jkbnE3z34hxaw0dPdnO8J7R1IC/mzXz5jRfN4aAHf8PMfwuxy/U1+/nQ8U6+f2mO51+d5L2H2+mN+ClqzUo6x2Iix3wiw0Iiy0IySzZ/Z7ZGn9tJV9hLT8TPoc4g7UHP9qUUnS6z/jP6A5g4a/anDzwOzUPmZye1AItXYemq2QCwkTdokn71PtyQ60eqJkl0NnHmzBl97tz+rAuYL5Z48eoCb07E6Ar7+PgD3TQH7kqJq7VZ4EnOmm/I+YtmF0r7YTjy0Yb8ZhP729W5BD+8PE86V7znMbdT0R700hHy0h704nU7KJVMHqNoOs9sPMNiMovW0B708MhQC8e6w/ffzrt0zZySTs7feb/TbUrZtR9Zq5GrTM3c2fNrOdaVyRPTcdwE+VLRDKpKeVOGMdBmPuqQfVQp9arW+symj0lAr51socjofJJXxpaJreZ5rMfFk+FFXMkZyMShkDHz48Ws+eZZ/9o4XNB2wORfaR6o75sQogrZQpHriylWUnkcCpoDHjpCXpr9bhz3Cc7pXIGrc0nemoqxmMgS9rt5YqSV4z33Cexam1F5YsYEZn8LRAa23tabXjYnpGfeND+PW/FFzGi+7xGzy6xGqg7oSqmngT8BnMCXtNb/812Pq7XHPw6kgU9rrbfNabmbArrWmkJJ43Ko7f+U20K2UGQ+nmVubRSxmi+SyZfIF0vki5pC8fa/HaU8B9Uk7wrM0JafMQ14Q+agj8tndqo4PebD5TN7a0M9UsdTiDVaa64vpnh5bJm5eIbmgJsnRto41h267y+FHSnmIb0EuZRZYHW4zc+lLpq594XLsHzdjPZ7HzIpCPxb1/q1SlUBXSnlBN4BPgJMAmeBZ7XWFzc85+PAH2AC+hPAn2itn9iu3VoH9GyhSCpbJJHJE03nia3mia7miaVzxFbz5IsmoIf9btqDXtqDHlqbPIT9bkI+Fw6lUAoyuRKJbJ6FRJb5RJb5eIalVO7WYDrkc9HkdeF1OXA7HXhUAX8xSVN2nuDqFN35CcLuEsrfCt2nbh9bFkLsiNaascUUL11bYiGRpTng5khXiMHWAG1BD363s6IB2o4k5kyKjPnL5vO2g2YbcbgX/K33P9xXge0CejlXexwY1VqPrTX2FeAZ4OKG5zwD/Adtfju8rJRqVkr1aK1nquz7PSZHLzB57lu3piM0gNZoNOa/O+8HTUlDqbRhoUWD0wE9bicH3A78bidup4N8sUhmpUQqWyCTLzADzNxqzFAbPok4Ff0+FyGfi6DXRcjrwo2CLJDV5vrFDafk3H7oPQ1dp8xeb6n+I0TFlFIc7AhyoL2J0fkkb07GOHdjhZ9fN+kGHErhcztwOR0ozI+bWnudtT96p3E7DtAeO0/zzCU8hbO3Hik63JQcbkpqPdQqNIrIwcc49tiHrewEUF5A7wMmNnw+iRmF3+85fazFw3VKqc8AnwEYHBzcaV8BcHl9uMNdt744qLUvEIAyXzhzLfOgUuZ+r8uB3+3A53bS5HXjcztQ5sUbe3jrX/mSJpUrkM4VyeQ1pbVfDi6nE7/bScjnwu9xodZfs0U7eALgDUOwC5raJYgLYTGlFIe7QhzuCrGaKzKfyLCYzJHJF8nkixRKem38ZwZ39iwbesk1f5B5/QE82WU82SU8uSiOUhZnMYdTFzEjQzPw9Pjt2dRQTkDfLALd/b+knOegtf4i8EUwUy5lXPse3QOH6B6wv/K8G2he+xBCNAa/x8lQWxNDbVscyquJ3rpduZwd9JPAxq0V/cB0Bc8RQghho3IC+lngsFJqRCnlAT4FfPOu53wT+MfKeBKI2TF/LoQQYmv3nXLRWheUUp8FvoPZtvjnWuu3lVLPrT3+BeAFzA6XUcy2xX9qX5eFEEJspqw9NVrrFzBBe+N9X9jwbw38M2u7JoQQYicaNwuNEEKIO0hAF0KIPUICuhBC7BES0IUQYo+oW7ZFpdQCcLOKJtqBRYu6sxvJ+2ts8v4a225+f0Na647NHqhbQK+WUurcVglq9gJ5f41N3l9ja9T3J1MuQgixR0hAF0KIPaKRA/oX690Bm8n7a2zy/hpbQ76/hp1DF0IIcadGHqELIYTYQAK6EELsEXsioCul/kgppZVS7fXui5WUUv+LUuqyUuotpdTXlVJ7ot6GUupppdQVpdSoUupz9e6PlZRSA0qpHymlLiml3lZK/fN698kOSimnUup1pdS36t0Xq62V0Hx+7WfvklLqXfXuU7kaPqArpQYwBazH690XG3wPOKW1fhBTqPu/rXN/qrZWdPzzwMeAE8CzSqkT9e2VpQrAf6W1Pg48CfyzPfb+1v1z4FK9O2GTPwH+Xmt9DDhNA73Phg/owP8B/DdsUvKu0Wmtv6u1Lqx9+jKmElSju1V0XGudA9aLju8JWusZrfVra/9OYIJBX317ZS2lVD/wj4Av1bsvVlNKhYH3Af8eQGud01pH69ur8jV0QFdKfRKY0lq/We++1MDvAt+udycssFVB8T1HKTUMPAy8Ut+eWO7fYAZRpXp3xAYHgAXg/1qbUvqSUqqeBUp3pKwCF/WklPo+0L3JQ/8K+JfAL9a2R9ba7v1prb+x9px/hflT/i9r2TeblFVQvNEppYLA14A/1FrH690fqyilPgHMa61fVUp9oN79sYELeAT4A631K0qpPwE+B/z39e1WeXZ9QNdaf3iz+5VSDwAjwJtKKTDTEa8ppR7XWs/WsItV2er9rVNK/RPgE8CH9N44NLDnC4orpdyYYP6XWuu/rnd/LPYe4JNKqY8DPiCslPoLrfV/Wud+WWUSmNRar/9V9TwmoDeEPXOwSCl1Azijtd6tGdJ2TCn1NPC/A+/XWi/Uuz9WUEq5MAu8HwKmMEXIf0tr/XZdO2YRZUYXXwaWtdZ/WO/+2GlthP5HWutP1LsvVlJKvQj8vtb6ilLqfwSatNb/dZ27VZZdP0Lf5/4t4AW+t/ZXyMta6+fq26XqbFV0vM7dstJ7gN8Bziul3li771+u1eUVjeEPgL9USnmAMRqo6P2eGaELIcR+19C7XIQQQtwmAV0IIfYICehCCLFHSEAXQog9QgK6EELsERLQhRBij5CALoQQe8T/D9dHC7JclBF8AAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"f = plt.figure()\n",
"plt.plot(y_d, g_αn, alpha=0.5, label='Never takers')\n",
"plt.plot(y_d, g_αa, alpha=0.5, label='Always takers')\n",
"plt.plot(y_d, g_\\alph, alpha=0.5, label='Compliers')\n",
"plt.legend()\n",
"f.savefig(\"alpha.pdf\", bbox_inches='tight')"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}