{
"cells": [
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"nbsphinx": "hidden"
},
"outputs": [],
"source": [
"import statsmodels.formula.api as smf\n",
"import pandas as pd\n",
"import numpy as np\n",
"\n",
"from auxiliary import get_treatment_probability\n",
"from auxiliary import get_plot_probability\n",
"from auxiliary import plot_outcomes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Regression discontinuity design"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This following material is mostly based on the following review:\n",
"\n",
"* Lee, D. S., and Lemieux, T. (2010). [Regression discontinuity designs in economics](https://www.aeaweb.org/articles?id=10.1257/jel.48.2.281). *Journal of Economic Literature, 48*(2), 281–355.\n",
"\n",
"The idea of the authors is to throughout contrast RDD to its alternatives. They initially just mention selected features throughout the introduction but then also devote a whole section to it. This clearly is a core strength of the article. I hope to maintain this focus in my lecture. Also, their main selling point for RDD as the close cousin to standard randomized controlled trial is that the behavioral assumption of imprecise control about the assignment variable translates\n",
"into the statistical assumptions of a randomized experiment.\n",
"\n",
"**Original application**\n",
"\n",
"In the initial application of RD designs, Thistlethwaite & Campell (1960) analyzed the impact of merit rewards on future academic outcomes. The awards were allocated based on the observed test score. The main idea behind the research design was that individuals with scores just below the cutoff (who did not get the award) were good comparisons to those just above the cutoff (who did receive the award).\n",
"\n",
"\n",
"**Causal graph**\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Intuition"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Key points of RD design**\n",
"\n",
"- RD designs can be invalid if individuals can precisely manipulate the assignment variable - discontinuity rules might generate incentives\n",
"\n",
"- If individuals - even while having some influence - are unable to precisely manipulate the assignment variable, a consequence of this is that the variation in treatment near the threshold is randomized as though from a randomized experiment - contrast to IV assumption\n",
"\n",
"- RD designs can be analyzed - and tested - like randomized experiments.\n",
"\n",
"- Graphical representation of an RD design is helpful and informative, but the visual presentation should not be tilted toward either finding an effect or finding no effect.\n",
"\n",
"- Nonparametric estimation does not represent a \"solution\" to functional form issues raised by RD designs. It is therefore helpful to view it as a complement to - rather than a substitute for - parametric estimation.\n",
"\n",
"- Goodness-of-fit and other statistical tests can help rule out overly restrictive specifications."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Baseline**\n",
"\n",
"A simple way to estimating the treatment effect $\\tau$ is to run the following linear regression.\n",
"\n",
"\\begin{align*}\n",
"Y = \\alpha + D \\tau + X \\beta + \\epsilon,\n",
"\\end{align*}\n",
"\n",
"where $D \\in [0, 1]$ and we have $D = 1$ if $X \\geq c$ and $D=0$ otherwise."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Baseline setup**\n",
"\n",
"\n",
"\n",
"* \"all other factors\" determining $Y$ must be evolving \"smoothly\" (continously) with respect to $X$.\n",
"\n",
"* the estimate will depend on the functional form"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Potential outcome framework**\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Potential outcome framework**\n",
"\n",
"\n",
"Suppose $D = 1$ if $X \\geq c$, and $D=0$ otherwise\n",
"\\begin{align*}\n",
"\\Rightarrow\\begin{cases}\n",
"E(Y \\mid X = c) = E(Y_0 \\mid X = c) & \\text{for}\\quad X < c \\\\\n",
"E(Y \\mid X = c) = E(Y_1 \\mid X = c) & \\text{for}\\quad X \\geq c\n",
"\\end{cases}\n",
"\\end{align*}\n",
"\n",
"Suppose $E(Y_1\\mid X = c), E(Y_0\\mid X = c)$ are continuous in $x$.\n",
"\\begin{align*}\n",
"\\Rightarrow\\begin{cases}\n",
"\\lim_{\\epsilon \\searrow 0} E(Y_0\\mid X = c - \\epsilon) = E(Y_0\\mid X = c) \\\\\n",
"\\lim_{\\epsilon \\searrow 0} E(Y_1\\mid X = c + \\epsilon) = E(Y_1\\mid X = c) \\\\\n",
"\\end{cases}\n",
"\\end{align*}\n",
"\n",
"\\begin{align*}\n",
"&\\lim_{\\epsilon \\searrow 0} E(Y\\mid X = c + \\epsilon) - \\lim_{\\epsilon \\searrow 0} E(Y\\mid X = c - \\epsilon) \\\\\n",
"&\\qquad= \\lim_{\\epsilon \\searrow 0} E(Y_1\\mid X = c + \\epsilon) - \\lim_{\\epsilon \\searrow 0}E(Y_0\\mid X = c - \\epsilon) \\\\\n",
"&\\qquad=E(Y_1\\mid X = c) - E(Y_0\\mid X = c) \\\\\n",
"&\\qquad=E(Y_1 - Y_0\\mid X = c)\n",
"\\end{align*}\n",
"\n",
"$\\Rightarrow$ average treatment effect at the cutoff"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Sharp and Fuzzy design"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEiCAYAAAAF7Y7qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfQElEQVR4nO3de5gdVZnv8e8vNwiCoAQUcjGg4RIQBFpEQQkCGjgHcLxwQBkEGTijgp5RUI8XhoHRUZhRjoqXeEBEjwIyPBolwEgwRFGQyCWaQCDEAAlBEiARCbl08p4/qtqubHbvXrvZXVWd/D7PU8+uXbWq6u3Vu/fbq1bVKkUEZmZmKYZVHYCZmQ0dThpmZpbMScPMzJI5aZiZWTInDTMzS+akYWZmyZw0bFBIWiwpJJ1WdSx1Jum0vJ4WVx0LQB5LSJpSdSxWTyOqDsDqS5KAdwPvBQ4EdgY2AH8GlgG/A34FzIyIv1QVp5mVx0nDmpK0A/AT4PDC4m5gNTAB2B04FPgn4HTgylIDtMGyIH9dXWkUVltOGtaXq8gSxgbgUuDbwMMRsVHSCGAyMJWsFWKbiYjYq+oYrN6cNOwFJE0CjsvffjYivlhcHxHdwNx8uljS6JJDNLOKuCPcmnldYf6n/RWOiOdbrZc0StJ5ku6T9JykVZJulTS1xTa7SfqkpJskPZhv91dJ8yVdKmlCi21n5Z25F0gaKenjkuZIWlns5G3shJZ0tKQbJS2X9LykeZI+K2nr/uqgn5//EEk/kbQi3+8CSZ+XtG3i9jtJ+ldJ9+R1t0bSIkmXS9qnxXbjJH0l/zmek7RW0uOSfp8vf32TbVp2hEsak2+7KI9jmaQfSzqw1faSpvSsy9+/RtIVkh7L41oi6TuSxqbUiVUoIjx52mQC3gNEPh09wH0szrc/G7gjn18HPFvY90bgA31sP6tQbi2wguxUWc+ylcBh/Wz7ReD2fH498HR+zCl5udPydYuBD+XrAngmL99zrLuBlw2wHj7QJO61+fz9ZH1CASzuY/uj8nh6tl8H/LWhbk5tst3++c/bU6678PP3LLuyyXY966Y0WbcHsLRQZg2wqhDHcX1tD0wprDui8Dn4S0NdLwXGVv034KnFZ7rqADzVbwImFr5c5gJ7DGAfPUnjaWAJcAIwMl+3J/DbfP2zwPZNtr80/yKfBAzLl40ADgZuLHzBjG6y7azCvp/Nk8PofN2OwMvz+Z6k8Vz+ZXwtMD5fNxr4x/yLMYDrB1AHBxa+EH8J7JUvHwmclCeDnoSwuMn2ryXrkA5gGrA3MDxfNwG4jN6E2NWw7S35ut8DhwDKl4/K6/TjwHlNjtnXl/7I/LMQwHLg7wqx7AXMZNMk1bj9lMK6p8lasHsVYjqRLIEEcFXVfwOeWnyuqw7AUz2n/Euq2CK4O/+S+gCwb8+XUIvtF9P73+heTdbvBDyfl3lfm7ENB+7Ltz2lyfpZhdiPa7Gf0wrlZpEnp4YyZxTKvL7NOGfk2y2geXJ7e2Hfi5usn5mv+0KLY/yfvMxPGpb3JJs3thlzX1/6pxQ+C29ust3WZC2nlKRxax91fU6+fjUwouq/AU/NJ/dpWF8+BFxE9l+4gAPyZZcDfwCekPRlSa/oZz/XRcQDjQsjYjlZawNgv3YCi4gNwE3528NaFJ0XET9L3O2/RsTGJsu/S9ZSgqx1kCS/ZPnt+dtLokm/T0TcTG8dNG4/EXgr2Wmlf29xqKvy16MkDS8sX5m/7pIacz/ek7/OjohfNa6MiDXAJYn7+kIfdd3TfzaarDVkNeSkYU1FRHdEnA+MBf4e+L9k/92vy4vsTHY+/o+SDm6xqztbrHs8f315s5WS3izpSkkP5J3gUehM/URebFyL/d/eYl1RN9lNii+Qf7nNyt92Je4PslNTPX9ft7Yo19e6Q/PXYcB8SU80m+hNni8hO/XW4+f56/ck/YekwyVt00b8jQ7MX29rUWZW4r76+kw8Xphv+pmw6vmSW2spIlYBP8gn8iuJDgM+QtbxOQb4T0mT8v82Gz3bYvfd+evIxhWSvkRvYoCsM/kZepPWtmRflC9psf8nW6wrWhERa1usX5q/7py4v8ayS/ss1duKabRr/joM6K8116OYFD4BvIas0/lj+bRB0r3ADcC0iGgVV6Od8tfHW5RJ2l9ENP1MRER3NggB0OQzYfXgloa1JSLWRMQtEXE88L188TiyG/06QtLR9CaMb5B1CG8VES+PiFdGxCuBr/QUb7GrDZ2KqQI9p5r+HBFKnBb3bBwRKyPircCbgYvJWl3dwEHA+cBDkk4eQFzx4n4sG+qcNOzFmFaY37OD++3pO7g5Ij4cEX/M+zGKXtnB442RNKrF+p57B1JbLo1lW9170Ne6J/LXMZJataZaiohfR8QnI+IwYAeyq9j+QNZvcEVCn1SP5fnrri3K+B6LLYCThr0Yfy3Mtzq9067x+es9zVbmAym+tYPHG0H2H3lfxzo8fzunjX3eTXalEWSniPrS18/R0x8zHDimjeP2KW8lTgfemS/qOdWY4u78dUqLMq3W2WbCScNeIL8be4+Eou8vzN/dZ6n2rcpf9+9j/T+SDZjYSZ+R1Ozv4f30JrFrUncWESuB/8rfntvsrnJJRwFv6mP7h+jtWP68pO1bHU/SywvzI/r4WXoUr+RqdhVTM9flr2+RdGjjSklbAecm7suGMCcNa2Yf4H5JN0g6Nb/8E4B8WI4DJH2XrHMVsiHSf93B4/dcEXSMpM/1nJ6RtIOkTwNfA57q4PFWk/3H/UNJ4/JjbS3pLOCbeZmfRsTv2tzv58j6VfYCbpC0Z77vEZJOJLuZcGWL7c8ha83tAdwh6YRi8pE0VtLfS5oJfKmw3TiyPovP5r+rEYVt9iO/qIHscupWV0MVXQPMI+tDuj6PZXi+zz3Jrtbq5ClDq6uqbxTxVL+JTW86Kw5X8RSbDkPRc8fxrk32sThff1qL41xJk+EsyK6cmV04xkayu4h7huP4Odk9JAHMarLfWfm6C/r5OU+jdxiRDxd+tqfJrtLqOf69wI4DrMuzGupsJb13macMI3Io2bNLerbvJhtSZTWb/h6+U9hmYsO67vx3t5ZNf5/vbnK8pjfn5ev2aohlTf7z9Mz/98K6Qxq2ndKzrp/66vP4nuoxuaVhLxDZTWeTgI8CPyb7cltL1pG6GniI7L/kk8jukm51GeZAjr8eeBvwL8CDZMNkiKxF80HgeDp8ZVREXEaWLG8i+5LfCDxAdqXRGyNiQC2biJhG9sX/M7JktBXwCPBvZEOiPNPP9reTtTTOJUukK8l+DxvIfi8/AN4H/K/CZkvJ6ugrZON+LSO7RLkbmE92Z/++EXEdbYjsJs39gK+SJVqRJYtryYYqKd4Xs7KdfdvQ0TMejdkWR9mjaL8LPBIRE6uNZujLL5X+L7JE8tI8+dtmxi0NM3vR8qvMPpm/vdUJY/PlpGFmSSQdoexZJl3KH7ylzEFkp9+OJOuPuLjKOG1weRgRM0u1PVk/10cBJD1DdpNgzxVdAZwbEalXZNkQ5KRhZqnuILuM+Eiy+2R6xqNaRDbg49cjop0bIG0IGtId4VOnTo2bbrqp/4JmZlbUasy2loZ0n8aKFSuqDsHMbIsypJOGmZmVy0nDzMySOWmYmVkyJw0zM0vmpGFmZsmcNMzMLJmThpmZJXPSMDOzZE4aZmaWzEnDzMySOWmYmVkyJw0zM0vmpGFmZslKSRqSrpD0pKQ/9rFekr4qaaGkuZIOLCMuMzNrT1kPYboS+DpwVR/rjwEm5dMbgG/mr1aSlavXsXTl81WHYWYl2GfX7Qe8bSlJIyJmS5rYosgJwFWRPRHqDkk7SNolIpaVEd+W7rYHl3PmVXNY172x6lDMrASLv/jfBrxtXR73OhZ4rPB+Sb7sBUlD0lnAWQATJkwoJbjN3fR7H2cyj/3t03Bv99hqAzKz2qpL0kgWEdOAaQBdXV1D91m1NdK9cSMHjOzNz2t32rvCaMyszuqSNJYC4wvvx+XLrASNj4m/8aNvriYQM6u9ulxyOx04Nb+K6hBglfszzMzqp5SWhqQfAVOAMZKWAP8MjASIiG8BM4BjgYXAauD0MuKyjM/xmVmqsq6eOrmf9QF8uIxYzMxs4OpyesrMzIYAJw0jGnvCzcz64KRhZmbJ6nLJrVUogEc3ZMMK7Psihhcws82fk4YBMHPdJADe8cYDKo7EzOrMp6fM19yaWTInDduEVHUEZlZnThpmZpbMfRpGELxuRDbU1+P3b4T9dq04IjOrKycNA/jbKLfLHlgGHFdtMGZWWz49ZS8Y5dbMrC9OGmZmlsxJw9zSMLNkThpmZpbMScPMzJI5aRjhW8LNLJGThpmZJXPSMHeEm1kyJw0zM0vmO8KNABZ0jwHgsNeMqTYYM6s1Jw0D4DfrJwLw/q6Dqg3EzGrNp6fMzCyZk4Zt0hHux2mYWStOGmZmlsx9GgYEbxq5GIBFv38O9jmx2nDMrLacNAyAPUesAODPf1pRcSRmVmc+PWW+uc/MkjlpmJlZMicNMzNL5qRhHuPWzJI5aZiZWbLSkoakqZIWSFoo6VNN1k+Q9EtJ90iaK+nYsmLb0oV7ws0sUSlJQ9Jw4DLgGGAycLKkyQ3FPgtcGxEHACcB3ygjNjMzS1dWS+NgYGFELIqIdcDVwAkNZQJ4aT6/PfB4SbFt8dzOMLNUZSWNscBjhfdL8mVFFwCnSFoCzADOabYjSWdJmiNpzvLlywcjVjMz60OdOsJPBq6MiHHAscD3Jb0gvoiYFhFdEdG10047lR7k5uqe9btwz/pdmLBPV9WhmFmNlZU0lgLjC+/H5cuKzgCuBYiI3wJbA34iUAki4N7usdzbPZaJr3191eGYWY2VlTTuAiZJ2k3SKLKO7ukNZR4FjgSQtDdZ0vD5JzOzGiklaUREN3A2cDNwP9lVUvMkXSjp+LzYx4EzJd0H/Ag4LXwtaCmKleznaZhZK6WNchsRM8g6uIvLzi/MzwcOLSseMzNrX506wq0iEcGRox7iyFEP8YfZM/rfwMy2WH6ehgEwYfgqAJ5+fFXFkZhZnbmlYWZmyZw0zMwsmZOGmZklc9IwP+7VzJI5aZiZWTInDSM8zq2ZJXLSMDOzZE4aZmaWzEnD3BFuZsl8R7gBcPu6VwFwxmG7VRyJmdWZk4YRAQ9uyB5oNX7SPhVHY2Z15tNTZmaWzEnDfMmtmSVz0rBN+SlMZtaCk4YBcNxW8zluq/ncMePaqkMxsxpL7giXtGNEPDWYwVg1ImDMsNUAPPv06oqjMbM6a6el8aikn0p6t6RRgxaRmZnVVjtJYyIwE/gk8ISkaZIOG5SorFTuBjezVMlJIyKWR8RXI+L1wBuBJ4HvS1ok6UJJrxq0KM3MrBYG2hH+ynx6KfAwMBa4R9KnOhWYlchNDTNL1E5H+D7AKcB7geeA7wH7R8SSfP1FwFzgi4MQp5mZ1UA7w4jMBn4EvCcifte4MiIWS7q0U4GZmVn9tJM0/i4iZjculHRwTxKJiPM7FpmVxneEm1mqdvo0ft7H8ps6EYiZmdVfvy0NScPIBpeQJLHpQBOvBroHKTYriZ+nYWapUk5PddN7fU1jgtgIfL6jEVklbln7GgA+ecxeFUdiZnWWkjR2I2td3Aa8pbA8gOUR8fxgBGblCeCxjTsA8IrxfgiTmfWt36QREY/ks755bwvgQW7NrJWWSUPStIg4K5+/qq9yEXFqpwMzM7P66e/qqT8V5h9uMfVL0lRJCyQt7OvOcUknSpovaZ6kH6bs1168KPSEy00NM2uhZUsjIv6tMP8vAz2IpOHAZcDRwBLgLknTI2J+ocwk4H8Dh0bEM5J2HujxrH3/Y+v7APjFNfdz0CfOrTgaM6ur/k5PvTVlJxFxaz9FDgYWRsSifL9XAycA8wtlzgQui4hn8n0+mXJse/EC2EbrAVj7/PpqgzGzWuuvI/zyhH0EsHs/ZcYCjxXeLwHe0FBmDwBJtwPDgQsi4gU3Dko6CzgLYMKECQnhmZlZp/R3eqrM6y9HAJOAKcA4YLak10bEyoaYpgHTALq6unxbWgf45j4zS1XWM8KXAuML78fly4qWANMjYn1E/Al4kCyJmJlZTbRMGpLuL8w/JunRZlPCce4CJknaLX9U7EnA9IYyPyFrZSBpDNnpqkXpP4qZmQ22/vo0zizMnzLQg0REt6SzgZvJ+iuuiIh5ki4E5kTE9Hzd2yTNBzYA50XEUwM9pqXz2SkzS9Vfn8avC/O3vZgDRcQMYEbDsvML8wF8LJ/MzKyGkvs0JI3KnwX+kKTn8teLJG09mAFaCdwTbmaJ2nkI0zeBPYGPAI+QjUX1abLLaT/Q+dDMzKxu2kka7wBeXbgEdr6kO4GFOGkMaW5nmFmqdpLGE8A2wMrCstHAsk4GZNWYvmZvAL584v4VR2JmddbOMCLfB26S9DWyeyrGAx8G+hz91oaOp+IlAOww5hUVR2JmdTaQYUQ+3fD+fwJf6kw4VoViP7hHuTWzVuo0jIiZmdVcWcOIWI1FoStcfnafmbXQzn0aL5X0ZUm/l/RIm8OIWM2dPnoOp4+ew/WXX1p1KGZWY+20NL4BHAhcCLwcOAd4FPjKIMRlJfK9fWaWqp1Lbt8G7B0RT0naEBE/lTQH+BlOHGZmW4R2WhrDgFX5/F8lbU92j8ZrOh6VmZnVUjstjfuAw4GZwK/ITlf9ley5FzaE+fSUmaVqp6VxJrA4n/8osAbYATi1syGZmVldJbc0ImJRYf5J4IxBichK54aGmaVq6z4NSR+Q9AtJ8/LXMyTfQ2xmtqVIbmlIuhg4AbiU3qHRzyUbLv0TgxGclSPcqWFmidrpCD8NODAilvQskPRz4G6cNMzMtgjtnJ56Np8al/2lc+GYmVmd9Tc0+u6Ft5cC10v6Ir1Do5+Hb+zbLFz9/H4A/L9/eEPFkZhZnfV3emoh2cU1xc7uIxrKvBX4eieDsvI9zygARr9k24ojMbM6629odI+CuwXw8zTMLFU7HeEASJoAjAWWRMRjnQ/JzMzqqp2h0XeRdBvZKavrgYclzZa066BFZ6UIgtGsYzTrWPPcc1WHY2Y11s7pp2+SjT/1sojYBXgZcA/wrcEIzMp10ui5nDR6Lj/9wbSqQzGzGmvn9NRhwC4RsR4gIp6T9Alg6aBEZmZmtdNOS+MZYHLDsj2BlR2LxirhG8LNLFU7LY2LgVskXU7vMCKnA58bjMDMzKx+2hnl9juSHgbeC+wHPA68NyJmDlZwVg43NMwsVVLSkDSc7GFLkyPi1sENyczM6iqpTyMiNgAbgK0HNxyrgke5NbNU7XSEXwpcK+lwSa+WtHvPlLKxpKmSFkhaKOlTLcq9S1JI6mojNjMzK0E7HeE940sd3bA8gOGtNsxPb12Wb7sEuEvS9IiY31BuO7JHyd7ZRlxmZlaSflsakraR9AXgBuAiYJuIGFaYWiaM3MHAwohYFBHrgKvJHujU6CLgS2TPH7eS+OSUmaVKOT11GXAccD/wLuCSARxnLFAcp2pJvuxvJB0IjI+IG1rtSNJZkuZImrN8+fIBhGJmZgOVcnpqKtkT+5ZJ+howGzink0FIGgZ8mezpgC1FxDRgGkBXV5f/Se6EgO8+n3Uhzfz44RUHY2Z1ltLSeElELAPIR7XdfgDHWUr20KYe49h0+JHtgH2BWZIWA4cA090ZbmZWLyktjRGSjqD3QUyN70m4d+MuYJKk3ciSxUlkNwn2bL8KGNPzXtIs4NyImJPyQ9iLU2yu+XEaZtZKStJ4Erii8P6phvcBtLzsNiK6JZ0N3Ex2pdUVETFP0oXAnIiY3l7YZmZWhX6TRkRM7MSBImIGMKNh2fl9lJ3SiWNauh2VPUfjqSf/zO47+ZGvZtacH+dqRATHb30/x299PzOu+0HV4ZhZjTlpmJlZMicN8819ZpbMScPMzJI5aZif3GdmyZw0zMwsmZOGmZklc9Iwwl3hZpbIScPMzJI5aZg7ws0sWTtP7rPN2OoYCcCO246qOBIzqzMnDSMCrlmzPwC/+sgRFUdjZnXm01NmZpbMScPMzJI5adgm5KcwmVkLThoGwPhhKxk/bCWPLlpYdShmVmPuCDcigqO2ypLFL36+kDcdtF/FEZlZXbmlYWZmyZw0zIOImFkyJw0zM0vmpGFmZsmcNMxjT5lZMicNMzNL5qRhfp6GmSVz0jAzs2ROGuY+DTNL5jvCDYAVG7cBYK9XbldxJGZWZ04aBsDP1k4G4LPvO7LiSMysznx6yjbpBvcgt2bWipOGmZklc9KwTTvC3dQwsxZKSxqSpkpaIGmhpE81Wf8xSfMlzZU0U9KryorNYI/hy9lj+HLmz7236lDMrMZKSRqShgOXAccAk4GTJU1uKHYP0BUR+wHXAReXEZsBBIeOeoRDRz3CbbfcXHUwZlZjZbU0DgYWRsSiiFgHXA2cUCwQEb+MiNX52zuAcSXFZmZmicpKGmOBxwrvl+TL+nIGcGOzFZLOkjRH0pzly5d3MEQzM+tP7TrCJZ0CdAGXNFsfEdMioisiunbaaadyg9tM+Y5wM0tV1s19S4Hxhffj8mWbkHQU8Bng8IhYW1JsZmaWqKyWxl3AJEm7SRoFnARMLxaQdADwbeD4iHiypLgMP+7VzNKVkjQiohs4G7gZuB+4NiLmSbpQ0vF5sUuAbYEfS7pX0vQ+dmdmZhUpbeypiJgBzGhYdn5h/qiyYrFNhTs1zCxR7TrCzcysvpw0zMwsmYdGNwJ4dMP2ALxlki9jNrO+OWkYADPXTQLg3995dMWRmFmd+fSUbXJznzzKrZm14KRhZmbJnDRsk0tu5QdqmFkLThoGwOtGLOV1I5by29tnVx2KmdWYk4YBcMDIZRwwchl3/ubXVYdiZjXmpGEee8rMkjlpmJlZMicNc1PDzJI5aZiZWTInDXNDw8ySOWmYmVkyJw3z8zTMLJmThpmZJXPSMDOzZB4a3QhgQfcYAN5z0LhqgzGzWnPSMAB+s34iAN855u3VBmJmtebTU+bnaZhZMicNMzNL5qRhhG/vM7NEThoGwJtGLuZNIxdz840zqg7FzGrMScMA2HPECvYcsYK5995TdShmVmNOGoZvCDezVE4aZmaWzEnD3A1uZsmcNMzMLJmThrmpYWbJnDTMzCyZk4aZmSUb0gMWLn7qOc648q6qwxjy1m/cWHUIZjZEDOmk8eyabmY+8GTVYZiZbTE0lB/1KelZYEHVcdTEGGBF1UHUhOuil+uil+ui19YRse9ANhzSLQ1gQUR0VR1EHUia47rIuC56uS56uS56SZoz0G3dEW5mZsmcNMzMLNlQTxrTqg6gRlwXvVwXvVwXvVwXvQZcF0O6I9zMzMo11FsaZmZWIicNMzNLNiSShqSpkhZIWijpU03WbyXpmnz9nZImVhBmKRLq4mOS5kuaK2mmpFdVEWcZ+quLQrl3SQpJm+3llil1IenE/LMxT9IPy46xLAl/IxMk/VLSPfnfybFVxDnYJF0h6UlJf+xjvSR9Na+nuZIOTNpxRNR6AoYDDwO7A6OA+4DJDWU+BHwrnz8JuKbquCusiyOAbfL5D27JdZGX2w6YDdwBdFUdd4Wfi0nAPcDL8vc7Vx13hXUxDfhgPj8ZWFx13INUF28BDgT+2Mf6Y4EbAQGHAHem7HcotDQOBhZGxKKIWAdcDZzQUOYE4Hv5/HXAkZJUYoxl6bcuIuKXEbE6f3sHMK7kGMuS8rkAuAj4ErCmzOBKllIXZwKXRcQzABGxuY6/k1IXAbw0n98eeLzE+EoTEbOBp1sUOQG4KjJ3ADtI2qW//Q6FpDEWeKzwfkm+rGmZiOgGVgE7lhJduVLqougMsv8kNkf91kXe3B4fETeUGVgFUj4XewB7SLpd0h2SppYWXblS6uIC4BRJS4AZwDnlhFY77X6fAEN/GBHrg6RTgC7g8KpjqYKkYcCXgdMqDqUuRpCdoppC1vqcLem1EbGyyqAqcjJwZUT8h6Q3At+XtG9EeLjnBEOhpbEUGF94Py5f1rSMpBFkTc6nSomuXCl1gaSjgM8Ax0fE2pJiK1t/dbEdsC8wS9JisnO20zfTzvCUz8USYHpErI+IPwEPkiWRzU1KXZwBXAsQEb8FtiYbzHBLk/R90mgoJI27gEmSdpM0iqyje3pDmenA+/P5dwO3Rt7Ts5npty4kHQB8myxhbK7nraGfuoiIVRExJiImRsREsv6d4yNiwAO11VjK38hPyFoZSBpDdrpqUYkxliWlLh4FjgSQtDdZ0lheapT1MB04Nb+K6hBgVUQs62+j2p+eiohuSWcDN5NdGXFFRMyTdCEwJyKmA5eTNTEXknX8nFRdxIMnsS4uAbYFfpxfC/BoRBxfWdCDJLEutgiJdXEz8DZJ84ENwHkRsdm1xhPr4uPAdyT9E1mn+Gmb4z+Zkn5E9o/CmLz/5p+BkQAR8S2y/pxjgYXAauD0pP1uhnVlZmaDZCicnjIzs5pw0jAzs2ROGmZmlsxJw8zMkjlpmJlZMicNMzNL5qRh1iGStpW0WNL7Csu2k/SopHdXGZtZp/g+DbMOkvR24Adkw3Evl/RN4BUR8c6KQzPrCCcNsw6TdCWwFdlwLv8J7BMRT1QalFmHOGmYdZiklwHzyYZsOC8ivltxSGYd4z4Nsw7LH3Q0D9gGuL7icMw6yknDrMPyZ5lMBG4he2qg2WbDp6fMOkjSzmStjBOBB/L5EyLiV5UGZtYhThpmHSTpWrLnEpyZv/8H4Fxg/834gVi2BXHSMOsQSe8AvkF2ue3KwvJbgd9GxGcqCs2sY5w0zMwsmTvCzcwsmZOGmZklc9IwM7NkThpmZpbMScPMzJI5aZiZWTInDTMzS+akYWZmyf4/Gllqz0QFY3YAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"for version in [\"sharp\", \"fuzzy\"]:\n",
" plot_outcomes(version, grid)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Alternatives**\n",
"\n",
"Consider the standard assumptions for matching:\n",
"\n",
"- ignorability - trivially satisfied by research design as there is no variation left in $D$ conditional on $X$\n",
"- common support - cannot be satisfied and replaced by continuity\n",
"\n",
"Lee and Lemieux (2010) emphasize the close connection of RDD to randomized experiments.\n",
"- How does the graph in the potential outcome framework change?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Continuity, the key assumption of RDD, is a consequence of the research design (e.g. randomization) and not simply imposed."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Identification\n",
"\n",
"Ad-hoc $\\times$ vs. thoughtful answers $\\checkmark$. Both are true, but only thoughtful consideration clarifies the strength of the regression discontinuity design as opposed to, for example, an instrumental variables approach."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Question**\n",
"\n",
"How do I know whether an RD design is appropriate for my context? When are the identification assumptions plausable or implausable?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Answers**\n",
"\n",
"$\\times$ An RD design will be appropriate if it is plausible that all other unobservable factors are \"continuously\" related to the assignment variable.\n",
"\n",
"$\\checkmark$ When there is a continuously distributed stochastic error component to the assignment variable - which can occur when optimizing agents do not have \\textit{precise} control over the assignment variable - then the variation in the treatment will be as good as randomized in a neighborhood around the discontinuity threshold."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Question**\n",
"\n",
"Is there any way I can test those assumptions?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Answers**\n",
"\n",
"$\\times$ No, the continuity assumption is necessary so there are no tests for the validity of the design.\n",
"\n",
"$\\checkmark$ Yes. As in randomized experiment, the distribution of observed baseline covariates should not change discontinuously around the threshold."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Simplified setup**\n",
"\n",
"\\begin{align*}\n",
"Y & = D \\tau + W \\delta_1 + U \\\\\n",
"D & = I [X \\geq c] \\\\\n",
"X & = W \\delta_2 + V\n",
"\\end{align*}\n",
"\n",
"- $W$ is the vector of all predetermined and observable characteristics.\n",
"\n",
"What are the source of heterogeneity in the outcome and assignment variable?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The setup for an RD design is more flexible than other estimation strategies.\n",
"- We allow for $W$ to be endogenously determined as long as it is determined prior to $V$. This ensures some random variation around the threshold.\n",
"- We take no stance as to whether some elements $\\delta_1$ and $\\delta_2$ are zero (exclusion restrictions)\n",
"- We make no assumptions about the correlations between $W$, $U$, and $V$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Local randomization**\n",
"\n",
"We say individuals have imprecise control over $X$ when conditional on $W = w$ and $U = u$ the density of $V$ (and hence $X$) is continuous."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Applying Baye's rule**\n",
"\n",
"\\begin{align*}\n",
"& \\Pr[W = w, U = u \\mid X = x] \\\\\n",
"&\\qquad\\qquad = f(x \\mid W = w, U = u) \\quad\\frac{\\Pr[W = w, U = u]}{f(x)}\n",
"\\end{align*}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Local randomization:** If individuals have imprecise control over $X$ as defined above, then $\\Pr[W =w, U = u \\mid X = x]$ is continuous in $x$: the treatment is \"as good as\" randomly assigned around the cutoff.\n",
"\n",
"$\\Rightarrow$ the behavioral assumption of imprecise control of $X$ around the threshold has the prediction that treatment is locally randmized."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Consequences**\n",
"\n",
"- testing prediction that $\\Pr[W =w, U = u \\mid X = x]$ is continuous in $X$ by at least looking at $\\Pr[W =w\\mid X = x]$\n",
"- irrelevance of including baseline covariates"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Interpretation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Questions**\n",
"\n",
"To what extent are results from RD designs generalizable?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Answers**\n",
"\n",
"$\\times$ The RD estimate of the treatment effect is only applicable to the subpopulation of individuals at the discontinuity threshold and uninformative about the effect everywhere else.\n",
"\n",
"$\\checkmark$ The RD estimand can be interpreted as a weighted average treatment effect, where the weights are relative ex ante probability that the value of an individual's assignment variable will be in the neighborhood of the threshold."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Alternative evaluation strategies\n",
"\n",
"- randomized experiment\n",
"- regression discontinuity design\n",
"- matching on observables\n",
"- instrumental variables\n",
"\n",
"How do the (assumed) relationships between treatment, observables, and unobservable differ across research designs?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Endogenous dummy variable**\n",
"\n",
"\\begin{align*}\n",
"Y & = D \\tau + W \\delta_1 + U \\\\\n",
"D & = I[X \\geq c] \\\\\n",
"X & = W \\delta_2 + V\n",
"\\end{align*}\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"* By construction $X$ is not related to any other observable or unoservable characteristic."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"* $W$ and $D$ might be systematically related to $X$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"* The crucial assumptions is that the two lines in the left graph are actually superimposed of each other.\n",
"\n",
"* The plot in the middle is missing as all variables are used for estimation are not available to test the validity of identifying assumptions."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"* The instrument must affect treatment probablity.\n",
"* A proper instructment requires the line in the right graph to be flat."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Nonlinear expectation\n",
"\n",
"A nonlinear conditional expectation can easily lead to misleading result if the estimated model is based on an local linear regression. The example below, including the simulation code, is adopted from Cunningham (2021). This example is set up closely aligned with the potential outcome framework. "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"df = pd.DataFrame(columns=[\"Y\", \"Y1\", \"Y0\", \"X\", \"X2\"], dtype=float)\n",
"\n",
"# We simulate a running variable, truncate it at\n",
"# zero and restrict it below 240.\n",
"df[\"X\"] = np.random.normal(100, 50, 1000)\n",
"df.loc[df[\"X\"] < 0, \"X\"] = 0\n",
"df = df[df[\"X\"] < 280]\n",
"\n",
"df[\"X2\"] = df[\"X\"] ** 2\n",
"\n",
"df[\"D\"] = 0\n",
"df.loc[df[\"X\"] > 140, \"D\"] = 1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now simulate the potential outcomes and record the observed outcome. Note that there is no effect of treatment."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def get_outcomes(x, d):\n",
"\n",
" level = 10000 - 100 * x + x ** 2\n",
" eps = np.random.normal(0, 1000, 2)\n",
" y1, y0 = level + eps\n",
" y = d * y1 + (1 - d) * y0\n",
"\n",
" return y, y1, y0\n",
"\n",
"\n",
"for idx, row in df.iterrows():\n",
" df.loc[idx, [\"Y\", \"Y1\", \"Y0\"]] = get_outcomes(row[\"X\"], row[\"D\"])\n",
"\n",
"df = df.astype(float)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What about the difference in average outcomes by treatment status. Where does the difference come from?"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"D\n",
"0.0 9836.848389\n",
"1.0 21643.244614\n",
"Name: Y, dtype: float64"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df.groupby(\"D\")[\"Y\"].mean()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we are ready for a proper RDD setup."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: Y R-squared: 0.783\n",
"Model: OLS Adj. R-squared: 0.782\n",
"Method: Least Squares F-statistic: 1795.\n",
"Date: Wed, 07 Jul 2021 Prob (F-statistic): 0.00\n",
"Time: 08:59:13 Log-Likelihood: -9407.4\n",
"No. Observations: 1000 AIC: 1.882e+04\n",
"Df Residuals: 997 BIC: 1.884e+04\n",
"Df Model: 2 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"Intercept 4377.3226 241.423 18.131 0.000 3903.568 4851.077\n",
"D 5889.2320 319.611 18.426 0.000 5262.044 6516.420\n",
"X 68.1328 2.699 25.247 0.000 62.837 73.429\n",
"==============================================================================\n",
"Omnibus: 799.728 Durbin-Watson: 2.046\n",
"Prob(Omnibus): 0.000 Jarque-Bera (JB): 27622.480\n",
"Skew: 3.372 Prob(JB): 0.00\n",
"Kurtosis: 27.849 Cond. No. 430.\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: Y R-squared: 0.973\n",
"Model: OLS Adj. R-squared: 0.973\n",
"Method: Least Squares F-statistic: 1.215e+04\n",
"Date: Wed, 07 Jul 2021 Prob (F-statistic): 0.00\n",
"Time: 08:59:13 Log-Likelihood: -8357.0\n",
"No. Observations: 1000 AIC: 1.672e+04\n",
"Df Residuals: 996 BIC: 1.674e+04\n",
"Df Model: 3 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"Intercept 1.005e+04 107.870 93.125 0.000 9833.768 1.03e+04\n",
"D 2.2077 131.768 0.017 0.987 -256.368 260.783\n",
"X -99.9678 2.202 -45.405 0.000 -104.288 -95.647\n",
"X2 0.9959 0.012 84.522 0.000 0.973 1.019\n",
"==============================================================================\n",
"Omnibus: 0.261 Durbin-Watson: 2.002\n",
"Prob(Omnibus): 0.878 Jarque-Bera (JB): 0.164\n",
"Skew: -0.004 Prob(JB): 0.921\n",
"Kurtosis: 3.062 Cond. No. 6.81e+04\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"[2] The condition number is large, 6.81e+04. This might indicate that there are\n",
"strong multicollinearity or other numerical problems.\n"
]
}
],
"source": [
"for ext_ in [\"X\", \"X + X2 \"]:\n",
" rslt = smf.ols(formula=f\"Y ~ D + {ext_}\", data=df).fit()\n",
" print(rslt.summary())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In a nutsheel, the misspecification of the model for the conditional mean functions results in flawed inference."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Estimation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Lee (2008)\n",
"\n",
"The author studies the \"incumbency advantage\", i.e. the overall causal impact of being the current incumbent party in a district on the votes obtained in the district's election.\n",
"\n",
"* Lee, David S. (2008). Randomized experiments from non-random selection in U.S. House elections. Journal of Econometrics."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
vote_last
\n",
"
vote_next
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
0.1049
\n",
"
0.5810
\n",
"
\n",
"
\n",
"
1
\n",
"
0.1393
\n",
"
0.4611
\n",
"
\n",
"
\n",
"
2
\n",
"
-0.0736
\n",
"
0.5434
\n",
"
\n",
"
\n",
"
3
\n",
"
0.0868
\n",
"
0.5846
\n",
"
\n",
"
\n",
"
4
\n",
"
0.3994
\n",
"
0.5803
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" vote_last vote_next\n",
"0 0.1049 0.5810\n",
"1 0.1393 0.4611\n",
"2 -0.0736 0.5434\n",
"3 0.0868 0.5846\n",
"4 0.3994 0.5803"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df_base = pd.read_csv(\"../../datasets/processed/msc/house.csv\")\n",
"df_base.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's put in some effort to ease the flow of our coming analysis."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"df_base.rename(columns={\"vote_last\": \"last\", \"vote_next\": \"next\"}, inplace=True)\n",
"\n",
"df_base[\"incumbent_last\"] = np.where(df_base[\"last\"] > 0.0, \"democratic\", \"republican\")\n",
"df_base[\"incumbent_next\"] = np.where(df_base[\"next\"] > 0.5, \"democratic\", \"republican\")\n",
"\n",
"df_base[\"D\"] = df_base[\"last\"] > 0\n",
"\n",
"for level in range(2, 5):\n",
" label = \"last_{:}\".format(level)\n",
" df_base.loc[:, label] = df_base[\"last\"] ** level"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The column `vote_last` refers to the Democrat's winning margin and is thus bounded between $-1$ and $1$. So a positive number indicates a Democrat as the incumbent."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### What are the basic characteristics of the dataset?"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEJCAYAAAB7UTvrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAByFklEQVR4nO29eXwV1f3//zozd8keYsIWAggEiglCVBQQxAK1RQS0HxGtWNqP29d+xE8/reJSP4ror4uiXRQ+ban106JYi/ipbFqXggsISNQEE0SIbFmUJYZAtrvMnN8fc+dmljP3zr2Ze3OTnOfj0Upm5s6cmTlz3ue8V0IpBYfD4XA4ACB0dwM4HA6HkzpwocDhcDicMFwocDgcDicMFwocDofDCcOFAofD4XDCuLq7AV1h9uzZ9J///Gd3N4PD4XB6GsRqR49eKZw6daq7m8DhcDi9ih4tFDgcDofjLFwocDgcDicMFwocDofDCcOFAofD4XDCcKHA4XA4nDBJcUklhDwHYC6AE5TScYz9BMDvAMwB0Abgh5TSjxPVnsYWH+qa2lGUl478LG+iLsPpJuy+30T2A/XcmR4RrX4p/N9AUMKRxjaUDe2HvEwPqhuaUd/UDl9QQr90Nw6dasOEoly4XUJ4+9C8DDS1BVA2tB+KB2ab2n74ZAve2HccWR4RuRkeDM1LR21TO4KSDF9QQnN7EPsamuESBPTL9CDL60J+hguf1DUjL92Nsz4Jueku+IMyGlv9mHRuHtqCFAOyvRicm46vmtux/eApnDzrw+B+aQjKMhpb/OiX7sG04nycbA3g5Jk2VH/ZgjSXgH6Zbkwcfg7a/RLe/uw4OgISxhf1w1XjB2PKqAJUHGvCxsovMbIgA21+CdtrTsEflJCf6cXkUfkAgPIjTRBFYFRBFnLT3TjZ4kf/LA9yMzzw+YP44PDXGJaXgVEDstAv3YXtNadwpj2I0QOzMHpAFqq/PAMCgsOnWtHul3D5N/rjgqH9UNXQjIKsNIwdlI1Wv4RdNSfx9ucncfX4wRg7OAevVjTALRIM6ZeOpjYfKAhKBudgUE4adh/+GqdafEj3iBial4EheRkozE0Lv9+G5nYABIW5adj/1VkcbWxBfqYXFAi/x9qmdnhdArLT3MhJd6O0MAf5WV6UH27ES3tq0RGQMKRfOmQAw/PScbSpDZkeF/ySDALALQqobjiDtz47Ee5rR351laN9lyQjSyohZDqAFgBrLITCHAB3QREKkwD8jlI6Kdp5J06cSMvLy2Nqy4aKetz3yl64BQEBWcYT147H/LIhMZ2Dk7rYfb+J7AfquQGgIyBDJIBEEf6vCgEQ69e3eMowXDT8nHDbW3zBmM/BcRbje40FlwCM6p+Jz4+3dqkNcQgGyziFpAgFACCEnAtgs4VQ+COAdyilfwv9/TmAb1JKv4x0zliFQmOLD1Mf34qOgBzeluYWsOO+mXzF0Auw+34T2Q9Y53Yajwj4pYSdntNDiVEwpHzw2hAAtZq/60LbTBBCbieElBNCyk+ePBnTReqa2uEW9LfsFgTUNbXH2FxOKmL3/SayH7DO7TQkZT5bTm+kx/UuSulqSulESunE/v37x/Tborx0dAT1U6yOoISivHQnm8jpJory0hGQ9TP0gCyb3q/d45xqg9NQJPb8nL5NqgiFegBDNX8XhbY5jlFdxivP9R7ys7x44trxSHMLyPa6kOYW8MS1400qIbvHdbUNaW7l8xJDC3XRcsFun8VThuHJ68rCbXfglJwu4sR77Srr/99kx86VKgnxNgJYQgh5CYqhuTmaPSEe6prake524awvGN6W7nahrqmd2xR6CfPLhmBqcUFUryK7x3W1DSzvo//7pB6vfFwPl0gQCMr4twuHYOyg7Ji8j7RtN3ofNZxuw3M7jkIUCCSZYvroAviDUkzeR41tQby05xiCNhYlJPQ/7aEugeDGi4dixxenbHkfNbcH0HC6Q3cOAmDu+YMQlCn+WX0cbkEx6M4ZNwhnfEGm9xEFxbb9JwECSDIwrjAbuemeuLyP9n91Bhsrv4RfY0V2iwTzJwzGxOHnhL2PKmpPY9nGKrT6NTYqF8GwczJw4ESnAXnmmP6Y/o3+Ye8jgOLulyvhC3Z9Ytovw9Plc6gky/vobwC+CaAAwHEAywC4AYBS+oeQS+pKALOhuKT+O6U0qgWZG5o5PZVEu0U74Za7saIe94a8nPyShCUzRuPGScOwo+ZUeLvquQVAc6yMay8cgpunjggLMTvtZX2bm5dMw9yV22P6ZuN5tqzfOOG40NTqR0XtaZ1A115v9+FGPPnmAbgEgrYueA88Mq8EP5w6IpafdL/3USKIxyVV29G5SyqnL2PHLddqgLUaRB969VO8VnU8fNziKcPw6NXn22oP69scnp+JhX/8QDebTnML+PvtUzBhaD9bbY12TKTnYHe8iGVcMbose0QCCiAQr18rgAfnjMVt00fF8hMuFLTw4DVObyHevpyIVXPN8bP41m/eM21/+yfTY1oxaO/nj+9+gV++vj/qOe0IONYxU4sLoj4HJ4Mhu+Ky7BKAsYOyUdVw1rTvP745EvfOPi+W01kKhVSxKSSV/CwvFwacHo92kNOqd+z0bdV1VqvFV91y4/02KmpPW263KxS032Zjiw9PvXXAdIxHAFo1qpbGFh/ue2UvOgJy+H7uebkSJYNzdBHgxmPufWUvVn//oqjPwe54wTqu5vhZbK85hYIsD6aMKuiS27NACM7Nz2AKhUyPGPd5jfRJocDhdDddXa2yBrmn3jqAldtqsGLB+KhG9ES45ZYZ1DnRtqv3YdXOuqZ2eEQCf9DwI0J07WQJOL9EMeeZ7XhygbJisBKCAHH0OWjv53dvH8CaXcfC+wQC3H/l2LgDG/0SxZZPjzP3jR2UE9c5WXChwOEkGSdSbLAGOQDwBWX8dF0FREGAR7Q+v+o6a9SDd2UFXTwwG4unDMOanZ0D4eIpwyxXCWt3HcXyzfvgEQmCMjW1sygvHUHZrN5eNq9U106r2BB/UMa9r+zF1OICSyFYWpjj2HPQr9xk+AyuWzIFnvinWRUWC1bK/v1fncGskkFdOrdKn7QpcDjdBUun7HUJ+OD+2HT5seimI9kKEmFfqzl+Vudxo70GoAi03Yca8QuDrcDjEvDaXdN0QkQ14IoCQUCiWDavBIsmDTddc2NFPe55uVLnPgoA2V4XXrh1EiYM7cc0Bhtdh+P1WAJg6324BNhy842Vf790OJbNN2UQigS3KXA4qQBrhu8Lynhx9zHcNWu07fOoM/2l66P7uUeyFdjRlxsFRzRBUjwwOzywa2fPHUEJlFJ4XaLOJqDiD8qY8/T7ePK6CeEVQyxxJyWDczDnme3wa0ZdrSpIPUYVWNVfnsHUx7fqhITRo8mKDRX1uHf9XggECMoUt04bwVy5sRAAWzHpLgGglECyMXE/1thm44z24EKBw0kgxgG0KC8dfsk8JKzcdtC2kVhFHTBf3H0MK7fVwCMqBmeZ6t0bu6IjVwd1kRAEJBnXXFCIjZVf2lJ9seweABCM4I/vl2hY5QMgPIu3Q/HAbDy5wFoVpAzklRCJgKAsgUJZfWgNz1OLC2wJyXtertQ949+/ewiiED20+aLheThvUA7+svNo1GO9bhEdNmMXGk5zocDhpDxWtoMlM4pNXjUuMT7Pn/wsL+6apXgdqcKHFVwW62pA3aYO6irrypXsM3YGUiu7RzTcgoC1u4/hf96pUa4VkOEVCYhAotpfrFYWjS0+3L2uIqS6UQda/QzcrvdVdUMzM6ZAYtg/jOw+3ITdh5tspU0PSBRpbvaqysj0MQOiHmMXLhQ4nARg5QI5tbgAN04ahpXbDurUPq0+CVX1zbbVF0a0aqBYU3hYCa+6pnaIJPLsVyAE1Q3NzEHJTnLAdLeAdoMe3i9JWLWtRmeo9UkUYKwiWPfHUolVN5yJqsv3SbKtVckHXzRGPSYaFErOJCFkKzHidQl4eF4JHtu8z9b5LhlxTpfbpJIqCfE4nB5JY4sPlbWn0dji022PlJ47P8uLh+eVms712JZ9pvPES36WFxOG9rO1QlCF11lfEB0BGUvX78V7B04i0yMiwFB1aWnzS7htTTk2VpjzV2qTA2Z6zYOtSIA/fv8i/GzOWHhcyjFpbgFLZoyGR2QPTeoqYurjW3HTs7sx9fGtzGubiT6LJ5Ri7srtludrbPHhvQMn8Nz2wzauFx23S8Cff3Ax7r5iDFya2xUJ8PBcxaAefn5RhNWuw10XVCpcKHB6FVaDdCLYUFFvOThFiwMYV5gLj+HrozJNem2PuqZ2uAy6cF9Qxi1/+RBzV27HNRcURj2HL6jM4FnPfH7ZEOy4byaWzys1DWwZHhcqapvx67cOwC0oyQEfuqoEV44bhDZTcIJCmz+IVdtqdELM6tpaSgtz4Y6SztQn0ZBQrDSdT33Xd7zwscnDKV4CQcUl9sZJw3T2CIl2ThDU5/fibZPxszljLc81siDLkTYBXChwehGRBmmnYc2wtYNTfpYXCycW6X6zcGJReOae6RHhN0zCfRJ1NDLVqt1aoVlV34wWn1lnHZAVXf7Gyi/xsyv1M/kfXT4SGW59OyMVKcrP8mLG2AEmLxpVTdQRkNHql+CXKJZtrMKcp9+3LG9JCDEJMZEQbNt/IqJgyM/y4qnrJsDrEuCNJhyCFC/uPhZ+VjXHz4bfdSxJ60SBIMMjWvp+Ugo0tfpDQXrWz1Nd9d0+fRT+a1Yx81wXn8vVRxyOjmiDtNNEq97W2OLDi7uP6farAw2gpGlQ6y2opLkFW0bFeDEKzbW7juKxLZF11m5BwKSR+dh5/0y8eOtk7LhvJm69bCRkgzpGuwpirdZYNSyWzBhtGuCDMiLOxD0uAT6DSqvVL+GRTdVRJwLzy4bgg/tn4oE50XME/ebtA7j0V8qzmvPMdlAbRmQjIqH40eWjUP7f38KPLh9p2i8DmPPMdlQ1NNuOqi4ewA4E3P+VOfVFvHChwOkVJKvUqjrgZXrEiB8yy7AZlJXtACxdRBNVBZAlNJdvqjYNykbUezLaKO78ZjG8LmIqUhRptaaqQl64dRJ23DcTHlfsQrDNL0EIrTi8rs62t/gkWxOB/CwvphUXRC2MI1NFjXbWF4Q/KCuG7hjxS8CqkAfVrZeNZF7TH5Tx6KZqPHRVia2iT0cbW5jXstoeD9z7iNMrSGSJTRWjl87CiUVYV15n4fppNYgo2xORZiISzNw/omA5K9e6gGrbpH0GAMHt00eG4ysieVwZk8ut3XWUmf3ULbK9cbSoajdJpuFI5PDvo7iVhtvvEiAF5FCAGCzVVSzsuJOqUJmiuuEMSgtzIAgEEuNCviDF161+7LhvZlSPMa+LrV602h4PXChwegWJHmRZA9668jpsXjKNmR5BNWwGDFW7Sgtzw38nsvqbtt1qAJhRaEpUSRvx2OZ9cAsCfEEJ1100FN+9oBBulxgecBtbfJaD/qp3anDjpGEA2IKH5bLa2OLD8k3VprZmeET86PJRWPVODURCoq4ivC4RAdlsp2hu94fbbHwWxriLoAykuQRIMeSeiGXN4JMobltTjiUzipHmEhGQ2Ab0lduU5xjNJbmpzR/T9njgQoHTa0jkIGuVZbPVLzE/ZNWwuXT93nBZzBUL2PWi482SGu0+7axs5pcNwezSQVi7+xhWbavBxsoGvPJJHfPY4fmZEdNMs1ZrqsvqigWdqSvqmtpDqxT9oB+UKW6cNAw3ThqGbftP4JFN1UwjuIpRqLUHgpAp8B9rFQ8hNU9SY4sP1Q3NOHSylaku60hEMiINvqCMldtqIEeI2fDYDF7MsHBEsNoeD1wocHoViaqVEY96KlFCym7FtFhWNmqwmBowpmY6VX/703UV+ONNF0V8Bp35mPbqA8+CVKdGKspLZ+bzWTavJNymGWMH4L83VOn2u0UCgQAeUTQJteqGM7htTTl8QRmBkLB58B9VqDzWhH9UNHSpqpkjUIpItmq7qs4h/djHWG2PB25o5nBswPKesaOeshtEZpdIwWZaA6uV4V1d2Wjbs3b3MVOaZyNBGbhlzUe4eHhexGcwv2wI/rR4YkSXVV1Qm0eERyT4+TXjdNlPWc/7qesm4IP7Z4UN1aogzM/yIjfdzYxDWPdRfdIFgkhgjv2QaES7hdZdORLVX7K9jKy2xwNfKXA4NkmGDSAaVllW73j+I8jorElgd2XT2OLDqm0HbV///ZpG/HnxRSjITrN8BqWFOSaX1fZAUBeDYedZWh3DOjbTI0YVbInGJQA/njUGN04apss/5QtKEAQSMa32uvI6/HjWmKh9akA2e7/V9njgKwUOJwacnvnHijL4mfXsbQG9S6bdlQ0rcCoat7/wEY42tlo+A/Xa2rlyUAa+/dv3dC6qrGdpjHFQ1U11Te2oOX6WGa2+oaIec1du77bB7NoLh2DNzRdj98++hbtmjUZ+llfnfvvaf14W9Rx23adPnmW721ptjwe+UuBwkky8hW1UW4IgEECiTPdNreHXzmzcTtI6I5LMzo6qva+SwTkmLx2ZKrWTWVlVG1t8IWP3QZ3NgAK4d/1eSLKMoKzEJhBC8NDcEowrzMWpsx2mNNbJxCMCV5cVorQwN6ITgdYzLp705oqx/Ax2WiTj+/z4GQfuRoELBQ4nicRbipPlTkmgDEpaz03j4BLN8G505fVLcngAjoQxHsBYWvPaC4qYvyOAyctGLVijqn98QcVtc+n6StPgqWSWpXjwH1UQSWzxBU4jEoCC4M61n0Q0+Nc1tWNqcYEuDiGW9OYbKuo1ab/ZlA3Jtd4ZI31SKCSiBCGHEw07wV1WsGwJXpeI26ePxKp3aroUm6GuKKobmnGmPYjqL5vx3PbD8IgiOoISghI1zfq18QD/rPoKD76qeAqpeexe/qiWeS2fRHUpwtVnYmUPiLQC6E6B4BYAQRBC3k7KTRvfZaQJgF37VGOLD/eur4wqpOubOxy7tz4nFJwoms7hxINVrIMd/3Qrw7Hq1x/vJEedIFXVN2PZxqrw4OMWO6OVAYSquymqHTUe4M61n8AvSfAzyoH6JYpxhdmoajB7xTy2ZR8mjTgHrX4Jze0By0I80cqMdgdelwCAYsGFQ7Ghsh4+TSya9l3GEt0dCaWmhYDOwkBsTp3lwWtx0ZWZWm+Gr5y6zgsfHMaGvV/i6vGDcdOlI0z7G1t8aG4PmEpx2vVPjxaxHc97UydILoGYgsQCEsUzW2tw5bhBKB6YHa7upo8HYEfnqlQ1nIWbAAHD2C7LFFc+/T48ooCAJNuqWJYq+IMyPC4Br1bUmyKute+yKxMALUpMhw2bD3HuGfYpoeDUi+pN8JVT15nwyD/R3KEMEHuONGHFm5+j8pHZ4f3aZyzJMtwiUVIexKjumVpcgNXfvwgAQWlhTpf6rKqWiDQb90syrnz6fTx1nRKNrMYDeETBtvunUSAo51U2BjQRzXaL2QsEEYPAEg2F4gKs9fXJ9IqQZKp7l07l4srP8mLFggn4aRSbAmj0+tB26VNCIRlJ03oSfOXUdV744HBYIKg0d0h44YPDuOnSEcxn7HUBqxZdwPRY0aJdwW2vOWUS3lqdNGBdnpKFErAWfXQNSBR3azyGivLS0R6IvEKIB0EgEChFpCaJRKlR4FSRm2iIhDAjr7VkeATcMnUE5k8oRPHAzrTWZgO+hDu/ya6FEI1Om88Z/M/Wg9h1pMl0THaac0N5nxIKyc5MmerwlZM1dlVqG/Z+abn9pktHMJ+xRxSRm+6JfF7N6kLrwqhNOyEKAtwiQbtfCY7Srj4irfYaW3x4Zqv9gLWARHVJ7QjR5wkVBQICGtUYGgm3SNDOWlaEIACuLB2Ef31+IikW5gy3iMWXDsfqdw+F35xLUO5da/xu88v48/bDWP3+IdNzVwdzNa/U6vcOYdU7NXGtxvOzvJg+pj82VdQzhUKam+c+iptUiEpNFbq6ckq0LaK7bB2xqNSuHj8Yexgf6dXjBwNgP2O/JKO5PcDM5AmwV3BGgjIQlOWwoVOSqKUXjJHqhmamV88PpwxH/2wvVrx5gPErRT1R19QOr0vQqX7S3AL8jIC6WAjKFG4RCFichgLYXPVVl64RCwFJwh/fPaTzuqIUWH51KR5+tUonl1TbgvG5q4n4Vm07CF+QhlVu97xciZLBObqVhV0OnGDHI1htj4c+GdHc3VGpqUK8+XyAxJe+dOL88dRrjrWC202XjkBumn6Wlpsmho3N6jPW5uXxBWXc/nw5Lv3Vv/DMvw6azq14nMSvI44eHcs+98zzBuKGS4bpisgDygy5tDAHALt8Z0CipqhoY/1pI+pur4tADNkJrARCdxCQzSmyJQoc/Oqs5UJF+9y1NZ2Najq/RDHnme26Pm23r3otos+ttsdD0lYKhJDZAH4HQATwLKX0V4b9wwD8FUC/0DH3U0pfS1b7+irxrJwSbYtw4vzxGtDjUalVPjI7ovfR1OICGDM2q0FoT711ACu31WDFgs72VdU3mzxbXAIgCgI8oZTTkqzMrllEW+2VFubAJUCn7lEH/vwsL26cNCycJRWArogOq3znPVeMwa/f1q8ujPWnjQgCwbrbJuHt/Sfwh3cPKdPwHsALu49a7gvIMjI9It47cFIXjMfCH5TDfVprL/JLEpbMGB1+5kbOdrDtORke54bypKwUCCEigFUArgRQAuB7hJASw2H/DWAdpfQCADcA+J9ktI0T+8op0aUvu3r+rtRrjleldtOlI/DyHZcy3VGj5RfyBTvbZzXwLp8/Dh/cr+TS+eD+WVh+dSnzXF5X5NWeqpJbPn8cvC6CDLcIr4vg1wvLkJ/lRc3xs3jpQ33g2bryuvDvjO8l0yNi7ODscHnOjGhLhBBBmWL9x3X40/uHbB2fChCocQpmvC6ChRcVYe7K7bjj+Y9seWe5BQHVDc26vuoLUjz1llIf2rg6bmzx4bPj7LKbHot2xUOyVgqXAKihlB4CAELISwCuBqDt/RRATujfuQAaktQ2Towk2ourq+fvigE9Ec4IdvILaYWese2ZXhHjhuTqgp0WTRoOUGD5pmq4RQFBWQ7PMAGgsvZ0eOXXmTfnFJ7bcUSJD5BlPDyvFOMKc8PHbaiox9L1e03ePWrbWPfREZRw61/3wOtyQaaAP9oSQcNLe+psH5sKiIbVlcq3xvbH/Veeh7krt0fMhGpEeZaEqSr0aVYS6juva2qHi4DpoTXEQQ/KZAmFIQC00486AJMMxzwC4E1CyF0AMgF8i3UiQsjtAG4HgGHDhjneUE50Eu3F1dXzd1WoxKJS05a7ZBWvMd4PlSmzCLy2faaymTJltn3R5OGYPW6Qrp2samsv7j6mG8zUWexjm/dhx30zdRG4fsao55c678t4H5KsxNr6owSy9QbS3S7cNHkYfv+ufnWz/YtGNDSbJyJWaOMaar9usyw7apzIKIFs7HPmprtju5kIpJL30fcA/IVS+hQhZAqA5wkh4yjVh/NRSlcDWA0AEydO7BmKyF5Ior24unJ+J4SWnRQEaiI3mVIEJIo0t7KEZ9kv5pcNQcngHFTUnsa5+RnYeehrrNxWE561a9sXS9u1s8imVr/JFqO1DRhR1Re56R40t/stBzWZAjtqTmF+2ZDwfcx5Zjtiq1bcs0gTAUP4CQKyjCmjCrBm51HdQK6o1IitbLNeF8HyeaWYMVZx7536+FbLY1nJDaeOysd2RqbUUy09L81FPYChmr+LQtu03AJgNgBQSncSQtIAFAA4kZQWcmLGzsDZXedPtNBqbPGZUjarqgOWUZxl+P7g/pnM9mkT1KnRy8Zrs4La1GIudukIKvWTPaIYjoVgEZAolq6vRL8MN0oLc9Hql+AVBeaqoqfjEoDlV4/DuMJcVDU0h+s/q++stDDHFNAWkGWUFubgiWvH4yd/r4gYRuELUpSF7HeVtaeZgjjDI0KmFA/NLdFVqwOAksIcplDoiTWa9wAYTQgZAUUY3ADgRsMxxwDMAvAXQsh5ANIAnExS+zi9kEQKLStff8C87Lfyptpx38xwtlAjrAjm+WVDDEFtnWmuwwNLDIFdavZTNVW1SBRDqkgI2gz+ob4gxR0vfBwerFizYoEoxtjuLofcFQghmHTuOSgemI0JQ/thdukgk+C2WslNLS6ASxQgRRCWXpGEVxlFeenoCJq9zP5w04Wo/brdJJDmlw3B163sFYHV9nhIivcRpTQIYAmANwB8BsXLqJoQ8ighZH7osLsB3EYIqQTwNwA/pLSH+Klx+iDWM3Ljst+ON5XWT93Ke6rm+FmDp0r0ugeMssVhWH74sizjB5cOZ3rZtPmV6m6Pbd6Hb5cMNO13iwSTR+RHblCKE5AovvO7zgpxRs+8xhYfhudnYvOSaaZa0YqXWeQhlQhE1zdYQ1xhbjoe27KP6T3nsnihVtvjIWk2hVDMwWuGbQ9r/r0PwNRktYfD6QosX3+A7RIaLarZuCq485vFTO+pCgt1QyRinbUHZOB/PziCh+eV4NFN+wAK+AyZXQkItjDSe/iCFDsOsSuD9SQkGVi6ProK8KGrOr3q1bxQxvfsFgkEAl01Oa0dKN3twllN/u10t4v5ntVJRIZFOgur7fGQSoZmDqfHkJ/lxa8XlmHp+kqIREBQlnDXzDHMoCPV8H23xgbhC8q47a97AEJMKqCV22pgnMf7JUlJJyHFpsePpzqZAKCxxY+gJDN/255KoccJQhRIVBXgg69WIcsrIhjyJJpfNoSpWlJtW6qHmprexMpLrmxoP0vvuSNftzLba7U9HrhQ4HDiJBZjNiuqWXFNNY+6HlHA7dNHYuW2gyAg8AdlyBR44P+qdKm3O4ISZKq4hVoRj36/LSDjt28diGE90nPxiARBiZru1egGzIp9ARBO+bF0/V70y/CYym6qfcLKRsQSIsUDs5nbAWA/o2gRABw+2ebQE+mjQoEXleE4hV1jthrVrBp1IxGQZZyT6QnlA1IGIVmT8M7rAm6aPAzP7TiSsIL1fUEgiIISyX34VCt++/aBsAB1iwQrFkRXAWrpXPkBKxZM0LkkR0rbYjWxMG7fXnMKUx/fahkcl5vWO+MUkgIvKsPpDuxENav89IoxeHTzPssB3yUI+POOI73SJTSZSDLFT/7+CdLcLogCwTXjB2PG2AGYMqoAgD4qXBv7IhLCDDhTgxJ/uq7CFIkcKcJeexzQ6X6q7tMKFSv6ZXChEBe8qEzfIpVWhOqgcs/LlZCpdTK7TI+II6faIubOCUhKSUh/7w8ito1bJBAFgmnFBXj7M3Nok5VtJSADgZCh95VPGrBp71e45oJCbKxsMCWo087eqxqa8eimfcz3FJSB6oYzmD6mP4DoEfbRJqpWqistrT7nOkOfSp2d6ERunNQh0am9jdhJfVx+5Gv4JWopEAClRsL6j2ot93tEAcvmlUY8R19k+DnpWP39i3D/7LFwMQL47GrZ/JKMdeV1lgnqVBfVRZOG40+LJ8JtOYJ2XpCVol4NTDO6GbOSN9pZZZ7pCNi7QRv0qZUCL8fZN0j2ijDSTE9drQSCEtbssk45ofKtsQPx7sGT8EtsD5/X/nMaigdmIzvNhaVR0jP3JWpOtuG2v+6BDJIQgalNUAcoE8z9X56BlUbnTHtQV0RJt8qo74yU9kkyiCFWwRj8qFVdWamQcng5zvjIz/Ji4UVFuo9z4cSiblctcJwlmWVGIwkgrcdJu83KZG/s+8oyLE4kBA3NHcjLVLxcHrhyLH71+n50cMEAAFAcgSILBLcIfOu8gVhwYRHuWPtxTIZ6tyBg7e5j+J93auASiKnYkJZ71++FDKqbIKh97/rVOyNW1WNNVFWhct0fPsChU2ZPo36ZzvXrPiUUGlt8WPeRPl3vuvI6/HjWGC4YehHJXBFaCaDqhjNRS2qyCMrAkFwv6pvNaiiJUtzylw9BoaRj8IgCFwgxEpCA16uO483q44i1uJ1fkrBqW42t1ZmaJoSV/trYX9LcAmSZwusyB7gZYQkEAPA5GDvSp4QCL1TfN0h0am+tAdtKAAE05uhjFZZA6Dy3+i+qq5PMiQ2LEBFL3CLBkhmjsfq9QzGp7ESiD4Kzsg+89p+XWaZeV1m721r9eOgku/hOPPQpocBtCj2bWLyJEpUllWU/YAmg0sJc2y6onNRAJEqZULcooM3gckpAMaIgI+aI8la/hKqG5nDiQ6Nra0BS0mUUD8y2PIdaJGnVthrLY5zsaaQn55ybOHEiLS8vj+k3GyvqTR8wj1NIfVIhvqSxxWcKIEpzC9hx30wAMAmgjaFKZqKgfPyJCjTjOEOaW8DmJdNQUXsayzZWm2IRlPrYne/eblZYtY9oJyZrdx3F8s37lIhqmTKdE7Sp0QWYM9dqGVGQjm33zIzldi2VZ31qpQDoi52UDe0XUUJzUoN4vIkSEaMQSf3IqnFNQ/8vywCVKURCTLn4OcmBQIlgVj2TXAJw46RhWFdeZ0oxkZfpwX9vqDKdw7hKUJ2cyopy8dlXZ0IR6xIIIToVEyuV+mNb9sEflMOxJiznBG1q9Gh4hJ5XozllSIUZJyc2YrUFJeodx6J+VAWZL6hVXlO4hMS4TPZFYkn253ERbLnrMjQ0t0MtXJSf5cWPZ40xTR60Kp5oM3QAqKhrxi+/Ow5D8jJQmJuGuSu36/ZbpVJ3yjkBAPplOhfR3KeC16zy1EcKOOJ0P/EMxol4x6wgJCsDdnXDGaYhsyera1MJr0vA//77JfjZlWPBKP2gQwzlIyoemI3pYwZg+pj+OgHAWuXNLxuCHffNxB++fxHcNmoVPLKxGreuKcfuI19H7SPRnBPiIeigarJPrRS491HPJD/Li4euKsHyTdVwiwIkSi0H40S/YysDtlEPfK9FYBk3K3QNtVSlWhozN92NP950EW57/iNmOVGBAH+/fTLcLlEXTGaH/CwvSgtzbAlyNePtg/+ows+vGcfMlKo9r5VzgrESmygQuASEVVNBmTLvM8fLcx/FBfc+6plsqKjHY1v2Kfl+JIpl80os1UHJeMfGzKj6EplSKLspH/2dRiRKqcrSwlxsrzmFS3/1L1AK+CVqqfKgFFj05926IjeRVIlGWxSrEE40lm+qxuxxgyxLrQLsyUVji88kgAgottw1Ha1+CZkeEVf85j3m+b4602G7fdHoU+qjWJb/nNRAqw5q8UnwB5VykFbqoGS/Y6O6yhekXCAkCIkClbXNAIC711XAF6Twh561lQZeqUFNbakSWfmyWJMMr0vAL787Dku/PYapunKL9vKpGVVXqgDSXcutFOaZEHKKKR2UxTxXYT/nJj19aqUAJM5/nZMY4lEH2X3HNcfPdtkLzU4Gy0gQxBRD1edZue0gRhRk2PLIYWHVd6w83HbcNxNPXDs+7FoclGQsmVGMb5cOQn6WF/0yPXjwH3pPJYnSuFamLAHU6pNQVd8Z51A8KAdVX5kD1bLTufqoS9gtjMLpfuJVB0V7xw+/+qkuB9biKcPw6NXnx9y+TI9oqmEsEDD1viy4QIgNlyjgVIvfcr9aE9klCOyaB8Egs+9EmnyorsWBIEVAplj1Tg1WvVODJ64dj0WThgMUtuxd0cjP8uKhuSUmIfPo5n0Yek4GSgtzcKadnQ3Vans89EmhwEkN7MQSJCJlRc3xs6aMpWt2HsPiyeciL9NjexWp2hLULJdpbkHRcfN8RAmjNaRCdIvEpKbzughWLJiAqcUFeOrNz/Hih+YU5BIFdtScMtkVrCYfmR5R41qsoAYvqrEFiyYPx+xxgxzRPowrzA3XclbxBWXc8fxHkEExI1SjwciEoty4r2mkTwqFVCq+0leJJZbAaZVfRe1p5vbndhzG/31Sb6tNrGpYskyxfH4pHviHOfCJ4xxPvfU5Hplfikc37YMoEEgyxZIZxbhx0jDkZ3lRc/ws1n/Erp8hyezAR6vJR6tfslQPamuxROqbapoKgKK0MDdi/y3KS2cGOKqxEv+sPs78nZMrzj4nFHjwWvcTT4RyvCo/1gSgzMIr5JWP6+EL6ttUMjiHmaiMpW7wuhS3QU5i8QUpquubseWuaaZ3s6GiHktfrgwboFlo7Qra/mHlEWSVwyogy6iqb8b1q3dajicbKupx97qKsA3ELRI8dd0EyzFH634tCgTthvoJhCgeVUb2f3U20iOLiT7lfcSD11KDuqZ2UIPSncrU8Qp4VtXXigdmY/GUYbpj55w/EB7R/DnMefp9ZvU2lrrBF5QwrjDXVrATp2u8+GEtrnrmfVQ1NKOuqR2NLb7w9x1JIACdNilW/zB6BGm92byh95rmFsLV0x7bss9yPGls8eHe9ZU6o3hAoli6PrIHlOp+LckUxq5kZavK9vIiO3HBg9dSA8U4q+/dPoki0yM6do1oq5FHrz4fiyefG/Y+ysv0YOrjW3XnUFVDfilo+r1W3UBlCp9EIQgENz33Ib53yVD8fY9St4NXRkscvqASLJbpESFRiju/WcxU9SjlOSnS3a7wbB6A7dWqdgWh6vuL8tKjjid1Te0QiQDAHJAWzQMqfD6RwEUQjrMYUZCJz740rwoaW52b2PYpocCD11KDVr+ENLdgyjbK8haJFzsTgOKB2TpXVK1OWS2TqBVeokCwbf8JzBg7APlZ3nByxTnPbAdAw/ezrrwOW+6ahobmdtz61z1w8LY4DNR+8/TWgyAG7bpHJHjtPy8zORBU1p62zD+Um+42qQut1JeRxhPFPsCIapfZLqvMAjwuEasWXYDcdA+K8tLx5Bv7mUJhYE4a89nEQ59SH6mzO6+LIMMtwusifTZ4zU6h+URhJYSdFM7xTADUfDcv3DoJr901DcRQAL7VJ2HZxmqdKqnVL8FrUDu5Q+6Q08cMwA8uHeHQHXGiEZAoZKqkhVADF5+8Tsl5ZFQLsfpHeyCI29aUM9WFLKIFSuZnebFiwQRdgJtbJFixgD3mWPXZ0sLccNtzLeIRrLbHQ59aKQCqlT6UCJ32Td2vMS3Dkhmjw54bySCam2k83mEsD494XFm1M8JwMRSBoDVUj1edlS5dX4mpxQURXRmvXrkdlXXNsTwajk2MK02VoEzhdQlYtejCcCZUFsb+oaYn8QXlsMovmvMDENkzrrHFh+H5mfjnj6ejobkDdryP7vxmMVZuq4FHZPfZI43scpxW2+OhTwmFznTGnZ3JzovvTbB07U+9dQArt9VgxYLkeWJZfUzxeIdF8vCIlLxOqx9mvX/19yyfd1+Q4sXdx3DXrNE64RGQKOZPKMTs376HII9M6zIeESb1m9clYPX3J6L26zYs31RtMix7RAG56e6okwxt/2hu9+POtZ8gIHXmOGLZG1nnYqmWYu3H2uMBitunj2RO1AZYjFNW2+OhTwkFK++WvmRotkrL4AvKSReQxo+JJbDuebkSJYNzLNNQNLb4sPTlCqaHh9YorKJ+fIBiSPaKBEQgET/aVz5mqxFWbqvBjZOGYX7ZEJztCIaiWgnWldfF9Bw4ZpZ+ezRcooin3joAr0vpn9p3NT0UxDVpxDmY88x2XcCgVk0YbXBW+wfL9dSobrQ70Mfqcs06ftU7St8y0hpgJ+az2h4PSbMpEEJmE0I+J4TUEELutzhmISFkHyGkmhDyotNtyPSIpiVnR0B21Osl1bEqHA7og3G6A8VbQ6/S80sUc57ZbqnfXbv7GNOQq3p4aNF9fKF+4JOozpXQaGupa2pnuqoCyoxUdYd8bMs++CWKVj/3NnKCczK9eOqtA/Br1DmUELxw8yUYnp8Zfj95mR7cNaMYXpeADI8Ir6tTr89yQb9n/V7UHDcbaqPZB2JxZ1cnXloifVuxHH/4lDnvUaTt8ZCUlQIhRASwCsAVAOoA7CGEbKSU7tMcMxrAAwCmUkqbCCEDnG6HYhQkOo8Sr0gc9XpJddTOv3R9pS50H+iaJ5YTUeJV9c3Md+G3WMU0tvgsi5mzPDwiJa9zCwLW7j6G/3mnRjcTnFpcYClEO4LWromc+CEAHt5YbUpjIRDgxmd3w+tS3DMXXlSEdR91uv66BUDQOAdUN5yBYChF7A/KmPP0+3iSEUAWyT4Qizt7rE4OrOP9ksQ8/lgjW7BYbY+HZK0ULgFQQyk9RCn1A3gJwNWGY24DsIpS2gQAlNITTjeiKC/d5FFCBNLnXFLnlw3BB/fPwt1XjIHX1fUU0xsq6nHpr7bie3/ahUt/Fd1rg4U627aCNXOymsW7BLaHR6RVkl+SsGpbjWkmCCDssWZEzX0f6byc2BAAuBh5jQBlVe+XOtNgr9l1TLfqC8iKrefeV/Zi7a6juG1NObOUpl+ilrN8q0pssQz0saZvV4/XeinJoRxNRpLhuZcsoTAEgNZSVxfapmUMgDGEkB2EkF2EkNmsExFCbieElBNCyk+ePBlTI3g9hU7ys7y4a9ZofHC/4oK5476ZcRmZG1t8uOflSviCMtr8EnxBGXe/XBmzqytrCa2F9QGyPlSXAPz2eiUpmrGddU3teGhuCbwuISxMvCJBmlvAkhmjTQJGFUTzy4bgT4snIsOgZkx3u8IzRbVf9SVVZCIQRQKBmAWwWyThiOKo5yAEyzfvixg4GKuqNNaxY37ZELxw8yX496nn4oWbL4n6bU0tLoCo6f8BC8F11fhC5u+ttsdDKhmaXQBGA/gmgCIA7xFCzqeUntYeRCldDWA1AEycODFm/w5eT0FPV9OIVzc0m2Z1AYmiuqEZ08dE1gBqVU5Ws201WpX1ARrdCtsDQRBC8MD/VekMgVoDYUdQgixTeF3K4L1k5uiwQW/VO3pVlFYQlRbmQpKt1W1qv6puaMa//+8eXnYzTtwCQRvD1fQPiy7EnX/7xFY904CkZFG1TrAN+CX7qlK1n04tLohYZlOLNjX701troqZmV1e9WkHGUk9lWaSzsNoeD8kSCvUAhmr+Lgpt01IHYDelNADgMCHkABQhscfpxvB6Ck5iNXuLPKtjeXIY4woemluCoXnpAAhKC3OY5+kcjM/gtjXl8AXlsFuhmtDO6NkBdGadXLntYFgoRPIR315zCpJGaLlFgju/WaxrixJc5IFLIJC4VIiLYCjfj/bxuQSgIDvNNAFgLQRcAnDNBUOieoAtmVFsawyIx0U6Ump2Ky86u+qp2iZ2PILV9nhIllDYA2A0IWQEFGFwA4AbDce8CuB7AP6XEFIARZ10KEnt67N01UBcWpgDlwDdB+oSYDmIq9e0qnKlnYltrzmF25//KOoHqUZ6GmdaIiGoYKQz0OILUjz0ahW2fn7C0kdcba/R7fWP730RLraitouV16kvoqa0ju03wI8uH4XfbdWv2IKy8ly1q/w3qxuw6p3DpnP8v+kj8OcdRyNexyOC6e5pJJ5svoB1avaK2tOWQsFusKWVetJJtWVShAKlNEgIWQLgDQAigOcopdWEkEcBlFNKN4b2fZsQsg9KBqmllNLGRLSH11NQcCKNeH6WF79eWIal6yshEgESlbFiwYSIz7WuqT2UpKwTdamsGvli/SCZpQz9Er5u9Uc1Ar9W9RUAhK+jlnzMSXejtDDX0ruoxdcZ3dwvQzmWldepLxKrQFB+A/z+3S9MKwWth2B+lhf/rPoKq98/wjzHeYP7wSUcY+5ToSDMQjtG4k2gaZWa3Wq7ih3V9uFT7BVBbloPTHNBKX0NwGuGbQ9r/k0B/DT0v4TRF+spsIRgvLMgFrHaaarqm8MDqopxqVzdcMakglHTa7POn5/lxU+vGINfvLZft/3Jtw7gnm+Pwa/fOgC3IKAtEIQUZbz2BSmW/K0CgKImemR+KfyStduyL0hxxwsfQ6YUD11VEvnkHLgEAlmmzLUbK+21tubx2l1H8eCr7CJGi6cMw5RR+dFTZ4eMuNH6ejz5sxpbfGj1S1g4sUinwlo8ZZitOuCRVNuNLT78n0Ug5acNZ6Ke2y6pZGhOOE4OhD0FKyHodBpxu3YaK9fTh+aWmAqlBGT76bU3VNTjyTc+N233B2U8+cbnWDavFOOG5GL3oUb84vX9jDOwCUgUj27ah5unjsDv37XWZraFZrKPbq7G0Lx01Jx0Tsfb+6C4/8qxtt8DCXkjNbb4sHxTtWm/1yXgqesmYO4ExQNn2TxznWMjdvq6tuCNnfrLxm/tZ1eOxTmZHpQN7WdLIERDMUYDHYz5yZSR53T5/Cp9KktqrJGGPZ1IUZjdlUac9Q4yPSLGFebq2sya7Vml145WXMUvUTy2ZR8yPSJ+/fYB0/6FE4uQ5haQ4WYLHIEAA3O8YIQqmPAFKRcIUQjKiofQ3PMH2jo+zSWirqld6TuMuBQKYMqo/PDfiyYNx8+vGQePSJDhYQ9xdryPtAVvAjLFQ3NLIpZnvXd9pe5b+/XbBzBj7ABHBAKgrFysvGyNFdq6Qp8SCn2tnkIkIZiomI1oKblZ70CrHogUr0Ap0NweMJ07WowDoNy3anTWkukRsWjScOy4byaeWHA+c+BvD8j45euf8wR3DEQCW8LSyIo3D2Dzp+x6w8ZwBPUbtapfvGxeianfLpo8HDsfmIVH549jri6jeR9pJ1QtPgn+oIzHNu+z7Ndrdx8zZQiINOGMJ3V9fpYXFw7rx9xXzdVH8dGZ4mFv2DuiNwevRROCTsds6FNyy7h56rmYMqpAl8LYnLJY1rl2WsUruEUCSZZx59qP4ZdkLJlRjCkjz8GRxjacm58R1ZjslySUDe3HFEiZHhFrdx/D794+YOkGzyuosZGoszNLl0CwfH4pHtuyj+mFE85ISwgCkoxl80qxaNJw3Tm0NrSyof3QYYhqFkl076NY1KtW6VasUlWo34lLIPBLFMvmlZjugcWGinp8YpGKfVwEb79YsS0UCCHXUUpfZmxfQCld71iLEgxV/5+S8F+9FTtubl2J2dB+fIC5vOHv3z2E3797CC4B+PXCsvDSW61Y9tyOI1j/0TGsfu+QzrXTmOf+5mkj8Nz2w/BJwFmfEoPw1Ft6NdC4wdk4cKLFUoUkU2Dfl2dMz2PhxCLMefr9qMZJIwKxrpfb13BSXAZliq9b/ZZBYtEmMtqJiRrMaHy1ghB9aROLVoEVeAYAS2aMZgoQY8nNB/9RhdaOIG6/fJRle9TfsdJ/AED/7O5Jnf1nACahACW6uEcIhc56ChRq3dTebmhOVAS30ahmVR8XUHTIalGa/CwvNlTU4971nXUt/IZgs34ZHjy5YELIJTQHdU3teGHnMfiC1umBqxglCrUEJIp71u/Fa3dNw+Yl01BRexrn5mfghj/tstTTRoILhMTxzNYDmDA011SQRjsJmcBw72Q5krAmfqqNIpqR2W6RJpYA8boE5mqE5Y4NAL94fT8y01yWK4ZoSRd3Hf4as0oGWd5PLEQVCoSQkaF/CqHgM+0djQTQ4UhLkkBdUzuobN/NsbfgdAQ36+Nbue0gIkUxi6RTv2osdKQSkGR8+zfvhbu96g46NC89okuoXfxBGd/57XsgaiF0SY5LIHASi19Cp4vv3BKMK8xFVX2zSaVkNPrazVZr145onFABQGXtaWb9ZrsCJNMjWqoil2/ah9mlg2wLHi0jCzKj3o9d7KwUaqCIWwLgC8O+rwAsd6w1CYYVbRrJzZHDhvXxeUQRt08fiZXbapidXqJy1DTTxviBgETx4D+qkOkRHZuZSxQABYJy30mX3hNRXXwf/EcVMtxCOB9SJFdyO9lqvS6Ch+aWhCcodsuzRotvsrMiV89hZYNxiyRiHI5qD2V9Xxef65xLalShQCkVAIAQ8i6l9HLHrtwNsKJNrdwcOdZY6VtvnDQMN04ahhd3H8Pv/nUgPAt3CdBFOceaZjre9+O2SMHMSU1EwTwxAMBMkMcy+lolSExziWHnhHMyPXhsc+QVhxG78U3RAs+MtgQjrBogWuaXDUGrL4gHGDEY+78645jrayw2hX9jbSSEjKKUGlcQKUkycpGnColM5aEG9TyyqQouIkCG3ovrxknDMGFoP5xpD4TtAlrvo/kTChNWstItKoNAQJbx02+NiSlQjeMcHpFAFIgt/3mXQPDHmy6E2yWGkxpGw0oFxFL5aP899fGtMQevOhHoyV5dE1AAHpdg2xOyziLx3WdfnsHcCc5kZohFKHxKCLmFUvq6uoEQ8iMAjwEosP5Z6hCL7q8nk+hUHhsq6vHIpmoEJCAAWedXzrq29vmu3XU0oQLh9f+8DK1+xRVw237H6zT1eox5h+Ih2+vCqkUXAqChQb7zhGpfUa/hFgmeum5C2Ei6YkHI5VQgaPWZV4iZXjHqAGqcsav/rmQkR7QzuDsR38Q6hyAQbF4yLdxf7YxDVq/GyfVwLELhFgDPEkI2APg1gGcAFAKY6WB7Ek5vr6dgtdQtGZwTU+eLdP571+td4ySqeBex0lRrZ2KNLT48sjFy+oGuEJAoGprbMX3MANQcP4vjZ3qMD0TK4IS2rT0QxJn2AKaMyseKBRPCkzBfMAiZEgQ1BiKBQFcQSft9Go3LD11VgnFDcm33YeNqOd7B3YnJpNU5YlX5DM3LiGl7PNgWCpTS1wkh5wN4B8D/A/A3AFeHymv2KJpa/Th4/CwyPWKvEwqsZSqVKeY8sx1esesrh7qmdogMlzqRCMw01dqZmPJbIaGlK8+0B/EfL5TjtSp2tCxHwS0SfPeCQmys/BKEAu0OumEFZWDJ3z4JrwJ23DcTa3cfwzNbD+oEAgCAgmkbUMtizh43KK4JnNVqOd7B3YnJpBPnyLbIhmq1PR5iCV7LAvAkgFwAvwFwM4AfIlQFraegrYgEAAsnDsETC8ps/bYnpNxmzYZ8IZcbf+jDX7p+L/pleHS6fkC5v+qGZqhFbbQZVbUV0lhpkSUqo2xoP5PrqHYmVpSXDprggMGfrKvgxuUozC4ZiHu+8w3kZXowd3whPq1rxoo3zTmhukpAoli6fi+23DUN//NODfO9+CSK2q9bLb+peNypIxmGuzIwO+Ha3RMKfMWiPtoLYAeA8ZTSZkLICwCeJ4TMp5TOTUzznIVVEWldeT2K+2ebogmNAqCnpNw2LlN9QQmCQHReD76gjDue/yhsIFZLVt7zcmX4w1WjkClgCslfsWA8frKuUiccZp03ALsPf63b5hKgm4nlZ3mxYsEE/OTvFWE1hdORwVwgROfNfcfx5mfHIQqdnjnxxvdH+50oRC90dO/6T3V9satEMwxH8xJK5YlfTjp7RWC1PR5iEQr3U0rXqX9QSisIIRcD+IVjrUkw22tOMbc/8cZ+XHtRkS51s1YAPHRVCR7bsq/LKbeT1eG0s6FMj4i5K7ebjlHLUar2BqOdICgD97xcAUL04fsP/qMKP7tyLFwG98HXPj2O1wwJzijV64sBZQBxiQI8RElpwMtWJh8ZACggSzRcutQOBIoQ174yj0sApTKsvIYlmTJzTmnR9kUnsgvEazvoCRO/Y6darLeP6e/INWznsqKUriOEuAkhlxFCrg9tFgE85EhLkkCBRWcTNdkMWemml2+qtqwUZpcNFfWY+vhW3PTsbkx9fCvW7j5qK0tiPNkUtcKneGB2OBtqBiNIj8oUFbWnmXYCAgGsNDFPvLEf0WowA8rgsfOLzuJ5nWlGZLQHZAQkdqEVIzZS1XASjEcU8Mz3LkCGRz+PdAkE1100DF4XgdeY3hTA9RcXhfsgY7cOUSDYtv9ETH2dRTwZgCOlmbc6Ptbv0gneOcie2Fptj4dYbArnA9gIwAegCMDfAVwO4AcAro/w05Rhyqh8iIQw0u+aUzfrlp6iYEqYFotLGkvHqUbqqoU7WDOSeGYuVr+ZWlyAZ/51EH/Zqa9f65Mozs3PYNoJKGRIsnne4Bbt+Z8Dio5fptSysI8deJ6h7sXrErBiwXhMGZXPLHm6obIeAMHN087Fn7cf0q0a/r6nFt86byCmFhdg9fcvwi1rPrK8TqtPwrKN1fjvDVVdnqWrSRcrak/bKnITSyxCd64oLh7WD29/Zna1vtgipXY8xJL19vcAHqaUjgUQCG17F8A0x1qTYPKzvPjN9RN0s35jtK1Vvv9l80rirj1gle+/1S9ZzkhYM5el6ysjzkyizXZe/NBcu9brEuB2iVixYDzcov65PHldGZbNN5eXDMQQYByQKO5+uRI1x8+iuT0AvyFklbVC4XQfLkHpE9leF7wugruvGIMP7p+J+WVDdDNwbWqYFp8EX1DGn7cfgUj0/VwtVTr18a1o9Uu4rDjfeElkejvPxfom4pmVb6iox9yV27F80z7MXbkdGyvYZSxV7KqcYl1ROI3XIiWP1fZ4iMWmUArghdC/KQBQSlsJIT0qHFidQWyvOYmCrDRMGZVvK7nV/LIhmF0an3tctJwsrBkJa+biC1K8uPsY7po1OrxNqyqKNNsB1BWPfkSnofZNGNrP+rlQhEsSBuVOLybTfYhKjQrjzD4gUXznt+/B7TKnF5YZy4Den9Q8ObgJEIjhQboE4MezxqBkcDYq65oxfXQBJo7QD+LqqnPb/hN4ZFO1rta2UeCrqHmM7n1lL3bcNxO/fetzrP2wFh5RAAVw9YRCbKho0KUzUfvt9ppTMc/K4ym7aycWobHFh237T1iqkqONCU7YFL0uC6FgsT0eYhEKRwBcBKBc3UAIuQRKwrweg52ln5XbWrzuZNoOx4rUZM1Iiiwyg67cVoMbJw1jekQ9dFVJxNlOpKpVkZ7LosnDw/7itV+3YcnfPjGd5+LhefjD9y/Czi8amfslCkgMlRNrzOICoet4XQJuuHgo/mpQF0aCEIKntx4MOxw8vbUGi6cMw6NXn286dkBOWsyeXm5BwL3rK/Gv/ScBAB2hCcIrH9fD+NYDsoxMjxhXTfV401JEcldVvw+REFMurmQasbPT2EO21fZ4iEV99BCALYSQ5QC8hJAHoNRR+G/HWpNgYln6qcEzTnkJzS8bgh33zcSLt07Gz68ZF1UVlZ/lxZIZo03n8YhK52bdy6Obq/HvU8+F10VM5zYu/T0iwc+vGYdFk4ZbPpea42fDy3b1eVhRUXcaADB2UDYI1wh1Kx6XgIfnleDFD+0LBEBZ0RkH+jU7j6H8cGO4H6gOE3eu/RiSLMMtEnhd5mEk3S3AY7As+yU5LBCMLJkx2vRNtPqluGqqdyUtBeu7134fWoGQ6RUTYsSOjNXH5dxHF0tE82ZCyHcA3A4lqnkYgGsopR871poE40Riq64Qa6TmjZOGmVJRq53bSr301w+OAiC4ffrI8IpCxWomFEsUtJU/tEsgWLv7GP7nnRpmFSpOcvCIBK/dpeTTERxSxF2/ehcyPC74JQkyVYSHtoiNl6G5oACWzS/VZST9twuG4MUPa03H+oIyrhw3CDdOGqbrm40tvm5LS6GF9X1kekQsn1eKGWMH2Ehi173jTqzE4n3kAXAJFJH0NYBMAP9FCAGldHGC2ucoTiS2cgo7qigl2Mu6c7PsFKr+dtU7NczKT6zr2omCVpftpYU5zKRpMgVWbTtoKl7OSR5eF8GKBRNQPDAbjS0+pkdZPEi0swwqC22NEm3COqMdrqnVzxQKblFRyRQPNCey6860FCpWzid2BILV7+Mdd3LS2UO21fZ4iEV99FcA/wXgDBQ7whea//UI4vFf7m5UtdMLt07CjvtmhvWQ2nvJcJunarHEURifi8joFdqZzW+uL9P5nLtFgiUziuER7Ru73E5We+dAAHDztBHhYMH8LC9umz4y8o80RIshsIM6ezb2U1UdUzwwGwsnmvXookAsB0ir/m8Hp1TAXR03nBx3SgtzdV6CgPL9lRbmxnwuKwhlGB+ZBxLSBGAEpfS0Y1fvIhMnTqTl5eXRDzSQ6qHssaDkKzpjykOf5haw476ZMUdcVzc049a/7jFFqHpdAj64v/N86nUBGu6Qaq56O1x7geJxwhcWZrqSvtobsieMK8xFpkfE7N+9b05CZ/oNwcNzS/Ho5uqoKz2RKBHpLtHsMGHV54zf2+p3v8CKNz+HWxQgR4jTSTW6Om44Ne5srKjH0vWK04okU6xYENfzs5wGxCIUKgF8m1KaMukn4xUKvZGNFfVMN9pYqaw9jZue3W1SF9x9xRjcNWt0uGNnekRTKu6NFfW4e10F7MgFtwBbx/U1fjhlOF75uD6iusYOmR4RQVnGrPMG4M3q4+EqeASAS9RXI1NtT+pgE8kedPcVY8K6/7W7jmDdR53+/yxPJa3XTkCSsWxeKRZNHt5tE7PeMiF04D4cEQp3A7gOwO8A6AQDpXRrPK3qKvEKhfLDjXjv4CmmH3ZPxokO39jiM834PSLBsz+4GLVft+GxLftAZQqfROF1CaCU4pZpI3DrZYqqYucXjbjrb59wt9I48boESLIM47jsFQBKCK4JpbtmuUZGOufNU8/FlFGKTQiAZT9pbPHhxd3H8MzWA6bVokckePHWSXC7xHBOLWNpW+1KgdWXAODn31W83pJNT8htlEQcEQqHLXZRSql95aWDxCMUbnp2F7bXdObjuaw4H8/fOtnppvVotKuOjqAESinS3KIuUImFSJQBiFVTV4sAxJjoom/hFgkEAnjEUFnRK8bgnAxPOF2DKvyrGprx6KboKh8gdnWiKhxWbqsBIUBHQIZLUBIlekNlJI3Zd7O9Lrxw66Sw63Jl7Wnc+KddJuHlcQnYeX9sqs2uwhJQ8ahYexGWQiEWl9QRzrSl+yg/3KgTCADwfk0jyg83xrRi6A1L0MYWH3Z+0YhTLR2YVtxflxtG9dyobmhWbBUSogoEQNGDRxIIqoMkFwiRSXOJWLXoAuSme0zVx9TZbdi1uXRQePBm6flVYnWBzM/y4q5Zo3HluEGY8/T7ABBevYS9jaLkAyvKS0eAEeXsFknS3TF7mltod+KcH1MP4D2LTILvHTxlWyik2hI0HgG1oaIe//VShUbF85lJH5yf5UVuusfR8GKuUrKHki6CINMjRk3Zrg7eqp6/qr6ZaTD2SzKa2wPhQES7tPoleF0i/BFSbHtFAiIQk0dNfpYXy+aV4sFX9SVYgxKNqy0q8fT5VHJHT3WS5hhICJlNCPmcEFJDCLk/wnHXEkIoIWSi022YProgpu1GrJLUvXeg6+l+48GYjjta0i9AuYelL1eYBug1O4+h5vhZ3bZMj6jzQVdx8YjlmGA9rkguoAFJxn+s/Sg8Q9cSzdV49rhB+OD+Wbj7ijHhxHZKTioZd6792HY/AZS+Uvt1a9RAREqUAvSsydGiycPx8++Og8clINMrxt0WlXj6PNAz3dEjkcjU3UlZKRBCRACrAFwBoA7AHkLIRkrpPsNx2QB+DGB3ItoxcUQ+xgzMxIHjreFt3xiYaXuVYBVFfMcLHyfdtS6epF+Acg/EQqtfUXtap0Zq9UtIcwsmQ2EMSVI5YK+Qvj9lOP7yATsNhUw16jobKdtZq1d19aBVAaoeTXb6ibESH6CofVj5jryiENHovWjScMwuHRR3W1Ti7fMqTga0dSeJ1lYka6VwCYAaSukhSqkfwEsArmYc9xiAxwF0JKIRjS0+HD7Vptt26FSbbWlrle20LZTud+n65KXQZaXjtpsXxqq0TZkht1GmR2TWMrDpm8AJIUAxrmpzTs0cO8DWbz0igcfVObt9aG5JOPcVEDmvjqoCNAYVRusnjS0+UyU+BYo/L74IHkOuI62gsprBxtsWLfH2eWM7nMxplmySkbo7WTaFIQC08e11ACZpDyCEXAhgKKV0CyFkqdWJCCG3Q8m/hGHDzGkcIlHd0Gzq6AGJorqhGdPHRP9ItWH3Aki4jKCKLyibUlsnCjs6UpbuNT/LiyevK8OPX9KrkBZPGaZbJaizEcIlQJeRAbx4yyVwu0RdXh+rmbcWCuBvod9W1Tfrcgk9ce14DM/PjGhAjUeXXtfUzqxz4RZEFGSn4UmL1CvRZrBd1etzu4DybqhhpkZl6qjBPCUMzYQQAcCvAfww2rGU0tUAVgOKS2os1znTHmBu37r/hC2hAGg9c8xRxACwcttBUyK6RBAtL4z2A/VLEm6eOiLsp67eg5X3kXY2wnGGtoCE6Ro1ZX6WF09dNwE/XVcZMeJYAHDTcx9a1gnfvGQaI2+VHC6CE0/+oKK8dGbeJInK4dobRjWMHdVOVxPVOZ3orifCsvP5JKoretRVkiUU6gEM1fxdFNqmkg1gHIB3iJJ3eRCAjYSQ+ZRSB0OW2da9F3YdxV0zR8fUOaeP6Y8lM4rx1FsHdPs8opg0NzcrHSnrA/39u4fw+3cPwSUAv15YhvllQzB3QiHzvPGWzQSUesq8fCYLc99T39+z2w/h2fcPQyAEvqCsi/b2SRSQaLjIkRa3oOjy1YESUOIJCKWYu3J7eKZuV5euXVmuWDAed2tsCsYKhcbEinZdPruq13fCLtCTXcpZdr40d2SbTqwkSyjsATCaEDICijC4AcCN6k5KaTOAsAsQIeQdAPc4KxCsEYX4/KaV1Nb6zKDJXs6ysp5GGtSDMrB0fSXTOKdNYcGynbgIouYq4gLBjCIOKNMFMz/Li/tmn4dbp40MP/uK2tOmqmag1LK4i1o1T/VWUgWJOlMHFNUpYJ14jqX62fXArPDvSgtzoq4u7Kp27GQIjkRXfp9qLuWxYvX+nBxzkmJoppQGASwB8AaAzwCso5RWE0IeJYTMT0YbIiFJNK6HqqS2npBybm7Ryn+KxGyc07r6zV25HfMZqwievC4+CAF+9MLHmPLLf2HtLrbHkWoALR6YjRljB5hUSn7G63xobkm4r6nxBFrcgoC1u49h8i//hcXP7cHi5z7EpF+8bXLjtDJeAsD0MQMwfUz/sIpINSIbDco9weWzu+srO0EynnPSbAqU0tcAvGbY9rDFsd9MRBusco7fNn1k3A/VznI22ctVteMsXV/JTIGg6oa17TOqm16taECGO3rKCk50ZIrwLP/BV6sAgoi5f4y6c19QMqWUyPSIGKdJl8yaCPglCau21eiM2UEZuPvlSpQMzgnbkVgrS4EQnQOGdobdHgiCECWpnna2neounz0xqpk1diT6OfeprPasugMAMMume6AVkdzc4g226Srzy4aEg5i0+deNumGA7eonEnYRH07XWb5pX9TZqbaOwGv/eZlpv0T1q1vWDHLJjNFML6KARDHn6ffDfZElUNr8Em5bU46NFfWmGXZQVs7RGcC5V7diSFWXz57mvdRdY0dKeB8li6qGZsvticiW2tVgm66iTYEQSTfM+ljaAzJEYj/JHcc+dnP/aHXndrxu5pcNQcngHFTUnsa5+RloC8iW1df8GptD58pSnzbbF1SOWf39iyI6Hli5YqeaQTce76XuTPFtNXZsrzmVULtInxIKVj7h0XzF4yVVlquKt5T1ashqUFAfS1AGrp4wGBsqv0x0U/sEkhy7DUurMlBrWRgN16qKB1C8kLwigURpuBiLEW1fnF82BP0yPLjj+Y908TfKCpJEXTUaXbFT1aAbi+qlO+/BauyobjiT8Ilmn1IfGcvYRdvOIpacIz1puTq1uAAPXDkWaYzERn5JxuvVx5HBa2jGxTVlg+EWCTK9YlyGQbXPAcCRxlbMXbndpFLQzSxDqzqfRBGUFZXhr747zlQC1S9Jur5YWpgD2ZCUIyDLKC3M0ZdrZXwuHlFEdUMzKmtPo+b42ZQ26NpRcXW3Udpq7ABol6O6o9GnVgrTivtDcX5ibY+OOnNwCQR+iWLZvJKYDIapGmyj3pcsU/gtVk0iQdTEaBw2r336FQRCcMf0UZaBjVZqCmMQokyVla1xlhjJDdkjijivMBdPLSzDT9dVhFNgyxTYUXPKVPeb1V+NK5Wrnnlf58TQEVRsEB5RDBvGtaS6QddId6/yrd5FaWFuwieafUoo5GV6YtquhRXl++A/qgCqZIK0ItU9MuxGL7dzm0LcKIKW4umtiorFiJWagqVXNqIOVJHckNVBoygvHaIgIBg6LmCwKwCR+6vWxnH9xUOxZuex8L6gREEB+IKhFNs2EvmlMmxvruTeg9W7SPREs0/pA5RC8/a3a6lraoeL4cmxfFN11CWlEx4ZTqbKVc9Vc/wstu0/AZHwXNjJQM2zpUVNPsdSU7C8wkznDA22Wu+jtJCeyCsSnbqqrqkdHkZUtFH1EK2/Nrb4sK68TrfNuL70GhL5peIKORLq89SqliVZxo4adk2WRLbD+C60nmk77pvpuJ2jT60UGpraYtqupSgvnalacYuJX1I6afBau+solm/eBwJFHeQVCbNmAidR6AXw2t3HTGq5SLN/lwCIggCPaJ4lsozR2hmmUzYuO2lQiECwZck0Uxt6ElOLC6CdBwbl5HoPRqKrUeGR6FNCodZi8LfariU/y4tl80oUlZEGo7+40zjp1rp211FTFSytQPCIxNKmwOk6LkEx5qo0tviwattB03GqmsJKrxxJHRlpsHDKxsUSLsa60k9cO16XZLEnoqysxE6VGHqebSQe+pRQyPCwb9dqu5FFk4YDFOHkZFKosE4iO4hTBq/GFh+Wb6q23J/pFbH029/AL1/fzw3KDpPmFkApZQYNGgcdAFgyoziqjj+REfjRiEdY9UR6kvegk/QpoTCkX1pM21ksmjwcs8cNSlrnj9Yx7QbX1DW1wy0qHiwsJJliZP9MeESBCwUHyfSIWD6/FDPGDrAVNOh1EZMx2mlVgRPnc1pYpSI9xXvQafqUUCg/etpy+zUXDmXuY2H1USUi+jFSx4zF1lCUlw7JomCOS4Cluxuna0iUMgUC0PMHnUTqtVOFVPceTAR9SihIFgPeiTMdzLTGsZDI6EdWx4zV1hApSZ4oCOHfPXRVCZZvqoZACDr4iiFuMj1iVPViY4sPw/MzsbmHG2R7O6ko/BKZfqNPCYWivAzm9ncOnMDUx7fiobklGFeYa7sQiXpMMnIcGTtmPLYGJZWBG7f+9SP4pc7fUark229qC+CxLfvgcSkqJJdAIlYF45jxiMCy+eOi9iPWJGKCoUY2h8Mi0ek3+pRQmDTiHOZ2vwRAkvHgP6p0Mzzjg7Z6Gd0R/WjXCGYUYoW56TqBACjeLrf+tRwypeCLg+iIxBSbFebXCy+wrGin0t2JEjk9l2T0nT4lFOxk+lTz3hsfdKSX0R1eCnb00foUCTKuvXAILh2Vz4xN4K6o9on2qCprT0dcJXR3CgVOzyUZfadPCYV6G/EIKsYHHellTBjaLy6DYVf1gpGMYCwh9uKHtXjxw1owArM5DkCgFLDRBpax3k9fdXXkdJ1k9J0+JRSsgtRcAkxqE+ODjvYyYvVScEovaGUEixR1ys0E0fnWef2xdf/JmJ6VSyTwBeWwS+9P11WYoo/nlw2J2eso1eoScLqPZHis9SmhkGkRpHbDxcNw3uAcPLZln+WDtvMy7HopJEMvGK1Os5ZMr4hAUAYhpM/FKDw4ZywOn2rD3/ccC6uFXAIw4xsD8f7BRlvPwyUACy4qwsbKLxHQxIEEZSAodwoJ7Tu2O4nozpz+XBilJtpiSmWhut5O0qeEwthB7Ic34xv9UTYsD0PPSYdVdTLAOZ9lJ6OUI6U7YBXOMZLpEbF8Xim+bvXjF6/vj/1mejjZXhde+bhWZycIysCyjVW4adJw/GXn0ajnkGTgpT11UY8zvuNok4juNEinapEcDvc+cpTNexuY23//Tg2qvjzLfMjGgdf4IRv31xw/G1WCd0UvqF6vqr7ZtLIxdgxViN33yl68/dkJ5vkCEkVTmx+Pv/F51Gv3RjqCMkQiANBHegdl4PndxyJ6GqmwdntEAolCV/GsI9hZ1MbOLLy7DNLcOyp14d5HDlNzopW5vfyYks441lqoRol98fA8vF/TGN6/eMowPHr1+abrxasXVK8nEhL2korWMfKzvLhr5mi8+/kJsJyv/JKMn7/W91YIgKL2mVZcgOWBfcz9kkzhFglExOadleEW8cSC8fivv3+i205DEeV2Z3rdZZDm3lGpSzLeTZ+qpzChKIe5XTT8bayFyirHxyrXpxUIALBm5zHUHD/LvGasOdG111MFgrHNViX5Mj0iUyD0ZUQCLJw4FEcbW5kzfRUCgILCFcOXEpBl5KS7kO7Wz7nS3S5UNzTbLvOorZGQzLoE3DsqdeHeRw7T4mePjMZBQVsL1Uoi28kpDwAVtact1UixhM9Hu16kJHmvVX1l6xo9HVFQ9Pt2kGini24k4onfoJSiMNeqxi4xvUdRIJYzve7IvdPTczL1Zrj3kcOcm89Oc/Gd0gHYduBUTLVQ7Xr3lDmUuiDa9R6aW8JMkvfQVSVYta3GkTakOlYCwSsChAgJyeV05zdH4n8/OIo2zeot3e1Cq1+yqLGbY3qPrT4JVfXNlmkuuiP3Tl9MBNdTSPS76VNCIcvLvt2Lhp+D/++742OqhWqU2O2BICRZv+q4rDjfMXcx9Xr3rN8Lv2Fwy/SIGFeYyzRCLd9UHcplpD+fohbpIxCSsHs9b3AuZEP2WXXyMGFoP+bH+9Bcc7Gmx7bsw+xxg1Jq8E3FRHAchUS+mz4lFI40soPXjjS2MR+ylURW1TNTiwuw476ZqG44g9vWlIcLoqt8eOTrLmdfNbanZHAO5jz9vk6toVZ/YxqhRMGkAnEJBIQoao7uCEsQBaVdHQkydKh1dT0uAZKs5LEClGjjQAR1EAEgWHgbZbgFU5oUt0gwZVR+1MmD8f2PK8wNl8sMn4sbchMGj7eIjT4lFIZbqI+stgPmj5rlOTI8PxMiI3eESJz/0IsHZuPJ6yZYDkJG1YREKZbNK8Fjm/eFcyBJsowAu9ZOQnELBL+5fgLGDsrB3JXbTftdgqJ66QhIkCmN6gpqxSPzSk2FkBpbfFHTe1CwBcKPLh+JKaPysfOLRjy34whEgUCSKVYsGB9TIJoKq7YFN+RGJt6BncdbxE6fEgrfGMT2PrLabsTKR3jzkmk6f3QViSbmQ49U9co4a1XTgb9w8yU40tiGNLeIB/7vU5w16pMSiFsECAiWzSvF3AnKB8lSoVAK3DR5GP70/mHT4OwWlAE7WtqJDI+IcUNymanGWaUv7VDf1Ibbnz8CtyAAoPjR5cW4cdKwuIV9oo2FvW1mHO/AzuMt4qNPCYXCXHbZTavtRqx8hBua27FkRjGe3nowrJ5wCTDV5I1GLB+zlU5RKzCq6pvx2OZ9oDKFT6KhWsHWxYYSRUBS6gw8tmUfstNcmF82BEMZtS0kCjy7/TCzhkNAVuI+/r6nFoB1Og45pEozYtcxQABM/l0b9yreW+p7X/VOja5kZjyDVqKMhb1tZtyVgZ3HW8RH0uIUCCGzCSGfE0JqCCH3M/b/lBCyjxCylxDyL0LIcKfb0OqXkObW33KaW2D6/bNgDSwdQQm3rSnH6vcOQSCKqmHl98rw3A8vwdTiAttt21BRj6mPb8VNz+7G1Me3YmNFve3fGsnP8qIoLx2PbdmHjoAcTpPdEVDz8CQ/TapfogaffPaUP5LO/6UPj2HLXZfhT4svgtdlvgevy9qPPz9LqSoXCZcA/OSKMchwR/4stDEhrHgVq7gDVpsmDO3n6Aoh3rakKurAriVSTI4WHm8RH0lZKRBCRACrAFwBoA7AHkLIRkqpNpT0EwATKaVthJAfAXgCwPVOtsOqM9jtJMZlv6qf90kIqyX+9P4hZmbMSCRimRsprkGMsaKakxXYBAB//eAIpo8ugFskEYWAEQJFgE8fMwArFkzQvAcJS2aMNql01JWXatQdek4GsrwiWnydk4BMj4il3xmDkf2zUVqoqBGf2XowYju0A0tdUzuo4dlQmXbLbLQ3zoy7MrDnZ3mx8KIirNl1LLxt4cSiHvsskkWy1EeXAKihlB4CAELISwCuBhAWCpTSbZrjdwG4yelG5Gd5sXBiEdbsjL+TaJf9ze0B3Ln2Y51+PlJmTCsS8THHkiVVxSrPjygQyDKNEqZnj7aAjKe31uDprTW4rDgfuw9/bTtAjMJ+qnJVjQIoKySvSABCTKoziVLMC9k51HMtm1eKB1/V2zsAdt3lTI9oKljkkygyPcY4+cTTG2fGXbG/NLb4sO4jfaLCdeV1+PGsMVwwRCBZ6qMhALSho3WhbVbcAuB11g5CyO2EkHJCSPnJkydjakRjiw/rys2dxLi8bmzxobL2tOWyOz/Li0yPiEMnW0ylLY3YWeom4mPWpkjwhlw009wC0twCViwYj8nn5pl+I1Fg/oTBEDWaGZdA8PC8Eoy3SBFiF9YY+X5NI34w5VxbvycAnryuzJSqnKV+0a28Qm6kPonCF0oP7nURXdqI7TWndKq77DQXfv7dcfC4BGR6RaS5Bfz8u+Pw4m2TTSlJuqqSdJLuSouRaGJNCaPSFdVTXyblDM2EkJsATARwOWs/pXQ1gNUAMHHixJh0GnZm5NqkcwFJxrJ5pVg0WW/eePjVT3VLUlEgyHCL8EsSZKrXi9sZ3BPljaKdTasqFLUt94Zm0Ube3Hccb/zXdKz/uA5/3n4YHlHA8o3VcaV7SHMJuP/KsZhWXICNlQ14eqs5svp/Pzgc9TwigDd+Mt0UCGhlmI+kOktziVi16ALkpnvCz2Lq41tNqrsd983E7NJBUQ3BXVVJOk1vjUSOJ1irN66ckkGyhEI9gKGav4tC23QQQr4F4EEAl1NKHbeOFeWloyOon8EZ0xmrM0yVB1+tAgiwaJIiGGqOn9UJBEDJpvnLfzsfU0blY0fNqbgGd7sfc6zuhqyPqbL2tOWgqXhTdeAvHxxBQKK6ojExQ4B5EwoVX/4JhUyh4BEFU9CfEQnA61Vf4S6NUIjkZRNJdRaQZZQW5oafCetZaMus2vECS7U8QTwSWSEV301PIFlCYQ+A0YSQEVCEwQ0AbtQeQAi5AMAfAcymlLKT/zsANQQNaf+ua2qHSMxeLcs37cPsUiUFQUXtaeZ5OwJSXIFMWqJ9zE65G0YaNBV1mDkZYCQuK87HnqNNADr190Qgug+weGA2Fk8ZZrLnbKzU17hwC0osgnFhsnJbTdiQHM0wrx0MIrXJ6lnEOpvsrbPz3gB/N7GTFKFAKQ0SQpYAeAOKNuA5Smk1IeRRAOWU0o0AVgDIAvAyUQbmY5TS+U62o66pHelul84wnO52hdVHRXnpCDBsBG6xM4ulVYI77fZEzNQaW3y4d30lfEHaZQ8l7aCpxjCoSLKM2q/bTQOlQNiBY16XgN/ecAEAmNRUxnY9evX5WDz5XF0RomnFBbr8USBEMXgbpIJHjJyh1qgGtFKdGdvk1GySz85TF/5uYiNpNgVK6WsAXjNse1jz728lug3RZoX5WV6m54kkdwZEsWa8i6cMc7xOqpG1u4/BF9QPlHY9lBpbfKhuaIa21Kg6aFY3nMGtf90TthkEZSXIbP6EQp1R/qbJw3DJuefg7pf15T3Vwdquv33xwGzds9K247Y15ZZBadEy1LJm93YHAz6b5HA6STlDcyKxMytcNHk4QBSVkVsk4YRq2mNYM95E0tjiY6a/9ktSVDXHhop63KNJBOcSgF8vLMP8siHIz/IiN90Nr0uEX+pcPYkCwauf6E0+68rrsHjyuTBq1wKyjEyPiMra03EPqGo7PKJgEgoZHhGywQ00EbpiPpvkcBT6lFAA7M0KF00aHtXzxDjjTSRK3h7zgLlkxuiIA5mictqr84YKysDS9ZVhtRNz1i3RUHZVfRZPVo2AhROLMHfl9oTYObwugj/cdKHOMKzCZ/ccTmLoc0LBLrHMHBOdgIw9YAq6/Dss6praLbO3VjecQW66G0V56cwkeo9t1tctZtUIyPSImLtyuyOR2Faz/+ljBkT8DRcGHI6z9Dmh4HTCsGQkIItXXVKUl87M3uoL5WvSpuLYcd9MnWDL9rqi1giI5M4Zz2DNZ/8cTvdDjC6aPYmJEyfS8vJy28c3tvjCgUoqaW4BO+6bGdcA5PT57Fwv1gFzY0W9rriMSABB0OccsmpztOvFe/+9LbUzh9MDscyK2adWCk7nGEp2ArJ41CWd3j2K9xFAcefaTxDQGJat2hztevGsYHpbamcOp7fRp4SC02HvPSWMPj/LG9bNN7b4HG1zLCofXvSEw0l9klZPIRVwOmFYIhOQRUvKFy+JaLNVYjojPEEZh5P69KmVAuC8MTMRxtFEq1i6y6DbU1ZWHE5fpk8ZmnsCyTZeJ5uNFfUmGwS3KXA4SYcbmnsKvbF6lhbudpoacA8wjhVcKKQYfUHFwoPOuhfuAcaJRJ8yNPcEUrF6lh2jd6IM4xxn0XqAnfUF0RGQce8re/l744ThK4UUJJVULHZmlXzm2XPo7epJTtfhK4UUxa6bZyKxM6vkM8+eRV9QT3K6BhcKHEvsxBXw2IOeRSqqJzmpBVcf9XK64mViZ1bJZ549j1RST3JSD75S6MVsqKjH1Me34qZnd2Pq41uxsaI++o802JlV8plnzyQV1JOc1IQHr/VSnAyCs7Pa4H7vPQv+vvo8PHitr+Gkl4mduAIee9Bz4N5inEhw9VEvhev6OSy4txgnGlwo9FK4rt85elNgHvcW40SDq496MdzLpOv0NlVLUV462gNB3bb2QJCvIDlhuFDo5XBdf/z01qJAhCgV+PR/czgKXH3E4VjQG1UtdU3tSHOJum1pLrFH3xPHWbhQ4HAs6I3G+t54Txxn4UKBA6B3GVOdojca63vjPXGchQevcXqdMdVpemOgV2+8J05M8OA1Dpveakx1kt5orO+N98RxBq4+6uP0RmMqh8OJHy4U+jjc8MjhcLQkTX1ECJkN4HcARADPUkp/ZdjvBbAGwEUAGgFcTyk9koi2jF+2BWd8QI4X2Lv8qvB2Kz1rNP2rcb/6d6ZHRKtfCv9XHWhZ56o5fhbba07B6xJAANQ2teLQqTY0nO5AyeBsTBpxDsqPnoZHJCjslw5fUMY5mR4QAO8dPIGDJ1pRPCALw/IycPBkC/qludERDOJYYzty0twYkJuGNl8QXpeAgCSj+stm5HrdCMoUJGRXUpWMA7M8uHH1TrT7JTS2+ZGb5lbUDZkeZHgFtAdkTCjqh3a/hI+Ofo0BOV4MykkHBZDhEdEvw4OmNj86/BLys7yYNOIc7Kg5hW37T2JAjgfTR/fH/uMtaPUFIAgEYwflYHBuOr5qbsf+r87CLRK0ByQ0nvWjINuDNLeIlo4g+md70S/DjeNn/RiQ5UFuhgfjCnPgdokIBCVUNZxBQJLhD8rwugQ0tfnQ5peR7hExNC8DQ/LSUVqYi6ZWPypqT+Pc/AwcPNGCj442wS0SFOVl4Dulg1A8MBuNLT5UN5wBQFFamBt+V+p7KsjyYMqogvD71L5j7XstP9yI9w6ewvTRBZg4It92n4rUr6LZAez2V6t+adVHtee1OsZ4bCAo4UhjG8qG9kPxwOyY2hkPfcVWksj7TIqhmRAiAjgA4AoAdQD2APgepXSf5pj/ADCeUnoHIeQGAN+llF4f6bzxGJrPvX+LaduRX11laWyNZoQ17l94URHWfVQHKlP4JAqXAARlJUOpJFNQSpHudunO9fCrn2LNrmMx3QcncVxWnI+dhxoRDC2g3CLBU9dNQPmRr03vyS0SiAJBR0CGVyQgAgm/15ue3YXtNY268z5/6+S4SpwunFiEdeV1UZ0B7PZXAOE2S5SCEII0l4iOoMTso9rztgeC4eMjXSMoyeFnCACLpwzDo1efb6ud8dBXHCYcuk9LQ3OyhMIUAI9QSr8T+vsBAKCU/lJzzBuhY3YSQlwAvgLQn0ZoYKxCQV0hGMnyAEEqmNJMb14yDXNXbrdMP81KTx0LaW4BL9x8CRb8cVdcv+ckD7cA2H3NaW4Bq753AW5Z85Fp358XX4Q7//ZJxJTmdvoVKw16tHTp8fRXq+/Aqi3RrvH2T6YjL9PjWFp3FSdTxacyDt6npVBIlk1hCIBazd91oW3MYyilQQDNAPINx4AQcjshpJwQUn7y5MmYGsESCADQ4gfT2FpRezqiEZZlpI0FtyDgvYOn4v49J4nEkArCLQh4c99x5r439x2Pq8Qp6xpGZ4BoTgPx9Fer7yDea1TUnk6Ic0NfcZhIxn32OEMzpXQ1pXQipXRi//79Y/ptjoUgzfKAaWwtG9ovohGWZaSNhYAsY/rogrh/z0kiMayoA7KMb5cMZO77dsnAuEqcsq5hdAaI5jQQT3+1+g7ivUbZ0H4JcW7oKw4TybjPZAmFegBDNX8XhbYxjwmpj3KhGJwdQ2tU1lL16FXMKM/igdkRoz9Z0aGLpwxDmluAV1Rmlq7QE05zC3CLBC4BunNNHJGPxVOGOXmbnC5yWXF++L0BIZvCwjLme3KLBGlu5WBv6N9PXDses0oG4bJi/UL3suJ8zCoZFFeJU7VfRYpCjhatrN2vbbNLUO4j2+ti9lHjd6A9PtI1XIbRZfGUYSgemJ2QqOq+EqmdjPtMlk3BBcXQPAvK4L8HwI2U0mrNMXcCOF9jaP43SunCSOeNN6KZex/pvY++OtOBNI+INJcL6W4Cr1tAu1/m3kfc+4h7H6UoDtxn9xqaAYAQMgfAb6G4pD5HKf05IeRRAOWU0o2EkDQAzwO4AMDXAG6glB6KdE6e5oLD4XDiovuFQiLgQoHD4XDiotu9jzgcDofTA+BCgcPhcDhhuFDgcDgcThguFDgcDocTpkcbmgkhJwEcjfPnBQBSMZyYtys2eLtiJ1XbxtsVG11p1ylK6WzWjh4tFLoCIaScUjqxu9thhLcrNni7YidV28bbFRuJahdXH3E4HA4nDBcKHA6HwwnTl4XC6u5ugAW8XbHB2xU7qdo23q7YSEi7+qxNgcPhcDhm+vJKgcPhcDgGuFDgcDgcTpheLRQIIdcRQqoJITIhxNJ1ixAymxDyOSGkhhByv2b7CELI7tD2vxNCPA616xxCyFuEkIOh/+YxjplBCKnQ/K+DEHJNaN9fCCGHNfvKktWu0HGS5tobNdu783mVEUJ2ht73XkLI9Zp9jj4vq/6i2e8N3X9N6Hmcq9n3QGj754SQ73SlHXG066eEkH2h5/MvQshwzT7mO01Su35ICDmpuf6tmn0/CL33g4SQHyS5Xb/RtOkAIeS0Zl8in9dzhJAThJAqi/2EEPJ0qN17CSEXavZ1/XlRSnvt/wCcB+AbAN4BMNHiGBHAFwBGAvAAqARQEtq3DkoKbwD4A4AfOdSuJwDcH/r3/QAej3L8OVDSiWeE/v4LgAUJeF622gWgxWJ7tz0vAGMAjA79uxDAlwD6Of28IvUXzTH/AeAPoX/fAODvoX+XhI73AhgROo+YxHbN0PShH6ntivROk9SuHwJYyfjtOQAOhf6bF/p3XrLaZTj+Ligp/xP6vELnng7gQgBVFvvnAHgdSqbTyQB2O/m8evVKgVL6GaX08yiHXQKghlJ6iFLqB/ASgKsJIQTATADrQ8f9FcA1DjXt6tD57J53AYDXKaVtDl3filjbFaa7nxel9ACl9GDo3w0ATgCIrV6rPZj9JUJ71wOYFXo+VwN4iVLqo5QeBlATOl9S2kUp3abpQ7ugVEBMNHaelxXfAfAWpfRrSmkTgLcAMKNwk9Cu7wH4m0PXjgil9D0ok0ArrgawhirsAtCPEDIYDj2vXi0UbDIEQK3m77rQtnwApymlQcN2JxhIKf0y9O+vALAL+nZyA8wd8uehpeNvCCFOlZiy2640Qkg5IWSXqtJCCj0vQsglUGZ/X2g2O/W8rPoL85jQ82iG8nzs/DaR7dJyC5TZpgrrnSazXdeG3s96QohaujclnldIzTYCwFbN5kQ9LztYtd2R5+XqUtNSAELI2wAGMXY9SCndkOz2qERql/YPSiklhFj6BYdmAOcDeEOz+QEog6MHiq/yfQAeTWK7hlNK6wkhIwFsJYR8CmXgixuHn9fzAH5AKVUrnMf9vHojhJCbAEwEcLlms+mdUkq/YJ/BcTYB+Bul1EcI+X9QVlkzk3RtO9wAYD2lVNJs687nlVB6vFCglH6ri6eoBzBU83dRaFsjlGWZKzTbU7d3uV2EkOOEkMGU0i9Dg9iJCKdaCOAflNKA5tzqrNlHCPlfAPcks12U0vrQfw8RQt6BUkL1FXTz8yKE5ADYAmVCsEtz7rifFwOr/sI6po4o9clzofQnO79NZLtACPkWFEF7OaXUp263eKdODHJR20UpbdT8+SwUG5L6228afvuOA22y1S4NNwC4U7shgc/LDlZtd+R5cfURsAfAaKJ4znigdICNVLHcbIOizweAHwBwauWxMXQ+O+c16TJDA6Oqx78GANNLIRHtIoTkqeoXQkgBgKkA9nX38wq9u39A0bWuN+xz8nkx+0uE9i4AsDX0fDYCuIEo3kkjAIwG8GEX2hJTuwghFwD4I4D5lNITmu3Md5rEdg3W/DkfwGehf78B4Nuh9uUB+Db0K+aEtivUtrFQjLY7NdsS+bzssBHA4pAX0mQAzaGJjzPPK1EW9FT4H4DvQtGr+QAcB/BGaHshgNc0x80BcACKpH9Qs30klI+2BsDLALwOtSsfwL8AHATwNoBzQtsnAnhWc9y5UKS/YPj9VgCfQhncXgCQlax2Abg0dO3K0H9vSYXnBeAmAAEAFZr/lSXiebH6CxR11PzQv9NC918Teh4jNb99MPS7zwFc6XB/j9aut0Pfgfp8NkZ7p0lq1y8BVIeuvw3AWM1vbw49xxoA/57MdoX+fgTArwy/S/Tz+hsU77kAlPHrFgB3ALgjtJ8AWBVq96fQeFY68bx4mgsOh8PhhOHqIw6Hw+GE4UKBw+FwOGG4UOBwOBxOGC4UOBwOhxOGCwUOh8PhhOFCgcOJEULIkVAQGIfT6+BCgcNJMkRJ5f3/dXc7OBwWXChwOBwOJwwXChxOnBBCLiFKYZ/ThJAvCSErQykT1EIovyFKsZQzhJBPCSHjCCG3A1gE4F5CSAshZFP33gWHo6fHJ8TjcLoRCcBPAJRDST72OpQCO7+FkndmOpTiP80AxkJJLb6aEHIpgDpK6X93R6M5nEjwlQKHEyeU0o8opbsopUFK6REoyebUdNQBANlQhAGhSsGnLy1OxeGkDFwocDhxQggZQwjZTAj5ihByBsAvABQAAKV0K4CVUBKXnSCErA6l9uZwUhouFDic+Pk9gP1QakPnAPgZlAyWAABK6dOU0oug1GYeA2CpuivZDeVw7MKFAocTP9kAzgBoCeXd/5G6gxByMSFkEiHEDaAVQAcAtRLccShpxjmclIMLBQ4nfu4BcCOAswD+BODvmn05oW1NAI5Cqby2IrTvzwBKQl5LryattRyODXg9BQ6Hw+GE4SsFDofD4YThQoHD4XA4YbhQ4HA4HE4YLhQ4HA6HE4YLBQ6Hw+GE4UKBw+FwOGG4UOBwOBxOGC4UOBwOhxPm/wdfG7oiU3XYwwAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"df_base.plot.scatter(x=\"last\", y=\"next\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What is going on at the boundary? What is the re-election rate?"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Re-election rate: 90.93%\n"
]
}
],
"source": [
"info = pd.crosstab(df_base[\"incumbent_last\"], df_base[\"incumbent_next\"], normalize=True)\n",
"stat = info.to_numpy().diagonal().sum() * 100\n",
"print(f\"Re-election rate: {stat:5.2f}%\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Regression discontinuity design"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"How does the average vote in the next election look like as we move along last year's election."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAENCAYAAADgwHn9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9JElEQVR4nO3dd5icZbn48e+9s71la/pukk1vhLIsIfQmAZWqQFARRVEEPceCP6zHdkQ9R/F4jiCgFOmIoEGaFEMJSUhIJQlJNn2zKbvZ3nd2nt8fMzv7vjPvlC0z2+7PdeVi3jrPzi7vPc/9NDHGoJRSSgEkDHYBlFJKDR0aFJRSSvlpUFBKKeWnQUEppZSfBgWllFJ+GhSUUkr5xSUoiMgDInJMRD4IcVxE5HciUi4im0Xk5HiUSymllF28agoPAUvDHL8EmOn7dzNwTxzKpJRSKkBcgoIx5i2gJswplwN/Nl6rgRwRmRCPsimllOqRONgF8JkEHLRsV/j2HQ530dKlS83LL78cy3IppdRIJKEODJWgEDURuRlvioni4uJBLo1SSo0sQ6X30SGgyLI92bcviDHmPmNMqTGmtLCwMC6FU0qp0WKoBIXlwA2+XkiLgXpjTNjUkVJKqYEXl/SRiDwBnAsUiEgF8B9AEoAx5g/Ai8ClQDnQAnwuHuVSSillF5egYIxZFuG4AW6NR1mUUkqFNlTSR0oppYYADQpKKaX8NCgopdQQZIxh19FGPJ74ro6pQUEppYagbzy9iYvueovPP7w2ru+rQUEppYYYj8fw3AbvUK0VO6qobe6I23trUFBKqSGmsc1t2z6uQUEppUav2hZ7EDje1B6399agoJRSQ0xNQFCobtKaglJKjVqBbQjVWlNQSqnRq7al07Z9qK6VFTuOUR+wPxY0KCil1BATWFO476093PjgWq68eyVdMR63MOzWU1BKqZGosq6VdftrqWpsZ9vhBsdz9lQ3s6eqiZnjsmJWDg0KSik1yN7eVcWND66NqhZQE+PuqZo+UkqpQXbvm3uiTgsdbYxto7MGBaWUGiDWeYo8HsOXHlnHkjtf551d1SGvqaxrZeXu0McDHWto61cZI9GgoJRSA+CZ9ytY9ON/8vWnNmKM4f0Dtbyy9SiV9W18+k9rQk5s97eNhzC9aDs+pjUFpZQa+r71l000trt5bsMh3t9fy/7jLbbjq/YcD7rGGMOz63uWo8/LSI74Pke1pqCUUsPL9iONNLXZxxQ8vuZA0Hm1LZ2UH2sCICUxgZvPLol471gHBe19pJQa1ZZvqmTNnuPcfHYJU/IzONbQxs9e2E5akotry4o4uTi31/esbmzHBOSEXtl6hKrGdgqzUvz7Glp7AsfY7BTmjI/c1fRog6aPlFIqJvZVN/PvT27gsTUH+P7fPgDgnjd3s3xTJU+tO8hVd7/Lr17+MOJ9AtsLqprag2Y2dXsMqwNSSE3tPbOhZiQnMiU/I+J77a1u5qO/e5tfvvwhu6uaIp7fWxoUlFKj1j+3HaH7ef72rmpaO7rYcKDOds6Taw9GvE9zh32q64M1LRx3mMRuw4E6nlp7wN8bqdkSFDJTEpmUkxZVubdWNnDPit3sq26O6vze0PSRUmrU2n640bb9/v5adhyx72to7cQYg4j493V2efjFSx/S1Obmux+dS2tHl+2a/cdbGJ+dGvR+D6zc63/9qdOKOXtWoX87IyWR5MTov6cnuYTFJflRnx8tDQpKqVFr3f4a2/bT6w7S2ml/wLs9hna3h9Qkl3/fC5sP86d3vA/4zNRElpUV2645VNdKghDWY2sO8PeNlf7tzJTgx7EILJ6WT3OHm2MN7RyxNDKXTskjw+Ga/tL0kVJqVDpS38bBmlbbvuWbKh3PXb3nOHe+uJ3tvjmJfvfGLv+xP72z15YGAujyGPYFdEl1YmtTSHEFHU9NdPHEzYtZftuZnDLF3uB91qyCiPfvCw0KSqkRo741+qmlA2sJ4dz44FrufWsPX3h4HQA5aUm244FBIdDJxTkR36P7W/9d1y7y7/vZFQv8rwsy7WMYzp5ZSCxo+kgpNSL8aPlWHnp3H1efPJlfX7Mo4vnr9tX2+j0O1bXS0uHGFZAbKg/TC2hMWhInTM5hfUADdqDu9NFliybR1Oamy2O47MSJ/uN1AQFv3oTsXpY+OlpTUEoNe/WtnTz07j4A/rq+goa2yDUGa02hpCC4K2hWqvN35urGDg7X2weQvbc3dK0jPyOZGWMzI5anu6bgShA+c/pUbjxjGkmunkf0pQsn+F+fO7uQhEiNFn2kQUEpNeytDXgoV0WYH6ip3c22Sm/7gAg89LmyoCkmFk3Ocbz2aGNb0KjicLWOvIxklkzv6SV0UnGOY9dTp4ZmqwvmjOXTi4u5aN44fnn1CWHP7Q8NCkqpYS9wUFikoLDhQK1/fMKc8dkU56fzp8+WkuTyfvs+f85YxgS0G3T78HADnV32wWpHwkw9kZeRTElhJj+9fD6XLBjPnVctdLx3pKCQ6ErgZ1cs5P4bShnn0N11oGibglJq2Fu9t3dBYa3lm/2pU729ek4qzmX5bWeysryayxZN5Nf/3Ol47eaK+l6VLd/XQPyZ06fymdOnApCdFvzojUX30r4YGqVQSqk+qm/tZGulffnKSNNLv29pTyidmud/PXdCNnN9DbihHtK9DQpOM59mpwbXFJy6pA4GTR8ppYa1tXtrgtYjCFdT6Ozy2KayKJ3iPOFdZoiH9I6jjY77u00vtDda52WkBJ2T3Yf0UbxoUFBKDWubKuqC9oULCtsPN9Dim5ZiUk4aE0PMN5QZovdRJOfOHmvbzneoKTi1KQyV9FHcgoKILBWRHSJSLiJ3OBwvFpF/icgGEdksIpfGq2xKqeHrYE3wyOGqpnY63B6+9sQGrrl3FXstE8dZ2xNKp4aeFruvD+myaXm27cAxDeCcPhpVNQURcQG/By4B5gHLRGRewGnfB542xpwEXAfcHY+yKaWGt4ra1qB9xxrauHtFOcs3VfLe3hruf3uP/5itPSFE6giie0g7DSAbk5bExfPHAZCYIJw6NS/onKHc0ByvmkIZUG6M2WOM6QCeBC4POMcA3Z/wGMB5EhKllLJwCgrVTe389rWe+Ym6Vz0zxgTUFIIf2N2iCQqLinIcr/vJ5Qu49bzp3HfDKYwfE9x91KmmkJ40NBqa4xWaJgHWSckrgNMCzvkR8E8R+SqQAVzodCMRuRm4GaC4uNjpFKXUCPHhkQZu/8tmphZkcNc1i0h02b/Htru7ONoYPEagOmAtg6yURG56aC1v7aryjzHISk1k1rjQK51F+uae5BIWThrDEw7XjctO5faL54S8NrBNISPZFbMRyr01lBqalwEPGWMmA5cCj4hIUPmMMfcZY0qNMaWFhbGZEEopNTT89B/b2HKonuc3VfLchkNBxyvr2vw9jyblpDk26gI0trt5/cNjtkFnp0zJdcz3d4tUUyiblsek3OBG6mi6lgb2PhoqqSOIX1A4BBRZtif79lndBDwNYIxZBaQCsZkbVik1LKws7xmU9tf1FQC8tOUwtz62nnfLq6mo7WlknpSbZlv/OJJw7QkQOSicN3sshZnB75eRHPkBH9imMFQamSF+6aO1wEwRmYY3GFwHXB9wzgHgAuAhEZmLNyhUxal8SqkhJrBbaWObm0N1rXz1iQ24PYZXth7hdMucQpNz06hKTODDI+HHEXQrmxZ+1bJI397PnT2WnHT7N34RSE+OoqaQOnRrCnEpiTHGLSK3Aa8ALuABY8xWEfkJsM4Ysxz4JnC/iHwdb6PzjcYEDklRSo0WHxyyjxzecaSRB97Zi9s3aZHbY3jbt9YxQFFuetT3XlZW7J/eIpRI396nF2YEDZpLTBDbsp2hBKePhkYjM8RxmgtjzIvAiwH7fmh5vQ04I17lUUoNPR8cqmdTRR3TCjKCFsFxe4x/CUwnk3PTaGwLv9gNwJ1XLQxaPtNJalLo7PqpU3MREQKf/4ET5YWSEVCbSEkchUFBKaXCeX9/DdfdtzrqB2ugybnptFnWV05MEHIzkoPSULnpzrOfBnL6xj+tIANjDHdde6J/X1qSK2hd597e2zOEkiIaFJRSg2rX0UZW7TnOD/++Neprbr94Nlsr63lxyxH/vqK8NGaNy+TZDYdoae/id8tO4pt/2RgUFHLSnXsoRZKalMAb3zwn6IGekdL7oBCoy6NBQSmlqKht4WP/+w7tbk/U15xUnMPnzpjKzqNN/qAgAuOzU0l0JfDcV3qy0E6DxHL7GBQSxLm9IC2KhuVIhlJQGErjFJRSo8zbu6qDAsKvP7mIX1690L997uxClpUVkZuexFfOnc5TN59OenIiJxbl8JnFU0hMEG49d0bQwDZwnngu2vRRIFeIBuRouqBGoukjpZQC9h1vtm1//6NzuerkSYgIWalJvLe3hpvOnEZRXjo/v3Jh0Df1n16xgB9dNj/kIDSnmkJf00ehOhXNHJcVdTfYULSmoJRSwD7L7KX/c92JfOGsEv+D/9KFE/jRZfMpyvN2NQ3V1TPcqOQx6cHTSSQn9u2xF2oaiu9dOpe0JBci8IdPnxL1/c6e1TMjw2WLJvapTLGgNQWl1KDZf7xnRPKU/IwwZ/ZNYPqor7UECJ0+Gj8mldXfvYCG1k5/AIvGnVct5PvPbSE/M4VPLZ7S53INNA0KSqkB09bZxW9e3Ulqkouvne+c5+/m8Rhb+mhqfvQP1GhlByyUk5vRu/aE3PQkals6ASgpDB20xqQlObZfhDMpJ40HP1fWq2viQdNHSqkB84c3d3PfW3v43eu7+Mfmw2HPPdbYTlunt5E5Jz2pX9/iQwkcOdzbnkf3fPoUEsQ7I+qdV50wkEUbsrSmoJQaMNY1DB5YuZcrTprEWzur+POq/XzilMksXTDef9xaS4hF6gj6nz5aXJLPyjvOJ9mVQL7D5HcjkQYFpdSAqKyzL3aT5EqguqmdLz/6Pi0dXazec5xzZxeS6ltMxtrIHIvUETjVFHrfHXXCGOc1nEcqTR8ppcLq8hi+9Mg6ltz5Ou9YJqALtGKHfVLjI/Vt/GHFblo6vKN9m9rdtrWS91kamacO0ZrCaKRBQSkV1rPrK3hl61Eq69u449nNIc/7145jtu1Dda38MWACu51He/rz22oKBbGpKQQGhb4OXBtNNCgopcL657aj/tdO6yEDtHS4WVkeuhbRrfxYE+BdK3mHJUDEqk0hcPBa4PoHKpi2KSilwjpY02Lbfnd3Ne/vq+XasiLGZnkXpX92/SF/miiclz84QnN7F+PHpPhTSWlJLmaHWSu5PwIHqiW7hs4U1UOVBgWlVEjGmKCg8NkH3qOzy7DtcAP3fPoUPB7Dgyt70kTJrgQ6upwnuNt1rIldvtpCt0sXTojbymNhBj8rH00fKaVCOljTSnNADaB7vYOXPjhCW2cX97+9h91V3m/9WSmJfHvpbNv5p0RYC/ma0skDWOJgF8wZC3inuFgyQ5d9j0RrCkqpkDZV1IU9fsYv3uB4c4d/+5pTi1hUlGM75+MnTOBoQ5tje8TU/HTKpuUNRFFD+ukVC1iw9iBnzCjo9ajj0UiDglKj3MGaFlo6upg1LjNo0rlNB+vCXmsNCFPy0/nyOdODUjQXzB3Hy1uPOAaFW8+bEdWaxv0xMSeNr180K6bvMZJoUFBqFNtwoJZl96+mrdPDrHGZ/Ojj820plnX7a6O6z8cXTeSXVy8k3be2wNmzCnlrZxUfmTeOorz0oOkl3v72eTS0dTJ/4piB+2HUgNCgoNQo9ptXd/rnH9p5tIkbH1rLG988h8m56VTWtbIxQk2h25fOLvEHBIAHbzyVXccamVGYCcD1pxXz0gfeVdL+4+PzejWbqIovDQpKjVLbDzfwdsAI5Q63h7tX7ObnVy7khQgT2nXLSklk/sRs2z5XgjBnfM++M2cU8McbSmnucPPxE4bO2gEqmPY+UmqU+uPbPd1IrQ2wf1l3kEN1rfxjc6V/X0lB6MFlnywtitguICJcOG8cl584KeRiNWpo0KCg1AixuaKOR1bto6GtM+jYwZoWlv72La68eyV1LR20u7tsD/0HP3cqpb6uo51dho/85k02VdQD3mmjrysrcnzPSTlp3Hre9Bj8NGqwaPpIqRFg//Fmrrl3FW2dHt7fX8tvrzvJdvyXL3/oX0f4wZX7OH16Pu1ub1vCtIIMTi7O5RsXzeL6P64BsI1NOGdWIdMKMm33S05MYO33LiQj2RV2IR01/OhvU6kR4Kf/2OZvMP7bxkrqWjpYWV5Nu7sLd5fHtuDNPSt2s2r3cf/24pJ8AJbMKOC282bY7pufkcwdl8wlL2DFssLMFMakJWlAGIG0pqDUMHesoY03PrTPUHr+r9+kprmDjy+ayKdOK7Ydc3s8vLa9Z5K706fn+19/46JZHKpr5bkNh5iUk8afbypjemGmbcprgLHZo2PBmdFIg4JSw9xD7+7DY+z7anyDyp7fVBk0d5HHwNbKBv/2YsuI4oQE4TfXLOKWc6dTnJfuXxAnL2CcwdgsDQojlQYFpYa5wFpCoHBjDaYXZjA2O9W2T0SYFTBraXaa/VGRk6aL1YxUmhBUahhr6XAHzTraG9bUUTiBXU4TXdqtdKTSoKDUMLa1soGuwNyRg/RkFz+/ciFZqT3f+DNTErnu1OIwV4WWpA3MI5amj5QaxqwT1o3NSuFYY7vt+I1LpjIxJ5VPnFJEXkYyxXnpPLu+ghOLc/jYCRPJy+hbGigwvaRGjriFexFZKiI7RKRcRO4Icc41IrJNRLaKyOPxKptSQ427y8Mjq/Zx94pyOtzOC9YAbPYNMAO44qRJtmOFWSn86LL53Hz2dP/D/8yZBfzm2hO54fSpvQ4I937mFNKSXCyaPIZPxngNBDV44lJTEBEX8HvgIqACWCsiy40x2yznzAS+A5xhjKkVkbHxKJtSQ9HPXtjOQ+/uAyApIYEvnl1iO26MYXdVMy9vPeLfd/H88dz31h7/dlFu2oCW6eL549nww4tISUyI+XTXavDEq6ZQBpQbY/YYYzqAJ4HLA875IvB7Y0wtgDEmfJcKpUaox9cc8AcEgFctYwq6/fqfO7nwN2/6axFJLmHBJPukdGnJA78ecWqSSwPCCBevoDAJOGjZrvDts5oFzBKRlSKyWkSWxqlsSg0JHo/h+3/bwnef22LbHzhwbOfRRu5eUW7bN3dCNimJLmZbcv0XzR0Xu8KqEWsodSFIBGYC5wLLgPtFJCfwJBG5WUTWici6qqqq+JZQqRh6Zn0Fj64+ELS/qrGdqsZ2Wjrc/G3DIT5y11u2wWpJLuELZ3nTS7+4eiGTc9M4a2YBy07rW88iNbrFq/fRIcA6zeJk3z6rCmCNMaYT2CsiO/EGibXWk4wx9wH3AZSWlkbui6fUEObu8vDYmgNkpCTy6raeNNGFc8exp7qJPVXeWsL2ww08te6gbY0DEXj+tjOZOS6TlERvquik4lze+X/nx/eHUCNKvILCWmCmiEzDGwyuA64POOdveGsID4pIAd500h6UGsEeXLmP/3xxe9D+Oy6ZwwMr9/qDwsryal7+4IjtnGVlxSyYpMtZqoEVl6BgjHGLyG3AK4ALeMAYs1VEfgKsM8Ys9x37iIhsA7qA240xx0PfVanhzykg5GckM70wg3kTehqO733L/v3ozqsWctXJgc1ySvVf3AavGWNeBF4M2PdDy2sDfMP3T6kRqb61k3fLq6msbyM9RO+gsml5iAjzApa47Pb9j85lWZm2F6jY0BHNSsVJa0cXH/vftzlY0xr2vDLfrKVzxmchAiag5ezi+eNjVUSlhlTvI6VGtOc3V0YMCNATFNKTEzmpKMd2bOGkMRTlpceieEoBGhSUipsn3wvubhpobFYKc8b3pI3u+fQpXHnSJLrHi924ZGqMSqeUl6aPlIqDHUcaWX+gDvCOK/jbrWfwxYfXUVnfxiULxvOZ06fw9w2VXHPqZFwJPSOGx2Wncte1J3L7xbNpbnczUyeiUzGmQUGpOHh8zX7/64/MH8/8iWP45zfOYfPBOk4qziUt2cWS6QUhr5+YM7DzGCkVigYFpWKsoa2TZ96v8G9f7+s5lJmSyJIZoQOBUoNB2xSUioF3y6v549t7aG5389R7B2nu6AJg1rhMlkS52plSg0FrCkoNsIraFj79pzV4DKw/UMumgz1rHnz+jGk6y6ga0jQoKDXANh6s809Y9+KWnqkp8jKSgxbCUWqo0fSRUgOsodXtuP/zZ0wlNWng1zhQaiBFHRRE5JMh9n9i4Iqj1PBX09wetC83PYnPnzltEEqjVO/0pqbwpxD77xuIgig1UlQ3dQTtu/3iOaQna7ZWDX0R/0pFpHtx2ATf1NfWVrISoC0WBVNquKpptgeFL541jWtPLQpxtlJDSzRfXcoBgzcY7A44dgT48UAXSqnh7LglffTw58s4Z1bhIJZGqd6JGBSMMQkAIvKmMeac2BdJqeHtuCV9lJ+RPIglUar3etOmcJXTThGZPkBlUWpEOG5JH+VpUFDDTG+CwhYRucS6Q0RuAdYMbJGUGr6MMdRqUFDDWG+Cwk3AH0XkbhGZISIvAV8GdJVwpXwaWt24fSPXMlMSdVyCGnaiDgrGmJeAhcCZwA7gOHCqMWZzjMo2ojyyej+ffeA93t9fM9hFUTFUbWlkzs/UWoIafnozeC0T+G9gDHAXcClwY2yKNbIcrGnhB3/7gDd3VnH1PasGuzhqgLV1dvHmziqqm9pt3VE1daSGo96MptkMrAROMMbUi8ijwCMicpkx5mOxKd7IsP1ww2AXQcXQd5/bwrPrDzFhTCrfXjrbvz8/I2UQS6VU3/SmTeEOY8xnjDH1AMaYjcCpwM5YFGwk8ZjI56jhqba5g2fXHwLgcH0bz2867D+m3VHVcNSbNoWnRSRJRM4SkWt9u13AD2JTtJHDYzQqjFSvbjtq237jw2P+19qmoIaj3rQpLMRbK7ifnnmQzgEeiEG5RpQurSqMWC9sORzymLYpqOGoN+mje4AfGmPmAJ2+fW/i7Y2kwtCawshU19LByvLqkMcLMrVNQQ0/vQkK84FHfa8NgDGmGdAVxSPQmDAyvbmzyj8mwYnWFNRw1JugsA84xbpDRMrwTpinwtD00ci0tTJ8r7LivPQ4lUSpgdOboPAD4AUR+TGQIiLfAZ4Bvh+Tko0gmj4aGXYdbeTuFeXsOtoIwDZLUPj8GT0L6CQnJvC9S+cytSAj7mVUqr+iHqdgjPmHiFwM3AysAIqBK4wx62NUthFDY8LwZozhnjd3c9erO+nsMjzwzl7evP082/iTzy6ZwsxxmeytbuZTpxUzJV8Dghqeog4KIpIMlOFdV6EGyAD+XUQwxtwQo/KNCF0BUcEYg4iEOFsNNS9uOcKvXt7h365u6uB/3yj3z4aamZJIUW46y8o0EKjhrzcjmh8GFgHP411cR0Wps8tj2+7yGBJdGhQGS1tnFy9sPszMcZmcMDnHv9/d5eGlD45QWddKfmYKF80bx5i0JJ7fVBl0jz+82bPe1JzxWSQk6O9TjQy9CQpLgWnGmLoYlWXE6nDbg4LbY0jUyTMHzS9f/pAHV+4jySX89ZYl/sDw0Lv7+NkL2/3nLZmez6M3ncbqvcfD3m/exOxYFlepuOpNQ/MBoM8dr0VkqYjsEJFyEbkjzHlXi4gRkdK+vtdQ09llTx9pb6TB9fiaA4D393LTw+swvvTe3zfaawTv7j7O85srqWvxDsspyEzmaxfMDLrf3AkaFNTI0Zug8Gfg7yKyTETOt/6LdKGIuIDfA5cA84BlIjLP4bws4N8YYQv3ONUU1OBod3fRbvl9VDW287eNh6ht7uCDyvqg82//S8/M8KeV5HPz2SXMGZ9lO2eeBgU1gvQmfXSb778/D9hvgJII15YB5caYPQAi8iRwObAt4LyfAr8Ebu9FuYY8pzYFNTj2VjcH7XvgnX0ku1yOvcQ6LL+700vyyUxJ5OHPl3HV3e9yqK6ViWNSmR0QJJQaznrTJXVa5LNCmgQctGxXAKdZTxCRk4EiY8wLIjKig4Lb4wlxpoq1HUcag/Z9UFnPC1t6UkdfOXc6T609aFtrGeD06fkAjMtO5YWvncmLW46wuCRPV1dTI0pv0kcxIyIJwG+Ab0Zx7s0isk5E1lVVVcW+cAOg3a01hcFQ39rJk+8doPxYk3/fzqPBQcEYb7fTbufOHstVJ0+ynTNzbCYllsFoOenJXH9aMSWFmTEouVKDpzfpo/44BBRZtif79nXLAhYAK3z998cDy30L+Kyz3sgYcx9wH0BpaemweLoG1RS6hkWxh71vPr2R17YfIyPZxUv/djbF+ensPNoTICaMSeVwfZvtmvRkFycW5bBgUjZtnR6qm9qZOTaTa8uKdWyJGhXiFRTWAjNFZBreYHAdcH33Qd/CPQXd2yKyAvhWYEAYrrRNIXbe3uWtLZ41sxCAleXV3LNiN6dMyeW17d61DZo7uvjx81u5/4ZSW03hxiVTufOlD233u3j+eJITE0gmgZ9esSBOP4VSQ0dcgoIxxi0itwGv4F2Y5wFjzFYR+QmwzhizPB7lGCyBXVK199HA+N/Xd/HrV70L//3iqoV84pTJfO2JDRxv7uCdgCmtX//wGCXffdG/nSBw7alFQUHBqcupUqNJvGoKGGNeBF4M2PfDEOeeG48y9UaXx/Dwu/to6XBz05klpCVH37gY2CX15Q8OU36siRvPmMaJRTkDXNLR4aUth/0BAbwD0sakJQU1DocytSCDnPRkSqfksm5/LQCXLhzPNJ3ETo1ycQsKw90rW4/wk394e9CKCLeeNyPqazsC0kf//U/vw2zXsSZe+NpZA1fIUaLD7eG7z22x7att6eSWx5znZszPSA4KFkt8PYluOnMaGw/WMTYrhf/4+PzYFFipYUSDQpT+0zL9wX+9sqNXQSGwTaHbgeMt/S7XaFF+rJG05EQm5aSxqaKO2pbOyBcBZ84o4JGbymhqd+NKEN7aWUVNcydXnDQRgEsWTuD96QWkp7hIcg2JznhKDSoNClFKTer7AyMwfdSttbOrz/ccTf6+8RD/9uRGMlMSeepLi1m1u2cuomtKJ1NZ1xbUhlBSmEFVYzvfung2IkJWahIASxdMCLr/mPSk2P4ASg0jGhSilNKPGexC1RTcHkNnl2fUf0Nt6+wiJTHBscunMYYfP+9N2zW1u/nR8q24LDOSLplewHlzxnLtvav40DcwbVlZMT+/coF2IVWqD0b306gX+lVTCDMuYbTXFp547wDzfvgyn/zDKvZWN/O1Jzbwhzd34/H10NpcUU+NpT1g7b5aVu+p8W+fPj2fMWlJ/PmmMj4ybxxnzijgGxfN0oCgVB9pTSFK/aophEgfAbR1dJGdOnLTF10eY/tmb3WssY3vPOttMF63v5bz/nsFAMs3VbLjSCO7q5rYXBE8SV236YUZjMtOBWBsVir33TBiJtZVatBoTSFKKQE1BdOLNTYDex9ZjeSawv+8tov5//EyP39xu+Pxe1bsdtwP8NyGQ2EDAvTMRaSUGjgaFKIUGAPaOqOf1C5UmwJAS8fIDArN7W7uem0nbZ0e7ntrD8eb2m3Hdx1t5DHfugbRmJybxufPmEair9aRmCBcU1oU4SqlVG9p+ihKbQHf6GtbOkhLTovq2nDpo5FaU1h/oNa2vfFgHRfMHQdARW0Ln/nTeyF7ZVldOHcs7W4PXz1/JmXT8vjGR2axuaKOcdmpTNfJ6JQacBoUohQYFOpaOpmYE11QCJc+ahuhNYU1lsZggA0HvEHB4zHc9vgGjjR4J6JLT3ZRlJvODsucRHdetZB9x5u54fSpTAr4jDNTElkyvQClVGxoUIhS4Df6utboplOA0OMUnO47nLR2dIWc7uO9vfag0F1z+MeWw2w8WAdAkku4/4ZSROD6+72L7Z0xI59lZcWxK7RSKiwNClEKCgpRjqiF4Anxwt13ODDG8IWH1/Hmziq+99G5fO4M+/pLbZ1d/gd/t00H62jpcPNLywR0N51ZwhkzvN/6b794Nlsq6vn20tkxL79SKjQNClFq7ehPUAhTUxiG6aP1B2p5/UPvtNQ/fn5bUFDYdLAuKGXW3NHFPSt2c6iuFYDc9CS+ct50//HeTBuilIod7X0UpcCHd21LdOkjj8eEnSo7sK1iONgfMGdT9/oQze1uOrs8rNpz3OkyHrf0NrqmtGhEj89QarjSmkIUjDFBaZ761uhqCuEamWF4dkkNXK2sqrGd8mNN3PLo+yS67APV5k/MZmtlA4BtptJZ43Sxe6WGolFdUzDGRDUIraPLQ+CX/booawrhUkcwPNsU9lQ127a3H27gK4+9T2O7m9qWTv8MpgkCt5w73ekWlBTqugVKDUWjNiis3VdD2c9f59p7V0dM4bR1BD/Yo526OVJf/OEYFPZWN9m2v/To+zS0uYPOW1SUw+klzqOOdcF7pYamURsUPnX/Gqoa23lvXw1v+BpNwZsff2nLYd7aWeXf5/Tgro8yKITreQTDc5zCvoA2hVCB79xZY8nPTGFcdoptf0FmCmPStD1BqaFoVAaFXUcbbbl+a/fJv66v4JbH1nPDA+/xzi7vHP0tHcHfgqNtaB5p6aO6lg7brKXhnDu7EIB5E7Jt+6dr6kipIWtUBoVn1lfYtq3f+n+8fKv/dfeSj04P7roBamhu7cUcSkPB3urmkMeuPGmSvwYwKSeNhZPGADA3ICho6kipoWvU9T5yd3l4bv0h274DNT3pkGZLOqd7v1ObQ31LJ8aYiPP2R2xTcKiFDGXhgsKpU/O44fQpvLD5MFeePIkE3+R18yZqTUGp4WLUBYW3y6s51mifsbP74W+MIUGw9TRqbnfT6tDQ3NHloand7V/mMZSRlj4KFxTmT8xmUVEOJxXn2vYH1hR0Ijulhq5Rlz5q7ehicq59krXK+lY63B5qmjuCup5uOlgX8sH9+vZjjvutIgaFYdbQ3L3kpZPZ453HHkzNt9cMphZoTUGpoWrUBYVLF07grdvP44kvLvbvM8Y7nfO+48HfgtcfqA0ZFJ5cG3k9gPaIXVKHT5tCW2eXv/E9UHJiAqlJzpPjuRKEr/jGK1w4dxxT89NjVkalVP+MuqAAkJAgnD49n8Ulef59+2ta2FvdEnTu+gN1trz/WTML/MtLrt5T40+n/GNzJWX/+Rq3Pr7eNiAuYpfUYZQ+entXtT9ABg4+y00Pn0b79tI5bPjBRdx/wym6frJSQ9ioDArdpuT1PNgO1rQEDcoCb3dVa4pnan4G580e699+cu0Bdh1t5JtPb+JYYzsvbD7M9sM9KZZwC+zA8EofvbL1iP/1xfPHU5jVM/4g1CA1q9yMZA0ISg1xozooFFvSGPuPt7DPoaZQ19JBi+XbfFqyi2VlPctA/vX9Cv79qY22NNGxxp65gUZKQ3OH28Pr24/6tz8ybxy/uvoEXAlCVmoi/++SOYNYOqXUQBl1vY+sivN6gsKBmhYO1bYGneMx0GSZwiEtycU5swoZn53KkYY2qps6qG6yD+ayTqsdcZzCMKgpGGP4zrNb/FN7jMtOYdHkHBIShFXfOZ/MlETSk0f1n5JSI8borilYg8Jx54ZmsM+ImpbsItGVwDWlk0Pe1zraOdI4hY4uD+4QgcMYw+NrDvDb13bS2Bb9+g0D7YGV+/irZcDfF88q8Y9BGJuVqgFBqRFkVAeFKZb00a5jjf5prLNSE8mxNJxaRy+n+XrYfLK0CGt6PMHy2jpZXqSGZoC2EIFjzd4avvvcFn772i4efndfxPvEgsdjeOCdvf7ta0onc9OZ08JcoZQazkZ1UBiTlkRWqvdbrnV8wrSCDFISez6aBoegUJSXzpUnTgKgMCuFL5xV4j+n3lJTiNSmAKFTSLssi9nvPBrcCB6NDreHJ987wPJNlVFNEx7o/QO1/tXSctKT+NkVC7WxWKkRbFTX+0WEKfnpfHCowbZ/an6GLQVkbSNItSxU//OrFnLVyZOZOyGLNy2zqlprCpHSRxC6W6p1yg2nSfmi8ez6Cu541juH05i0JM6ZVdir65dvrPS/vmTBBJITR/X3CKVGvFH/f7i1XaHb1IIMUhJ7Hv51rT0BIs0yQCs1ycWZMwvIz0yxpZtsbQrR1BRCBIUWW1DoW4N0d0AAuP0vm3p1rbvLw4tbDvu3L1s0sU9lUEoNH3ELCiKyVER2iEi5iNzhcPwbIrJNRDaLyOsiMiUe5SrOC55yYVpBOsmuno/GWlNIT3YetZuTnux4fn/SRy3tPbWD5gHopRRNT6e2zi5/mmnjwTr/EprjslMom5YX7lKl1AgQl6AgIi7g98AlwDxgmYjMCzhtA1BqjDkBeAb4VTzKNsVhyoWp+RmkJPV8NI2WLqmhpnLItQaF1uh7H0HoWoB1fMRAzKba7vawt7qZ7z63hfN/vcLWgAzwyKp9LPzRK9zwwHs0t7s5WNszbqN0Sp5/JLdSauSKV02hDCg3xuwxxnQATwKXW08wxvzLGNP9FFoNhO7zOYCc0kfTCjJsNQWrtJBBwdJbqbl3NYVQbQq2mkJ7/2sKHV0eLv7tWzy+5gB7qpq586XtNFve476399DZZXh7VzVfeWw9lXU9g/AmjEnt9/srpYa+eAWFScBBy3aFb18oNwEvxbREPoFBISc9iZz0ZFJCPPzTQqSPslKT/F1UG9vd/mAQTZfU6NoU+lZTSAz4dm+tuXR2GXb6ejh1eQwHa3oG7725s4r/eW2Xf3u8BgWlRoUh19AsIp8GSoH/CnH8ZhFZJyLrqqqqnE7plQljUm0Pzmm+aZ17W1NwJYht3eHudoWoGppDpY8GoKE5VLqrW/dU2FUBa0yAvewaFJQaHeLVJfUQUGTZnuzbZyMiFwLfA84xxgQ/pQBjzH3AfQClpaW973gfINGVwOTcNP9i9NN8c/9b2xSsQgUF8LYrdAeD+tYOCrNSompTCF1T6KkdtLs9dHnMgOf1//p+Bev31zIuO/xDX9NHSo0O8QoKa4GZIjINbzC4DrjeeoKInATcCyw1xkRevWYAFeWl+4NC9wIwKaFqCiHSR0BAt1RvcOhXm0JA7aClI/JKb1ZdHkNTe/i007r9tazbXxvxXuPHpEU8Ryk1/MUlfWSMcQO3Aa8A24GnjTFbReQnInKZ77T/AjKBv4jIRhFZHo+yAZwypWf5yFLfa6eagitBSHKF/qaeY0kf1fq6ckZVU4gifeS0HUmkgBAtERhrmSZbKTVyxW1EszHmReDFgH0/tLy+MF5lCfSFs0pwiTA2O4XTp3vXBXBqU0hLcoWd4sHeLdVbU7B2Zw2lJYr0EWDrKdStvrWTDN8kfYFCTaI3KSeNhtZOGqMMGoWZKSSFqDkppUaWUT3NRbfMlES+esFM2z6n3keRGm3tA9i8NYXuwV/h9LWm8NyGCr759CZKCjN54Wtn2kZhQ+iANCk3jYLMZDZV1Dsez01Psk3VoY3MSo0e+vUvBMeaQnL4jyvXoU3heJNje7mNU5uCx2MiBoWvP7UJj4HyY028svUogUIFhck5aUEBxGrJ9ALb9vgIjdBKqZFDg0IIKQ4Tv4XreQT2hua6lg6MMbZ5kEJx6n3U5g7eZ00nBTZgtzikgkKljybnpvGpxcUhy7Nkhn1pTe15pNTooemjEJxmA40cFOzzHzW0uaMbvOaQPnIawWytKewKmErbqbzh0kcfO2EiB4630NTh5t4399iOnxYwx5H2PFJq9NCaQgiONYUw3VEhsKbQSY2lPSHUYDhwrik4BYo9VU1sOFCLMYatlfb2AKeeSaEbmtNxJQhfvWAm37lkrq3cAFPy7ZMEBh5XSo1cGhRCSHbIuUeqKWSm9FS8Wjrc1DT3tCcUhunS6VhTcJjW4r//uZMr736XR1fvZ2ulfQ0Ip2kwGsLUFKy6FxrqFtjTKPC4Umrk0qAQQl9qChmWoNDU7uZ4U09NYWx2mKDgUFMIN9fRD/6+NcqagvM9AtsITp3aky7Ky/CmwL56/gzA28h84dxxIcuilBpZ9CtgCE6D1yJ1SbWutdDS0WVLH43LCt1Y6xwUwg9U2xZQU3CqbTiljy6YMzbo5/j2xXN4ddtRWjq6+PUnFwHwjYtmcfH88UwtyIj4cyulRg4NCiGEGrwWTmZgTaE5uppCW5QNzbbjAddY0007jjRyy2Pvs6eq2b/vF1ctZMbYTE6YnBN0r/FjUlnz3Qtoancz1he8RIQFk8aELYNSauTRoBCC0+C1SEEhPdnaptBlTx8FtCm4EoQuj7dnkmNDc2fvpqiw1ix++PcPbAEBID8zhdKpeYGX2cpuLb9SanTSNoUQnGoKoZbi9F+TmOC/rstjOFzfsz7B2ID0kbVW4RQUeruoTnf6qMPtYc3emqDj2lislIqGPilCcGxTiBAUANJTXHS0eAeWWZezDOx9lJmSSL1vfqS2Tg8ejyHBMi12NOspW204UMfV97xLZV2r43ENCkqpaOiTIoS+tCkAZCQn+tdUOHC8Jyh09+rp5koQUpMSaOv0BpA2d5ctfePUJTWcIw1tHGloC3k8uxdTbiulRi9NH4WQ6lBTiCoopPScYx0nkJ9pDwoi9vsF1gx6W1OIRGsKSqloaFAIIdnl0NAcTfooRGNtfoY9fSQEBIXO0L2JBoK1DUMppULRoBBCX8YpgPPDNy3JFRRQRMTWRhE4U2pf12TuFtiG4bTeglJKBdInRQh96X0U6pzA9gRwqCl0BM562veg8NgXTuOlfzvLv71QxxsopaKkOYUQnGoK0bQpONUUAtsTABB7AAlMH4Vajc1qUk4ahxx6Gy0uyceVIPz582W8svUInzptSsR7KaUUaFAIyammEE36KD0l+Jz8EDWF1DBtCk7rIwQqykvjaEMbbk/P9NxZKYm4fF1bz55VyNmzCiPeRymlumn6KIREV4L/4dotmobmDIeG5ryM4CkuRCQgfWQPAtG0KYzPTg0qU3aadj1VSvWdBoUwAmsL0XVJjS59JNiDzJcfXc9l//cOxxq9Yw3CzZLabfyYtKAgNEaDglKqHzQohBHYrhBNUIi6oVmC77e5op6n1x4EoCmKhuYxaUlB76dBQSnVHxoUwgiqKUSRPnJqaHbufSSObRQHalqob+mkusm7QE+SS/j5lQvJSknk6pMn285NcolD+kibiZRSfadPkDCsNQUR54V3AqU7pY9C1BScahVHGtptC+jMHp/F9acVc92pRSQkCC99cNjf3nBScQ7/3HrUdr3WFJRS/aE1hTCsNYW0JBciEuZsr0yn3keZkRuaux2tb+MDS1CYP8E7xqB7srwHbjyVBZOyufnsEk6ZkhdUU9CgoJTqD60phJFiWac5mvYEcJ7mIlSXVKd01JGGNtv6ywsmZduOLy7J5x9f7RmYpm0KSqmBpDWFMJIt6aJol6R07pLqnD5yumd9ayfr9tX6t+dHGI0cGIS0S6pSqj80KIRhbUOIppEZ7LOkdt/Dqe3AqfdRt+5RygkCc8dnO57TTWsKSqmBpEEhDGtNIdr0UeA4hfyMZMe2CCG451Cg6YWZEc8JDApaU1BK9YcGhTBsbQpR1xQCgoJDIzP4agoR7rkgionsAtNHWlNQSvWHBoUwUvpQU0gPOM+pPQFC9z6yOrEoJ/L7afpIKTWANCiE0ZegkJAgtge1U88jCJ4628nikvyI7xc0eE2X3VRK9UPcgoKILBWRHSJSLiJ3OBxPEZGnfMfXiMjUeJUtlOQ+NDSDPaUTuqYQ/p55GcnMHJsZxXtpTUEpNXDiEhRExAX8HrgEmAcsE5F5AafdBNQaY2YAdwG/jEfZwknpQ5dUsPdAynNaS4HINYWyqXn+AWvhGGPfTo5i1LVSSoUSrydIGVBujNljjOkAngQuDzjncuBh3+tngAskmiHEMdSX3kdgH6sQMn0kznMfdVtckhfVe3V2eSKfpJRSUYpXUJgEHLRsV/j2OZ5jjHED9UDkpHoM2XsfRf9RWWsK+Q5rKYA37RMufXRaFO0JQNCaD0op1R/DLtcgIjeLyDoRWVdVVRXT97I29C6ZXhD1daf7zk1PdnHylFz//l9dfYL/9Y8vm09Gsss/jcXpJfmUTfPWDuaMz2L2uKyo3uvShRP87QifO2Nq1GVUSiknYgKT0rF4E5HTgR8ZYy72bX8HwBhzp+WcV3znrBKRROAIUGjCFLC0tNSsW7cuZuU2xrDxYB3g7R4abTbL3eXhrV1VzCjMojg/3Xa/VbuPU5CVwizfQ7++pZNVe46zZIY3AL2zq5qyaXkUhBjf4KSyrpVtlQ2cNavAVrtRSqkQQj7M4hUUEoGdwAXAIWAtcL0xZqvlnFuBhcaYL4vIdcBVxphrwt031kFBKaVGqJBBIS6zpBpj3CJyG/AK4AIeMMZsFZGfAOuMMcuBPwGPiEg5UANcF4+yKaWU6hGXmkKsaE1BKaX6JGRNYdg1NCullIodDQpKKaX8NCgopZTy06CglFLKT4OCUkopv2Hd+0hEqoD9fby8AKgewOIMpKFaNi1X72i5em+olm2klavaGLPU6cCwDgr9ISLrjDGlg10OJ0O1bFqu3tFy9d5QLdtoKpemj5RSSvlpUFBKKeU3moPCfYNdgDCGatm0XL2j5eq9oVq2UVOuUdumoJRSKthorikopZQKMKKDgoh8UkS2iohHREK20IvIUhHZISLlInKHZf80EVnj2/+UiDivrdn7cuWJyKsissv331yHc84TkY2Wf20icoXv2EMistdy7MSBKFe0ZfOd12V5/+WW/YP5mZ0oIqt8v/PNInKt5diAfmah/mYsx1N8P3+57/OYajn2Hd/+HSJycX/K0YdyfUNEtvk+n9dFZIrlmOPvNE7lulFEqizv/wXLsc/6fu+7ROSzcS7XXZYy7RSROsuxWH5eD4jIMRH5IMRxEZHf+cq9WUROthzr3+dljBmx/4C5wGxgBVAa4hwXsBsoAZKBTcA837Gnget8r/8A3DJA5foVcIfv9R3ALyOcn4d3OvF03/ZDwCdi9JlFVTagKcT+QfvMgFnATN/ricBhIGegP7NwfzOWc74C/MH3+jrgKd/reb7zU4Bpvvu44liu8yx/R7d0lyvc7zRO5boR+D+Ha/OAPb7/5vpe58arXAHnfxXvtP8x/bx89z4bOBn4IMTxS4GX8M52uhhYM1Cf14iuKRhjthtjdkQ4rQwoN8bsMcZ0AE8Cl4uIAOcDz/jOexi4YoCKdrnvftHe9xPAS8aYlgF6/3B6Wza/wf7MjDE7jTG7fK8rgWNA4QC9v5Xj30yY8j4DXOD7fC4HnjTGtBtj9gLlvvvFpVzGmH9Z/o5WA5MH6L37Va4wLgZeNcbUGGNqgVcBx0FXcSjXMuCJAXrvsIwxb+H9IhjK5cCfjddqIEdEJjAAn9eIDgpRmgQctGxX+PblA3XGGHfA/oEwzhhz2Pf6CDAuwvnXEfzH+J++auNdIhL92p0DV7ZU8a6Vvbo7rcUQ+sxEpAzvt7/dlt0D9ZmF+ptxPMf3edTj/XyiuTaW5bK6Ce+3zW5Ov9N4lutq3+/nGREp6uW1sSwXvjTbNOANy+5YfV7RCFX2fn9ecVl5LZZE5DVgvMOh7xlj/h7v8nQLVy7rhjHGiEjILmC+6L8Q76p13b6D98GYjLdL2v8DfhLnsk0xxhwSkRLgDRHZgvfB12cD/Jk9AnzWGOPx7e7XZzbSiMingVLgHMvuoN+pMWa38x0G3PPAE8aYdhH5Et5a1vlxeu9oXAc8Y4zpsuwbzM8rZoZ9UDDGXNjPWxwCiizbk337juOtkiX6vul17+93uUTkqIhMMMYc9j3AjoW51TXAc8aYTsu9u78xt4vIg8C3oi3XQJXNGHPI9989IrICOAn4K4P8mYlINvAC3i8Fqy337tdnFiDU34zTORXiXaN8DN6/qWiujWW5EJEL8Qbac4wx7d37Q/xOB+IhF7Fcxpjjls0/4m1D6r723IBrVwxAmaIql8V1wK3WHTH8vKIRquz9/rw0fQRrgZni7TWTjPeXv9x4W23+hTefD/BZYKBqHst994vmvkF5TN9DsTuHfwXg2EMhVmUTkdzu9IuIFABnANsG+zPz/f6ew5trfSbg2EB+Zo5/M2HK+wngDd/nsxy4Try9k6YBM4H3+lGWXpVLRE4C7gUuM8Ycs+x3/J3GsVwTLJuXAdt9r18BPuIrXy7wEey15piWy1e2OXgbbVdZ9sXy84rGcuAGXy+kxUC974tP/z+vWLWeD4V/wJV4c2rtwFHgFd/+icCLlvMuBXbijfLfs+wvwfs/bDnwFyBlgMqVD7wO7AJeA/J8+0uBP1rOm4o38icEXP8GsAXvg+1RIHMAP7OIZQOW+N5/k++/Nw2Fzwz4NNAJbLT8OzEWn5nT3wzedNRlvtepvp+/3Pd5lFiu/Z7vuh3AJQP8Nx+pXK/5/l/o/nyWR/qdxqlcdwJbfe//L2CO5drP+z7HcuBz8SyXb/tHwC8Crov15/UE3t5znXifYTcBXwa+7DsuwO995d6CpXdlfz8vHdGslFLKT9NHSiml/DQoKKWU8tOgoJRSyk+DglJKKT8NCkoppfw0KCjVSyKyzzcATKkRR4OCUnEm3mm8fzbY5VDKiQYFpZRSfhoUlOojESkT76I+dSJyWET+zzddQvciKHeJd6GUBhHZIiILRORm4FPAt0WkSUSeH9yfQim7YT8hnlKDqAv4OrAO78RjL+FdXOe3eOecORvvwj/1wBy804rfJyJLgApjzPcHo9BKhaM1BaX6yBjzvjFmtTHGbYzZh3eiue6pqDuBLLzBQIx3wafDIW6l1JChQUGpPhKRWSLyDxE5IiINwM+BAgBjzBvA/+GdtOyYiNznm9ZbqSFNg4JSfXcP8CHedaGzge/inb0SAGPM74wxp+Bdl3kWcHv3oXgXVKloaVBQqu+ygAagyTfn/i3dB0TkVBE5TUSSgGagDeheBe4o3inGlRpyNCgo1XffAq4HGoH7gacsx7J9+2qB/XhXXfsv37E/AfN8vZb+FrfSKhUFXU9BKaWUn9YUlFJK+WlQUEop5adBQSmllJ8GBaWUUn4aFJRSSvlpUFBKKeWnQUEppZSfBgWllFJ+GhSUUkr5/X+tl9gtRLcn7wAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"df_base[\"bin\"] = pd.cut(df_base[\"last\"], 200, labels=False) / 100 - 1\n",
"df_base.groupby(\"bin\")[\"next\"].mean().plot(xlabel=\"last\", ylabel=\"next\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now compute the difference at the cutoffs to get an estimate for the treatment effect."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Treatment Effect: 0.096\n"
]
}
],
"source": [
"h = 0.05\n",
"df_subset = df_base[df_base[\"last\"].between(-h, h)]\n",
"stat = np.abs(df_subset.groupby(\"incumbent_last\")[\"next\"].mean().diff()[1])\n",
"print(f\"Treatment Effect: {stat:5.3f}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"How does the effect depend on the size subset under consideration?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Regression approach"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we turn to an explicit model of the conditional mean. We first set up explicit models on both sides of the cutoff and then aggreagte the model into single regression estimations."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"\n",
" Republican\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: next R-squared: 0.271\n",
"Model: OLS Adj. R-squared: 0.270\n",
"Method: Least Squares F-statistic: 339.2\n",
"Date: Wed, 30 Jun 2021 Prob (F-statistic): 3.05e-187\n",
"Time: 12:05:34 Log-Likelihood: 1749.4\n",
"No. Observations: 2740 AIC: -3491.\n",
"Df Residuals: 2736 BIC: -3467.\n",
"Df Model: 3 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"Intercept 0.4278 0.007 57.880 0.000 0.413 0.442\n",
"last -0.0971 0.077 -1.264 0.206 -0.248 0.054\n",
"last_2 -1.7177 0.205 -8.359 0.000 -2.121 -1.315\n",
"last_3 -1.4636 0.142 -10.338 0.000 -1.741 -1.186\n",
"==============================================================================\n",
"Omnibus: 203.681 Durbin-Watson: 1.866\n",
"Prob(Omnibus): 0.000 Jarque-Bera (JB): 1087.416\n",
"Skew: -0.022 Prob(JB): 7.42e-237\n",
"Kurtosis: 6.086 Cond. No. 113.\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\n",
"\n",
" Democratic\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: next R-squared: 0.379\n",
"Model: OLS Adj. R-squared: 0.379\n",
"Method: Least Squares F-statistic: 776.5\n",
"Date: Wed, 30 Jun 2021 Prob (F-statistic): 0.00\n",
"Time: 12:05:34 Log-Likelihood: 2055.2\n",
"No. Observations: 3818 AIC: -4102.\n",
"Df Residuals: 3814 BIC: -4077.\n",
"Df Model: 3 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"Intercept 0.5393 0.007 71.995 0.000 0.525 0.554\n",
"last 0.3553 0.071 4.998 0.000 0.216 0.495\n",
"last_2 0.1932 0.174 1.107 0.268 -0.149 0.535\n",
"last_3 -0.2111 0.114 -1.856 0.064 -0.434 0.012\n",
"==============================================================================\n",
"Omnibus: 439.976 Durbin-Watson: 2.136\n",
"Prob(Omnibus): 0.000 Jarque-Bera (JB): 1993.314\n",
"Skew: -0.477 Prob(JB): 0.00\n",
"Kurtosis: 6.409 Cond. No. 114.\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
]
}
],
"source": [
"def fit_regression(incumbent, df, level=4):\n",
"\n",
" df_incumbent = df[df[\"incumbent_last\"] == incumbent].copy()\n",
"\n",
" formula = \"next ~ last\"\n",
" for level in range(2, level + 1):\n",
" label = \"last_{:}\".format(level)\n",
" formula += f\" + {label}\"\n",
"\n",
" rslt = smf.ols(formula=formula, data=df_incumbent).fit()\n",
" return rslt\n",
"\n",
"\n",
"rslt = dict()\n",
"for incumbent in [\"republican\", \"democratic\"]:\n",
" rslt = fit_regression(incumbent, df_base, level=3)\n",
" title = \"\\n\\n {:}\\n\".format(incumbent.capitalize())\n",
" print(title, rslt.summary())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"How does the predictions look like?"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"dfs = list()\n",
"\n",
"for incumbent in [\"republican\", \"democratic\"]:\n",
" rslt = fit_regression(incumbent, df_base, level=4)\n",
"\n",
" # For our predictions, we need to set up a grid for the evaluation.\n",
" if incumbent == \"republican\":\n",
" grid = np.linspace(-0.5, 0.0, 51)\n",
" else:\n",
" grid = np.linspace(+0.0, 0.5, 51)\n",
"\n",
" df_grid = pd.DataFrame(grid, columns=[\"last\"])\n",
"\n",
" for level in range(2, 5):\n",
" label = \"last_{:}\".format(level)\n",
" df_grid.loc[:, label] = df_grid[\"last\"] ** level\n",
"\n",
" tmp = pd.DataFrame(rslt.predict(df_grid), columns=[\"prediction\"])\n",
" tmp.index = df_grid[\"last\"]\n",
" dfs.append(tmp)\n",
"\n",
"rslts = pd.concat(dfs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's have a look at the estimated conditional mean fuctions."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEJCAYAAACE39xMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlL0lEQVR4nO3deXiU1d3/8fc3CRAIARTCmoQEkrDITgARFdksVouKVFFsrV2sAq51a+3Tp7VPW+ziVqktdUGtS5WfC1qLioDIKhFBdhJCSMK+QwJkPb8/EkImCWSQSSYz83ldF5e5zznMfG8In9yeOfe5zTmHiIgEvjB/FyAiIr6hQBcRCRIKdBGRIKFAFxEJEgp0EZEgEeGvNx47dqybM2eOv95eRCRQ2ek6/HaFvm/fPn+9tYhIUNKUi4hIkFCgi4gECQW6iEiQUKCLiAQJBbqISJBQoIuIBAm/rUMXEQk1JaWOd77ajgHXDYz1+esr0EVE6lhpqWPOul08/slmMvbk0TqqMWN7tSeqiW8jWFMudWjBggVcddVVAMyePZtp06adduyhQ4f429/+VnG8Y8cOJkyYUOc1ikjdcc4xf9Mexk1fxORXV5KxJw+A/fmFvLx0m8/fT1fo30BJSQnh4eFn9XvGjRvHuHHjTtt/MtAnT54MQMeOHZk1a9Y51Ski/rMscz9//mgTadsOerRHN4ngx5d04eYL433+nrpCryIrK4vu3bszadIkevTowYQJEzh27BgJCQk89NBDDBgwgLfeeouPP/6YoUOHMmDAAL773e+Sl1f2k3fOnDl0796dAQMG8Pbbb1e87syZM5k6dSoAu3fv5tprr6Vv37707duXJUuW8PDDD7Nlyxb69evHAw88QFZWFr169QLgxIkT3HrrrfTu3Zv+/fszf/78itccP348Y8eOJTk5mQcffLCe/7REpKpVOYf43vPLmThjmUeYRzYK46fDu7DwwRHcPTqZ6MhGPn/vBnuFnvDwf+rstbOmXXnG/k2bNvH8888zbNgwfvjDH1ZMhbRu3ZqVK1eyb98+xo8fz9y5c4mKiuKxxx7j8ccf58EHH+QnP/kJ8+bNIykpiRtuuKHG17/rrrsYPnw477zzDiUlJeTl5TFt2jTWrl3LqlWrymrMyqoYP336dMyMNWvWsHHjRi6//HI2b94MwKpVq/jqq69o0qQJ3bp148477yQuLu7c/5BE5Kxs3HWEv3y8mU/W7/ZobxRu3Dg4nqkjkmjbIrJOa2iwge5PcXFxDBs2DICbb76Zp59+GqAioJctW8b69esrxhQWFjJ06FA2btxIYmIiycnJFb93xowZ1V5/3rx5vPzyywCEh4fTsmVLDh48WG3cSYsWLeLOO+8EoHv37nTu3Lki0EeNGkXLli0B6NmzJ9u2bVOgi9SjLXvzeHJuOh98vYPKj2gOMxg/IJa7RyUTd36zeqlFgV4DM6vxOCoqCij7oGPMmDG8/vrrHuNOXl3XpyZNmlR8HR4eTnFxcb3XIBKKcg4c46lP03l7ZS6lzrPvqj4duGd0Ckltm9drTQ020GubFqlL2dnZLF26lKFDh/Laa69x8cUX89VXX1X0X3jhhUyZMoWMjAySkpLIz89n+/btdO/enaysLLZs2ULXrl2rBf5Jo0aN4tlnn+Wee+6pmHKJjo7m6NGjNY6/5JJLePXVVxk5ciSbN28mOzubbt26sXLlyjo5fxE5vZ2Hj/PMvAz+vSKH4ipJPrpHO+4bk0LPji38Ups+FK1Bt27dmD59Oj169ODgwYPccccdHv0xMTHMnDmTG2+8kT59+lRMt0RGRjJjxgyuvPJKBgwYQNu2bWt8/aeeeor58+fTu3dvBg4cyPr162ndujXDhg2jV69ePPDAAx7jJ0+eTGlpKb179+aGG25g5syZHlfmIlL39hw9wa9nr2P4nxbw6vJsjzC/JLkN704ZxnO3pPotzAHMOVf7qDqQmprq0tLS/PLeZ5KVlcVVV13F2rVr/V2KiDQA+/MKmLEwk5eWZnGiqNSjb3DC+fzs8hSGdGldnyWd9olFDXbKRUTEnw4fK2LG51uYuTiL/MISj76+ca342ZgULkluU+0zN39SoFeRkJCgq3OREHbkRBEvLNrK859v5WiB5yKDCzq24L4xKYzs3rZBBflJXgW6mY0FngLCgeecc9Oq9D8BjCg/bAa0dc618mGdIiJ1Kq+gmJmLtzJjYSZHTngGebd20dw7JpnLe7YnLKzhBflJtQa6mYUD04ExQC6wwsxmO+fWnxzjnLu30vg7gf51UKuIiM/lFxTz0tIsZizM5NCxIo++LjFR3DM6hat6d2jQQX6SN1fog4EM51wmgJm9AVwNrD/N+BuB//VNeSIideNYYTEvL93GjIWZHMgv9OhLaN2Mu0cnM65vJ8IDIMhP8ibQOwE5lY5zgSE1DTSzzkAiMO80/bcBtwHEx/t+YxoRkdocLyzhX8u28Y+FW9iX5xnkcec35c6RyYzv34mI8MBb1e3rD0UnArOccyU1dTrnZgAzoGzZoo/fW0TktI4XlvDq8m38/bPqQd6pVVPuHJnEdQNjaRSAQX6SN4G+Hai8OUhseVtNJgJTzrUoERFfORXkmezLK/Do69gykikjk/juwDgaRwRukJ/kTaCvAJLNLJGyIJ8I3FR1kJl1B84Dlvq0QhGRb+BMQd6hZSSTRyRxfWosTSLO7tkGDVmtge6cKzazqcBHlC1bfME5t87MHgXSnHOzy4dOBN5w/rr1VESEMwd5+xaRTBnRlesHxQVVkJ+kW/9FJCgcKyzm1WXZ/GNh0Ae5bv0XkeB0rLCYfy0rW35Y9cPO9i0imTyiKzcER5DXSoEuIgEpv6BsHfk/P6++jjzIrsi9pkAXkYCSV1DMS0uyeO7zTA5WubOzY8tI7gjCDzu9pUAXkYBw5EQRLy3O4vnFW6vdot+pVVMmj+jKhIGhGeQnKdBFpEE7fKyIFxZv5cXFW6ttmhV7XlOmjkhi/IDYoFhHfq4U6CLSIB3ML+T5RVt5aUlWtW1sO7duxpQRSVzbv1NA39npawp0EWlQ9uUV8NznW3llafUHSyS2iWLKiCSu6dcxIPdaqWsKdBFpEPYcOcGMhZn8a/m2ao966xoTxV2jkrmqT8eA2v2wvinQRcSvdh4+zj8+y+T1L7IpKPYM8m7torlzVBJX9OqgIPeCAl1E/CL34DGeXbCFt9JyKSzxDPKeHVpw16ikBv+EoIZGgS4i9SprXz5/W5DB2yu3U1zqufVIn9iW3DUymVE9GuYzOxs6BbqI1Iste/OYPi+D91bvoKRKkA+Ib8Vdo5IZnhKjID8HCnQRqVObdh3lmfkZfPD1DqruBTgk8XzuGpXMRV1bK8h9QIEuInVi7fbDPDMvgznrdlXruzipDXeOTGJIl9Z+qCx4KdBFxKdW5xzir/PSmbthT7W+Ed1imDoymYGdz/NDZcFPgS4iPvHltoM8/Wk6n23eW61vTM923DkyiT6xreq/sBCiQBeRc7I8cz9Pz0tnccb+an1X9GrPnSOT6dmxhR8qCz0KdBE5a845lmbu56m56SzfesCjzwy+06cjU0cmkdIu2k8VhiYFuoh4zTnH4oz9PPXpZlZkHfToCw8zru7XkSkjkuga09xPFYY2BbqI1Mo5x2eb9/L0p+mszD7k0RcRZlzbvxNTRiSR0CbKPwUKoEAXkTNwzrFg016e/DSd1TmHPPoahRsTBsYx+bKuxJ3fzD8FigcFuohU45xj/qY9PDk3na9zD3v0NQ4P4/pBsdxxWRKdWjX1U4VSEwW6iFQ4Y5BHhHHjoDhuv6wrHVoqyBsiBbqI1BrkNw2O547LutKuRaSfKhRvKNBFQljFHPnczayuEuRNIsK4aUg8tw9XkAcKBbpICHLOsTB9H098splVVT7sbBIRxqQhnbl9eBfaKsgDigJdJIQ451iyZT+Pf7KZL7d5riNvHBHGpCHx3DG8q4I8QCnQRULEF1sP8JePN1W7s1Nz5MFDgS4S5FblHOIvH2/i8/R9Hu2Nwo2Jg+KZPEKrVoKFAl0kSK3fcYTHP9lUbRvbiDDju6lxTB2pdeTBRoEuEmS27M3jiU8288HXOz3awwyu7R/L3aOSiW+tOzuDkVeBbmZjgaeAcOA559y0GsZcD/wacMBq59xNPqxTRGqx49BxnpqbzqyVuR7P7Dy5++Hdo5O1aVaQqzXQzSwcmA6MAXKBFWY22zm3vtKYZODnwDDn3EEza1tXBYuIpwP5hUyfn8ErS7dRWFLq0TemZzt+dnkK3dtrP/JQ4M0V+mAgwzmXCWBmbwBXA+srjfkJMN05dxDAOVf92VMi4lP5BcU89/lW/vl5JnkFxR59w5Jac//l3egfr0e9hRJvAr0TkFPpOBcYUmVMCoCZLaZsWubXzrk5VV/IzG4DbgOIj4//JvWKhLzC4lJe/yKbv85LZ19eoUdf37hWPPStblyU1MZP1Yk/+epD0QggGbgMiAUWmllv59yhyoOcczOAGQCpqakOEfFaaanjgzU7+fNHm8g+cMyjL6ltc+6/vBvfuqAdZuanCsXfvAn07UBcpePY8rbKcoHlzrkiYKuZbaYs4Ff4pEqRELd0y37+8N8N1TbO6tgyknvHpDB+QCzhYQryUOdNoK8Aks0skbIgnwhUXcHyLnAj8KKZtaFsCibTh3WKhKT03Uf5w383Mm+j58dSLZs2YuqIJL43tDORjcL9VJ00NLUGunOu2MymAh9RNj/+gnNunZk9CqQ552aX911uZuuBEuAB51z1R4CLiFf2Hi3gibmbeeOLbCqtQKRJRBi3Dkvkjsu60rJpI/8VKA2SOeefqezU1FSXlpbml/cWaahOFJXw/KKt/G1+BvmFJRXtZnDdgFjuG5NCR93dGepOO7emO0VFGgDnHLNX7+CPczax/dBxj76Lk9rwi2/3oGdHrSWXM1Ogi/jZ6pxD/Ob9dazMPuTRnty2Ob+4sgeXpcRo5Yp4RYEu4id7jpxg2pyNvL3Sc9HY+VGNuW9MChMHxRERHuan6iQQKdBF6llBcQkvLs7ir5+me8yTNwo3bh2WyNSRSbSI1AeecvYU6CL1aMGmPfzm/fVs3Zfv0T6mZzse+XYPEtpE+akyCQYKdJF6kHPgGL/9YD0fr9/t0Z7ctjm/+k5PLkmO8VNlEkwU6CJ1qKC4hBmfZfLM/AwKik/thBgdGcG9o1P43tDONNI8ufiIAl2kjnyevpdfvbeu2vTK9amxPDi2O22aN/FTZRKsFOgiPrbn6Al++8EG3l+9w6O9V6cWPHp1LwZoS1upIwp0ER8pLXW8+kU2f5yzkaMnTu1PHh0ZwQPf6sakIZ21gZbUKQW6iA9s2nWUh9/+mq+q3Bx0Tb+O/OLKHrSNjvRPYRJSFOgi5+BEUQnT52fw98+2UFRyal+khNbN+L9renNxsh40IfVHgS7yDX257QAPzvqaLXtPfejZKNy4Y3hXJo9I0ra2Uu8U6CJnKb+gmD99tImXlmZRebPSgZ3P4w/je5PSLtp/xUlIU6CLnIUlW/bx4KyvyT14akfEqMbhPHxFdyYN6UyYPvQUP1Kgi3ghv6CYaf/dyCvLtnm0X9Ytht9d25tO2qNcGgAFukgtlmfu5/5Zq8k5cOqqvGXTRvzvd3pybf9O2tpWGgwFushpnCgq4c8fbeL5xVs95spH92jL76/tTdsWWoooDYsCXaQG63Yc5p43VpG+J6+iLToygt+Mu0BX5dJgKdBFKikpdcxYmMnjn2zyWFd+aUoMj13Xmw4tNVcuDZcCXaSSe/+9itmV9mBp2iicR67swaQh8boqlwZPgS5SrriklPe/PhXmfeNa8cT1fekS09yPVYl4T4EuUq7UUfHhZ0SYMev2odqrXAKKvltFamCGwlwCjr5jRUSChAJdRCRIKNBFRIKEAl1EJEgo0EVEgoQCXUQkSHgV6GY21sw2mVmGmT1cQ/8PzGyvma0q//Vj35cqIiJnUuuNRWYWDkwHxgC5wAozm+2cW19l6L+dc1ProEYREfGCN1fog4EM51ymc64QeAO4um7LEhGRs+VNoHcCciod55a3VXWdmX1tZrPMLM4n1YmIiNd89aHo+0CCc64P8AnwUk2DzOw2M0szs7S9e/f66K1FRAS8C/TtQOUr7tjytgrOuf3OuYLyw+eAgTW9kHNuhnMu1TmXGhMT803qFRGR0/Am0FcAyWaWaGaNgYnA7MoDzKxDpcNxwAbflSgiIt6odZWLc67YzKYCHwHhwAvOuXVm9iiQ5pybDdxlZuOAYuAA8IM6rFlERGrg1X7ozrkPgQ+rtP2q0tc/B37u29JERORs6E5REZEgoUAXEQkSCnQRkSChQBcRCRIKdBGRIKFAFxEJEgp0EZEgoUAXEQkSCnQRkSChQBcRCRIKdBGRIKFAFxEJEgp0EWDv0QLufuOrimPD/FiNyDfj1W6LIsHKOcfs1Tv49ex1HDxWVNE+KPE8P1Yl8s0o0CVkZe8/xi/fW8vCzZ6PQ7xxcBy/+HYPP1Ul8s0p0CXkFBSX8NznW/nrvHROFJVWtHdsGcm06/pwaYoejyiBSYEuIWXBpj385v31bN2XX9FmBrcMTeBnl6cQHdnIj9WJnBsFuoSEzL15/P7DjczdsNujvWeHFvxhfG/6xrXyT2EiPqRAl6B26FghT32azitLt1Fc6iraoyMj+NmYFG6+sDMR4VrsJcFBgS5B6XhhCTOXZPHsggyOnCj26PvuwFgeuqI7bZo38VN1InVDgS5BpaiklFlf5vLk3M3sPlLg0Tc48Xx+eWUP+sS28k9xInVMgS5BobiklHdX7eDpT9PJPnDMoy+xTRQPje3Gty5oj5luGJLgpUCXgFZUUsp7q3YwfX6Gx8oVgJjoJtwzOpnrU+NopHlyCQEKdAlIBcUlvL1yO39bkEHOgeMefa2aNeKnl3bllos606yxvsUldOi7XQJKXkExry3fxnOfb2XPUc858ujICH5ySRduHZag9eQSkhToEhB2HznBzCVZvLpsW7VVK62aNeLHFyfy/YsSaKEglxCmQJcGbf2OIzy/aCuzV2+nqMR59LWNbsKPLk5k0oWdad5E38oi+lcgDU5JqWPexj28sGgrSzP3V+tPaN2Mnw7vyvgBnWgSEe6HCkUaJgW6NBiHjxfxVloOLy/dVm3pIcDAzudx26VdGN2jHeFhWn4oUpUCXfxuw84jvLx0G+9+tZ3jRSUefeFhxthe7fnhsEQGdtYe5SJnokAXvygoLmHO2l28uiybL7IOVOtv2bQREwfF8f2LEujUqqkfKhQJPF4FupmNBZ4CwoHnnHPTTjPuOmAWMMg5l+azKiVo5Bw4xutfZPNmWg778gqr9XdvH80tFyVwTb9ONG2s+XGRs1FroJtZODAdGAPkAivMbLZzbn2VcdHA3cDyuihUAldxSSkLNu3l1eXbWLB5L85zsQoRYca3erXnlqEJDEo4T7fni3xD3lyhDwYynHOZAGb2BnA1sL7KuN8CjwEP+LRCCVi7Dp/gjRXZ/HtFDjsPn6jW375FJDcNiWfioDjatoj0Q4UiwcWbQO8E5FQ6zgWGVB5gZgOAOOfcf8zstIFuZrcBtwHEx8effbXS4JWUOham7+W15dnM27iHklJXbcylKTHcNDie0T3aai9yER865w9FzSwMeBz4QW1jnXMzgBkAqamp1f+lS8Dac+QEb6bl8PoXOWw/dLxaf+uoxkxIjWXS4M7Et27mhwpFgp83gb4diKt0HFvedlI00AtYUD732R6YbWbj9MFocCstdSzZsp9Xl2/jk/W7PZ4IdNKFXc5n0pDOXH5BO90EJFLHvAn0FUCymSVSFuQTgZtOdjrnDgNtTh6b2QLgfoV58DqYX8hbX+bw2vJssvZXvwGoVbNGTBgQy41D4uka09wPFYqEploD3TlXbGZTgY8oW7b4gnNunZk9CqQ552bXdZHif845VmYf5F/LsvnPmp0UFpdWGzMo4TwmDenM2F7tiWykq3GR+mau6hqyepKamurS0nQR39DlFxTz3qodvLJsGxt2HqnWHx0ZwXUDYrlpSDwp7aL9UKFIyDntul7dKSo1ytybxyvLtjErLZejBcXV+vvGtmTSkM58p29H3QAk0kAo0KVCaaljweY9zFyyjYWb91brj2wUxtV9O3HzhZ3pHdvSDxWKyJko0IW8gmLeSsvhpSVZNX7Imdgmipsv7MyEgbG0bKoHSIg0VAr0EJZz4Bgzl2Tx5oqcatMqZjCqe1tuuSiBYV3bEKbtakUaPAV6CPoq+yD//DyTOWt3UXXpeHRkBBMHxfG9CxN0A5BIgFGgh4jSUsfcDbv55+eZrMg6WK2/S0wUt16UwPgBsUTpcW4iAUn/coNcYXEp767azoyFmWTsyavWf3FSG350cSLDU2I0rSIS4BToQepYYTGvf5HDPxdmsuuI506HjcKNcX078eNLEunRoYWfKhQRX1OgB5kjJ4p4eUkWzy/aysFjRR590U0iuOnCeG69KJH2LbVdrUiwUaAHicPHinhh8VZeXLyVIyc8V6y0ad6EH12cyKQL42kRqWWHIsFKgR7gDh8v4vlFW3lx0dZqSw9jz2vKT4d35bsDY7W3ikgIUKAHqKMninhhURbPLcrkaJUr8sQ2UUwZkcTV/TrSSA+QEAkZCvQAc7ywhJeXZvH3z7ZUmyPvGhPFXaOSuapPR8K1YkUk5CjQA0RRSSlvpuXw1Nx09hwt8OjrEhPF3QpykZCnQG/gnHPMWbuLP320icx9+R59cec35Z5RKVzTv5OCXEQU6A3Zl9sO8Lv/bGBl9iGP9rbRTbhrVDLXp8bROEJz5CJSRoHeAGXvP8a0ORv4cM0uj/boyAhuH96VHw5L1B7kIlKNAr0BOXKiiOnzMnhxcRaFJace8dY4PIzvD+3MlBFJnBfV2I8VikhDpkBvAEpLHbO+zOWPH21kX16hR993+nbkwW91I+587XwoImemQPezldkH+fXsdXyde9ijvX98K/7nqp4MiD/PT5WJSKBRoPvJ/rwCHpuzkTfTcj3aO7SM5OErujOub0fMtHJFRLynQK9nJaWO177I5k9zNnrsudI4IozbL+3C7Zd1pVlj/bWIyNlTctSjtdsP88g7a1hdZXplTM92/OqqnponF5FzokCvB/kFxTz+yWZeXLzV45FvnVs349ffuYAR3dv6rzgRCRoK9Do2f9MefvnOWrYfOl7R1jg8jNsv68rky7pqF0QR8RkFeh3Zn1fAox+s571VOzzaL+ramv+7phddYpr7qTIRCVYKdB9zzvH+1zv59ex1HMg/tab8vGaN+J+renJt/05avSIidUKB7kN7jpzgkXfX8sn63R7t1/TryP9c1ZPWzZv4qTIRCQUKdB9wzvHeqh387+x1HD5+ao/yDi0j+f343ozopg89RaTuKdDP0b68Ah55Zw0frfO8Kp80JJ6Hr+hOtJ7hKSL1RIF+Duas3cUj76xhf6W58tjzmvLH6/pwUVIbP1YmIqHIq0A3s7HAU0A48JxzblqV/tuBKUAJkAfc5pxb7+NaG4wjJ4r4zez1/L+VnrftTxoSz8+/3YPmTfRzUkTqX63JY2bhwHRgDJALrDCz2VUC+zXn3N/Lx48DHgfG1kG9fvfF1gPc++9VHuvK27eI5I8T+nBpSowfKxORUOfNpeRgIMM5lwlgZm8AVwMVge6cO1JpfBTgCDJFJaU8OXczzy7Y4nG359X9OvLouF60bKa5chHxL28CvROQU+k4FxhSdZCZTQHuAxoDI2t6ITO7DbgNID4+/mxr9Zusffnc/cZXHnuwtIiM4HfX9uY7fTv6sTIRkVN89kBK59x051xX4CHgl6cZM8M5l+qcS42JafjTE845/t+XuVz59OceYX5R19Z8dO+lCnMRaVC8uULfDsRVOo4tbzudN4Bnz6WohiCvoJhfvrOGdyvdut8o3Lj/8m785JIuhIXpbk8RaVi8CfQVQLKZJVIW5BOBmyoPMLNk51x6+eGVQDoBbO32w0x9bSVZ+49VtCW2ieLpif3pHdvSj5WJiJxerYHunCs2s6nAR5QtW3zBObfOzB4F0pxzs4GpZjYaKAIOArfUZdF1xTnHS0uy+N2HGygqOfXJ54SBsfxm3AVEaTmiiDRg5px/FqSkpqa6tLQ0v7x3TQ4fL+KhWV8zZ92uiraoxuH87treXNO/kx8rExHxcNr5Xl1yAmtyDzPltZVkHzg1xXJBxxY8c9MAEttE+bEyERHvhXSgO+d4dXk2j76/nsKS0or27w/tzC++3UMPnxCRgBKygZ5fUMwjVVaxNG8SwWPX9eHKPh38WJmIyDcTkoG+ZW8ed/zrSzbvzqto69GhBc9OGkCCplhEJECFXKDPWbuL+99aTV5BcUXbDalx/ObqCzTFIiIBLWQCvbiklD9/vJm/f7aloq1JRBi/vaYX16fGneF3iogEhpAI9AP5hdz1+lcsythX0RZ3flOenTSQXp10o5CIBIegD/Q1uYe5/V9femx3O6JbDE/e0F87JIpIUAnqQH97ZS4/f3sNBcWnliTePSqZu0clay8WEQk6QRnoRSWl/P7DDby4OKuiLToygidv6MeoHu38V5iISB0KukA/kF/IlFdXsjRzf0Vbctvm/ON7A+kS09yPlYmI1K2gCvQNO4/wk5fTyD14ar78Wxe04y/X99NzPkUk6AVNyn24Zic/e3M1x4tKKtruG5PC1BFJmi8XkZAQ8IFeWup48tN0nv701BbszZtE8MQN/RjTU/PlIhI6AjrQ8wuK+dmbqz22vE1o3Yx/fj+V5HbRfqxMRKT+BWyg5x48xo9fSmPjrqMVbZckt+GZGwdofbmIhKSADPQVWQe4/ZUv2Z9fWNH2o4sT+fkV3YkI99lzr0VEAkrABfqbaTk88s6aikfENQo3fndNb64fpP1YRCS0BVSgPzMvnT9/vLniuHVUY/7+vYEMSjjfj1WJiDQMATU/cXFyDI0jykru3j6a96YOU5iLiJQLqCv0fnGt+NOEPny4ZiePX9+PKN0sJCJSIeAS8ep+nRjXtyNmullIRKSygJpyOUlhLiJSXUAGuoiIVKdAFxEJEgp0EZEgoUAXEQkSCnQRkSChQBcRCRLmnPPPG5vtBbb55c3PTRtgn7+LqGehds6hdr6gcw4k+5xzY2vq8FugByozS3POpfq7jvoUauccaucLOudgoSkXEZEgoUAXEQkSCvSzN8PfBfhBqJ1zqJ0v6JyDgubQRUSChK7QRUSChAJdRCRIKNBrYWbnm9knZpZe/t/zzjC2hZnlmtkz9Vmjr3lzzmbWz8yWmtk6M/vazG7wR63nwszGmtkmM8sws4dr6G9iZv8u719uZgl+KNOnvDjn+8xsffnf6adm1tkfdfpSbedcadx1ZubMLGCXMirQa/cw8KlzLhn4tPz4dH4LLKyXquqWN+d8DPi+c+4CYCzwpJm1qr8Sz42ZhQPTgSuAnsCNZtazyrAfAQedc0nAE8Bj9Vulb3l5zl8Bqc65PsAs4I/1W6VveXnOmFk0cDewvH4r9C0Feu2uBl4q//ol4JqaBpnZQKAd8HH9lFWnaj1n59xm51x6+dc7gD1ATH0V6AODgQznXKZzrhB4g7Lzrqzyn8MsYJQF9tNVaj1n59x859yx8sNlQGw91+hr3vw9Q9nF2GPAifosztcU6LVr55zbWf71LspC24OZhQF/Ae6vz8LqUK3nXJmZDQYaA1vqujAf6gTkVDrOLW+rcYxzrhg4DLSul+rqhjfnXNmPgP/WaUV1r9ZzNrMBQJxz7j/1WVhdCLhnitYFM5sLtK+h65HKB845Z2Y1rfOcDHzonMsNlAs4H5zzydfpALwC3OKcK/VtleIvZnYzkAoM93ctdan8Yuxx4Ad+LsUnFOiAc2706frMbLeZdXDO7SwPrz01DBsKXGJmk4HmQGMzy3POnWm+3a98cM6YWQvgP8AjzrlldVRqXdkOxFU6ji1vq2lMrplFAC2B/fVTXp3w5pwxs9GU/WAf7pwrqKfa6kpt5xwN9AIWlF+MtQdmm9k451xavVXpI5pyqd1s4Jbyr28B3qs6wDk3yTkX75xLoGza5eWGHOZeqPWczawx8A5l5zqrHmvzlRVAspkllp/LRMrOu7LKfw4TgHkusO/Eq/Wczaw/8A9gnHOuxh/kAeaM5+ycO+yca+OcSyj/97uMsnMPuDAHBbo3pgFjzCwdGF1+jJmlmtlzfq2s7nhzztcDlwI/MLNV5b/6+aXab6B8Tnwq8BGwAXjTObfOzB41s3Hlw54HWptZBnAfZ17h1OB5ec5/ouz/Mt8q/zut+kMuoHh5zkFDt/6LiAQJXaGLiAQJBbqISJBQoIuIBAkFuohIkFCgi4gECQW6hBQzyyq/cUYk6CjQRc6Cmc00s//zdx0iNVGgi4gECQW6hCQzG1z+gI5DZrbTzJ4pvzUcK/OEme0xsyNmtsbMepnZbcAk4EEzyzOz9/17FiKetDmXhKoS4F4gjbINm/5L2a6ZTwKXU7atQQplW+Z2Bw4552aY2UVArnPul/4oWuRMdIUuIck596Vzbplzrtg5l0XZhlQnt4otomwXvu6UbY+xodL+8CINlgJdQpKZpZjZB2a2y8yOAL8H2gA45+YBz1D26LI9ZjajfKtgkQZNgS6h6llgI5DsnGsB/AKoeDqJc+5p59xAyp5DmQI8cLKrvgsV8ZYCXUJVNHAEyDOz7sAdJzvMbJCZDTGzRkA+Zc+ZPPk0pt1Al/ouVsQbCnQJVfcDNwFHgX8C/67U16K87SCwjbKnFP2pvO95oGf56ph3661aES9oP3QRkSChK3QRkSChQBcRCRIKdBGRIKFAFxEJEgp0EZEgoUAXEQkSCnQRkSChQBcRCRL/H94eH2CGIBTRAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"rslts.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Regression "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are several alternatives to estimate the conditional mean functions.\n",
"\n",
"* pooled regressions\n",
"\n",
"* local linear regressions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Pooled regression"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We estimate the conditinal mean using the whole function.\n",
"\n",
"\\begin{align*}\n",
"Y = \\alpha + \\tau D + \\beta X + \\epsilon\n",
"\\end{align*}\n",
"\n",
"This allows for a difference in levels but not slope."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"
OLS Regression Results
\n",
"
\n",
"
Dep. Variable:
next
R-squared:
0.670
\n",
"
\n",
"
\n",
"
Model:
OLS
Adj. R-squared:
0.670
\n",
"
\n",
"
\n",
"
Method:
Least Squares
F-statistic:
6658.
\n",
"
\n",
"
\n",
"
Date:
Wed, 30 Jun 2021
Prob (F-statistic):
0.00
\n",
"
\n",
"
\n",
"
Time:
12:05:34
Log-Likelihood:
3661.9
\n",
"
\n",
"
\n",
"
No. Observations:
6558
AIC:
-7318.
\n",
"
\n",
"
\n",
"
Df Residuals:
6555
BIC:
-7298.
\n",
"
\n",
"
\n",
"
Df Model:
2
\n",
"
\n",
"
\n",
"
Covariance Type:
nonrobust
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
coef
std err
t
P>|t|
[0.025
0.975]
\n",
"
\n",
"
\n",
"
Intercept
0.4427
0.003
139.745
0.000
0.437
0.449
\n",
"
\n",
"
\n",
"
D[T.True]
0.1137
0.006
20.572
0.000
0.103
0.125
\n",
"
\n",
"
\n",
"
last
0.3305
0.006
55.186
0.000
0.319
0.342
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
Omnibus:
595.910
Durbin-Watson:
2.143
\n",
"
\n",
"
\n",
"
Prob(Omnibus):
0.000
Jarque-Bera (JB):
3444.243
\n",
"
\n",
"
\n",
"
Skew:
-0.225
Prob(JB):
0.00
\n",
"
\n",
"
\n",
"
Kurtosis:
6.522
Cond. No.
5.69
\n",
"
\n",
"
Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: next R-squared: 0.670\n",
"Model: OLS Adj. R-squared: 0.670\n",
"Method: Least Squares F-statistic: 6658.\n",
"Date: Wed, 30 Jun 2021 Prob (F-statistic): 0.00\n",
"Time: 12:05:34 Log-Likelihood: 3661.9\n",
"No. Observations: 6558 AIC: -7318.\n",
"Df Residuals: 6555 BIC: -7298.\n",
"Df Model: 2 \n",
"Covariance Type: nonrobust \n",
"==============================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"Intercept 0.4427 0.003 139.745 0.000 0.437 0.449\n",
"D[T.True] 0.1137 0.006 20.572 0.000 0.103 0.125\n",
"last 0.3305 0.006 55.186 0.000 0.319 0.342\n",
"==============================================================================\n",
"Omnibus: 595.910 Durbin-Watson: 2.143\n",
"Prob(Omnibus): 0.000 Jarque-Bera (JB): 3444.243\n",
"Skew: -0.225 Prob(JB): 0.00\n",
"Kurtosis: 6.522 Cond. No. 5.69\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\"\"\""
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"smf.ols(formula=\"next ~ last + D\", data=df_base).fit().summary()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Local linear regression"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now turn to local regressions by restricting the estimation to observations close to the cutoff.\n",
"\n",
"\\begin{align*}\n",
"Y = \\alpha + \\tau D + \\beta X + \\gamma X D + \\epsilon,\n",
"\\end{align*}\n",
"\n",
"where $-h \\geq X \\geq h$. This allows for a difference in levels and slope."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Bandwidth: 0.3 Effect 8.318% pvalue 0.000\n",
" Bandwidth: 0.2 Effect 7.818% pvalue 0.000\n",
" Bandwidth: 0.1 Effect 6.058% pvalue 0.000\n",
" Bandwidth: 0.05 Effect 4.870% pvalue 0.010\n",
" Bandwidth: 0.01 Effect 9.585% pvalue 0.001\n"
]
}
],
"source": [
"for h in [0.3, 0.2, 0.1, 0.05, 0.01]:\n",
" # We restrict the sample to observations close\n",
" # to the cutoff.\n",
" df = df_base[df_base[\"last\"].between(-h, h)]\n",
"\n",
" formula = \"next ~ D + last + D * last\"\n",
" rslt = smf.ols(formula=formula, data=df).fit()\n",
" info = [h, rslt.params[1] * 100, rslt.pvalues[1]]\n",
" print(\" Bandwidth: {:>4} Effect {:5.3f}% pvalue {:5.3f}\".format(*info))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There exists some work that can guide the choice of the bandwidth. Now, let's summarize the key issues and some review best practices."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Checklist"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Recommendations:**\n",
"- To assess the possibility of manipulations of the assignment variable, show its distribution.\n",
"- Present the main RD graph using binned local averages.\n",
"- Graph a benchmark polynomial specification\n",
"- Explore the sensitivity of the results to a range of bandwidth, and a range of orders to the polynomial.\n",
"- Conduct a parallel RD analysis on the baseline covariates.\n",
"- Explore the sensitivity of the results to the inclusion of baseline covariates."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n",
"\n",
"* **Cattaneo, M.D., Idrobo, N., & Titiunik, R. (2019)**. [A practical Introduction to Regression Discontinuity Designs: Foundations](https://www.amazon.com/Introduction-Regression-Discontinuity-Quantitative-Computational/dp/1108710204?asin=1108710204&revisionId=&format=4&depth=1), Cambridge University Press.\n",
"\n",
"\n",
"* **Cunningham, S. (2021)**. [Causal Inference: The Mixtape](https://www.scunning.com/mixtape.html#:~:text=Causal%20Inference%3A%20The%20Mixtape.%20An%20accessible%2C%20contemporary%20introduction,allow%20social%20scientists%20to%20determine%20what%20causes%20what.). *Yale University Press*\n",
"\n",
"\n",
"* **Hahn, J., Todd, P. E., and van der Klaauw, W. (2001)**. [Identification and estimation of treatment effects with a regression-discontinuity design](https://www.jstor.org/stable/2692190). *Econometrica*, 69(1), 201–209.\n",
"\n",
"\n",
"* **Imbens, G., & Lemieux, G. (2007)**. [Regression discontinuity designs: A guide to practice](https://scholar.harvard.edu/files/imbens/files/regression_discontinuity_designs_a_guide_to_practice.pdf). *Journal of Econometrics*, 142 (2) :615-635.\n",
"\n",
"\n",
"* **Lee, D. S. (2008)**. [Randomized experiments from nonrandom selection in US House elections](https://www.sciencedirect.com/science/article/abs/pii/S0304407607001121). *Journal of Econometrics*, 142(2), 675–697.\n",
"\n",
"\n",
"* **Lee, D. S., and Lemieux, T. (2010)**. [Regression discontinuity designs in economics](https://www.aeaweb.org/articles?id=10.1257/jel.48.2.281). *Journal of Economic Literature*, 48(2), 281–355.\n",
"\n",
"\n",
"* **Thistlethwaite, D. L., and Campbell, D. T. (1960)**. [Regression-discontinuity analysis: An alternative to the ex-post facto experiment](https://psycnet.apa.org/record/1962-00061-001). *Journal of Educational Psychology*, 51(6), 309–317.\n",
"\n",
"\n"
]
}
],
"metadata": {
"celltoolbar": "Slideshow",
"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.8"
}
},
"nbformat": 4,
"nbformat_minor": 4
}