{ "cells": [ { "cell_type": "markdown", "source": [ "# Linear Econometric Models\n", "\n", "So far in our Scientific Python adventures we have covered:\n", "\n", "* NumPy\n", "* Pandas\n", "* Matplotlib\n", "\n", "For our next adventure into Python's Scientific stack we are going to look into estimating simple linear econometric models. This is just the tip of the iceberg as far as statistical modelling packages and capabilities in Python - but it is enough to get you hopefully confident enough to go looking for further modelling libraries on your own.\n", "\n", "This notebook is structured a little bit like how a simple project might look. We do all the tasks needed from downloading the data, some data visualization and the econometrics together. This should show you how far we have come on our Python Adventure. \n", "\n", "## Loading some data\n", "\n", "The plan is to use the data from Acemoglu, Johnson and Robinson's 2001 AER paper (hereafter AJR) on linking institutions to economic development as an example on how to run OLS and IV style regression models. That means we have to load the data.\n" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "### Importing into Pandas\n", "\n", "AJR provide multiple datasets on their replication page that we have put the 'data' directory. \n", "Each dataset allows us to construct estimates from one of the tables in the paper.\n", "\n", "We'll start by working with table 1's data, `maketable1.dta`. \n", "It's a Stata datafile, so we will need to use the function `read_stata` from the pandas library:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 1, "source": [ "import pandas as pd\n" ], "outputs": [], "metadata": { "collapsed": true } }, { "cell_type": "code", "execution_count": 2, "source": [ "data = pd.read_stata('data/maketable1.dta')" ], "outputs": [], "metadata": { "collapsed": false } }, { "cell_type": "code", "execution_count": 3, "source": [ "data.head()" ], "outputs": [ { "output_type": "execute_result", "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
shortnameuro1900excolonyavexprlogpgp95cons1cons90democ00acons00aextmort4logem4loghjyplbaseco
0AFG0.0000001.0NaNNaN1.02.01.01.093.6999974.540098NaNNaN
1AGO8.0000001.05.3636367.7706453.03.00.01.0280.0000005.634789-3.4112481.0
2ARE0.0000001.07.1818189.804219NaNNaNNaNNaNNaNNaNNaNNaN
3ARG60.0000041.06.3863649.1334591.06.03.03.068.9000024.232656-0.8722741.0
4ARM0.0000000.0NaN7.682482NaNNaNNaNNaNNaNNaNNaNNaN
\n", "
" ], "text/plain": [ " shortnam euro1900 excolony avexpr logpgp95 cons1 cons90 democ00a \\\n", "0 AFG 0.000000 1.0 NaN NaN 1.0 2.0 1.0 \n", "1 AGO 8.000000 1.0 5.363636 7.770645 3.0 3.0 0.0 \n", "2 ARE 0.000000 1.0 7.181818 9.804219 NaN NaN NaN \n", "3 ARG 60.000004 1.0 6.386364 9.133459 1.0 6.0 3.0 \n", "4 ARM 0.000000 0.0 NaN 7.682482 NaN NaN NaN \n", "\n", " cons00a extmort4 logem4 loghjypl baseco \n", "0 1.0 93.699997 4.540098 NaN NaN \n", "1 1.0 280.000000 5.634789 -3.411248 1.0 \n", "2 NaN NaN NaN NaN NaN \n", "3 3.0 68.900002 4.232656 -0.872274 1.0 \n", "4 NaN NaN NaN NaN NaN " ] }, "metadata": {}, "execution_count": 3 } ], "metadata": { "collapsed": false } }, { "cell_type": "markdown", "source": [ "## Data Visualization\n", "\n", "AJR's main hypothesis is to investigate the relationship between institutional differences and economic outcomes.\n", "\n", "To do this the paper operationalizes the variables of interest as follows:\n", "\n", "* economic outcomes are proxied by log GDP per capita in 1995, adjusted for exchange rates, `logpgp95`\n", "* institutional differences are proxied by an index of protection against expropriation on average over 1985-95, constructed by the Political Risk Services Group, `avexpr`\n", "\n", "Before we jump on to running regressions, let's check out whether there is any relationship between these variables graphically:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 5, "source": [ "import matplotlib.pyplot as plt\n", "plt.style.use('seaborn-notebook')" ], "outputs": [], "metadata": { "collapsed": false } }, { "cell_type": "code", "execution_count": 7, "source": [ "fig, ax = plt.subplots()\n", "\n", "ax.scatter(data['avexpr'], data['logpgp95'])\n", "ax.set_xlabel('Expropriation Risk (Institutions)')\n", "ax.set_ylabel('Income')" ], "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "Text(0, 0.5, 'Income')" ] }, "metadata": {}, "execution_count": 7 }, { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoPElEQVR4nO3dfZwdVZkn8N/TN28d0pAXSNKdCCGIecNZ1IyDIEIkMZoRZJXhIyoj7KDiB0WDrLDLaJJ14sqaEdeZ0VFR4iKgEhGJk5WEt6gguElEoBPiC5CY1455IZ3kdjrpfvaPqpvcvqmqW/fWy6lT9ft+PvnQ99bte09XN/XUec45zxFVBRERFVeL6QYQEZFZDARERAXHQEBEVHAMBEREBcdAQERUcAwEREQFN8h0Axp16qmn6qRJk0w3g4jIKmvXrv2Lqp7mdcy6QDBp0iSsWbPGdDOIiKwiIpv8jjE1RERUcAwEREQFx0BARFRwDARERAXHQEBEZIFybx827T6Icm9f7O9t3awhIqIiOdrXj8UrNmDV+p3Ytq+MjpGtmDN9HG6bNw2DSvHcyzMQEBFl2OIVG3DXk68ce7xlb/nY4wWXzojlM5gaIiLKqHJvH1at3+l5bNX6nbGliRgIiIgyqqu7B1v3lj2PbdtbRld3Tyyfw0BARJRRY9uG4aShJc9jw4eWMLZtWCyfw0BARBRCkrN2TONgMRFRgDRm7fjp6u7BIZ/AU+7tQ1d3D84Yc1Lkz2GPgIgoQGXWzpa9ZfTr8Vk7i1dsSPyzx7YNQ8fIVs9j7SNbmRoiIkpaWrN2/LQOKWH2tLGex2ZPG4vWId7jB41iICAi8tHV3YNt+7xn7WzfF9+snSDa4PPNYCAgotjkbUA1rdSMn3JvHx7d0OV57NENXbGdZw4WE1FkJgdUk9Q6pIQ508cNWNlbMWf6uKZTM5WB3rFtwwLfI2gdQaVHEsdgMQMBEUWWRhkEU26bNw2AMyawfV8Z7VVBrlGVgLmycwe27+tB+8hheMeM8Z4B82hfP+785ctoaQH6+k98rzh7JIkFAhH5LoB3A+hS1XPc50YD+CGASQBeAXClqu5Nqg1ElLx6A6qfnTs1tkFNEwaVWrDg0hn47Nypoe7igyx8qBPff2bzscdb9/Xgridfgapi4WXnDHjt4hUbcPfTvrtLRuqR1Eqyz7YUwDtrnrsVwKOqejaAR93HRGSxLAyopqF1SAlnjDmpqYvv0b5+fO7B5wcEgWrL1m4ZkO8PCq6lFuDq805vqkfiJ7FAoKq/ALCn5un3APie+/X3AFye1OcTUTpMD6jawLm79w4CAHDgcB827zl07HFQcIUC1104Odaxl7RHccap6nb36x0AxoX5JhFZKCIqIrpt27bkWkdEDasMqHqJM31hq6C7+2o9R473CNIOrsaG81VVEXIqrKouVFVRVeno6Ei4ZUTUqNvmTcO1F0zCxFGtKAkwcVQrrr1gUqzpC9P8psbWmzIbeHdf5XBVIAgKrrOmxLeQrCLtWUM7RaRdVbeLSDsA7wmyRGSVOAdUs8Zvauwtc6fg9oc31p0yW7m73+IzDbRiz6HDAx5Xz1basreMkjt76LEXd2LQcol1am7aPYKHAHzY/frDAH6a8ucTUYKiDKhmlV+tofd+4ynP5xct7xzw/UF399X+47kdOFo1T7QSXGdNOQ3A8SmklZlGcdY6SiwQiMh9AH4NYIqIbBGRfwDwJQBzROQPAGa7j4mIMikov79xZ7fn8/f+ZjM+9+ALAy7qt82bhqvPOz3ws5Y/t/2Ei3u5tw+Pb9zl+fo4ax0llhpS1at8Dl2S1GcSEcUpKL/vtcir8vzdT2/CoJIcW0w3qNSCL1z+evT1K+79zZ99P6923UWYqbksQ02UA80OQlLygmbv1EvPe92xf/7d0/Ha0/wv3LXrLtKaPcQSE0SGRB2EzLuw9XiSFFRraMq4Nqzf7p0eArzv2G9/eCP+uOug7/fUXtyTqnVUi4GAyBC/+jzPvLR7wAUmT3V7wshaATu/WkO3zJ2CxSs24N7fbA5VC6jc24eVnTsCP8vr4h5nrSM/4kznt8fMmTN1zZo1pptBFEm5tw9z7ljtOaWw5FNkbOKoVqyaf1GuZuR4WbS80/MO+NoLJhkNhH49lM89+IJnTaDa9m7afRAXffkJ3/e/4o0T8KX3/ZVvsIvaQxKRtao60+tYcfqZRBnSzCBknur2+DG9I1gQv6mxCy6dHmoxXdvQwYHv/9/nTQ/s8SQ5NZepISIDghYZ+fUIkigt0MxdZpK5+7RmycQp7GK6XQeCg/iuAz0YPWJIUs0MxEBAZEAzg5BxDg42k4dPI3cfFCCzXsCucsfuT+q8Q73jyWFqiMgQv/o8D3z8/MTr9vitlg1ardrM9zQqzwXsTh89HCOGerd/xNASTh89POUWHcceAZEhQSmFJOv2NLORTJqbz6QxS8aE1iElXH7uBM89Ca5400SjQY6BgMgwv5RC/VRDc5rJw6eZu89jAbtKWu3xjU6dzRYB+hXoGDkMc92tKk1iaoioYJpZrWpi85k8FbCrpNW27nMGjPvdWfuXTB2LBZfOML5QkIGAqGCaycPnOXeftKC02uMbd2WihAhTQ0QF1EwePq+5+6TZMCWWK4uJCixr6wjyKGgVeZqrxbmymIg8NZOHz1PuPg02pNWYGiIiqpJEjyfraTWmhoiIkM7KaZNptaDUEHsERBnHnHw6/MqCA/GV/05qbUhUDAREGZW1uvymJRkQ01w5nUUMBEQZlcYdqpes9UDSCIg2TPFMEgMBUQaZuEPNag8kjYBoc9XTOBSvf0lkgTB3qHFLo7poo9LaqMaGKZ5JYiAgyqC0a/tkdWewNAOiX1nwrEzxTBJTQ0QZFLRxTRJ3qFnNkaeZsslj1dOw2CMgyqgod6jl3j5s2n0w9J382LZhGO5z0WsdUjKWIzeRsiniymn2CIgyqpk71KwO+EaR9VW5ecBAQJRxjSxC8pthc7SvH9ddONk3mHR19+DgYe/ew6HDfUanTxY5ZZMWBgKinAga8L33N5txzzObfXsIY9uGYcIo71x8x6hsTJ/M6qrcPLCzr0iUY43m9yuCBnz7+hE4JTRv0yebPYdFxR4BUUZEze8HzbCp5bUoLQ+5+DyOkaSBgYAoI6KuoA2aclrLa0poHnLxpspy2I4hkgjmUwlxLeiqnnLaAsDvJjhoDr6t0yezuijOBuwRUKFlJZUQ14Ku2rv6O3/5Mu5+etMJr7Mx719PVhfF2YCBgAotK6mEuFfQVu7qF1w6HYNKYnXeP6yiF46LgoGACitLNeiTKimRh7x/WGmX5cgTjhFQYZmo8BnklrlTML297Vhev9QCTG9vwy1zp0R+b1vz/o0qcuG4KNgjoMIymUrw2vzl9oc3Yv327mOv6esH1m/vxu0Pb+SMl5CK1AOKEwMBFZaJVILf4PRNs1+XmTRVHnAVcmOMBAIR+RSAjwAQAN9W1a+aaAdR2ouo/Aan95ePcMYLGZN6IBCRc+AEgTcD6AXwcxH5mar+Me22EKWZSgganP71S7vRfsowbN134rgEZ7xQ0kwMFk8D8IyqHlLVowBWA3ivgXYQHZPGYGrQ4PTOV3vwlsljPI81k6YyvUCO7GIiNfQCgMUiMgZAGcA8AGsMtIMoVfUGpxdcOgNtrYMjpamyskCO7JJ6IFDVDSJyO4CVAA4CeBZA4G2LiCwEsAAA2tvbE24hUTLqDU63tQ6OnKaKY4Gc14wmyjdRVbMNEPkigC2q+vUwr585c6auWcMOBNmh9qJafcdee9cf9Y693NuHOXes9uxxTBzVilXzLwq8sNfrTTBA2E1E1qrqTK9jpmYNjVXVLhE5Hc74wHkm2kGUlKCLalKD01Fr7fj1JvpV0SLCdFOOmVpH8GN3jOAIgBtUdZ+hdpBFbLojrZeiSWKeezML5CrntG3oYN8ZTT9euwUHqraxZGnn/DESCFT1QhOfS3aybQDUVA2jRhbI1Z7TsW1DsWP/Yc/3PeCzlzEXuuUHVxZT5mWlQmhYjaZoGu3pBL0+7AK52nPqFwSC2LjQzaZeZZoYCCjTslQhNKywKZpGezphXl9vgVy5tw+b9xzCw507Qv88I4aWPHsFNi10s61XmTYGAso0GzcbCZuiabSnE+b11Xe81edlwIVwbxn9Ae0ff/JQ7Oo+fKw3oapY+pTdm9vY1qtMGwMBZZqtm43US9E02tOp9/qbZr8OX3nk9753vLUXQj8TR7XioRveiu7DRwZMeRWxd3MbG3uVaWMgoEyzdbOReimaRns69V6/aHknlq3beuy56jvez86d6nshrDVn+jiMHjEEo0cMOfbckT7FNedPwidnnT0gQNjCxl5l2hgIKPPSrhAaJ79poo32dIJeP+6UYXjqpd2en79q/U68/69f43shBJyCYx2jTjynQXl1m9jaq0wTAwFlXh43G2m0pxP0+rdMHoOf/HbrCc8Dzh0vIL4XwgkjW/Hda2bi9NEnFtzLS17d1l5lmjhcTtbI23aLjW6r6Pf6hZfOQMfIVs/vaR/ZitNHD8ec6eM8j7996lhMGX/yCee0Xl7dtqqm3MIymPFaQ41irSHKmzjWESxa3ul5x3vtBZOw4NIZA9I8W/aWUWpxtsKcMHIY3jFj/AnTKDftPohZS55Av8floSTAYzdfbGVevcjrCIJqDbFHkFOsR2+PRns6Xq+vt/F9Jb02a8ppAJwgAABb9/XgridfweIVGwZ8RiWv7sXmvHreepVxYSDImaN9/Vi0vBNz7liNWUuewJw7VmPR8k4c7QuaOU62q2x8X/k1V298X1Hu7cPjG3d5fn9tuqeSV/fCvHr+MBDkTGWAb8veMvr1+ABf7R0f5UfYfH6YaZTVmFcvDs4ayhEunCmmsPPkG51GmcfZWuSNPYIcafSOj/IhbD6/2XQP8+r5x0CQI3kd4KNgjVzgme4hLw2lhkTkNFX1Hm0i47hwJh42TjEMu/qa6R7yEmodgYj8DYAfAWhR1deIyEwAH1XVjybdwFpcRxAsyT1x886vpMJNs1+HPYd6rbho2hjEksTzcVzQOoKwgeBJAB8BcI+qvsF9rlNVU19nzkAQDv8HaJzfoqwRQ0s41NvHGvYW4f4DJ4pj8/ohqrpeRKqf643cMkpMEnvi5lnQjKvKpiy21toporzUSUpL2NB4WERGAFAAEJHpADgFhRJhYlV00IyrWjbW2rFBXL/3vNVJSkPYHsFiACsBdIjIUgDvBPChpBpFxWSyOx80x76WqRr2eU33xf175/4DjQsVCFT1/4rIRgBzAQiAf1LVPybaMiqcJLrzYS+eQTOuaqU9FTfv+e64f+/cf6Bxof+KVPUlVf2Gqn6dQYDiFnd3vpmaS5U59iOGBt9tpz0VN89lQ5JI47BOUuNCBQIReauI/FJEtolIl4jsEpGupBtHxRH3quhmLp6DSi347NypOLl1sOfxUgtw9XlnpLr4Ku/57qRWw3PhXGPCjhF8F8BtANYCsPsvj0JLMycdZ3c+Ss2lru4e7HjV++Kj/cB1F56Zajom7/nupNI4XDjXmLB/0XtV9X43PbSp8i/RlpExJkpZx9mdj3KXGVSmo2NU+vnlvJcNSTqNwzpJ4YQNBPeKyPUiMlpEhlf+JdoyMsZUTjqu7nyUi2fW8stZa08SmMYxL+zK4qsAfBtA5f8uAaCqmvpfIVcWJ6vc24c5d6z27KpPHNWKVfMvSvziE0dKqt7WjUGyVqYja+1JSl6nx2ZFHCUmXgFwBYB1qmp0qysGgmTlZa/aOC6eSVyYorwnL5QURRwlJrapKq++BZCXOdhxDBbGWaYjzFqAehd6lg2hpIQNBI+KyO0Afoiq0hKquj6RVpExeStlnfbF0+9iHrRo6rZ503K9YIyyL2wgqJSTuLLqOQUwOd7mUBaErW1PxwXd8R/p08DprEf7+nH305uPPccCaZS2UGMEWcIxgvQwJx1e0OD0NedPChx3Oa1tKHbsP3zCsbQG56kYgsYIQvc7RWS6iNzg/psaX/Moq6LMwTZRQdSUegvY2oYO9p3OelrbUOz0CAIA95mm9IQtMXE1gFUAznX/PSIiH0yuWWQrE4vRTKu3gK378JGAtQDjMWFUfheMkR3CjhHcDOBNqroDAERkPICHAdyTVMPITkXcECTMTKugcZdBJUl8cJ5pPgoSevP6ShCofF2zWxlRpBo/Ngs708pvOmslSDzcuQM79vVg/MhhmDtjfCyD83kvYU3xCBsI/iQiiwB80338EQAvJdMkslXcBdJsuosNO9MqaDqr1Pw3DkXsoVHjwgaC6wF8DcBzcKaNPgLgY81+qIjMB3Cd+17PA7hWVTkqFkEWLppxLUaz8S42ygK22ov11n09sVysi9pDo8aF3aGsC8D74/hAEZkA4EYA01W1LCI/ct97aRzvXzRZumjGtRjN5rvYRhewJXmxznsJa4pP2FlDt4jI6KrHY0Tkv0b43EEAWkVkEIDhALZFeK9Ci6NSaJxTPaNWkrR9I5ZGz2VSG7MA+S9hTfEJmxq6SlVvrzxQ1d0i8gEAX270A1V1q4gsAbAZQBnASlVd2ej7UPS7ySR6E42kSLzSWbbexTZ7LpOs7ZS3ciGUnLCBwGv8KvSMowFvJDIKwHsAnAlgH4D7ReRDqvr9gO9ZCGABALS3tzfzsbkU9aKZZAomKEUSdNEMujCe1jYUbUO9t5E0rdlzmfTFmuVCKIywt31/EJGbxNEiIp8B0OwG9rMBvKyqu1T1CIAHAJwf9A2qulBVRVWlo6OjyY/Nnyhdf5MpmKB0VtBGLDv2H8Zl//arzC1Qi3ouk9yYpdJDWzX/Ijx288VYNf8iLLh0hmcvpUirwWmgsHf1NwL4PoAvwpnp8xSAq5v8zM0AznN3OCsDuAQAiwc1IcrdpKkUTJh0VvVdbG3PIIsDx1HPZRr76zbbQ8vqLC2KV6jfsqpuU9W3AxgD4FRVvURVmxrgVdVnACwDsA7O1NEWAN9q5r2o+btJUwOJYS6alQvjQze8FeNPHur52iwNHMd1Lk3tr2tqa1LKjtB5fhE5C8BZAAZVVhWr6opmPlRVF8DN+VM0zd5NmhpIbGRwtPvwEXR1Bxdky8LAsc2DslxrQEDIQCAi/xPOArANACq3YQqgqUBA8WtmAxYTA4mNXDRt2i3N1kFZW2dpUbzC9gj+DsBZqro/ycZQutLITXtppByDLXfaps5lVDYFW0pO2ECwnUEgv+LYzrGREheNXDRtu9O2bV9hm4ItJSfUDmUi8r8AvAbA/Ri4Z3HqqSHuUJYtac04yUItpbyq/h2eWCKbs4byImiHsrCB4HGPp9WdSZQqBoJsCdqiMc3pnQwU0fEc5ltQIAhMDYnIdPfLG2JvFVkvCzNOustHsHB5J3790m7seLWHc+AjsC2tRfGpN0bwH+5/a7sN4j43OfYWkTW6unuw1WOQEUh+xkklnXH/mj/jwOHj6wmyuOCMKOsCA4GqnplWQ8guR/v6cecvX0ZLC+BV7SHpGSe1tX1qcQ48UXjsO1NTFq/YgLuf3uQZBIBkZ5wEpaQqopZwJiqSpiqIUrEFXYhLLcAH3nx6otM7gxZBVXAOPFF47BFQwwIvxApcd+HkRAdqg2r7VHAOPFF4DATUMNM7XwWVqh4xtBRbCWeiomBqiBqWhdWotSuOx5/SirdMHo0Fl85AW2s2N68hyqpQC8qyhAvKzKlecDS4JJlYjcpFUEThNL2gjAgILiNhusgaF0ERRcdAQHXV2483bxfi2l4Gex2UdwwEFCgLZSTSUtvzaT9lGE5pHYxXy0ex/VVu4Uj5xUBAgYq0cUltz2frvh5s3Xd8UVqc5SvYy6AsYSDIiaQuLEXZuCTMauWKKD0hbhRPWcRAYLmkLyxZmCoaJK4AGGa1ckWUnlC98RYiExgILJfGhSWLu4TFHQCDej61mu0JFWm8hezCQGCxtC4sWdyPt9kA6NeDCOr51Gq2J1Sk8RayCwOBxdK+sGRlzn5QAFzZ6R0Aw/QgTlyt7Mwa2t9zNJaeUFHGW8g+DAQWK+qFJSgAbt1Xxj8++Dxuf99fDUgRhelB+PV84hqHyPp4S5I4SyrbOE0hJuXePmzafRDl3r76L45JUPG1PF9Y6lUf/fG6rVi8YsOxx/VSaLW/s0rPp3L+ah9Hcdu8abj2gkmYOKoVJQEmjmptukieib+5Rh3t68ei5Z2Yc8dqzFryBObcsRqLlnfiqN9GFmQEewQRmZ4OmMWB3KSFyedXj5GETaGlcdcax3iL6b+5RnCWlB0YCCIy/YeexYHcNNw2bxq6y0ewbN1Wz+PVF/h6KbTRw4dg0fLOVC+sUcZbTP/NhVXu7cPKTs6SskG2bh8s02jKIUlxpi9sMKjUgi9c/np0jPQeB6keI6mXQvvKI7/HXU++gi17y+jX4xfW6vRSVmTpb66eru4ebPXpiW3Zy61Es4SBIIIwKQdKTuuQEubOGO95rHaMxC83f9Ps16VyYY0rn2/T31zb0MHw61CVWpzjlA1MDUVQ1Fk7WRJ2jMQvhbZp98FEp+CmufAta39z3YePwG9MuK/fOT56xJB0G0WeGAgiKPJ0wKxodIykNjef9IU17ny+TX9zY9uGYcLIYQMK91VMGDksU0Gr6Jga8tBINz7O6YDUvGbHSJKcgptUPt+Wv7nWISW8wyd1944Z4zMVtIqOW1VWidKN54IZe1X/3uPcdnPT7oOYteQJ9Hv8L1YS4LGbL46UdrLhby6pc0uNC9qqkoGgyqLlnZ5d7msvmJSpaXmUjLgvrOXePsy5Y7Vn2mniqFasmn9RZi/gcbMhaOVdUCBgSHbZNC2PkhH3FNyirvz2UrTpzbbhYLGLlSEpCUVc+U32YSBw2TQtj06U1dRDUVd+k10YCFw2Tcuj42ypu5OVEt5EXhgIqrAbbx9b6u4QZVnqgUBEpgD4YdVTkwF8XlW/mnZbarEbbxdu/UgUj9QDgapuBHAuAIhICcBWAD9Jux1BstKNz2reOys4wE8UD9OpoUsA/ElVNxluR6bYkvc2jQP8RPEwfVV5P4D76r1IRBaKiIqIbtu2LYVmmVXJe9tQFtkkztMnioexQCAiQwBcBuD+eq9V1YWqKqoqHR0dyTfOIC5sa4wtdXeIssxkauhdANapqvdVr6CylPe2YYyCA/xE0ZkMBFchRFooTjZc2LKQ97ZxjCIrA/xENjISCETkJABzAHwsjc+z6cKWhYVtnJtPVCxGroKqelBVx6jqq2l8nm2Dr5+ZfTZGDx+4jd/o4YPxmdlnJ/7ZSY1RxLVVYxbk6WchAsxPH02cjYuOrvzW09hz6MiA5/YcOoIrv/U0VnzqbYl+dtxjFDb1xurJ089CVC33f702bfYNAHsO9GLjzm7PYxt3dmPPgd5EP78yRuGlmTEK23pjQfL0sxBVy30giPvClrQXd+wP3PD7xR37E/38OOfm52kqbJ5+FqJauQ8Eti06mjr+ZPhlGUotzvG41ea845qbb1tvLEiefhaiWrkfIwDsqio6esQQTBnXhvXbT0wPTRnXhtEjhsT2WUE57zjm5mdhKmxc8vSzENUqRCCwbdHRAx8/H+/9xlPYuLMbff1OT2DKuDY88PHzY/2cetNEo87Nz8JU2Ljk6WchqlWIQFBhy6KjYUMGYcWn3oY9B3rx4o79mDr+5Fh7AkB6s6ls6o3Vk6efhaiaqKrpNjRk5syZumbNGtPNsN6m3Qcxa8kT6Pf49ZcEeOzmi2MNmlle1d1o27L8sxD5EZG1qjrT61ihegR0XNo57yz2xppdF5DFn4UoitzPGiJvts2mSgLXBRA5GAgKrMglnLkugOg4poZiYmPe2LbZVHHKUrlvItMYCCLKQ/2ZIua8uS6A6Dg7rlQZxjyznThGQnQcA0EEzDPbrchjJETVmBqKgHlmuxV5jISoGnsEEdhW2ZS8VcZIaoNAlA1ouHkN2YQ9gghYfyafokwAyMPkASoeBoKIWH8mf6Ls2cz9nslGDAQRMc+cL1GK8dm4LSoRwDGC2PjlmckuUTag4eY1ZCsGAqIqUSYAcPIA2YqBgKhKlIVmXKRGtuIYAVGNKBMAOHmAbMSNaYh8RCkkaGMRQso3bkxD1IQoxfiKWMiP7MUxAsqsrK7OzWq7iJrFHgFlTlZX52a1XURRMRBQ5mR1dW5W20UUFW9jKFOyWto7q+0iigMDAWVKVlfnZrVdRHFgIKBMyerq3Ky2iygODASUGZW597OmjPU8bnJ1LlcNU55xsJiMq52N037KMExvb8P+nqOZWp3LVcOUV1xZTMYtWt7pubnP1eedjusunJy51blcNUw2ClpZzNQQGRU0G+fxjbsyebFlyXHKGwYCMoqzcYjMYyAgozgbh8g8I4FAREaKyDIReVFENojIW0y0g8zjbBwi80zNGvrfAH6uqleIyBAAww21gzKAs3GIzEp91pCInALgWQCTtYkP56yh/OJsHKLkZG3W0JkAdgG4S0R+KyJ3iggLtxNn4xAZYiIQDALwRgDfUNU3ADgI4NagbxCRhSKiIqLbtm1Lo41ERIVhIhBsAbBFVZ9xHy+DExh8qepCVRVVlY6OjsQbSERUJKkHAlXdAeDPIjLFfeoSAOvTbgcRETlMzRr6JIB73BlDLwG41lA7iIgKz7paQyKyC8Am0+3w0QHA1kEMtj19trYbYNtNidL2M1T1NK8D1gWCLBMRVVUx3Y5msO3ps7XdANtuSlJtZ4kJIqKCYyAgIio4BoJ4LTLdgAjY9vTZ2m6AbTclkbZzjICIqODYIyAiKjgGAiKigmMgICIqOAYCIqKCYyAgIio4BoKYiEjJ3V/hZ6bb0ggReUVEnheRZ0XEqh1/bN3yVESmuOe78m+/iHzadLvCEpH5ItIpIi+IyH0iYs3G0iLyKbfdnVk/5yLyXRHpEpEXqp4bLSKrROQP7n9HxfFZDATx+RSADaYb0aRZqnqu3+5FGVbZ8nQqgP8ES86/qm50z/e5AN4E4BCAn5htVTgiMgHAjQBmquo5AEoA3m+2VeGIyDkAPgLgzXD+Xt4tIq8126pASwG8s+a5WwE8qqpnA3gUdfZyCYuBIAYiMhHA3wK403RbisLd8vRtAL4DAKraq6r7jDaqOZcA+JOqZrWQopdBAFpFZBCc/cZtKeA2DcAzqnpIVY8CWA3gvYbb5EtVfwFgT83T7wHwPffr7wG4PI7PYiCIx1cBfBZAv+F2NEMBrBSRtSLyUdONaUBetjx9P4D7TDciLFXdCmAJgM0AtgN4VVVXmm1VaC8AuFBExojIcADzALzGcJsaNU5Vt7tf7wAwLo43ZSCISETeDaBLVdeabkuT3qqqbwTwLgA3iMjbTDcopIa3PM0adz+OywDcb7otYbk56ffACcQdAE4SkQ+ZbVU4qroBwO0AVgL4OYBnAfSZbFMU6pSFiKU0BANBdBcAuExEXgHwAwBvF5Hvm21SeO4dHlS1C06e+s1mWxRaw1ueZtC7AKxT1Z2mG9KA2QBeVtVdqnoEwAMAzjfcptBU9Tuq+iZVfRuAvQB+b7pNDdopIu0A4P63K443ZSCISFX/m6pOVNVJcLr5j6mqFXdIInKSiLRVvgbwDjjd58zLyZanV8GitJBrM4DzRGS4iAic827FID0AiMhY97+nwxkfuNdsixr2EIAPu19/GMBP43hTU1tVUjaMA/AT5/9nDAJwr6r+3GyTGmLtlqdu4J0D4GOm29IIVX1GRJYBWAfgKIDfAviW2VY15MciMgbAEQA3ZHmCgYjcB+BiAKeKyBYACwB8CcCPROQf4OzUeGUsn8Xqo0RExcbUEBFRwTEQEBEVHAMBEVHBMRAQERUcAwERUcExEBSQW3H0xZoKmJNMt8uLWzriwhCv+3Rljrj7+HoRmR9zW1REnhOR37n/vayRzxORpSLyiZCfdaOI3Ox+fXGUyrAicq6IXFnz3LMi0up+PeDc1Xmvy0XkzVWPZ4rIPc22rc5n3S4iH0jivWkgTh8tIHcV9LtVNfbFYyLSAmf1e+Q/LBEpqWqoEgBJ/kxVn6EA2lT1gIi8C8CPAIxyC5iF+f6lANao6r/Wed1wAM8DOEdVyyJyMYAlzVaHFZFr4JybK3yOv4KQ5y7szxAHETkNwK8ATFNVG+t4WYM9AjpGRKaKyJ9F5Az38QIR+YH79UIR+ZGIPOb2Jn7sVgCtHLtfRFbCWd07UkT+Xpx9Dp4TkZ9Urei8xq2j/pCIrHffb0LVsUfc178A4PUi8oRbzwki8gERecYtMvdbEbnEff42OHVvlrl3utPdNi1xj5dEZIk4dehfcL8uuceWisi/u+34g4j8H3fFbD1PABgBYFTVOah83vkiss5tS6eIXOVxrme55+Ycj/d+H4BfqGrZ4/smichfRGSxew42ishb3WNj3fP3vPvvDnEWT/0PALPd9nzNfa2KyAifczeg51J5LCJz4dRGutV97d/X9lbq/N5XisgP3XPypIiMDzpfqroLzkLBS0L8PigKVeW/gv0D8AqAF+EU3XoWzh1e5djVAJ6GU25iI4CT3ecXwqk2Oc59/F04d6mVY5sBnOo+PgdOaeJ29/EXAPzQ/foaAGUAU9zHCwAsqzp2AMBZVe15As7dKgCMwfFe7BQ4tYaqf6Zzqh4vrGrfxwE8AmCI++9RAB93jy2Fc9c5zD3WCWCOz3lTACPcrz8Apy681+f9FMBV7tcCYGTVZ30CwAcB/D8AE3w+5zsArq96fHHldwRgktuOyjn5IIAn3a/nA/hm1feNqjqvywJ+ltpztxTAJ7weexyrblu93/teAK9xH38bwOKg8+U+/jyAL5n+fybv/1hioriuUI9UgKre7d5pPwjgQlXdX3X4Z3q8QNp3APxL1bEVqvoX9+tZ7uNKudxvAvhd1Wt/paob3a/vhJMGqT72J582nwXgPrcHcQTAeBEZr07doSCzASxV1V4AEJG7APxnAN9wjz+oqj3usXXu56zyea+nxKnPNB7A231e8ziAfxSRswCs0uOF8QCnDEYZwCU157baRABBO90dUNXK8acB/HPV1/NF5Mtwau0/HPAeSaj3e39SVf/sfv00nBIbQPD52gFn3wlKEFNDNIA4dXtmANiHxmqdH4ipCUHvcx+Ar6vqDDiVRo/CuZOPqqfq6z4E1+A6X1XPBPA5AD8Qj20aVfWrcFIouwD8i4j8U9Xh38EJItMCPqOM4J/rsFd7VfXXAN4AYC2cnt3jAe8R5CgGXhvi2orS8zzXOV/D4JwPShADAdX6MpwLyRwA/y7O7msVf+sO4AHOne1jPu/xOIB5lRwwnO0Bq++wLxCRs0O8T62RAF52v/4vAIZWHdsP4BSf73sEwIdFZLCIDIZTtdHvjj+sfwawE8D1tQdE5HWq+idV/Sac7TSrS3uvg1P18h4RucjnvZ+Hk/pqiIicCWC/qv4AwE0A3iTO4H3QuYHH8T8C+Gv3Pdvh3On7vbZavd+7X7uDztc0DOxVUAKYGiquZSJSfYd2HZyUxMUA/kZVe0RkEZxUTOVC8Es4d8ET4AwKf8brjVX1BRG5FcAqcWbavISBVTafBLDEDQY74Ny9hvFpAA+KyF44G4vsrjr2NTi7lR2Ck7+v9i0Ar4VTKRNwUibfDvmZnlRVxZne+QMR+WbN4Rvdc9YL5+79kzXf+5w7AL5cRD6hqrUpnAcAfB3OuEMjLgZwk4j0wbnJu15V+0XkUQA3i8jvAKxW1Rtrvq/23H0bzt/Hejj1+qtTNXcDWCoifwfgK3DGhio/V73fux/P8+UO2r8dwBcbOQnUOE4fpVBEZCGcwcWbI77PNQiYykgOEXkYwK2q+tu6L84pd5bSh1Q17I0CNYmpIaJs+gSAdtONMOxkALeYbkQRsEdARFRw7BEQERUcAwERUcExEBARFRwDARFRwTEQEBEVHAMBEVHB/X+eJZUNNAft0wAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ], "metadata": { "collapsed": false } }, { "cell_type": "markdown", "source": [ "So there definitely appears to be a positive correlation between instututions and income. " ], "metadata": {} }, { "cell_type": "markdown", "source": [ "## Univariate Regression\n", "\n", "OK, so now we have visualized data let's run some OLS regressions. We will start with the simplest regression we can think of to map:\n", "\n", "$$\n", "logpgp95_i = \\beta_0 + \\beta_1 avexpr + u_i\n", "$$\n" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "To run OLS regressions we'll need a new package, `statsmodels`.\n", "\n", "First we need to load the `statsmodels` package:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 10, "source": [ "import statsmodels as sm\n", "sm.__version__" ], "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "'0.12.0'" ] }, "metadata": {}, "execution_count": 10 } ], "metadata": { "collapsed": false } }, { "cell_type": "markdown", "source": [ "There's two ways to write out the specification for a regression, the easiest is by writing a formulaic description of the regression we want to estimate:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 11, "source": [ "import statsmodels.formula.api as smf" ], "outputs": [], "metadata": { "collapsed": false } }, { "cell_type": "code", "execution_count": 12, "source": [ "simple_reg = smf.ols('logpgp95 ~ avexpr', data=data)" ], "outputs": [], "metadata": { "collapsed": true } }, { "cell_type": "markdown", "source": [ "To estimate the regression specification, we use the `fit` method:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 13, "source": [ "results = simple_reg.fit()" ], "outputs": [], "metadata": { "collapsed": false } }, { "cell_type": "markdown", "source": [ "The `results` object now contains a bunch of information.\n", "Most of the 'standard output' we have come to expect from software like STATA are contained in the `.summary()`:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 14, "source": [ "print(results.summary())" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: logpgp95 R-squared: 0.611\n", "Model: OLS Adj. R-squared: 0.608\n", "Method: Least Squares F-statistic: 171.4\n", "Date: Mon, 13 Sep 2021 Prob (F-statistic): 4.16e-24\n", "Time: 12:08:11 Log-Likelihood: -119.71\n", "No. Observations: 111 AIC: 243.4\n", "Df Residuals: 109 BIC: 248.8\n", "Df Model: 1 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept 4.6261 0.301 15.391 0.000 4.030 5.222\n", "avexpr 0.5319 0.041 13.093 0.000 0.451 0.612\n", "==============================================================================\n", "Omnibus: 9.251 Durbin-Watson: 1.689\n", "Prob(Omnibus): 0.010 Jarque-Bera (JB): 9.170\n", "Skew: -0.680 Prob(JB): 0.0102\n", "Kurtosis: 3.362 Cond. No. 33.2\n", "==============================================================================\n", "\n", "Notes:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n" ] } ], "metadata": { "collapsed": false } }, { "cell_type": "markdown", "source": [ "An observation was mistakenly dropped from the results in the original paper (see the note located in maketable2.do from Acemoglu’s webpage), so our coefficients differ slightly from the ones in the original paper." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "We can get the predicted values from the regression using the `.predict()` attribute. Let's use these in a plot that has predicted versus actual values.\n", "To do this, we'll need to drop missing data from our observations because the `predict()` method doesn't do well with missing data:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 18, "source": [ "# Drop missing observations from whole sample\n", "\n", "df_plot = data.dropna(subset=['logpgp95', 'avexpr'])\n", "\n", "# Plot predicted values\n", "\n", "fig, ax = plt.subplots()\n", "\n", "ax.scatter(df_plot['avexpr'], df_plot['logpgp95'], alpha=0.5, label='observed')\n", "ax.scatter(df_plot['avexpr'], results.predict(), alpha=0.5, label='predicted')\n", "ax.set_xlabel('Expropriation Risk (Institutions)')\n", "ax.set_ylabel('Income')\n", "ax.set_title('Predicted vs actual values for simple OLS regression')\n", "ax.legend()" ], "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "" ] }, "metadata": {}, "execution_count": 18 }, { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEaCAYAAAAcz1CnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABTUUlEQVR4nO2deZxcVbWov3Vq6Lk7U6eTTgIJgSQdkhBCCDJHIEaQi6IIDvhELg5cR656xftUog+feOU6Kw4PRa8KAoLT5UrCJIPMEEggBEgImTtpkk7PXVXnrPfHPtVdXanqru6u6qrq2t/vl3Sdoc7ZZ6i99lp7DaKqWCwWi6V0cfLdAIvFYrHkFysILBaLpcSxgsBisVhKHCsILBaLpcSxgsBisVhKHCsILBaLpcSxgmCEiMhNInKt//l0Edk8RudVETl6LM6VT0RkjYj8JgfH7XtuY4WInCoir4hIh4i8I4fnydl7mKvnYQEReUFEVuazDeNaEIjINhHp9n+AzX4nUJ3t86jqQ6o6P4P2XCYiD2f7/IXAeL62LPA14IeqWq2qf8zVSTJ9D/OB/35sEJEuEdkrIjeIyISE7WkFjYicJiL/EJFDInJARB4RkRPHrPE5RlWPVdUH8tmGcS0IfP5JVauBZcBy4EvJO4hIcMxbZSkljgReGMkXx8O7KSKfBb4JfB6oA96EuSfrRCQ8xHdrgb8CPwAmATOArwK9GZ47a/dPDOOyzxyXF5UKVd0F/A+wCPpMLB8XkVeAV/x154vIehFp9UcgS+LfF5HjReQZEWkXkd8D5QnbVorIzoTlWSJyh4jsF5E3ROSHItIE/AQ42ddQWv19y0TkehHZ7mstPxGRioRjfV5E9ojIbhG5PN31icglIvJU0rqrROTP/ufzRORFv/27RORzaY4zV0Tu89vdIiK/TRq5DefaHhCRKxK+O0BrEJHvicgOEWkTkadF5PR015fUxk0icn7CctBvzzJ/+TZ/1HlIRB4UkWPTHOcwLUYSTG+DPRsRmSIif/XflQMi8lCqTkJEtgBHAX/x702ZiDSKyJ/9770qIh9O2H+NiNwuIr8RkTbgshTHTPksU7yH2/z353kR6RSRG0WkQUT+x//uPSIy0d93tn/tH/HftT3p3hF//zf5v5FWEXlO0pg2xHTkXwU+qap/U9Woqm4DLgZmA5emO4fPPABVvVlVXVXtVtW1qvp8mvMddv9EpM6/9j3+/bpWRAL+/gER+U//XX9NRD7h34egv/0BEfm6iDwCdAFHicgCEVnnP7/NInJxBs8m7fviP6dz/M9lIvJd/xns9j+XJT5fEfmsiOzzr+dDQ9y/zFDVcfsP2Aac43+ehRmV/R9/WYF1mFFGBXA8sA84CQgAH/S/XwaEgdeBq4AQcBEQBa71j7US2Ol/DgDPAd8BqjAC4zR/22XAw0lt/A7wZ78dNcBfgG/4294KNGOEVxXwO7/dR6e41kqgHTgmYd2TwHv8z3uA0/3PE4Flae7Z0cAq/7rrgQeB747w2h4ArkhYHrAPphOYDASBzwJ7gXJ/2xrgN2na+BXgtwnLbwM2JSxf7t/LMuC7wPqEbTclPLdUbe67v0M8m29ghF/I/3c6IEO9h/7yg8CP/fu3FNgPnJVw3VHgHZiBWkWK46V8liS8hwnnfQxowIyk9wHPYN71cuA+4Bp/39n+td/sP9vFfrvOSWjXb/zPM4A3gPP8Nq7yl+tTtPWtQAwIptj2K+DmwZ43UOsf+1fAucDEIX7zh90/4E7gp/51TQWeAD7q7/8x4EVgpn8v7/HvQzDhHd4OHIt5T+uAHcCH/OXjgRZg4RDPJu37wsB+6mv+M5uK+f39g/4+a6V/L7/mH+M8jHAa9J5k1FeO9gCF/M+/wR1AK6Yj/zH+D8t/2Gcl7HtD/IYnrNsMnAmcAewm4YfuP6BUguBkzA8o1Yt/GQM7QgE6gbkJ604GXvM//wK4LmHbPNIIAn/7b4Cv+J+PwQiGSn95O/BRoHaY9/AdwLPDvbaEH1FaQZDiGAeB4/zPa0gvCI5Ourbfxq87xb4T/HtW5y/fRAaCIINn8zXgT+meRYr3MHFA4gI1Cdu/AdyUcN0PDnG8lM+S1ILg/QnLfwBuSFj+JPBH//Ns/9oXJGz/D+DG5OcBfAH4r6Rz3w18MEVbLwX2prmO64B1GTzvJv+57cR0hH8GGtLsO+D+YYRgLwkCFXgvcL//+T58oeAvn8PhguBrCdsvAR5KOudP6Reo6Z5N2vcl6f3YApyXsG01sC3h+XaT8PvDCPc3DfUODvWvFExD71DVCap6pKr+i6p2J2zbkfD5SOCzvurWKsa8MQto9P/tUv/O+7ye5nyzgNdVNZZB2+oxI/mnE875N389/nkT25junHF+h3nJAd6H+ZF3+cvvwowgXheRv4vIyakO4JsObvHV2jaMcJkygmsbEhH5nBgzzyH/2usSzpUWVX0V2AT8k4hUAhdgrj2u6l8nIlv89m/zvzbkcZMY6tl8C3gVWCsiW0Xk6gyP2wgcUNX2hHWvY0bZcXYwOBk9S5/mhM/dKZaTnSeS37fGFMc8Enh30m/lNGB6in1bgCmS2lY/3d8+KKq6SVUvU9WZGO24EaPppSP5dx0C9iS09aeYETcc/htLde+Tj3dS0rW/H5jmb0/3bDJ9XxoZ+DtPfgZvJP3+ujj8GQ6bUhAEg5HYse8Avu4Ljfi/SlW9GaPuzRARSdj/iDTH3AEckebF16TlFsyP8diEc9apmdzGP++sDM4ZZx1QLyJLMQLhd30nVn1SVd+O+QH8Ebg1zTH+r9/OxapaixnRxa97ONcGZkRdmbAc/7EgZj7g3zC24omqOgE4lHCuobgZc41vB170hQMYAfh2zMiuDjPSJc1xB7RPRKYlbBv02ahqu6p+VlWPwgiifxWRszNo925gkojUJKw7AtiVsJzqXvZvzPxZjoTk9213in12YDSCxN9Klapel2LfRzEj8ncmrhTjvXcucO9wGqeqL2G0g0WD7ZbU1l5gSkJba1U1Pm+0B2MWipN4/emO9/eka69W1Sv99qV8NsN4X3ZjhE2cdM8gq5S6IEjk58DHROQkMVSJyNv8H+yjGJX0UyISEpF3AivSHOcJzMt1nX+MchE51d/WDMwU31NCVT3/vN8RkakAIjJDRFb7+9+Kmexa6I98rxnsAlQ1CtyGGX1MwggGRCQsIu8XkTp/nzbAS3OYGow57ZCIzMB4egz72nzWA+8UkUoxE7D/nHSeGL6pSUS+grEHZ8otwFuAK0kQeP5xezF25UqMYEvHc8CxIrJURMoxZgVg6GcjxrHgaH9wcAhj7kl3T/tQ1R0Ys+I3/Pu3BHNfMvLRH+azHAlf9p/XsRg7+O9T7PMbjDa22tfAyv2JzJnJO6rqIcxk8Q9E5K3+72c25t3eCfxXwu6Of6z4vzIxE7OfjR9bRGZhBgCPZXIxqroHWAv8p4jUiogjxiHiTH+XW4FP+892AsbsNRh/BeaJyAf8awmJyIki0jTYsxnG+3Iz8CURqReRKZj5sJzHb1hB4KOqTwEfBn6IsVW/iu+xoaoRzIjmMuAAxk54R5rjuMA/YezM2zEv+yX+5vswE9Z7RSSuEn/BP9djvinjHmC+f6z/wajA9/n73JfBpfwOMxq+LUmF/ACwzT/HxzDqbCq+inG1PQT8d+J1juDavgNEMELiVxhbfpy7MaaWlzHqbw9Dm0T68H/gjwKnMLCz+rV/vF2YScC0HYaqvoyx3d6D8RxLjoNI+2wwczD3YITmo8CPVfX+DJv/XoymshszkXmNqt6T4Xch82c5Ev6OueZ7getVdW3yDr4wezvw7xhBvgMzYEjZn6jqf/j7Xo/pHB/3v3O2qia6gb4Xo4XF/23BzAWdBDwuIp2Y57kR41yQKf8L4/DxIua3fTv9ZqyfYwTF88CzwF2YAYqb5lraMQOQ92Ce316Ma2yZv0u6Z5Pp+3It8JTfng2Yyf2cB0DGZ60tFksJ44/SXwNC2ZoDKkZE5FzgJ6p65JA7jyOsRmCxWEoWEakQ4/sf9E2h12C0tJLCCgKLxVLKCMYcehBjGtqEscuXFNY0ZLFYLCWO1QgsFoulxLGCwGKxWEqcostsOGXKFJ09e3a+m2GxWCxFxdNPP92iqvWpthWdIJg9ezZPPfXU0DtaLBaLpQ8RSZuixpqGLBaLpcSxgsBisVhKnKIzDaUjFovhedlMuVJ6OI5DMDhuXgmLxZIh40IjaG9vJxKJ5LsZRU8kEqG9vX3oHS0Wy7ii6Id/sViMQCBAZWXl0DtbBiUcDtPV1UUsFrOagcVSYERdj87eGFVlQUKB7I7hi/7X7nme7bSySCAQsCY2i6WA8Dxl3aZmNu1pozviUhEO0DS9llVNDThOpuU7BmdcmIYs2WNg7R2LxZJv1m1qZsPOQ6hCeSiAKmzYeYh1m5qH/nKGWEGQQ7Zt28aiRYMVUsoPK1eutLEYFksREHU9Nu1pI+AIrqf0RF1cTwk4wqY9bUTd7GjvJWtTyaW9LZdY+73FUjp09sbo6o2xq7Wblo5eYq4SDAhTqsuYMaGCzt4YEyrDQx9oCIqnB8wSnqfc/cJefnT/q9zwwBZ+dP+r3P3CXjxv9FlYv/3tb7No0SIWLVrEd7/7XcB03O9///tpamrioosuoqvL1JK/+uqrWbhwIUuWLOFzn/scAPv37+dd73oXJ554IieeeCKPPPIIAGvWrOEDH/gAp556Kh/4wAd405vexAsvvNB33vgIv7Ozk8svv5wVK1Zw/PHH86c//QmA7u5u3vOe99DU1MSFF15Id3f3qK/VYrHknqqyILsPdbOvvRcQggEHEPa197L7UDdVZdkZFJbc0DJubws4MsDeBrD62GlDfDs9Tz/9NL/85S95/PHHUVVOOukkzjzzTDZv3syNN97IqaeeyuWXX86Pf/xjPvShD3HnnXfy0ksvISK0trYC8OlPf5qrrrqK0047je3bt7N69Wo2bdoEwIsvvsjDDz9MRUUF3/nOd7j11lv56le/yp49e9izZw/Lly/n3//93znrrLP4xS9+QWtrKytWrOCcc87hpz/9KZWVlWzatInnn3+eZcuWjfo+WiylRv6sCOnm7bI3n1dSGkGivS2RbNjbHn74YS688EKqqqqorq7mne98Jw899BCzZs3i1FNNffdLL72Uhx9+mLq6OsrLy/nnf/5n7rjjjj7X13vuuYdPfOITLF26lAsuuIC2tjY6OjoAuOCCC6ioqADg4osv5vbbbwfg1ltv5aKLLgJg7dq1XHfddSxdupSVK1fS09PD9u3befDBB7n00ksBWLJkCUuWLBnxdVospUYurQhD0dkbo7GunIbacgBifh/VUFtOY105nb3ZqSpaUhpBZ2+M7ohLeShw2LbuiJs1e1siyV44IkIwGOSJJ57g3nvv5fbbb+eHP/wh9913H57n8dhjj1FeXn7Ycaqqqvo+z5gxg8mTJ/P888/z+9//np/85CcAqCp/+MMfmD9//mHft1gsIyNXVoRMqCoLUlkWZG59NUdMqqQrEqMybDQSEbJmGiopjaCqLEhF+HAhAFARDozqpp5++un88Y9/pKuri87OTu68805OP/10tm/fzqOPPgrA7373O0477TQ6Ojo4dOgQ5513Ht/5znd47rnnAHjLW97CD37wg75jrl+/Pu35LrnkEv7jP/6DQ4cO9Y3wV69ezQ9+8APiVeeeffZZAM444wx+97vfAbBx40aef/75EV+nxVJK5NKKkAmhgMOChhpe2dfO068fZMPOQzz9+kFe2dfOgoaarJmoSkoQhAIOTdNrcZNUOtdTmqbXjuqmLlu2jMsuu4wVK1Zw0kknccUVVzBx4kTmz5/Pj370I5qamjh48CBXXnkl7e3tnH/++SxZsoTTTjuNb3/72wB8//vf56mnnmLJkiUsXLiwb6SfiosuuohbbrmFiy++uG/dl7/8ZaLRKEuWLOHYY4/ly1/+MgBXXnklHR0dNDU18ZWvfIUTTjhhxNdpsZQScStCKuJWhFyjAijx/8xf9ddniaKrWbx8+XJN9IGP5xgKhzMz6YxFlF4xM9z7abEkUqxu2emIuh4/uv9VUnWTIvDxNx+d0+tMPL/rKVHXIxRwCDgy7POLyNOqujzVtpKaIwBwHGH1sdM4a8HUcfXCWiz5ZLwOsOJWhPgcQRzXUxbPrBtx35GpwOzsjdHREyPgCKGAM2B+M5vzmiUnCOKEAk7WJ4YtllIlnxOquWZVUwPAACG3eGZd3/rhEBeYG3e10tYdo7YiyKIZE1IKTM9T/rGlhQ27W4nF+gPJ5tZXIyKjntdMJGeCQER+AZwP7FPVRf66ScDvgdnANuBiVT2YqzZYLJbcM9SE6lkLpha11p1NK8JdG/bwp/W7aO2KEnE9yoIOW/Z3op7y1sXTB+y7blMzL+5up766nOa2HuKBZABzplSPSiNJJpdP5ybgrUnrrgbuVdVjgHv9ZYvFUsQUwoTqWBC3Ioyk8/U85a4Ne/jevS/z5LYDbNrTxmstnew42MXLze389fndAzyQEoXrUVOq+uIIPE9p6exlYWPtiDSSdORMI1DVB0VkdtLqtwMr/c+/Ah4AvpCrNlgsltwTd8tONaGaTfNFMbNuUzNPvnaA/e29dPa6xHzPxc7eGOWhGJ6ntHZFqK8p71sfj3kSEebWV9NUHaW282WaK47mlLmTszr3MtZPqEFV9/if9wIZiTQRWQNcAzB9+vTBd7ZYLGNKriZUxwvx0b0AHT0xYgkC01XoibrsPdTdFzUMA4WrE+ng3Rs/xsSenTi4eBIguG8BfOguCB0efDoS8vaE1PitZuS7qqprVFVUVRobG3PcsvzzwAMPcP755wPw5z//meuuuy7tvq2trfz4xz8e9jnWrFnD9ddfP+I2WiyJrGpqYPHMOkRMxybCiCdUC5Wo69HaFTksiCzd+jjx0b2q4qbo8VyFSEwHCIJQwKFpWjVzWu7nsmffzdTuVwlpBEEJihDa/xL88m1Zu7ax1giaRWS6qu4RkenAvjE+fz9uFHrboawGAqGxOaXrEgikjmxOxwUXXMAFF1yQdntcEPzLv/zLaJtnsYyY8eyWnc419uz5U7l3874hXWbjo/v27ih9sWEJBDAxCZ1J8yyrnKfZ5r1IZewQiqAooi7qQVRCBFpexulogeopo77GsX5SfwY+6H/+IPCnMT4/eB5s+is8dD08/F3zd9NfzfpRsG3bNhYsWHBYyunZs2fzhS98gWXLlnHbbbexdu1aTj75ZJYtW8a73/3uvqRyf/vb31iwYAHLli3jjjvu6DvuTTfdxCc+8QkAmpubufDCCznuuOM47rjj+Mc//sHVV1/Nli1bWLp0KZ///OcB+Na3vsWJJ57IkiVLuOaaa/qO9fWvf5158+Zx2mmnsXnz5lFdr8WSitFMqBYq6SqEXb9u84D1MVd58rUD/O2FvQO+HzedhYIOoaAQSDbtixGk67e39ieyc6M4+zZyVEUPDh6C4CDmr3q4nuLGotD8Atkgl+6jN2MmhqeIyE6Mjf864FYR+WfgdeDi9EfIEZvvgt3rwQlAqMKI593rzbam80d36BQppwEmT57MM888Q0tLC+985zu55557qKqq4pvf/Cbf/va3+bd/+zc+/OEPc99993H00UdzySWXpDz+pz71Kc4880zuvPNOXNelo6OD6667jo0bN/blJVq7di2vvPIKTzzxBKrKBRdcwIMPPkhVVRW33HIL69evJxaLsWzZMptqwmIZgnSusQDP7jjIitmTUVW27O/oKxyzYXcroLz12Ol9msGqpgZcT1n7wh4OdEbBU1QhEDCCYWJVmNcPdLNuU7OJu+hth0g3bsUUPJykhNNGP4jhwJQmsmHPyKXX0HvTbDo7V+ccEjcKzRuMEEjECZj181aPykyUnHL6+9//PkBfx/7YY4/x4osv9u0TiUQ4+eSTeemll5gzZw7HHHNM33d/9rOfHXb8++67j1//+teAKTJfV1fHwYMDwzDWrl3L2rVrOf744wHo6OjglVdeob29nQsvvLAv5fVg5iaLxWJIl7E46nr0RDwiMY+dB7vY196LI0IwIMRiHs+83krAcfqC6RxHOG/xdFzX4+cPbWHr/i5cVcocl6nhCPNmTiYcdPrjLspqIFxBNOLSG6yjKnaQfgOOgKe0VMyhJjSBCVm4ztLy6/KlLKGKw7dFus32ykkjPnyqlNPQn0JaVVm1ahU333zzgP0GyzI6XFSVL37xi3z0ox8dsD5eMc1SeKRLNzDe8vYUI+lcY0MBh/KwyfnT0mGEQJxgwKEiHEgZTHfuouk8s+MgoTc28SHuoI1qKkITmNo7gZaWebxYe1p/2oiGxYR2Pcv6hn+iqfkvVMYOIXi4BNhfdRR/WPwTPl7okcUFiS9lU/oqhSvM9lEQTzl98skn96WcjqeCBnjTm97Exz/+cV599VWOPvpoOjs72bVrFwsWLGDbtm1s2bKFuXPnHiYo4px99tnccMMNfOYzn+kzDdXU1NDe3t63z+rVq/nyl7/M+9//fqqrq9m1axehUIgzzjiDyy67jC9+8YvEYjH+8pe/HCYsLGPLaCchxzuFIAjTucYCHD9rIpGY11dHGMBTpaG2HEckZS6g+596hi8+/RaC9E8Mu11wr3MxDWz2aww0mQ3zzyMATDnwKJu9cwHoDVSwceq76CybmFXX3NISBIEQNCzunyOI47nQuHTU3kPxlNOXX345Cxcu5MorrxxQX6C+vp6bbrqJ9773vfT2mlDxa6+9lnnz5vGzn/2Mt73tbVRWVnL66acP6NzjfO973+MjH/kIN954I4FAgBtuuIGTTz6ZU089lUWLFnHuuefyrW99i02bNnHyyScDUF1dzW9+8xuWLVvGJZdcwnHHHcfUqVM58cQTR3WtltGTLj/Ps9sPEnSccZm3JxMKLYFdulxDZ8+fytpNzX4uII9gwKGhtpyjphgLwIBgOjdG9N6vsfIf3zus0w0CZ3fcyhN1V7DQ2U4IF3DAcaDpfI46+i1s37iVjfs9OmMmx9Di6dmNLC65NNR4npkwbt5gzEHhCiMc5p9nbvwI2bZtG+effz4bN24c8TEKAZuGemxIl97YU+WJbW9w4pGTDxuBjkXa40Lg7hf2pg1Oy6cgTKeh3LVhN8+83kpFONBnIjqsvevWENu0DufAxoETvwqI+fPqrEs4emYDzulXpTRRj1ZDsmmoE/GlLPNWj3kcgcUSJ90kZCRmJiGjrkcgyakhV+VUC4lCTmCXLmPxW4+dTsBx0mcnjXTBzieQSIr8mgmXOYcdOGVz0pqoc5kxufQEQZxAaFQTw8nMnj276LUBy9iRbhIyHDSTkKk6u1zk7RnJKDOXtvt81BUfLUMG07XvhWgPbtlkgh270h4nUncUoYbFeRmYlq4gsKREVQ/zfrJkn3STkKpmEjKZbOftGYkdfixs98WcwC6EywQ6gBoGxOrWTDM5gZwwLqk7XReIHXW2MVHngaI3NjqOQyw2PtLcFgKu6+KMYq7Ekjnp8vN8btX8nOftSRctu25Tc1a/M1xyWVc8Z8R6Yf3N8Pdvps5WEK6EmSsI4rFN5hPDzAnE//UQ4Eszf0/lce8Y1TzlaChc8ZohwWCQ7u5uurq6CAQCdjQ7QlQV13VxXZdgsOhfi6JgMJNCLvP2jMQOP5a2+2xWBMspcceTjX+A1h0QDEN1A0w+5vBsBWd9Cee+a5nY/TA72ufTE3GJSpBf1vwLkYZFXLCkMa9Cblz84mtqaojFYnijzBdUyogI4XDYCoE8kG4SMFeTgyOxw4+l7b5oEthtvgt2Pg0dzUYIgJkPAJgyb0C2Ak8CrGv8GC/HLmTPrm1sj9QQrKhmUWMti2dOyLuQGze/etuBWSyZMRI7fD5s9wVXVzzSZTr6mmlmQrd5A6hrUtfEJ3hFjGCYPBeikb5sBX0xI8EK6o9sYpKn9ERdFjbm1yU2ju09LZYSYySFZEq6+Iwbg/uuhZ1PQLTHTPw2LAYnZD4ne/m4MYhF+rIVpDKrBRyhqizIy83trFrYkPf7N46fnsViScdICsmUQvGZlNx3LWx/1HyOVwTb/QzsesJkKKhuAE0wSweCZr3vCloMNZ2tRmCxlCAjscMXje0+m/jBYIdlLA6EoOsARLvN5DAYk1AsAhNmwYwT+lxBi8ElNv8tsFgseWMkdviCs91nm8TqhX4wWMrawMEKmHAEdO2HullmXmDS0XDshRAq69utGMxqVhBYLBYL9LmDunueJ9LdSbiiikD9/PQF4sMVsPR9RjsYIl1NobvEWkFgsVgsgPfSf7N1wyO0dMSIuh6hQC9T9j/JUeUTcbpaDs9YfMTJJlgMhkxXU+hmtcJpicViSUnU9WjtihB1bZxMznCjbN34GPvaoyimuIwC+9qjbJUjYNYKs1+0x/w94mQ460vDPk2h1nS2GoHFUqAUWl7+fJOVZHeJ9v8EM060u42Dhw4hzkAzkIhwsL2d6PlfJRSu6I8jiGsC4wQrCCyWAiVd4RrIbYGaQqgMlkhWBOIQdUg6tZxuLUtZCL5by+jUciaEK2HyUVm9tkLBCgKLpQDJR17+QtVARiUQ4xrAaw/B3g3Gzh/yy9Um5AOqqqygtXYeUzs2o9J/X0U9WmvnU1WZos75OCL/4t5isRxGPoKQxiK76HAZSiCmnTfxPJMB9KHr4cH/hCd+Cge2MMCZ3wkYDcGNEgo4lC36J/ZWzTdzBG4PCuytmk/Zon8qCM0ol1iNwGIpQMY6CCmxw/VUicQ8wkEn75XBRpzsbtOfYceTZvTvBIxmkJgQLk6kuy8f0KqF01knb+epXQdwe9oIlNcyf8akgnHxzCVWEFgsBchYByF19sbo6o2x+1AP+9t7ibmmGHt9TRmNdeV5qww2bIHoebDpL/DEz/qTwVXVm7xAQn9COPEFi58PCArfxTOXlMZVWixFyGhy+wzX5bSqLMjuQz00txn3yKDfATa39bD7UE/e0iAMu1DN5rtgxxMm8VvcK6ijGTRm8gHFE8KBiQVIURqyUF08c4nVCCyWAmUkI9TRTfimGHYPun5syDgq140am3+4wiR+iyMOIEYz6GoxqaMFaFyat9KQhYYVBBZLgTOc3D7pPGxcTzll7uS0wqSzN0ZjXQWq0NLRS8xVggFhak0ZjXUVeS0an7FA7G03Nv9QhckI2r7X1AcA8GImH9Cid8KcMwZNB1GKWEFgsYwTUnnYqCqvtXTw5Otv8Mz2A1SXhVJqCFVlQSrLghw9tYY5U6r9FAtmsliEgsiQmbY4fJyyGqMNKAMzgsbNRLNOhKYL8lYXuJDJ/9O1WCwDGGlAVyoPmy37O9jX3ovnKQFx0vrgJ09OB/y8OgWRIXOIYLA+AiGzfvd6XHGI1s0lNPEoArEekyLi2Lfn7xoKHCsILJYCYbQBXckeNq6ntHT04ojgBIRw0HSa6VxCCzZD5ua7TPBXmmCwRLx557J+x0F6dzxnagWEKiibdRxLF5xvPWMGwQoCi6VAGG1KieRRfdT1iLmK40BDbTmO9AuTVD74Bec+6UZN8Zc96w8vDBMPBvOLw8dZ99J+NugKQjNPIOR2Eg1UEdUA+1/aXxC1gQsVKyQtFvKf4XPEEbRJJLqcup4SDAoNteUcNaVqwH6DBaXl3X0yOSr4tb9Dy8scFkwQDwbzGRAU5wTpDdXhOcFh38NSxGoElpKmUPLrjDiCNonkUf0/trTw4u52RAqzMtYA+vICPQh7N5pRf3mtCQZLFRWcEAwG2buHpYgVBJaSJl8ZPpPJdkqJ+Kj+rcdOJ+A4hWf3TyRxMri3A3Y/a9w/Jx/TXxy+fc/AqGDPNXEACWahYqgNXKjYO2MpWfKR4TMdifb9eNvi5x7N6L3g7P6pSJwMloBx92zfY7ZNmdfvCtq2C7rboGpiymCwYqgNXKhYQWApWQrNlHD2/Kk8u/0gz+44SE/EozzscPysiZw9f+qoj12wBefj0cDxyeBguD8quKMZJs0126bMM9rAio9AxcSirQ1cqFhBYClZ8mlKSBUrcO/mfQQdhxWzJ/dl/1Q164ve4yXSlbq6V2I0MBiNIB4V7MXAjYBT0W8Kqh5cKBaFBlSAWEFgKVnyYUpINzm9cl79ADNVn5Yi5DUN9KhxY3DftbDzCVPvN1QOM1eYer+B4MBo4DiJUcHeyPICFawGVKDkRRCIyKeBD2Me8c9V9bv5aIfFMtamhHST04Vmpsoa910L2x/1g8H8esDbHzXrV60ZEA3cZx4SMSahhe+AOafbvEBjwJgLAhFZhBECK4AI8DcR+auqvjrWbbFYxtKUMNjk9GstnYSDqd1Vi9bjJdJlNIFUwWA7nzDbw5X9I/3EFBJxDcDmBRoT8vF2NQGPq2oXgIj8HXgn8B95aIvFAoyNKWGwUX8k5rFgWg1b9ndmxUyVtwL08ViAshpj54+bgw5rYI/ZPvko09k3nW+ihOPftRrAmJIPQbAR+LqITAa6gfOAp/LQDotlTBlqcvq8xdN54OX9ozJT5S1ALlViuMnzUgsBMOtrkibAAyGonJS7NlrSMuaCQFU3icg3gbVAJ7AeSF2l20dE1gDXAEyfPj3HLbRYcsNQk9NlocCozVTZCJAbkTaRKjHcvk1QPgm69g80D3kuHHHyQO8hS17JiwFOVW9U1RNU9QzgIPDyEPuvUVVRVWlsbBybRlosWSA5h1Em5SdHmutntPmKPE+5+4W9/Oj+V7nhgS386P5XufuFvXh+mci+a4n0mmRwbtR8MTkWII4TgIYmkwIajDkIjBA460vDujZLbsmX19BUVd0nIkdg5gfelI92WCy5YjATTa4mp0freZROm/BUcUTYvLOFI/eupSG6k4YqOHJaPc60xaZjT4wFSCTWC6ddBcHy1HEEloIgX64If/DnCKLAx1W1NU/tsBQReZsAHQFDmWhyMTk9kgC5+D0NB5202sRdz+3kzfIMq1rWURvZi+uE6OqexFannKO99SbwKzkWIE48MVwgZCaGLQVJXgSBqp6ej/NaipNCyRCaKfnKYTScALnkeyoCW/d3sGBa7WGZSmc138+M8i1UR1vwnBACVEXeoHPfy7j1ywjs3wT1C2HvhsPnApISw1kKk8IeVlks9I+uVRkwul63qTnfTUtJ3ESTiriJJpHh1kIYbP9M5iDg8HsadBwOdkXZ2tLZv5N6zNl3D+d2/5lZnc9T17OLimirqQsgQkVvC9FI1JiFZp9qOn3BVAYbQTTwWJDvuhOFShFGqVhKiULKEJopmZpohqvpZLL/UAFy8Y5w466BWkPAEabWlNHc1sP8WmVCZAcTu7ZR3/kylQEPzylDxSHkdgHQHZpAEJeQxIzNv7yuoGMBik2rHGusILAUNMWYeiFTE81wXT0z2T9xHiXxviR2hAc7I2zc1ca0ujLm1lf3mYKOmeBwxc41zH5jG0GNEvR6iVVMJVI5ma6oSzRQQdjtJuR20R2spby8jIATNCki4p1+gcYCFErdiULFCgJLQVOsxUaGymE0XE1nqP1Xzqs/LBgtccSb2BHWVoQIBoR97b0EvR6W1HbR0LmZN73+M2p69xEOB9FAEPEEiexHpRspm0GXTEAVytxOagIxahuOgpknFJz5J5li1CrHmsL8FVksPsVabGQoE81wNZ2h9r9rw56+9BTJI96zFkwd0BE6IkyvVM7edQNzWrYyPdBKwI1SFXsDAmWICOLFzGRvIIT2dlLd0EBdpBXPDRIITMc5/v2w5N0QLMv2rcs6xahVjjVWEFgKnmIuNpLOTXS4ms5g+4eDDltbOggkJWiLj3iXHTGhryN03AjHtKzjoj13UB3ZRVQdwhoBCRAghqgAFSCCAjHXxXM9XjgUIlreRH21x+wlZ8KiC0dxV8aWYtUqxxJ7BywFz3gsNjJcTWew/Y+qr2Lz3nbKU8zNxr2XKkLCnJb7mbd/HXW9u5nYsw1PggScMFUSQwIOEnOMFuB54DjE1CFKEJEYGijHFYfnZT6vcBKrc3NbckKxapVjiRUElqJhvBUbGa6mk27/lfPq2X6gK+2Id0JlmDNjj6AHHqMmuh9BcVDEixIICo4TAHUhEAa3FzAeoq5CzCljb+1SHj/yw0QDVXhOENnbwVlNXlF1oMWsVY4FoqnengJm+fLl+tRTNlmpZfww3IjpVPvf/cLe1IXvG6tYvfuneC/eSaSnG6JddFFB2OvCEcERQQIhwsSQshroOQgSwlOXLi1jV+1x3DX/60ZI+PREXa5cObcohXIxRadnGxF5WlWXp9pmNYJxSim/8MXGcDWdVPsnF76vDUY4pSHG2TsfgB2P43gu5WXlRN0eKl0T8CUiKEovYRChDIGGY+HI03DLp/Dr9lPoDdcedv5itquPN60yWxTn07SkxQbOlCbxwvenNJZx6rbvMbVjG6FDEbzoDqioBQmgQJQgjkQBodepJOR1IZ5LW3giExadR+i0z0DFBEKBEHN9LcPa1cc/VhCMM2zgTOkRdT027TrA6a/fwPz9d1Mea0PFIeqU48ViaG87EgihgQoiUkbYgYAXoSs0kdaKJeyoO4H1E9/Kh89YPGC0bO3qpYMVBOMIGzhTYvhlITvdMk7c+mNmtj9HudsB4iBAWawLhxgQMvtWTkEibcQkTETKeG76u3lp6nl4gTAh4TBzz3j01rKkxgqCcYQNnCkRYr2w8Q44sAViEWpwmHfoEWJOOYKH+rkk1XHAE3BdcASnrAaXGrp7ethVezwvTn8HMLS5x9rVxz9WEIwjbODMOCdeF3jjH6B1BwTDUN1AoGIKVdpOl9svBOJoIIRU1JgMoZ7LhOoqDk5Zxj8mX05P1Jp7LIZh9QwiUq+q+3PVGMvosIEz2aFgPa42/Rlefwza9xghAKbqlxejIlyORCNEAlWE3U4EY9oJOgGomwWzToITr8CpmcbccCUfL9RrtOSFjASBiJwE3IqpXzBLRJYDH1HVj+SycZbhYyf4Rk46j6uV8+rpjrr56zRjvfD87bD+N+Zz+x6T5rlyMohA1xtI7XQqWndQNulI9NAunEiHyRdUVguzTzc1ggP9P/dSMfcUrFAvMDIKKBORR4APA79V1eP9dS+o6rE5bt9h2ICyzLA/gOFzd5K7pKry6v4OHIEZEyrH3hU3Phfw8v/AwR3Q2Qzhaoh0AQrhKqiaAm4MZq6AXU+ABMz3AmVQPx/e8jUoPzwWYLxj3agPJxsBZWFVfTGxhB0QGXXLLDmjVEZ82SKVx9XWlk72t/fiCMyeXD12rrg9HfD87+G1B+DQLujcB6FKEAcinZjiwALRLlDPjPSDITjlUzD3LFskHutGPVwyFQS9IlKNX55aRBYCPTlrlaWkyYc2k+xx5an6QkCIuR5R1yPgBHLritvTBmu/Atsfge5DpuRjqNKYfzRBACBGG+hth2gPTJoDM/y6AI5TtEXis/XcrRv18MlUEHwdWAs0ishNwFuBS3PVKEtpkk91PtnjKhLziLguDkLAkQEdR9ZdcXvaYO2X4eW10NMKGgMnbPr8SAegUDkJVfCClTjRDiRcbeYIll4Kiy+CUOHXBUhHtp+7daMePhkJAlX9HxHZDKzGvJ7XquqrOW2ZpeTIhTqf6Sgz0ePKEdhxsIt9bb24rkdNRYhtb3Ry1JQqRCQ7rriRLji0E57+Fbx4J3S8AV4EHL/zcnsBgVA56kbp6YkQVeVQaCISqCJatZgjl7wZZ/E7RteOAiDbz926UQ+fjO+Iqm4FbshhWywlTLbV+ZGMMuOeVX95bjd7WnuoDAfwNMDEyjDNbcYSOnty1ehccd0Y3Hct7HwCWl4xNn83AuKPUD3Xv/CgiRtQj6gKXRokgEsA5VBwCpudhbwsK4qqLkAqcmHGsW7UwydT99HTgG8Ac/3vCKCqOjWHbbOUENlW50cyynQc4awFU9m46xAzJlQQdIRtb3TS0tGL6yr7O3o4b/G00bni3vtV2PYIOEGI+dNsXgzwzGfzyzL/RFAniBuL0hOcQFvZNF6bfDqvTFmFFwgXZV2AZHJlxrFu1MMjU43gF8D/Bp4G3Nw1x1JIjOWkbTbV+dGMMjt7Y0RiXl/HdPTUGuZMqSbqeriecsrcKSObr4j2Gk+gjXcYTx9V4+bphPzJYAUcUyAGMR5CKJ5Txv7yWWxsfFefAIgzHuzduTLj2DxJwyPTu3xQVW/LaUssBUM+Jm2zqc6PZpSZqmMKOELACSApErMNSbQXXrjTxAK8sQW6D0Kw3PxTNR2/BMwEcSAMGjXmIQlA5US8+efzx8CluIHyww49HuzduTbjWDfqzMj0LfqdiHwME13c5zaqql05aZUlr+TLBztb6vxoRplZ65j6gsHuhkM7oKMZghWAY+YEAIJlxjwUKDOeoY4DhKByCsx7K6z6GqGKWuaP87oA1oyTfzIVBPuAnwM/8pcF8+oePuSyFDX59MHOljo/2s58VB1Tbwc893vY+gC0+cFgwXIzyo92mcCvWAS8qIn47faIu4cyaQ7UL4BV/8cUk8lGe4oAa8bJP5kKgm8AK4FnVNXLXXMs+aYQfLCzoc6PpvMcUccU9wbafBd0tUKs2wR9iR8B7PZCqBIvVIk6ZUikHceNGgGw4G2w/Aqom5EyGrhUOkprxskfmQqC3apqE/yUAOPFBzsbneewOqb7roXXHzGdfiAAMQ8i7YBAxQQU6IlEiboeh0INUNaAN7WJWcedg3PcRYA/Od8VSdtW21FackWmv+p7ReSbwO8ZOEfwYk5aZckb480He0w6z0iXiQsAPM8FFQQHEUxlMM+jxwvSRRkBenFUaS+v56XQ8bwUOIVVNkGaJc9kKgji6SQuTlinQHEmNbEMyni3SY+aSNfAxG7te/Ei3bRGhECvhyqEPSFEjKDjoMEyIr0uPaEa2quOYevE0wbEAri6lxd3t9kEaZa8kWmKiTm5boilcCgVm/SwSYwKjvZAqNykfz79XzkYDdAViVERqCQU6yTqlKNeD6ouEp7ALmazbfIZh8UCdPRGeX5nK2XBgXMyNkGaZSzJ2ODrZxx9s794r6q+lJsmWQqF0ZhVxlU9hLgG8OSNsOspkw8o5Pv1b3+U6IPX82pZEzN6nqc7WAdA0O3CdcIcCE+nfsUV3NNyHK5z+L0MOg4xV0k19TIeAsYsxUGmKSY+AFwH3OWv+qKIfEFVf5uzllmKknFVEKSnDdZ+CfZvNoFhB7ZAeQ1MONJ4AwE4AXTHkzw49Wucob9mesdGok45XcE69lXN497GK/lo07HMf+1AynmXJTPreLm5o+gn5y3FTaZv2eeAE1R1L4CITAPuBqwgsAxgXBQEiQuAl9dC7yGQIIQrTBRwTxu0vg4TZ/ftHoz1MNlp4x9zPkkw1kNldD9doXpiwXLK/GjkweZdAk5zzifnx5WGZsk6w8k+ujfxc1K1Moul+AuCRLrhnq/Cpj9Bd6uJ+nUCxh20u92kfwhVmIIwnudHAoMTrmDWEXM4sDcKwXLagrOAwzvzdPMucSGxcdch2rqj1FaEsjY5P640NEvOyPRXuUVEvioijf6/a4CtuWyYpfiIB6OlIm7vHg5R16O1K0LUzXEMY7QXnv0t/PI8eOEP0NXiJ39TkxnU7YFg0PjJea4RAvE0EZ4LM1dw9uI5LJ5Zhwj0RF1ESNmZx+ddUgtETfo7euIamioDNLR1m5qzdg5L8ZOpRvAx4PvA85i39B7goyM9qYhcBVzhH2sD8CFVtaUvR0EhqP7ZCkYbs1Gs58FLf4VHfwQHt0HXGybrpxszSd9Qsxzv/J2gKR4f7e6vG3DEyXDWl0bladVvTuufnM+GOa3oNTTLmJGp++g+4D3ZOKGIzAA+BSxU1W4RudU/9k3ZOH6pUUiqf7aC0cZsnmHzXbDhNlMpDL+9noepDeCYdar92wIhmDALZp4EK65IWSB+uJ5WueysCyFdiKU4yOgNE5EviMikhOXJIvL5UZw3CFSISBCoBHaP4lglTTZU/2yaYFY1NWRkIhmsLYN1jKNuoxuFrgPGJXTPepMVVH17vwg4fnF4dY0GIH4n6rnGY2j26fCWr5kC8SnyAg33XmbbnJZIXENLhfVIsiSS6ZvwXlX9ZnxBVd8QkfcB3xruCVV1l4hcD2wHuoG1qrp2uMexjH40mQttYjgmklTmrJyNYj3PaADNG8yksADNm3yTj3+uQJlfLCZIn33LCUHlBDjGpIVOzAo68PAju5e5zO003tKFWHJHpm9Zqjd5RG+oiEwE3g7MAVqB20TkUlX9zSDfWQNcAzB9+vSRnHZcMtpOM5cmmMFMJIN1muk6Rk8VEQgHh9l5uVHj5fPaQ7B3gx8MVmFG+L3tJjFcuMqkjy6vg55DJlOoBGHCTDjytEEFQJyR3stcd9Y2XYglEzLtzF8RkX8FvoMRClcBr47wnOcAr6nqfgARuQM4BUgrCFR1DbAGYPny5dlzqShyRjOazOdE4lCdZmLHqKpsbemkua2HiZVhfvbg1sy0lkQNoKcD9jwL1Q0w+RjfBBSA2unQfQDKJ5rvRDrNZHDZNCMAVl8LZdVDXs9o72UuO+vRamiW0iBTQfApTEf9fzGePv8APjDCc24H3iQilRjT0NmATXE9AkYzmszXRGImnWZix7hpTxsHu6JMrSljbn115lrL5rtg93rT4TsBoxm0+6EwU+aZv5OPMaYgERAPwjVQPRUWXwQL39EXJzAUo72XY5HbaaQamo01KA0y9RraDZwlIlX+cudIT6iqj4vI7cAzQAx4FvjZSI9X6ox0NJmvugOZdpqrj53G6cdM4fv3vkLQcQYIjpQj7e5DcGArTPIncZs3JNj+w8bjB8zk8OS5ZhJYBKYvhlM+DZEOMy9QOal/3wzJ1r3MV72BcRENbhkVw0k6NxeYCwTjUcWqetegX0qDql6Db/O3jI6RjibzNZE4nE4zEjMpnZO1B0gQGmHgtsth73ozug+WQX0TTJ4H5b5ZxwkYs1D7HmMyikX65wkalxrBkcIDKFOKeVLWxhpYIPOkc9/ABIBtAuK+bkp/EjpLnhnJaDIfE4nD6TQzEhq3vA92P228e4JlZuPeDXDwNVj0rv4vTD7G/O1oNq6hghEC88/LynUV66SsjTWwQOYawbuBuaralsvGWMaWfNUdyLTTHFRoNFYRWvu/Ycu9mDGJGBNQeR0EgtB90EwSx7UCEZg0Fxa+HeacAWU1wzYBDUax1nAYL6VJLaMj06e8xwqB8Us2bNPD8TgZTqcZFw6bdx1Auw9QHgww54iZnLP357DtIcAzaSDAmIZ6DkHlRMCBmgaIdZq4gXBFvwaQ4STwSCi2usLFbNayZI9MBcGjInIzcBsDaxZb01CJMxqPk0w6TQdltTzBOa134LbuJOgITsc02P+yGf1LQkfliEkG53kQKoMTPmhs/73tWdcAxhPFatayZI9MBcGJ/t9PJqyzcwSW3HqcRLpg/c2w5V4C3QcIhHyh0boDOvYaO3+4Gnrb0b6YR0W8KDSeBBWmWhiVk1Id3eJTrGYtS/YYVBD45SkBPj4GbbEUGTnzOIl2w7prYM8GOLjFmHbKqk1lMMeBYNiM+ns78CbPI9L8EsFoJ6B4CG/UHEvDRTdmnGPdYig2s5YlewylEfy3/zd5Kkn8dUdlvUWWoqGzN0ZHT4yAI4QCA339R+RxEuuFjXfAU78wo35xzDq8/spgk+aAE0TDlfT29LAr0k63ziQY9Kikl+bJJ/L3Of/K4pcPWh94iyVDBhUEqjpnrBpiKS48T/nHlhY27G4lFlOCAWFKtYn+FZHheZzE00Fs/AMc2AZvvArBEAR9X3/PNfb93g5TKMYJcjA8g06N0BOBML1EpJwng8fyQPXlHGV94C2WYWF9wywjYt2mZl7c3U59dTnNbT2AsK+9F4A5U6qH53Gy+S7Y+bTx8ReMn3/MV0KDYYjETNSveuBG8XDYTx0bpq/mD63HMEVaaXUmEXHKocvlSE/7vJisqcNiGRo7XLIMm8S5gaOmVNFQWw4YLaGls5eFjbWZeZy4UWhvNnUB1DXLgVB/+gc3arSCULm/XkAgWtXAhomreL3+TLxQJfuCjUYIADHXCAHrA2+xZI79pViGTWI0qogwt76aOVOqiMQ8XPU4Ze7kzDODdh6Evc9CTaOJDhbMxHBPm+n4VY1nUPU0aDwezvw8TriOnQ++jqNQX1NGc1sPjp/2JBgQHBGaptdas5DFkiH2l2IZNqkqXzli3Eery0KpR+LxymButD8zqALltUYAdOwDjRnzz4QjzXrEzAmEq2H2afDWb0D1VELhMpqm1+J6OkAjicZc6mvKWHrEBOsDb7EMA6sRWIbNsKJRkyuDBcPwxpb+3D+JCeEkANX10NkCtTNMvYA5K2HJJf2pInwSg6BmTKhgzpRKjppSzXmLp1OWIm+OxWJJjxUEloxJTCMxZDRqX2WwB2Hvxv7KYNFu4xoKA+sCALTtgqppZv2ko+HYC02EcApsEJTFkj2sILAMyWBpJA7riD0PNv230QB6O2B3UmWwYNj862g2SeAcf2J4yjxTJ2DFR6BiYsbpIGwQlMUyeuwQyjIk8TQSqgxII7FuU3NfR9w3Gk+0/0sA3Jgx+7zxitkuvikoFjF5geJ4LkxfaiqE5TknUNT1aO2KEHW9lMsWy3jDagSWQck4jUR8MnjP+v7KYMGwSQkNAzWAuCkoGDamosTMoHkkWfMpDzn0xjzKAg49Mc+WcLSMW6wgsAzKkIVLOjuZsPW/4cCr0N1mXEFrZ/qmoPhE8F7j/eNGwKkwnkGL3gXzVhdUZtDkBHqv7mtnb1sP0+sqhlcvOQNsoXhLIWEFwTghVx1L2sIl6nFs24PU/O1ROLTTjO6r6kFCA4vEJ1YG85IqgzlOwWQGTdZ8XE9p6egl6Djsb+9lzpQqHJFRJ9SzheIthYgVBEVOrjuWRFfREC4ht5NooIojDzzEsfIqgc59RgiAXwbSLxSTWCR+0lxY+A6Yc3rWR//ZEoDJmk/U9Yi5JodSzPWIxLy+baMp4WgLxVsKESsIipyx6FhWLainftc6enc8B9FuJFTOkbqbabOOho5Y/zyAOIAYzaBjrzEVVU3MSWWwbAvAZM0nFHAIBuLRyg7hYH/bR5q+whaKtxQqVhAUMTnvWCJd0L4XZ/ezLAu8jjtnEtGYR0h7CWzfD4dC/UIgjheDibOhft6wXUGHw0gFYDoNIjlILuCYbKrxOYJ4CovRlHC0heIthYoVBEVMzjoWNwb3XQs7nzDRwF37oXYGgTlnEAgFQMuMOaizxR/97zOxANCfNC7uCpoDUglAT02yuY27WlMKwEw0iOQguaOnVjNrUqXxGoqOvoSjLRRvKVTsm1fE5Kxjue9a2P6ocfUMBM0k76EdJkr4qJX93kCtO6DuiP45gVgEJsyCmSfk1BU0UQCqKltbOtnf3kvM9UBgzpTdvGPpjAEmokw0iHTRytmahyjlQvHWS6qwsYIgS+TjRc9JxxLpMppAPBbACfR/btvVnyo67g0UKoO6WWZieNJcWPROCKZOC5EtEgXg1pbOvuyjwYADKK80d7BuU3NfBz9cE1pytHI2o5ezWSi+GDpX6yVVHFhBMEry/aJnpWPx5wKomWb+RntMDQAwo/1QJUQ6jWYQ6TB2/zzGAsQF4HM7Wtnf3ttnv/dUmVpTRjjoDOjgMzWhjUXHmo0cSfl+54ZDXBOLWw5dT62XVAFiBcEoybc74Kg6lsS5gHjnP/14CCSNfisnm7/RLr9oDHmPBVjV1EBXb4zHtx4AjJvn1BpTKhMGdvBDmdAqQgHufmHvmHaso9Ey8v3OZUrU9Xhx1yG2vdFvugsGHOpryggI1kuqgLCCYBQUkjvgiDqWxLmAuAaw6ykjIAj2m4REjBZwzGo4/aqCiAR2HOH84xrZ2tJJJOYRCjgDnkPiHMlQJrQHXt5fFB0rFNY7NxSdvTE27W3jYFc0wXQHzW09RGLWS6qQKIw3pkiJmxxSER+RFizJcwFxnIDxCGo8wSxHe8zfI06Gc64xo/8CSAcBpoNfNKPuMCHgenpYhbJVTQ0snlmHCPREXURg8cw6Vs6rH7RjzVaiuWwlriumdy4cdDjUE+0z3cVxRDjUEx0Qm2HJL1YjGAVF7Q6YPBeQSKwXVlzRP2dQMw3ClWPfxgzIdI4knQmttSuSU9/+XAe+JVJo71wk5lFXHuZgV2SAMPBUmVgRJhLzsApBYVA4b00RUtTugDXTUgsBMOvjnf/ko8a2XcNkuHMkySa0XHes2bbnF9M7V1UWpGl6Da/u66Clo7cvZcfUmjKOnlpdUEKr1LFPIgXD8R7JpjtgVolXCEtnzw9XwswV/XMEcTzXmIEKVANIx0gnX3PZsebKnl+w71wSoYDDwsY6XA/mTKkm6np917uwsbCEVqljBUECI1HjC65kYnKN4HAFNCxOnevnrC8d7jV0xMlmfQmRq441V5HfBffODULivY26EAxI32/KUjiIptKJC5jly5frU089lZNj3/3C3rQjw0LzHjmMeCzA7mdh36bDR/mNS6Hp/MG/W8BzAWNBtuMIoq7Hj+5/NaXZSQQ+/uajC7YDzzbFEPw23hGRp1V1eaptViPwKSa3vAEMyAvUBV0tUNsIs8/s1wCcgNEQ5q1ObyYq8LmAsSDb9Y+LyZ6fa2xt6cKmdN7EISgmt7wBxGMBwASCeS607oRtfx+4X6TbzBlYxpR0bqvWNGIpJKxG4FNMbnl9DJoXaLcpDRmPEg5XmInjcUqhmh6KyZ5vKV0KsHfLD0Wpxh+WF0ggXAW9HaYuQG8nVIb75wgKJBAsmxRL3h1rGrEUMgXYu+WPolPjU8UCVEyGsmrT6TvOwLxA45C4n74qA/z0121qznfTLJaiYcw1AhGZD/w+YdVRwFdU9btj3ZZkCk6N7z4EB7bCpKOgou7w7aliAQQonwDHvAVOK4y8QLmiaCf4LZYCY8wFgapuBpYCiEgA2AXcOdbtGIy8q/GxCNx2Od6e9WisBwmW40xfCu/+RX+h+DiDxQIkl5EcZ9jSjxZLdsh3T3E2sEVVX89zOwoK77YP0fv6k0Q9wdMATixK6PUnKLvtQzjv/e3AnQNBWLWmJGMBinKC32IpQPKtN78HuHmonURkjYioiOju3bvHoFl5pPsQ3dufIeIJipn/VSDiCd3bnzHmolTEYwFKRAhA/wS/6w2UBKmyj1oslvTk7ZciImHgAuC2ofZV1TWqKqoqjY2NuW9cHom2bEGjvSm3abSXaMuWMW5RYVN0E/wWSwGST935XOAZVbXuHQl0Vs5CJZhSQkclSKxyFhPGqC2F6pufSMFN8FssRUg+BcF7ycAslE2KoWOrmjCZrTVNTG97EU1IEieex77ahRw1YXLO21AsvvmJ5H2C32IpYvIiCESkClgFfHQszldMHVso4LBt5Q/Q+z9JQ8cmAhrFlRDNNQt5feUPmD8GAqxYauJaLJbskBdBoKqdQO6Htj7F1rG9ecEMPvncl9nSuYtpsd3sDTUyd9IMfrBgRs7PnSvf/GLQxjJlPF2LxQL5dx/NOcUYdPTte1+mpaOXiZPq6famMNERWjp6+fa9L3P1uU05PXe2ffOLSRsbivF0LRZLIoXVA+aAYssq2hWJ8eyOgwQdB0eEYMD/6zg8u+MgXZHctjfum5+Kkfjmj6cUEOPpWiyWRMa9IMh2x5aWSBe8sdX8HQX723vpiXgpt/VEPPa3p3YtzRbZ9M0fShuLuqmvsxAZT9disSQz7k1DOc8qmlgYJp7iYeaKEad4qK8pozzsmCiyJMrDDvU1ZaNrbwqSbd7ZKt04nlJAjKdrsViSGfeCAHJc7DteGMYJ9GcC3f6oWb9qzbAPVxkOcvysiTy57QDBBPfRmOdx4uxJVIaz98gGs3lnwzd/PKWAGE/XYrEkUxJvb86CjpILw/SdMNBfOnIEKR8+t2o+16/bzLM7DtIT8SgPO5w4exKfWzV/9G1OYChvqtH65hdljYc0jKdrsViSKQlBECfrQUfJhWESifaY7SOoBRwMOlx9bhNdkRj723uprynLqiYAY+dNlVNtbIwZT9disSRSUoIg66QqDBMnVG62j4LKcJAjJ+fmEY2VzbsYUkBkGhdQDNdisYwEKwhGQ6rCMGBKQx5xckFnAh1rm3chpoAYaVxAIV6LxTIa7HAmFW4Uug6Yv0Nx1pdMpw/GHAT9hWEKGJvC2cYFWCxxrEaQiOfB5rugeQNEuiFcAQ2LTb1fJ03HWMSFYUrZ5l2MEecWS66wgiCRzXfB7vW+K2iF8eXfvd5sazp/0K9GA+V0VsykKhCkWCoEl7LN28YFWCz9WEEQx40aTSCVK2jzBpi3OmUR+PGQf6YUbd42LsBi6ac0hn+Z0NtuzEGpiHSb7SmwdubixM6RWCz92Lc9TlmNmRNIRbjCbE/C5p8pbmyZS4vFUFr6rxs1I/uymsPNPIGQmRiOzxHE8VxoXJrSLGTtzMVNKc+RWCyJlIYgyNQbaP555m/ifo1L+9cnYe3M44N0cySjKUBji9dYionS6Kky9QZyHLM8b3V6zSEBm39mfDIaB4Dx4DxgKT3Gf081lDdQqqCxQAgqJw0qBOJYO/P4YzQOANZ5wFKMjH+NIO4NFEoxERz3BqqcNOLDWzvz+GI0gWY2SM1SrIz/t3IE3kAjIW5ntj/04mY0pU2LrSyqxRJn/PdacW8gL+kH6rlmfQbmH0vpMJrSpmNWFtViyTLjXxCA8fppXAoCRLvN30G8gSyly2gCzWyQmqVYKY0hyjC9gSylzWiS8ZVyIj9L8SKaygm+gFm+fLk+9dRT+W6GpQSwcQSW8YSIPK2qy1NtKw2NwGIZAaNJxleKifwsxYsdqlgKlqjr0doVKbicTYXaLotlpFiNwFJwFGp0bqG2y2IZLVYjsBQchRqdW6jtslhGixUEloKiUFN7F2q7LJZsYAWBpaAo1OjcQm2XxZINrCCwFBSFGp1bqO2yWLKBFQSWgiHuez+vobrgonNt1LBlPGOHMZa8k+yNUx5yiHkeAXHoiXkFE51ro4Yt4xUrCCx5J+6NE3Ckr+xn0HGYP72WU+ZOLpjoXJty3DJesW+xJa8M5o3zcnN7QXa2NuW4Zbxh32RLXrHeOBZL/rGCwJJXrDeOxZJ/8iIIRGSCiNwuIi+JyCYROTkf7bDkH+uNY7Hkn3wNt74H/E1VLxKRMFCZp3ZYCgDrjWOx5JcxFwQiUgecAVwGoKoRIDLW7bAUDtYbx2LJL/n4tc0B9gO/FJFnReT/iUhVHtphKTCsN47Fkh/y8YsLAsuAG1T1eKATuHqwL4jIGhFREdHdu3ePRRstFoulZMiHINgJ7FTVx/3l2zGCIS2qukZVRVWlsbEx5w20WCyWUmLMBYGq7gV2iMh8f9XZwItj3Q6LxWKxGPLlNfRJ4Le+x9BW4EN5aofFYrGUPKKqQ+9VQIjIfuD1fLcjDY1AsU5i2LaPPcXabrBtzxejafuRqlqfakPRCYJCRkRUVYuyeK1t+9hTrO0G2/Z8kau2Wz89i8ViKXGsILBYLJYSxwqC7PLVfDdgFNi2jz3F2m6wbc8XOWm7nSOwWCyWEsdqBBaLxVLiWEFgsVgsJY4VBBaLxVLiWEFgsVgsJY4VBBaLxVLiWEGQJUQk4NdX+Gu+2zIcRGSbiGwQkfUi8lS+2zMcirXkqYjM9+93/F+biHwm3+3KFBG5SkReEJGNInKziJTnu02ZIiKf9tv9QqHfcxH5hYjsE5GNCesmicg6EXnF/zsxG+eygiB7fBrYlO9GjJA3q+pSVV2e74YMk3jJ0wXAcRTJ/VfVzf79XgqcAHQBd+a3VZkhIjOATwHLVXUREADek99WZYaILAI+DKzAvC/ni8jR+W3VoNwEvDVp3dXAvap6DHAvQ9RyyRQrCLKAiMwE3gb8v3y3pVRIKHl6I5iSp6ramtdGjYyzgS2qWqiJFFMRBCpEJIipN14sCdyagMdVtUtVY8DfgXfmuU1pUdUHgQNJq98O/Mr//CvgHdk4lxUE2eG7wL8BXp7bMRIUWCsiT4vIR/LdmGEwXkqevge4Od+NyBRV3QVcD2wH9gCHVHVtfluVMRuB00VksohUAucBs/LcpuHSoKp7/M97gYZsHNQKglEiIucD+1T16Xy3ZYScpqrLgHOBj4vIGfluUIYMu+RpoeHX47gAuC3fbckU3yb9dowgbgSqROTS/LYqM1R1E/BNYC3wN2A94OazTaNBTVqIrKSGsIJg9JwKXCAi24BbgLNE5Df5bVLm+CM8VHUfxk69Ir8typhhlzwtQM4FnlHV5nw3ZBicA7ymqvtVNQrcAZyS5zZljKreqKonqOoZwEHg5Xy3aZg0i8h0AP/vvmwc1AqCUaKqX1TVmao6G6Pm36eqRTFCEpEqEamJfwbeglGfC55xUvL0vRSRWchnO/AmEakUEcHc96KYpAcQkan+3yMw8wO/y2+Lhs2fgQ/6nz8I/CkbB81XqUpLYdAA3Gl+zwSB36nq3/LbpGFRtCVPfcG7CvhovtsyHFT1cRG5HXgGiAHPAj/Lb6uGxR9EZDIQBT5eyA4GInIzsBKYIiI7gWuA64BbReSfMZUaL87KuWz2UYvFYiltrGnIYrFYShwrCCwWi6XEsYLAYrFYShwrCCwWi6XEsYLAYrFYShwrCEoQP+PoS0kZMGfnu12p8FNHnJ7Bfp+J+4j7yx8Tkauy3BYVkedF5Dn/7wXDOZ+I3CQin8jwXJ8Skc/5n1eOJjOsiCwVkYuT1q0XkQr/84B7N8Sx3iEiKxKWl4vIb0fatiHO9U0ReV8ujm0ZiHUfLUH8KOjzVTXrwWMi4mCi30f9YolIQFUzSgGQy2tKOIcCNaraISLnArcCE/0EZpl8/ybgKVX94RD7VQIbgEWq2i0iK4HrR5odVkQuw9ybi9Js30aG9y7Ta8gGIlIPPAw0qWox5vEqGqxGYOlDRBaIyA4ROdJfvkZEbvE/rxGRW0XkPl+b+IOfATS+7TYRWYuJ7p0gIv9LTJ2D50XkzoSIzsv8POp/FpEX/ePNSNh2j7//RmCxiDzg53NCRN4nIo/7SeaeFZGz/fX/G5P35nZ/pLvQb9P1/vaAiFwvJg/9Rv9zwN92k4j8xG/HKyLyaz9idigeAKqBiQn3IH6+U0TkGb8tL4jIe1Pc6zf792ZRimO/C3hQVbtTfG+2iLSIyNf9e7BZRE7zt031798G/993xARPfQ04x2/P9/19VUSq09y7AZpLfFlEVmNyI13t7/u/krWVIZ77WhH5vX9PHhGRaYPdL1XdjwkUPDuD52EZDapq/5XYP2Ab8BIm6dZ6zAgvvu0DwGOYdBObgVp//RpMtskGf/kXmFFqfNt2YIq/vAiTmni6v/x/gN/7ny8DuoH5/vI1wO0J2zqAuQnteQAzWgWYTL8WOx+TayjxmhYlLK9JaN+VwD1A2P93L3Clv+0mzKiz3N/2ArAqzX1ToNr//D5MXvhU5/sT8F7/swATEs71CeD9wJPAjDTnuRH4WMLyyvgzAmb77Yjfk/cDj/ifrwJ+mvC9iQn39fZBriX53t0EfCLVcoptiW0b6rkfBGb5yz8Hvj7Y/fKXvwJcl+/fzHj/Z1NMlC4XaQpTgKr+lz/S/iNwuqq2JWz+q/YnSLsR+EHCtrtUtcX//GZ/OZ4u96fAcwn7Pqyqm/3P/w9jBknctiVNm+cCN/saRBSYJiLT1OQdGoxzgJtUNQIgIr8ELgRu8Lf/UVV7/G3P+OdZl+ZY/xCTn2kacFaafe4HviQic4F12p8YD0wajG7g7KR7m8hMYLBKdx2qGt/+GPCfCZ+vEpFvYXLt3z3IMXLBUM/9EVXd4X9+DJNiAwa/X3sxdScsOcSahiwDEJO351igleHlOu/IUhMGO87NwI9V9VhMptEYZiQ/WnoSPrsMnoPrFFWdA3wZuEVSlGlU1e9iTCj7gR+IyLUJm5/DCJGmQc7RzeDX1Zuqvar6KHA88DRGs7t/kGMMRoyBfUO2SlGmvM9D3K9yzP2w5BArCCzJfAvTkawCfiKm+lqct/kTeGBGtvelOcb9wHlxGzCmPGDiCPtUETkmg+MkMwF4zf98OVCWsK0NqEvzvXuAD4pISERCmKyN6Ub8mfKfQDPwseQNIjJPVbeo6k8x5TQTU3s/g8l6+VsROTPNsTdgTF/DQkTmAG2qegvwr8AJYibvB7s3pNj+KnCif8zpmJF+un0TGeq5p2v3YPeriYFahSUHWNNQ6XK7iCSO0K7AmCRWAiepao+IfBVjiol3BA9hRsEzMJPCn011YFXdKCJXA+vEeNpsZWCWzUeA631hsBczes2EzwB/FJGDmMIibyRs+z6mWlkXxn6fyM+AozGZMsGYTH6e4TlToqoqxr3zFhH5adLmT/n3LIIZvX8y6bvP+xPgfxGRT6hqsgnnDuDHmHmH4bAS+FcRcTGDvI+pqici9wKfE5HngL+r6qeSvpd8736OeT9exOTrTzTV/Bdwk4i8G/g2Zm4ofl1DPfd0pLxf/qT9WcD/Hc5NsAwf6z5qyQgRWYOZXPzcKI9zGYO4MloMInI3cLWqPjvkzuMU30vpUlXNdKBgGSHWNGSxFCafAKbnuxF5phb4Qr4bUQpYjcBisVhKHKsRWCwWS4ljBYHFYrGUOFYQWCwWS4ljBYHFYrGUOFYQWCwWS4ljBYHFYrGUOP8fwqIOZ2X1500AAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ], "metadata": { "collapsed": false } }, { "cell_type": "markdown", "source": [ "We mentioned that there is two ways to call a linear regression, and so far we showed how to use the formula interface which is probably the recommended way to do this in most applications.\n", "\n", "If we end up in a situation where we have to run a regression many times the formula interface actually slows down the code because python has to deconstruct it to get the exact information it needs. Examples where this may occur (and I have come across) are if a regression is part of a nested loop over structural parameters (BLP estimation for IO people), or if we are coding up our own indirect inference procedure.\n", "\n", "Here's another way to get the same result that is usually faster, although we have to import some additional functionality:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 21, "source": [ "from statsmodels.tools.tools import add_constant\n", "\n", "data = add_constant(data, has_constant='add')" ], "outputs": [], "metadata": { "collapsed": false } }, { "cell_type": "code", "execution_count": 22, "source": [ "data.columns" ], "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "Index(['const', 'shortnam', 'euro1900', 'excolony', 'avexpr', 'logpgp95',\n", " 'cons1', 'cons90', 'democ00a', 'cons00a', 'extmort4', 'logem4',\n", " 'loghjypl', 'baseco'],\n", " dtype='object')" ] }, "metadata": {}, "execution_count": 22 } ], "metadata": { "collapsed": false } }, { "cell_type": "code", "execution_count": 23, "source": [ "from statsmodels.regression.linear_model import OLS" ], "outputs": [], "metadata": { "collapsed": false } }, { "cell_type": "code", "execution_count": 24, "source": [ "simple_reg2 = OLS(endog=data['logpgp95'], \n", " exog=data[['const', 'avexpr']], \n", " missing='drop')\n", "type(simple_reg2)" ], "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "statsmodels.regression.linear_model.OLS" ] }, "metadata": {}, "execution_count": 24 } ], "metadata": { "collapsed": false } }, { "cell_type": "code", "execution_count": 25, "source": [ "print(simple_reg2.fit().summary())" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: logpgp95 R-squared: 0.611\n", "Model: OLS Adj. R-squared: 0.608\n", "Method: Least Squares F-statistic: 171.4\n", "Date: Thu, 17 Aug 2017 Prob (F-statistic): 4.16e-24\n", "Time: 16:32:59 Log-Likelihood: -119.71\n", "No. Observations: 111 AIC: 243.4\n", "Df Residuals: 109 BIC: 248.8\n", "Df Model: 1 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "const 4.6261 0.301 15.391 0.000 4.030 5.222\n", "avexpr 0.5319 0.041 13.093 0.000 0.451 0.612\n", "==============================================================================\n", "Omnibus: 9.251 Durbin-Watson: 1.689\n", "Prob(Omnibus): 0.010 Jarque-Bera (JB): 9.170\n", "Skew: -0.680 Prob(JB): 0.0102\n", "Kurtosis: 3.362 Cond. No. 33.2\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n" ] } ], "metadata": { "collapsed": false } }, { "cell_type": "markdown", "source": [ "Notice that this gives us exactly the same results." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "## Multivariate Regression\n", "\n", "Univariate regression is rarely enough in social science. At a bare minimum we are worried about other factors affecting GDP that are not included in our model that are also correlated with institutional quality. As reasonably good econometricians we know that this causes bias / inconsistency in our estimates for $\\hat{\\beta}$.\n", "\n", "So let's add more regressors! AJR argue that climate influences economic outcomes, and it is probably correlated to institutional quality (Discuss / think about why). They operationalize this by including latitude as a proxy. \n", "\n", "Notice that our current data, `data`, doesn't have that variable due to the 'interesting' way Acemoglu structures his data and do files on the web. We'll need the data `maketable2.dta`:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 19, "source": [ "# import new data by overwriting the old one\n", "data = pd.read_stata('data/maketable2.dta')" ], "outputs": [], "metadata": {} }, { "cell_type": "markdown", "source": [ "Now, using the formula interface we can readily add this variable:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 20, "source": [ "reg2 = smf.ols('logpgp95 ~ avexpr + lat_abst', data=data)\n", "results2 = reg2.fit() " ], "outputs": [], "metadata": { "collapsed": false } }, { "cell_type": "code", "execution_count": 21, "source": [ "print(results2.summary())" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: logpgp95 R-squared: 0.623\n", "Model: OLS Adj. R-squared: 0.616\n", "Method: Least Squares F-statistic: 89.05\n", "Date: Mon, 13 Sep 2021 Prob (F-statistic): 1.42e-23\n", "Time: 12:11:26 Log-Likelihood: -118.09\n", "No. Observations: 111 AIC: 242.2\n", "Df Residuals: 108 BIC: 250.3\n", "Df Model: 2 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept 4.8729 0.328 14.855 0.000 4.223 5.523\n", "avexpr 0.4635 0.055 8.352 0.000 0.353 0.573\n", "lat_abst 0.8722 0.488 1.788 0.077 -0.094 1.839\n", "==============================================================================\n", "Omnibus: 6.082 Durbin-Watson: 1.756\n", "Prob(Omnibus): 0.048 Jarque-Bera (JB): 5.950\n", "Skew: -0.567 Prob(JB): 0.0510\n", "Kurtosis: 3.037 Cond. No. 57.4\n", "==============================================================================\n", "\n", "Notes:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n" ] } ], "metadata": { "collapsed": false } }, { "cell_type": "markdown", "source": [ "This is nice - but its a little hard to compare coefficients and the like across regressions that are far apart in the notebook. What we really want is a regression table that summarizes all of our results. `statsmodels` allows us to do this too:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 22, "source": [ "from statsmodels.iolib.summary2 import summary_col\n", "\n", "info_dict={'R-squared' : lambda x: \"{:.2f}\".format(x.rsquared),\n", " 'No. observations' : lambda x: \"{0:d}\".format(int(x.nobs))}\n", "\n", "results_table = summary_col(results=[results,results2],\n", " float_format='%0.2f',\n", " stars = True,\n", " model_names=['Model 1',\n", " 'Model 2'],\n", " info_dict=info_dict,\n", " regressor_order=['avexpr',\n", " 'lat_abst',\n", " 'const'])\n", "\n", "results_table.add_title('Table - OLS Regressions')\n", "\n", "print(results_table)" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ " Table - OLS Regressions\n", "================================\n", " Model 1 Model 2\n", "--------------------------------\n", "avexpr 0.53*** 0.46***\n", " (0.04) (0.06) \n", "lat_abst 0.87* \n", " (0.49) \n", "Intercept 4.63*** 4.87***\n", " (0.30) (0.33) \n", "R-squared 0.61 0.62 \n", "R-squared Adj. 0.61 0.62 \n", "R-squared 0.61 0.62 \n", "No. observations 111 111 \n", "================================\n", "Standard errors in parentheses.\n", "* p<.1, ** p<.05, ***p<.01\n" ] } ], "metadata": { "collapsed": false } }, { "cell_type": "markdown", "source": [ "If we decide we want this table as latex output we can do that too:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 23, "source": [ "with open('ols_summary.tex', 'w') as file_handle:\n", " file_handle.write(results_table.as_latex())" ], "outputs": [], "metadata": { "collapsed": false } }, { "cell_type": "markdown", "source": [ "Notice we did not provide a file path, so the results are saved to our current working directory." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "### Challenge\n", "\n", "It's your turn!\n", "\n", "1. Add the following variables to the regression using the formula interface:\n", " * 'asia', 'africa', 'other'\n", "2. Add the same variables as (1); but use the alternative interface that requires you to specify 'endog' and 'exog' variables\n", "3. Export the results to a latex file 'out/tables/ols-regressions'. Give the columns meaningful names, and make sure you don't have stars in your output to satisfy the new AER policy." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "## Robustifying Standard Errors\n", "\n", "So far our regressions have imposed heteroskedastic standard errors.\n", "\n", "A simple test for Heteroskedasticity is the Breusch Pagan test:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 45, "source": [ "reg3 = smf.ols('logpgp95 ~ avexpr + lat_abst + asia + africa + other', data=data)\n", "results3 = reg3.fit() \n", "\n", "from statsmodels.compat import lzip\n", "import statsmodels.stats.api as sms\n" ], "outputs": [], "metadata": { "collapsed": true } }, { "cell_type": "code", "execution_count": 47, "source": [ "name = ['Lagrange multiplier statistic', 'p-value', \n", " 'f-value', 'f p-value']\n", "test = sms.het_breushpagan(results3.resid, results3.model.exog)\n", "lzip(name, test)" ], "outputs": [ { "output_type": "stream", "name": "stderr", "text": [ "/home/lachlan/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:3: DeprecationWarning: `het_breushpagan` is deprecated, use `het_breuschpagan` instead!\n", "Use het_breuschpagan, het_breushpagan will be removed in 0.9 \n", "(Note: misspelling missing 'c')\n", " app.launch_new_instance()\n" ] }, { "output_type": "execute_result", "data": { "text/plain": [ "[('Lagrange multiplier statistic', 24.434599912931393),\n", " ('p-value', 0.00017909452201843713),\n", " ('f-value', 5.9276177047116967),\n", " ('f p-value', 7.2348350639932351e-05)]" ] }, "metadata": {}, "execution_count": 47 } ], "metadata": { "collapsed": false } }, { "cell_type": "markdown", "source": [ "There is some evidence for heteroskedasticity. Let's robustify standard errors. Statsmodels has 'HC0' 'HC1', 'HC2' and 'HC3' standard errors. The STATA default is HC1, so for better or worse let's get a summary with those:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 50, "source": [ "robust_ols = reg3.fit(cov_type='HC1')\n", "print(robust_ols.summary())" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: logpgp95 R-squared: 0.715\n", "Model: OLS Adj. R-squared: 0.702\n", "Method: Least Squares F-statistic: 97.09\n", "Date: Thu, 17 Aug 2017 Prob (F-statistic): 9.39e-38\n", "Time: 17:09:20 Log-Likelihood: -102.45\n", "No. Observations: 111 AIC: 216.9\n", "Df Residuals: 105 BIC: 233.2\n", "Df Model: 5 \n", "Covariance Type: HC1 \n", "==============================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept 5.8511 0.293 19.936 0.000 5.276 6.426\n", "avexpr 0.3896 0.051 7.638 0.000 0.290 0.490\n", "lat_abst 0.3326 0.442 0.752 0.452 -0.535 1.200\n", "asia -0.1531 0.181 -0.848 0.396 -0.507 0.201\n", "africa -0.9164 0.154 -5.946 0.000 -1.218 -0.614\n", "other 0.3035 0.174 1.741 0.082 -0.038 0.645\n", "==============================================================================\n", "Omnibus: 4.342 Durbin-Watson: 1.865\n", "Prob(Omnibus): 0.114 Jarque-Bera (JB): 3.936\n", "Skew: -0.457 Prob(JB): 0.140\n", "Kurtosis: 3.126 Cond. No. 58.2\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors are heteroscedasticity robust (HC1)\n" ] } ], "metadata": { "collapsed": false } }, { "cell_type": "markdown", "source": [ "There are also options for clustered standard errors. A nice summary of what is available albiet on a different data set is available [here](http://www.vincentgregoire.com/standard-errors-in-python/)" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "## Instrumental Variables Estimation" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "AJR discuss that the OLS regression results may suffer from endogeneity due to:\n", "\n", "\n", "* richer countries may be able to afford or prefer better institutions\n", "* variables that affect income may also be correlated with institutional differences\n", "* the construction of the index may be biased; analysts may be biased towards seeing countries with higher income having better institutions\n", "\n", "They propose a 2SLS solution using settler mortality to instrument for institutional differences: They hypothesize that higher mortality rates of colonizers led to the establishment of institutions that were more extractive in nature (less protection against expropriation), and these institutions still persist today.\n", "\n", "Data for the IV regressions are contained in the file `maketable4.dta` from AJR's replication files. So we will start by loading that:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 24, "source": [ "data = pd.read_stata('data/maketable4.dta')\n", "data.head()" ], "outputs": [ { "output_type": "execute_result", "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
shortnamafricalat_abstrich4avexprlogpgp95logem4asialoghjyplbaseco
0AFG0.00.3666670.0NaNNaN4.5400981.0NaNNaN
1AGO1.00.1366670.05.3636367.7706455.6347890.0-3.4112481.0
2ARE0.00.2666670.07.1818189.804219NaN1.0NaNNaN
3ARG0.00.3777780.06.3863649.1334594.2326560.0-0.8722741.0
4ARM0.00.4444440.0NaN7.682482NaN1.0NaNNaN
\n", "
" ], "text/plain": [ " shortnam africa lat_abst rich4 avexpr logpgp95 logem4 asia \\\n", "0 AFG 0.0 0.366667 0.0 NaN NaN 4.540098 1.0 \n", "1 AGO 1.0 0.136667 0.0 5.363636 7.770645 5.634789 0.0 \n", "2 ARE 0.0 0.266667 0.0 7.181818 9.804219 NaN 1.0 \n", "3 ARG 0.0 0.377778 0.0 6.386364 9.133459 4.232656 0.0 \n", "4 ARM 0.0 0.444444 0.0 NaN 7.682482 NaN 1.0 \n", "\n", " loghjypl baseco \n", "0 NaN NaN \n", "1 -3.411248 1.0 \n", "2 NaN NaN \n", "3 -0.872274 1.0 \n", "4 NaN NaN " ] }, "metadata": {}, "execution_count": 24 } ], "metadata": { "collapsed": false } }, { "cell_type": "markdown", "source": [ "A scatter plot provides some evidence towards the necessary condition of a correlation between the instrument and the endogenous variable." ], "metadata": {} }, { "cell_type": "code", "execution_count": 30, "source": [ "# add numpy so we can get a linear prediction to add to plot\n", "import numpy as np \n", "\n", "# Dropping NA's is required to use numpy's polyfit\n", "df_plot = data.dropna(subset=['logem4', 'avexpr'])\n", "labels = df_plot['shortnam']\n", "\n", "fig, ax = plt.subplots()\n", "\n", "ax.scatter(df_plot['logem4'], df_plot['avexpr'], alpha=0.5)\n", "# add a linear prediction line on the fly\n", "ax.plot(np.unique(df_plot['logem4']),\n", " np.poly1d(np.polyfit(df_plot['logem4'], df_plot['avexpr'], 1))(np.unique(df_plot['logem4'])),\n", " color='red')\n", "# add labels to points\n", "for i, label in enumerate(labels):\n", " plt.annotate(label, (df_plot['logem4'].iloc[i], df_plot['avexpr'].iloc[i]))\n", "ax.set_xlabel('Log of Settler Mortality')\n", "ax.set_ylabel('Expropriation Risk (Institutions)')\n", "ax.set_title('Graphical Evidence of a First Stage Relationship')\n" ], "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "Text(0.5, 1.0, 'Graphical Evidence of a First Stage Relationship')" ] }, "metadata": {}, "execution_count": 30 }, { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEaCAYAAAAR0SDgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAACTQklEQVR4nOydeXhMVx+A3zPZExJL7LFGCJEFQSk+qtZaq7V1odoqovYW1ZK2lFaUttRS2mgRqkUtRW1RqlVbLBG7IAhii+zJzPn+uDPTmWSSDLLhvs9zn2Tucs7v3rlzzz2/VUgpUVFRUVFR0RS2ACoqKioqRQN1QFBRUVFRAdQBQUVFRUVFjzogqKioqKgA6oCgoqKioqJHHRBUVFRUVAB1QCgSCCGChRBLc9geKYRolZ99PEK784UQH+WwXQohauZ1v/mBEMJJCLFeCHFPCLEqn/poIYQ4lR9tP6kIIarp7yPbhzz+AyHEoryWK1Mf0UKI57PZ9th85+qAYAEhRB8hxD4hRKIQ4ob+/6FCCFEY8kgpfaSU4fnVvhCilRBCJ4RIyLQ0tUK2wVLKT/NLtgLmJaAcUFpK+fKjNGTyEDO9nkeklLullLUfss1cB3UhRHMhxF79oHZbCPGXEKKRftsAIcSeh+n7Ucl0j90XQpwSQryRT/3EmK6TUn4mpXwrr/uylkf5zguahxpxn2SEEGOA94EgYAuQAAQAY4HFQKqFY2yklNoCFDM/uCql9ChsIQqZqsBpKWVGHrZZwtr2hBC2j9K3EMIV2AAMAX4G7IEWWLhnC4mrUkoP/YtVR2CdEGKvlPKxeHt+KpBSqot+AdyARKBnLvuFAvOA3/X7Pw+8ABwG4oHLQLDJ/tUACQwCrgLXgLEm24NRfsA/AveBSCDQZHs08Lz+fxvgA+Ccft+DQGX9tq/0fcfr17fI1MfSbM6nFRCTzbbewIFM60YB60yuxRSTbe/pz+8qMFB/3jX12xyAEOAScB2YDziZygCMAW7o23jDpF0nYCZwEbgH7DE59hlgL3AXOAK0yuG7qwOE6/eNBLrq138MpAHpKC8Bb1o4tjHwt/7Ya8AcwD6bfgzfuW1O11r/3Y4DjqI8uG31n6/ov99TQBugQyb5jljoMxC4m8N5pwBa/fF39euzvW/121/XX/NbwEeY34saYDzKvXgL5R4uZe09pv+eX86trczXEngDiNJfn/PAO/r1LkAyoNOfYwJQkUz3PtBV/93f1d8LdTJ9H2P138c9YCXgqN/mjjLg3gVuA7sBjRXHWfrOJwAngDvAD4Z9C3spdAGK0qL/0WWQ6UdsYb9Q/Zf+rP5GdtR/6b76z34oD7zumW7oMP1N6wvcNPlhBaP8WDuhPPCnAf9kuoEM+74HHANqAwLwR1FxALwKlEZ5qIwBYk1uSrMfRabzMbthM21z1v/wvEzW7Qf6mFyLKSbX7zpQT3+eyzEfEGYB64BSQHFgPTDNRIYM4BPATn8tkoCS+u1zUX68lfTXqBnKAFMJ5QHSSX/t2+o/l7FwLnbAWZQB1R54Tn9utXO7RvrtDVEGH1v9dxoFjMxmX8N3bs2AEAFURhn0aqM8mCuatONppXyu+nNfgvIGXjLT9gHAHgvyZHff1kV5qDbXX68QlAHJcC+OAP4BPPTfxQIgLLd7TN9XV5QHd/3c2sp8LVEGMU+U+/9/+vukQXb3sul1A2qhvMS11d8P7+vvCXuT7+NflIGklP47HqzfNg3lJcZOv7QAhBXHWfrOj+u/81LAX5i8VBXqM7CwBShKC8oDNTbTOsObZzLQUr8uFPgxl7ZmA7My3dDeJtu/ABab3LDbTLbVBZIz3UCGH+EpoJuV53MH8DfpI6cBQac/T9PFRb99KTBJ/78XykPU2eRaGAaE74HpJu3W0p93TZQfbyL6h5t+e1PggokMyZg8QFHeIJ9BeYAkG84lk+zjgJ8yrdsC9LewbwuUQVJjsi4M/VtxTtcom+s2EliTzTbDd256Pcdi+eEw0ORzTf15Pw/YZWozV/lQZgKhKLOtDJQBuJx+2wAyDQi53LeTMHnAo7wcpPHfvRgFtDHZXgFlwMjyQpXpHktFmamMNNmebVtkM7ia7LsWGGHST04DwkfAzybbNCizsVYm38ermX6n8/X/fwL8hv4FJ1MfOR1n6TsfbPK5E3DO2vsuPxfVqGzOLcDd1JtBStlMSllCv830el02PVAI0UQIsVMIcVMIcQ8YjDLFJJtjLqK8TRiINfk/CXDMxquiMsq0OgtCiLFCiCi9QfEuigosswzZcVVKWSLTkqjfthzoq/+/H7BWSplkoY2KZD1HA2VQHigHhRB39fJt1q83cEua69CTgGL6c3DE8nlXBV42tKlvtznKA8WifFJKXSYZK1nYNwtCiFpCiA1CiFghRDzwGblfX3eT6xmSzT7GayalPIsy0AQDN4QQK4QQFbM5LgtSyigp5QCp2IPqoZzz7BzOKaf71uz71H/nt0wOrwqsMbnuUSgP+nLZdHdV/1tyBb5GmaE9cFtCiI5CiH/0RvO7KA9Ua+/zipjcl/p74TLm90Dm32Ix/f8zUGYTfwghzgshxmdqO7vjLJHTs6DQUAcEc/5GeXvpZsW+MtPn5ShvY5WllG4oU8vMXkmVTf6vgqJnf1Auo0yXzRBCtECZ/vZCURWUQFFr5YVn1FagjBAiAGVgWJ7NftfIeo4G4lDe8n1MHpBuUsqcfjSmx6Zg4bxRrsdPmQYyFynldAv7XgUqCyFM7/sqKG+I1jAPOImiPnNFUT3lxfU1u5eklMullM1RHpIS+NzSfrk2KuVJlNlCvRyOz+m+vYaiwgEUt1wUlaSBy0DHTNfeUUqZ4/WUUqaizOx8hRDdH6QtIYQD8CuK+qqc/j7/3UTm3K7RVZTramhPoNyzud4DUsr7UsoxUsoaKCqv0UKINrkdlw158SzIc9QBwQQp5V0U4+K3QoiXhBDFhRAa/YPQJZfDiwO3pZQpQojGKG/SmflICOEshPBBMYytfAgxFwGfCiG8hIKfEKK0vv8MFNuErRBiEsqb2CMjpUwHVqG8IZVCGSAs8TMwQAhRVwjhDEw2aUMHfAfMEkKUBRBCVBJCtLeifx2KOupLIURFIYSNEKKp/uGwFOgihGivX++odz205DG1D+XN7X0hhJ1QYju6ACusuQ4o1zgeSBBCeKN48+QpQojaQojn9OeWwn9GUlD0+9UyDWimx3oLIcYYzl0IURllAP/H5HgPIYR9pnPK7r79BeXaNtMfE4z5ADgfmCqEqKrvr4wQwpqXKaSUaShOApMesC17FBvDTSBDCNERaGey/TpQWgjhlk3XPwMvCCHaCCHsUGxtqSiq4RwRQnQWQtTUDyL3UGYwulwOy44gIYSHEKIUMJGHexbkOeqAkAkp5RfAaJS37ev6ZQHKG01ON81Q4BMhxH2Um/xnC/vsQplybgdCpJR/PISIX+rb/gPl4bQYxRi5BUUFcxplCppCJrVWLlQUWeMQeppsX46i114ls3GNlFJuQlFP7EA5zx2ZdhmnX/+PXuWyDcWIag1jUYzp+1E8PD5HsQVcRpnRfYDykLiMYnjPcm/rH0JdUAyuccC3wOv6N2lrZeiHYkP5jvz5ETsA0/XyxQJlUTxSQBmUAW4JIQ5ZOPY+0ATYJ4RIRBkIjqM89ED5PiKBWCFEnH5dtvetlDISeBdlwLyGYmC+wX9urF+hzC7+0B//j75/a/keqCKE6GJtW1LK+8BwvZx3UL6PdSbbT6LYhc7r1U8VMx1/CsVW+A3KNe4CdNHfG7nhhXLPJqBoE76VUu58gPM1ZTnKb/g8iip0ykO2k6cYLOQq+YgQohpwAcVImJc+7ioqBYYQohiKUdhLSnmhkMV5bBFCRANvSSm3FbYsmVFnCCoqKtkihOiiV3O6oOjtj6F4yag8gagDgoqKSk50QzF4XkVRmfSRqlrhiUVVGamoqKioAOoMQUVFRUVFjzogqKioqKgAj2G2U3d3d1mtWrXCFkNFRUXlseLgwYNxUsoyOe3z2A0I1apV48CBA4UthoqKispjhRDiYm77qCojlVyJjo6mXr16ZuuCg4MJCQnhn3/+oUmTJgQEBFCnTh2Cg4PN9hs5ciSVKlVCp3vYgE4VFZWC4rGbIagULfr378/PP/+Mv78/Wq2WU6f+q3Wi0+lYs2YNlStXZteuXbRu3boQJVVRUckNdYag8kjcuHGDChWUxKI2NjbUrVvXuC08PBwfHx+GDBlCWFhYYYmooqJiJeqAoPJIjBo1itq1a9OjRw8WLFhASkqKcVtYWBh9+/alR48ebNy4kfT09EKUVEVFJTfybUAQQnwvlAL1x03WlRJCbBVCnNH/LZlf/as8OulaHXeT0sjQWQ5eFEIwadIkDhw4QLt27Vi+fDkdOnQAIC0tjd9//53u3bvj6upKkyZN2LJlS0GKr6Ki8oDkpw0hFKXm7I8m68YD26WU0/XFJcajZMBUKULodJKtUdeJuhZPcpoWjTaF2Ju30OkkGo2S/fj27dtUr14dAE9PT4YMGcLbb79NmTJluHXrFnv37uXu3bv4+voCkJSUhJOTE507dy6081JRUcmZfJshSCn/RElTbEo3lHqv6P92z6/+VR6erVHXORZzDynB0c4Ge0cXnNxKMyP0V0AZDDZv3kzz5s3ZuHGjoQwgZ86cwcbGhhIlShAWFsaiRYuIjo4mOjqaCxcusHXrVpKSLBVaU1FRKQoUtA2hnJTymv7/WLIvtZevXL9+nX79+lGjRg0aNmxI06ZNWbNmDeHh4bi5uREQEICfnx/PP/88N27cACA0NJQyZcoQEBCAt7c3s2bNKgzR8510rY6oa/HYaMwLgfV7/wsWfT0D/4AAnnvuOSZPnoynpyc//fQTtWvXJiAggNdee41ly5aRmprK5s2beeGFF4zHu7i40Lx5c9avX1/Qp6SiomIlhWZU1mdMtCqznhAiWAghhRDy6tVHqzQnpaR79+60bNmS8+fPc/DgQVasWEHPnj156623jG+769ato1y5clSuXJmAgAAmTJiAh4cHERER/PXXX0ydOpUjR45gZ2fH/PnzH0mmR0UIwZgxY4yfQ0JCjPEA7du3JyAgwLhUrFiRJk2UuiMDBgzgl19+MWsrMTWD5DRtlj7KV63JwM9C2bX3XzZv3sz69evx9PTkzJkz1KxZk59//pkDBw4QFRVFqVKluHDhAq6uSsG28PBwhBC88cYb9O7dG4DOnTsTHh6eD1dDRUXlYSnoAeG6EKICgP7vDWsOklIGSymFlFJUrPhotah37NiBvb09gwcPNq6rWrUqzs7OLFq0iJYtWxIREUHVqlVJSkqiZs2aRERE8PHHHxMdHc1ff/1F6dKlqVmzJsuXL+eZZ54pdJdKBwcHVq9eTVxcXJZtW7ZsISIiwjiQubq6MmVK9sWZXBxscbK3sbjNyd4GZ3sbevToQatWrTh37hwHDx5k2rRpXL9+HVA8ixo1asTq1avNjvXw8GDq1KmPcJYqKir5TUEPCOuA/vr/+wO/FXD/REZG0qBBg2y37969m4CAAKpUqcLBgwfx8FBK89rb2+Pu7s6VK1e4dOkSKSkp/PXXX8ycOZMrV64QExNTUKeQBVtbWwYNGpSrGmvEiBF06tSJtm3bZruPnY2GOhVc0WbyLNLqJHUquLLnz13Y2dmZDaj+/v60aNGCc+fOkZCQwJQpU7IMkv7+/ri5ubF1a3blmFVUVAqb/HQ7DUOpO1pbCBEjhHgTpVZsWyHEGZT6vNPzq39LpGt1JKdloDN52AUFBeHv709iYqJRZVS9enUuX75Mx44dOXlSKbebmJhIdHQ0wcHB1KxZkz59+nD9+nUaN25Mr169WLmycGtkBwUFsWzZMu7du2dx++rVqzlw4ADTpk3Lta22dcrh6+GGEJCSrkUI8PVwo22dchw/fpyGDRtaPG7FihX06dOHFi1acOrUKeOswcDEiRNznJ2oqKgULvnmdiql7JvNpjb51Wd2mLpRnk4rwa6dK+gUGUvbOuWY/fU3XLoSi1c1DxYtWkRISAhr1qwBoFmzZixatAh/f3+ioqLw8/Pj4MGDHDhwgJYtW/L2228D0KdPHwYOHGimxy8I0rU6ElOVEs2urq68/vrrfP311zg5OZntd+XKFUaMGMGWLVtwcHDItV2NRtDepzzPeZclMTUDFwflNolPSc8yczAlLCyMNWvWoNFo6NmzJ6tWrWLYsGHG7S1btgRgz549D3yuKioq+c9TkcvI4EZpoxH4BD7Lth+/YsG8eRx+oQ+OdjZci7mMlPDvhduYFpA7duwY7u7uHDlyhBkzZjB58mQiIiIIDAykWLFifP/99/z2m6L1unr1KmfOnMHLyyvfzydznEC6VrIlMpbhw0cQGNiQN954w7ivlJL+/fszfvx4s7QS1mBno8HV0c6sr8tpJdi9+2+zmARQrtWZM2eM6qi0tDSqV69uNiDAf7MEW9un4tZTUXmseOJTV2R2oxRCMDB4LscO/sOnr7Xhy6CX+HX2B9jY2XHhZgLhf/5JQEAA/v7+bN26lTp16gBQpkwZGjRowOeff87p06cpXrw4Li4uHDt2jOjoaCZMmFBgxuXMcQIAx2LucfB6Or169WLx4sXGfUNCQnB0dCQoKChP+qoZ8Az3EpIZ9UmIcZ+jR48yfPhwgoODjXEHV69e5erVq1y8aJ5xt127dty5c4ejR48+lDwqKir5xxM/IFhyoyxWqgwNBgTT4ZNVDJm1kqAZP2FjY0ut+s8wdc0B9h88xJEjR/jqq68oVqwYoLhobt26lT///JOwsDBefvllYmNjKV68OAA9e/YskAEhuzgBG40g6lo8w0eOMvM2+vDDD4mKijJzPTXNOvrOO+/g4eGBh4cHTZs2zbUvZUCdw66dO/D09MTHx4cJEyYQHh5Ojx49zI7v0aMHK1asyHIOEydO5PLly490HVRUVPIeIaVVoQBFhsDAQPkgBXLStTrm7jxrpgpKSddyIPo2tjaCRtVKmz3wUtK1DGnlSQln+7wUO8+4m5TGvPBzxpmBKXkte0H2paKikr8IIQ5KKQNz2ueJnyFYcqO0t9Wg0QjcizlkedN2srcxGlGLIrnFCeSl7AXZl4qKSuHzxA8IkNWN0kYjaFi1JNVKu5jtZ/C1t7MpupcltziBvJS9IPtSUVEpfJ6KVzxLbpQ2Qph5zzjZ2xh97Ys6BhkLQvaC7EtFRaVweeJtCLlh8OV3cbB97N54C1L2x/k6qaio5JENQQjRRgjxmRDiRyHEAiHEUCFEpbwTs3Cxs9FQwtn+sXzImcq+du1ahBDGyOrw8PAstQdMk9lt2LCB+vXr4+/vT926dVmwYIHVfWVm6tSp+Pj44OfnR0BAAPv27SMjI4MPPvgALy8vo3eTaS4jGxsbAgICqFevHi+//LKaFltFpQiQ7VNQCNFXCHESGA0kA7uB40A9YLsQYokQonzBiKmSG2FhYTRv3twq19f09HQGDRrE+vXrOXLkCIcPH6ZVq1YP1e/ff//Nhg0bOHToEEePHmXbtm1UrlyZDz/8kKtXr3Ls2DEiIiLYvXu3WQlNJycnIiIiOH78OPb29oWeMVZFRSVnG4I/0FxKmTWFJiCEaAc8C/yaH4KpWE9CQgJ79uxh586ddOnShY8//jjH/e/fv09GRgalS5cGlGyptWvXfqi+r127hru7uzElhru7O0lJSXz33XdER0fj6OgIQPHixY0puTPTokULNVBNRaUIkO0MQUo5PrvBQL/9DymlOhgUIoaax7+uWUOHDh2oVasWpUuX5uDBgzkeV6pUKbp27UrVqlXp27cvy5YtQ6fTPVTfrds8z+XLl6lVqxZDhw5l165dnD17lipVqhiD9nIiIyODTZs2GUttqqioFB7W2BB6CyFc9f9/KoTYLISwnO5SpUDQ6ZTcRXN3nmVe+Dm+mPs9NZ9ph04n6dOnD2FhYQghLB5rWL9o0SK2b99O48aNCQkJYeDAgQ/V95L9sXz240bmz19AmTJl6N27d5bCNz/88AMBAQFUrlzZGKGcnJxMQEAAgYGBVKlShTfffPPhL4iKikqekKuXkRDimJTSVwjRGPgG+AoYJqVsVhACZiavvYweR7ZExhqT9SXG3+WTV/6Hi1sp7G012GuUh/7GjRsZPHgwf/31l/G4rl27MmbMGP73v/+ZtRcXF0f16tW5f//+A/VtQKuT+Hq40d6nPL/88gsLFizg0KFDREdHm80S6tWrx4YNG6hWrRrFihUjISEhD66GioqKNeRVpLLBEtgWWCSlXA44PqpwKg9H5vxCR3dvoWGbbkxaupMJS7ZzPvoi1atX5/bt21y9epWoqCgALl68yJEjRwgICCAhIcHsLd5QIe5B+wa4cfk8t69dJOpaPOlaHREREdSuXZs333yTYcOGkZKSAoBWqyUtLS0Pr4SKikpeY01gmhRC9Ab6AF3169QENoWEIVmfIb/QofANPNdLqcuQnKYlMTWDnj17smLFCpYuXcobb7xBSkoKdnZ2LFq0CDc3N+7fv88XX3zBO++8g5OTEy4uLoSGhj5w3wCpyUms+XYKiffv8Z2rM7VrebFw4ULc3Nz46KOPqFevHsWLF8fJyYn+/fvzqCVQVVRU8g9rVEZNgXHATinlV0IIL2C4lPLdghAwM0+7yshSsj4DQkBQ65r5FlNRmH2rqKg8GnmiMpJS/i2l7C6l/Er/+UxhDQYqhZtfSM1tpKLyZJOrykgIURZ4F/A03V9K2Ssf5VLJgcLML6TmNlJReXKxRmW0G4gC/gGMlWaklEvyVzTLPO0qI1MKM7+QmttIReXxwhqVkTVG5ZJSykF5JJNKHmLIL/S09a2iopI/WPNqd1wIkaeuIUKIEUKI40KISCHEyLxsW0VFRUXl4bBqhgAcE0L8BaQYVj6sDUEIUQ94G2gMpAGbhRAbpJRnH6Y9FRUVFZW8wZoBYbl+ySvqAPuklEkAQohdwIvAF3nYh4qKiorKA5LrgJAPxuPjwFQhRGmUtNqdANVKrKKiolLIWJPczl0IsUIIcVO/LBdClHnYDqWUUcDnwB/AZiACE++lbGQIFkJIIYS8evXqw3atoqKiopID1hiVFwCngQCgPnBGv+6hkVIullI2lFK2BO7o289p/2AppZBSiodOfZCUBF9+CffuPdzxKioqKk841gwInlLKSVLKK1LKGCnlZKDGo3SqD3ZDCFEFxX6QlzYKyyxZAmPGQOXKyl99GmYVFRUVFQVrBgSN4QEOxof5o0Yi/SqEOAGsB4KklHcfsb3c6dsXpk+HYsWUmUKNGvDaa3DkSL53raKiovI4YM2DPQQ4LIRYKIRYCBziET2CpJQtpJR1pZT+Usrtj9KW1ZQoAePGwYUL8P33UKsWLF0KAQHQrh1s3YrFrG0qKgVAsWLFjP///vvv1KpVi4sXLxITE0O3bt3w8vLC09OTESNGGNOIh4eH4+bmRkBAAN7e3owdO7awxFd5QrAmud2PQHvgqH5pL6Vcmt+C5RsODvDGG3DsGGzcCK1aKYNBu3ZQv74ySJgUg1dRKUi2b9/O8OHD2bRpE1WqVOHFF1+ke/funDlzhtOnT5OQkMDEiRON+7do0YKIiAgOHz7Mhg0bzAoiqag8KFapfqSUx6WUc/RLZH4LVSBoNNCpE+zcCfv3Q+/eyiDx2muKOmnmTIiPL2wpVZ4i/vzzT95++202bNiAp6cnO3bswNHRkTfeeAMAGxsbZs2axffff09SUpLZsU5OTgQEBHDlypXCEF3lCSHbAUEI8ZP+734hxL+Zl4ITsQAIDIQVK+DsWRg+HG7fhrFjFQP0uHGg/shU8ol0rY67SWmkpqbSvXt31q5di7e3NwCRkZE0bGhevtzV1ZUqVapw9qx5YP+dO3c4c+YMLVu2LDDZVZ48cpohzNb/HQu8Z2F58qheHb76SvFAmjIFnJzgiy+U9QMGKDMIFZU8QKeTbImMZe7Os8wLP4ewscXLtyGLFi16oHZ2796Nv78/lSpVon379pQvXz6fJFZ5Gsh2QJBSHtT/W1lKuct0ASoXjHiFRKlSMHEiREfDd98pKqQlS8DPDzp2hB07VAO0yiOxNeo6x2LuISU42tkghIZuo2ew9c+9fPbZZwDUrVuXgwcPmh0XHx/PpUuXqFmzJqDYEI4cOUJkZCSLFy8mIiKioE9F5QnCGhvCaCvXPXk4OsJbb8GJE7BuHbRoAZs3Q5s2ipopLAwyMgpbSpXHjHStjqhr8dhohNl6J2dnXp30LUuXLWPx4sW0adOGpKQkfvzxRwC0Wi1jxoxhwIABODs7mx1bvXp1xo8fz+eff15g56Hy5JGTDSFQCBEEuAshhposE4CnKxG+RgNdusCff8I//8BLL0FEBPTrB56eMHs23L9f2FKqPCYkpmaQnGY5W4vGsTg/r1nHlClTWL9+PWvWrGHVqlV4eXlRq1YtHB0djTOIzAwePJg///yT6OjofJRe5Ukm24ppQohuQHegK7DOZFM88JOUslAS0hWZimnnzsGsWUpMQ3KyEucweLBilK5QobClUynCpGt1zN151qLWUQgIal1TrUKnkudYUzHNmhKa7aSUf+SpZI9AkRkQDMTFwbffwpw5cPMm2NvDq68q6THq1i1s6VSKKFsiYzkWc89MbaTVSXw93GjvoxqGVfKevBoQhlpaL6X89hFke2iK3IBgIDkZfvxRiV84c0ZZ98IL8N570LKl8uqnoqJHp5NsjbpO1LV4ktO0ONnbUKeCK23rlEOjUe8VlbzHmgHBmnlpI5OlBTAJaPfo4j1hODnBO+/AyZOwZg00a/ZfJHSTJvDzz6oBughgSBERHR2NEIJvvvnGuG3YsGGEhoYCMGDAAKpXr46/vz+1atXi9ddfJyYmJs/k0GgE7X3KE9S6JkNaeRLUuibtfcqrg8FjxPXr1+nXrx81atSgYcOGNG3alDVr1hAeHo4QwsyFOCIiAiEEISEhwH/3lyHtyMcff1xYp2GGNakr3jBZ+qKkwNblv2iPKRoNdO8Of/2lLD16wIEDSiR0rVrwzTeQmFjYUqoAZcuW5auvvjLmBsrMjBkzOHLkCKdOnaJ+/fo899xz2e77sNjZaCjhbK/aDB4zpJR0796dli1bcv78eQ4ePMiKFSuMLw316tXj559/Nu4fFhaGv7+/WRszZswgIiKCiIgIlixZwoULFwr0HCzxwHehlPIaUCsfZHnyaNYMVq+GU6cUg/O1a4rRuXJl+PBDuH69sCV8qilTpgxt2rRhyZKciwIKIRg1ahTly5dn06ZNBSSdSlFmx44d2NvbM3jwYOO6qlWr8u677xr/T0lJ4fr160gp2bx5Mx07drTYVkqKUqrexcUl/wXPBWsqppm6nA4TQiwBbhSAbE8OXl4wbx5cvAiTJimziKlToWpVGDRIUTOp5CuGFBGZGTduHCEhIWi1ORbtA6BBgwacVL+rp550rY4Dh48QEFA/x/1eeuklVq1axd69e2nQoAEODg5m29977z0CAgLw8PCgT58+lC1bNpuWCo4HtSH4AyeAl/JTqCeWsmXh44/h0iWYOxc8PJRI6Dp1oFs32LOnSEZAm6ZmBggNDWXYsGEABAcHG/WiKSkptG3bluDgYCB7HWtBkjlFRLpW+azTKde5Ro0aNGnShOXLc6/RlJsDhsqTjem9tOdsHBExd433UlBQEP7+/jRq1Mi4f69evVi1ahVhYWH07ds3S3sGlVFsbCzbt29n7969BXk6FrFmQPjcxIbwtpTyc6Dwh7LHGWdnGDpUUSX98otidDZEQjdtCr/+Cla8sRYl0tLS6NmzJw0bNiQ4ODhXHWtBkTlFBMCxmHvsPhtn3OeDDz7g888/z/WBf/jwYerUqZOv8qoUXUzvJY8atYk5c4JjMffYGnWduXPnsn37dm7evGncv3z58tjZ2bF161batGmTbbvFihWjVatW7NmzpyBOI0esGRAsvTrlf8nLpwEbG+jZE/7+G3bvhq5dYd8+JRK6dm0lviFTmuOiSEZGBr1798bLy4vp06cDuetYC4LsUkTYaATnbtzH8Pj39vambt26rF+/3mI7Ukq+/vprrl27RocOHfJZapWiSOZ7ySvgGTLSUvlnYxhR1+JJ1+qypCQH+OSTT/j888+xsbHJtu2MjAz27duHp6dnvslvLTmlrnAXQtQFHIUQdYQQdfVLU6DwrR9PEkJA8+bw228QFQVvvw0xMRAUBFWqwOTJStBbAWPQuycnJxMQEGBcJk2aZLbfF198gb29PbNnzzaui4yMpEGDBgUssTk5pYhISdcZ1UYAEydOzDJ7ee+994xup/v372fnzp3Y2z9dWVtUFDLfS0IIBgbP5dzR/Xzy6nM806QJ/fv3z5JLqlmzZnTv3t1imwYbgp+fH76+vrz44ov5eQpWkVPqihHASKAicNVk0z3gGynl4nyXzgJFNjAtr7l+XYl+njsX7txREu0NGACjRytG6nwkc9BUcM9A1h04ZwyaCg0N5cCBA8yZM4fg4GAOHTrE4cOH2b59O7VqKQ5oX3/9NRcuXGDWrFkABAUFsWfPHuzt7dm/f3++ym9ATRGhklc8CffSIwWmSSm/klJWBz6WUlY3WQIKazB4qihXDj79VKnN8PXXUL48zJ+vqJJefFFRM+UT2endt0ZldZP9+OOPiYmJYfbs2XTs2JHLly9TpkwZli9fzu7du6lVqxbJyclGHWtkZCRhYWFcv36dzp074+/vT926denUqVOen4edjYY6FVzR6sx/xVqdpE4F1yL/A1YpOjwt91JOKiODj9RsIYRz5qWA5FNxcYF331XSYaxcCQ0b/hcJ/eyzsHYt6PIuTjAnvbtBV2qKnZ0dsbGxdOrUibFjx9KyZUvKly+Pu7s7NjY2VK9enalTpwKwZs0apJT07duXSZMm0bZtW44cOcKJEyeMtoe8pm2dcvh6uCEEpKRrEQJ8PdxoW6dcvvSn8uTyNNxLOQ1rhlfQBOC+/m+CyWeVgsTWFnr1gn//hfBwJU/S3r1KJLS3NyxYoORTekRy0rsnp2lJTM2afqNOnTps3LiRIUOG4OrqSkJCAjqdjrVr11K8eHE+//xz6tWrx4gRI5g2bRoA165dw8PDw9iGn5/fI8tuibxOEWFjY0NAQAD+/v40aNDA6CoYHR1trGtct25dXn/9ddLT043HZWRkUKZMGcaPH58n56VS8DwV6UaklAW+AKOASOA4EAY4Wntsw4YNpYqeyEgpBw6U0t5eSpCyTBkpP/5Yyps3H7rJtAytnLX1lPzyj6zLrK2nZFqG1mx/FxcXeeTIEdmzZ0+ZnJws/f395c6dO+ULL7xg3GfdunWyePHicvLkycZ1mzdvlm5ubrJVq1ZyypQp8sqVKw8tc0Hi4uJi/H/z5s2yZcuWUkopL1y4IH18fKSUUmZkZMjWrVvLpUuXGvf9/fffZbNmzWSNGjWkTqcrWKFVVKSUwAGZy/PVmkjl2dassxYhRCVgOBAopawH2AB9Hra9p5q6dWHxYqXU54QJkJ6ueCRVqQLDhik1Gx4Qa3WlppG/fn5+REdHExYWZtEW0KVLF0qUKMHQof8lzm3fvj3nz5/n7bff5uTJk9SvX9/Mh7soEhsbS0pKCp6enjRs2JAJEyZgZ2dHvXr1zPb79NNPkVJy5coVBgwYwC+//EJYWBgjRoygSpUq/J2P9p+8JrfgwpEjR1KpUiV0JmrL0NBQypQpQ0BAAD4+Prz00ksWXTJVih7WWEJaWlj3v0fs1xZwEkLYAs6YezGpPCgVKsBnnykR0LNnKxHRc+cqyfReflmJbXgActKVZhf526VLF8aOHWsxIhNAo9Gg0ZjfbqVKlaJfv3789NNPNGrUiD///PNhr0C+k5ahpWu37uh0OooXL05iYiKnT5+mT5+s7zIZGRlcunTJGLOQlpbGtm3b6NKlC3379iUsLKygxX8oZC7BhTqdjjVr1lC5cmV27dpldmzv3r2JiIggMjISe3t7Vq5cWRinoPKA5GRUflkIsQqoJoT42WTZAjz0cC+lvAKEAJeAa8A9WYQK8DzWFC8OI0bA2bOwfDn4+yuR0M88o9RkWL/eKgN0TrrS7DyQPJt3YfLkyfj6+lol6o4dO4xvjffv3+fcuXNUqVLl4c89nzAMgKNmLeVmYgZ2Dk58vmwzJ05EsXXrVjO/83PnzhEQEMCMGTNwdXU12kUOHjxI69atcXJyomfPnqxdu9Zi7qTY2Fj69OljnIF06tSJ06dPZ5mBmKYLGTBgAJUqVSI1NRWAuLg4qlWrlifnnltwYXh4OD4+PgwZMiTbQS4jI4PExERKliyZJzKp5C85zRBOAxtRDMgbTZb5gOW0fVYghCgJdAOqo8Q4uAghXs3lmGAhhBRCyKtX1clErtjaQt++cPAgbN8OHTv+Fwnt4wOLFoE+w2JOZE7NnJMHUpwsxpCgYVaLePDgQQIDA/Hz86Np06a89dZbZnlgigqGAfDqhTNUqaU8mA0uuE2bNuXOnTucO3fOTFXm5uZGTEwM69YplWf/+usvtm3bRrVq1WjYsCG3bt1ix44dZv1IKenRowetWrXi3LlzHDx4kGnTpnHdioy4NjY2fP/993l41gq5BRcacvT06NGDjRs3mhnRV65cSUBAAJUqVeL27dt06dIlz+VTyXtyikM4IqUMBepLKZeYLGuklPGP0OfzwAUp5U0pZTqwGmiW0wFSymAppZBSiooVKz5C108ZQsBzz8Hvv8OxY9C/v2JXePttqFZNybh6+7bVzVnyQJq+7jBg7oHUqlUrNmzYYLZfdHQ07u7uxs/vvfceJ06c4OjRoxw/fpwxY8Y85EnmH4YBECBDpzOmujC44B6LPIFWq8XT05Pff/8dT09PIiIiGDp0KJ06dWLatGmkp6cTFRXFpUuXiI6OJjo6mrlz52Z5o965cyd2dnZmb+P+/v5Urlw5VzlHjhzJrFmzyMijAkwG+1BmO5JpAre0tDR+//13unfvjqurK02aNGHLli3GfQ0qo9jYWHx9fZkxY0aeyKaSv1hjQ5gohHATQtgKIXYLIRJze6PPhUvAM/p4BgG0AaIeoT0Va6hXD0JD4cIFeP99xUX1ww8VA/SIEcr6XHBxsMXJ3nJOFid7G1wcbPNY6MLlfko6J67Gsz/6FnftyxF1LIL01BRCBnfji0Fd6dunDzNnzrR4bL169UhKSiIyMpJ69eqZpT7u1q0b69evJzU11fjwPXL0GA0bNrTYlkEVZVjmz59vtr1KlSo0b96cn3766ZHON7N96ExaCbbv+ceY4sM0gduWLVu4e/cuvr6+VKtWjT179lhUGwkh6NKlS5G2D6n8hzUDwvNSyntAe+AKSnGcsQ/boZRyH/ALcAg4ppdh4cO2p/KAVKoEn3+uREDPnAmlSimR0DVrQp8+ipopG56WaE0D/5y/xZ2kNEBQoW4jdBnp1O87lu6fLOX9hetY8uOPeOnTiFSrVo3jx48bjxVCcOTIEQICAhg1apRZu6VKleL69RuEn71jfPjuPnOTi7cSzfIrAdy5c8f4f2xsrNET68svv8TZ+b/40AkTJjBx4kTGjRtndnxAQIBFw7clMtuHagY8w937SYyY/IVxH4PdJywsjEWLFhlnPRcuXGDr1q0WvYn27NlTJBK3qeTOg/yCWwKr9UbhR0oML6WcLKX0llLWk1K+JqVMfZT2VB4CV1clL9K5c7B0Kfj6KpHQgYHQurWiZrJggH4aojVBUZucvp5AOVdHdFIihODZwdO4cfIAS0Z2ZdbgLkz6cCLly5fPta133nkHDw8PPDw8aNq0KZD14VuuqhcRhw9nSQ9SsmRJoypq8ODBjBo1isGDBzN69Ggzr6369etTpUoVjh07xqVLl6hfvz41atTg9OnTbNmyhUR92dbQ0FA0Gg1Hjx41HluvXj3OnDufxT4khODN4Ln8tedPqlevTuPGjenfvz8ff/wxmzdv5oUXXjDu6+LiQvPmzY0ZYw02BD8/Pw4fPsxHH330EN+CSkFjzRz/hhBiHoohebreVTT7XK4qjxd2dvDKK9CvH2zbBiEh8McfSjR03bowdqyyTa/yMHggPeddlsTUDFwcbJ+4mQH8Zy+p4a4k9r15PxW74qVpNmgKbk62TH/Rj1LFlGtiOjMAjAWCQHkAZ8aScd4r4Bl+//5LFi/6judCJmJnozy07927Z7XMXbp0Ydq0adja2nL48GEmTZpEamoq3333HaNHj2bBggUAeHh4MHXqVDNX0OQ05XwNnmMGXEuX5eX3QhjSypMSzv9leu3fv3+W/levXm38f8CAAVbLrVJ0sOaX3A84BfSRUt4BPIAv81UqlYJHCGjbFrZsgYgIeO01OH0aBg6E6tVh+nS4e9e4+5NeHN5gLxFC4FmmGI2rlyKwWikaVy9F3YpuFHeye+i2LRnnhRC8ETyHqAN/UcurJj4+PkyYMCHHGUhycjLr1q1j7NixJCcns2DBAkqVKmXcvnLlSoKCghg/fjwrVqwwru/cuTORkZGcOnXKuM7J/umyD6lYJtdfs94baLaU8h/952i995HKk4q/P/z4I5w/D2PGQEKCEglduTKMGqXUhn7CyWwv0QiBo50NY9t7s272BONAaMhR1LlzZ2U/jQZbW1ucnJxwdnbGy8uLEydOsG7dOhwdHUlLS8PFwZbEuCtMfLERfyz7FoDoqAh+CB7G7asXcXBw4OWXX2btuvWUqVSVw0eOmskWHBzM2LFjcXJy4vbt28Y8Sp988gmtWrWiVq1aHDhwAHd3d6pUqcKbb75JfHw8t/UeZRqNhvfff5/PPvss2/M18KTah1Qsk+uwL4RoBnwB1NDvLwAppVTLaD7pVK6sqJA++ggWLlSioGfPhm++gd69FXVS/ZwLjT/OGOwihroQTvY2ODo5czvmHMnJyTg5ObF161YqVaoEwK1btwAoUaIEtra2ZGRkcPbsWfr160diYiJpaWmUL1+eKlWqcD76Ip7+TUhJuM+3773OhciDFC9Zhi69X2XRzKmEbt7L3J1njf3WqeBqscRnulZnMeFgWFgYJ0+epFq1asa0Er/++it2dsrMpl+/fkydOpULJt5lls73SbQPqWSPNcP+YuBboDnQCAjU/1V5WnBzg/feU1xTQ0OhTh0lErpBA3j+eUXN9AQWoLcUsW2jEXTq1ImNGzcCmBVQL126NHZ2dvTp04fBgwfTsWNHunXrRkREBFu3bqV48eLY2trSr18/Ggc2xKduHf79YzUtXnoLe0dnvlr6GzfOHGPUJyHEO5Q3GpylVILhzt1MNMqm00m0OmmWQuTYFcXeIKXk559/5tixY0RHRxMaGkrNmjXN3EJtbW0ZM2aMWaT1U5HNUyVHrBkQkqWUy6WU56WUFw1LvkumUvSwt1eC244eVbyQnntOiYTu0OE/NVNaWmFLmedktpf06dOHFStWkJKSwtGjR2nSpIlx3/T0dObPn8+UKVNYtmwZUVFRJOvTkickJKDVahk3bhwxMTFsWhlKba8afDn6dd4bM5oRvduRoc3gh6+mo8swv442GkFcQqpRpbM16jo6iVkKkZjbSVy9m0JSUhKVKlWiYsWKREdHM3bsWMaPH8+JEye4a2IHGjBgANu2bcuSVPBJtw+pZI813/jvQoiHTlWh8gQihJIOY/t2JW6hb184cUIZLGrUgBkz4AG8Yx4HcsvuKqXkblKaMdr4ww8/JCgoCFtbWxwdHQFwdXXF09OTChUqMG3aNOrWrUtlDw9KONvzcfBkDhw4QMfOXUlNSmT+uDeyyPC/PkN5Z9gIo5eSIUrcwDMdetKs12BiY2NJTU2lTp069OrVi+HDh/Pmm28SGxtLiRIljPvb29szfPhwbty4kU9XTeVxw5oB4R1goxDinhDihhDiphBCvYNUFBo0UNRHZ8/CyJGKJ9L77yv2h7FjlQC4x5jcsrv27t2Hfy/c5kJcIrO2nkarg0u3k5BSUrVqVeLi4oxv4Onp6SQmJrJlyxbee+89s5TRAJ6engwfFoSDswvXok+TGH/HbLvB2yenFCJOpSpw7dY9Dh8+TFRUFP/++6+ZC+iAAQOYM2eO8fPw4cORUuZZQjyVxxtrBoRAlER0fqg2BJXsqFYNZs1SBoDPPlNKf86cqcwYXnsNjhwpbAkfCkvZXf86epbN4XsRQvBS31eZOWUScXfj2bv/MNqMdLb8toovZ3/N2rVrycjIoHTp0oSFhZGYmMhbb72Fr68v3bp14/bt28ZU0hs3bkRKyeWL0Tg6OKCxscHJxdUoh6m3z9OWQkSl4LDG7fSipaUghFN5DClZUnFRjY5Wivd4eSmR0AEB0L49bN362Bigs8vu+uOn7+Lh14xLV2Np/+EPeDTrjC4jnb8XjAck6SlJJCbEs3evMmjs27ePbdu2odFojMnqgoODuXnzJufOnWPbtm389NNPeHl5ERAQQDFnRz6Y8S02tjYWo8EtuYiO6VCHGYO7MXtIV/r16W1MISGEMEscGBISYhY4t3TpUvz8/PDx8cHf35+33nrLzM6g8nSRUz2E/UKIf7NbClJIlaKFEIJXX/0vv2FmX/zQ0FCGjRkDAweyZeZMAqpXJ8DFhWJ//EHtdu0IcHbm9RYtlApvRRhLqpk3P/4WGxtbGrTvxbW7yVy7m0rN5t2o/mwX3D396DV/Ly9/u4faHQYwbOQYnJycaNasGa+99hq1atUytlO8eHHu3r3L7t27mTJlCocPKyqfMWPGcPHiRSYO6pujt0/mFCJ29o78sG4n509FYW9vb0yA5+DgwOrVq4mLi8tyfps3b2bWrFls2rSJyMhIDh06RLNmzaxKua3yZJLT3PKhE9ipPNm4uLhw/Phxi774mWnfsSPtz58HoFXDhoSUKkXgjh2wZw94eip2h7ffVor7FDEMqhkpQSclaRk6rkafxsPLByd7G0VtIyQgiL96npJVahuPrdt5IB+92pCfflhEfLySQrtNmza89NJLZn34+voSHh5usX+Dt48lMqcQmWqjfAZo0aKFMVeRra0tgwYNYtasWUydOtWsjalTpxISEmL87mxsbBg4cOADXyeVJ4ec6iHsymkpSCFVih7Z+eLnSPHiMG2aYoB+9124dUuJhK5cGcaNgyJW/MjORoN3ueKcuXGffy/c5kD0bS7eSuJuchre5YpTprgjFdyc0GVSgemkpIKbU7YP87yW0bSfjIwMNm3aZFa5LigoiGXLlmXJi5RbARyVp4+cVEajhRAOOWz3F0J0yB+xVIo6Ofni50r16krK7UuXYMoUJXHeF18ohuk33oDIyHyT+0GRAn1uX+Wh71axOjcvnEQKaPd8GyrEn6JMcQdcK1Tn9sWTnN6+kn3z3uO7gc/i71uPxMREWrZsyV9//QUoabJ79uxpbP+XX355pERwBnfY5ORkAgICCAwMNKarMODq6srrr7/O119/nW07x44dIyAgAE9PT7X+8VNMTkbldCBSCLFYCDFICNFNCNFbCDFJCLEXmA6cKRgxVYoCufniPzClS8PEiUpupO++UzySQkOVYj6dOsHOnYVqgE7X6jgVex+vcsVpVK00gdVK8UL7dtih5YfFi+jVuzdn/vkDX6d7NA7wIf7ica7+s57ho8fiWaM6NWrUYObMmURERPDss88a2z148CAnTpx4JNkyu8Pa2jvy+bLNHDp0mG+++QZ7e/PZyciRI1m8eLExDTaAj48Phw4dAhTVVUREBB07djQG0qk8feSkMvoG8AV2Ao2BIUAfFLvDUCllRynluQKRUqVQyc0X3yp1UU44OsJbbynBbb/9Bs2bw6ZNSiR0o0awYgXkUXnIB8HUqGyjUZLb2dpojFlJp02bxpIloXw67FWCXmjEd/Pnc+/GFcYNfImTJ0/SqFEjhg1T6kx//fXXzJs3j5iYGO7evUtgYCCjR49m0KBBrFu3jrp162apOFasWDHj/yNHjqRSpUrG2IWtUdf5/vsfGNW2NpeO7QOU9BYfzwlFCMEvv/xi1lapUqXo1asXixcvNq6bMGECY8eONbq+Aupg8JSTo9uplDJZSrlUSvmWlLKDlLKHlHKSlDKigORTKQJY8sU/FnMPz+ZdmDx5spm++pHQaKBrV9i9G/7+G3r2hEOHlEjomjXhq6+UzKsFRHb+/m6ly/Hm5K85d/4CNjY22NraMunDibzzZn8cHR1xcnJCCMHHH3+Ms7Mzjo6OnDhxAi8vL8qXL8/333+PRqMhIyODzp0707VrV3777Tfeeecds0L1BnQ6HWvWrKFy5crs2rXL6A6r0QgqVK/F4XDFlmOjEaxbvQo/f3+L5zNmzBgzb6NOnToxfPhwOnbsSN26dWnWrBk2Nja0b98+j66gyuOGmqxEJUey88W30QjiZDGGBA2zeFxoaKixSpiHh4fZW6hVPPMM/PKLUpNh6FC4cUPxSKpcGT74AK5de8gzsh5L/v63Y2P4/O3O1CpXnMTUDGxtbbl27RoxMTFIKUlNTcXe3h4pJVWqVOHAgQO4ublRqlQpfHx8SExMZOrUqZQqVYrdu3cb2/Xy8sLZ2dmsZKaB8PBwfHx8GDJkiBLgZjJzqVEvkEunjjJ19b+kJidyI+YiPvX8jMcmmAyg5cqVIykpySwOoX///hw7dowTJ06wd+9eFi5cSIUKFfLyMuY5sbGx9OnTB09PTxo2bEinTp04ffq0sY516dKljZ5dBrp3767aRqxAHRBUciSnNAnJaVpj6uVWrVqxYcMGQEmPkJCQQExMjHHx8PAgPDycwMDABxOgZk2YO1cxQAcHg62t4qlUrRq8+SZERT3qKeZIZn9/xbgsOXH1HrO2niZDJ9FqtcTF3SIjI4Ny5crx7bff4uDgwLVr1/Dx8cHJyYlSpUpRo0YN7t+/j1arpV69ely6dMkYBHbo0CG8vLwoWdqdu0lppGv/S2th8OLq0aMHGzduxF4j/5u5CEGt+s04eWAPx/dux/fZNtjaPLnZSaWU9OjRg1atWnHu3DkOHjzItGnTjLETzs7OtG/fnjVr1hiPuXfvHnv27KFLly6FJfZjgzogqORIkUmT4O4OkycrA8P8+VC1Knz/vVLms0sX2LUrXwzQmVNCVytdjDStjv3Rdzh+5R7a9DQQgnPnzpKRkYGfn/J2npKSQlpaGt7e3sTGxnLu3DnOnTuHVqulS5cuaDQa/Pz8WL9+PWvWrKFJkyZ0eHWI0U4zd+dZtDpJSkoqv//+O927d8fV1ZUmTZqwY9tW6lRwRaefuQS0eoHD4Rs5FL6Rl3v1RiOe3AFh586dxgSCBvz9/Y0R4AB9+/Y1qxC3Zs0a2rdvj7Ozc4HK+jiS64AghKhrYV3b/BFHpahR5CppOTnBO+8oM4PVq6FpU9iwAVq1giZNYNWqfDFAG3II7TpzE61WmhWrcS5dEVAGgcycO3eO1NRUUlNTWbt2LQD//qsE+l+5cgUnJyd69OjBxFmL+GLiKNJSU401EHQSPl+8krt37+Lr60u1atXYs2cPYWFhtK1TDo9SzgignKcPsdGnkSn3eaNjszw/96LEvn37iIqKIiAggPLly1OpUiUCAgLo1KkTGRkZ9OvXj+HDh7Nt2zYaNWrEmjVrWLFiBX379mXPnj00btwYb29vvL29WbhwYWGfTpHDml/zciGEsbCrEKIlMCeH/XNECFFbCBFhssQLIUY+bHsq+U/bOuWIP/kXIYO7MeOdroQM7sa8ET3p6FuRTZs2ATB79mwcHR3Ngp/Cw8Nxc3MjICCAgIAAnn/++bwTysYGevSAvXvhr7+ge3c4cAB69YJatWDOHDBxsXxU0rU6Lt5K5Hp8ClopuXI3mat3FY+cjPQ0anfoj0flyoSHh/Pee+8BSvbSRo0aERQUhKurK6VLl2bWrFmEh4dz48YNMjIy+OSTT/hu8feU8G5K5Vr12P/HGrN+1/26ivkLFxIdHU10dDQXLlxg69atpKQk41vJDV8PN4a08mT+VzP5dtaMJ7qYTbpWh8bOgX79XiEiIoLBgwczatQoIiIi2LhxI5cuXaJly5ZcuHCBAQMG0KNHD06ePMnhw4fx9/enX79+zJ8/n5MnT7Jnzx4WLFhgDK5UUbBmQBgF/CaEcBFCNEapoNb1YTuUUp6SUgZIKQOAhkASsCbno1QKE41G8MmIgUSfjuTE8aNEn45k3OjhtGjRwuiREhYWRqNGjVi9erXZsS1atCAiIoKIiAi2bduWPwI2awZr1sDJk8rs4epVJRK6ShWl/Ocj5OYxdbld8lc015JtSUu8bxadrMtI5+6Na2i1Onr37k1wcDBCCGJiYoiOjub69eu8++67xMfH8+qrr2JnZ8fhw4eN8RsGO027V4LYtfoHdDodWm0GNra2RO3/k/+1+c/rx8XFhebNm7N+/XpAqfVcwtmeLp1foHXr1g99nkUZ0+/gdFoJNuz8iy2RsWazNEMiQYMqqW/fvmzdupWSJUvSrVs3Fi5cyIABA4yR2e7u7nzxxRdMnz69UM6pqGJNttOdwFfAJmAp0ENKeSqP+m8DnFOzpz4eGNIkXDh3lk8++YSffvoJjUbDuXPnSEhIYMqUKVl86QuUWrUU+8KlS8pAAEokdNWqykBx6sFvW1OX2+JOdtg4OqFxKUlM5H7iUzIAQUbyfeKO7+bmjets2rTJOCimpaWRmJjI9evXmTx5MlJK/v77bz755BNKlizJnDlzGDt2rNFOU7lWPSZ8vwWNRkNs9FncK1XlszX7qVCmlJlMq1evpnfv3llqGxgIDQ3NkjPpccb0O/AJfJaMtDS+W7jQWFL06NGj/Pvvv8ZCRKA4OZw5c4a5c+fSt29fIiMjadiwoVm7gYGBRBahqPiiQE6pK4YaFqAE4Az8CbTUr8sL+gCF+ARReVDS09Pp168fM2fOpEqVKgCsWLGCPn360KJFC06dOmWWLXP37t1GlVHm5Gr5Rtmy8MknysAwZw5UqgQLFyq1oLt3V1RMVmDqcjumQx1mD+3OqXlD0aUlE7vxa84tDAKguGdDHN1K4eJSDFtbW5o1a4aDgwM1a9akfv36SCmND24bGxsaNGhA48aNjf1kttPs3RDG0mmjaf/6iMKx0xQhTL8DrU6SmqHj9UnfcDbib1b/9B0zZsxgwoQJuLq6mh337rvvkpqayqlTp/jf//5XSNI/fuR0pzXKtBwDbPivSM4jIYSwR1E9rbJi32AhhBRCyKtFLAHak44hXYXBDfKjjz7Cx8eH3r17G/cJCwujT58+aDQaevbsyapV/32lpiqjiRMnFqzwLi4QFKTEMqxapUQ9GyKhmzVTjNJabbaHm7rc2tk7MvLbtdQfsZB6I77H693v8R02D6GxwcW1BB0m/cSV6zc4duwY6enpxgf/q6++ytdff80nn3zCa6+9lm2aD1P31gbtezF+8e+81O0FYw2Ep5XE1AySUjM4e+M++6NvcSD6Nqfv2/LsoCk07fIqw0aMYuPGjbRq1Yp69eoZj5s7dy5RUVFUrFgRjUZD3bp1OXjwoFnbBw8exMfHp6BPqUiTrc+glDJrUde8pSNwSEqZq4JXShkMBAMEBgY+HtVVHnN0OsnWqOtEXYsnOU2Lk70NaZeP8+uvvxrz34CSFO3MmTO0bas4nqWlpVG9enVjyoYigY0NvPSSEvm8Z49S83n9euVzzZowejQMGKB4MJlgmv4aFM8qB1uNMWLb3VEidRk0fmUMpVyLkZaho0Tx4gQHBxMSEsKsWbNo3rw569atIzU1NccYjMzprF0cbJ/qmYEBFwdbrt5L5sb9VDRCGGMsbtxP5X5qOva2ikv0c889xwcffMC8efMYMmQIgLFIECgZX5s0acKLL75IQEAAt27dYty4cUyaNKngT6oIY43baW8hhKv+/0+EEJuFEHmRM7cvqrqoyJI5XUVi/D2mTxjO0OBZFDepXRAWFkZwcLDRC+bq1atcvXqVixeLoFlICGjRAtatU/ImvfWWolYaOlQxQH/8MZikdjBV5aSnpTDn3Rc5/NUgjs15h/sndpEUd4USlb3wKFuaOuVdzWIyEhISqFChAiNHjqRp06bMnTvX+KAyDeLLjMFOow4GpmTnOfXfeiEEa9euZdeuXVSvXp3GjRvTv39/Pv/8cwAqVKjA0qVLefvtt/H29qZZs2YMHDhQDVbLhJC5BPMIIY5JKX31HkbfoBiYh0kpH9rhWQjhAlwCakgp7+W2vymBgYHywIEDD9u1ihWka3XM3XnWLM5rW9gCti2fh3ulqrgXdzD+FO/du8emTZvw9vY27jt69GjKlStHkyZNCAkJyfbhVySIjVXsDN9+C3fuKLOEAQOUWUPNmqSma/n92DX6PFub4F8PcOVuMjqdpEYZFy6dPcmvsz5g1Nw1+Hq4cfXfTXz11VfcunWLvXv3UrlyZXQ6HcWLFzfLMvo4cOvWLdq0aQMoqSJsbGwoU6YMUkrS09PNsqkePXqUjRs30rFjR0BxQR4/fjzXr1/Hzc3tkeS4m5TGtzvPcvVeCjfvp5Kh1WFro6FMcQcqujkytHXNAqk78SQghDgopcxR3W/NgHBIStlACDERuCGl/M6wLi+FtRZ1QMh/7ialMS/8nDGRnSkp6VqGtPJ88n6ECQnwww/w5ZcQHY0UghttOrKzS38u1PQluGcgYX+doqNPef48G0fUtXjuxSfwxYA2LN22n26BnsYYgHr16rFhwwaqVasGKFlLEwowKV9eExwcTLFixRg7NmsRxYULF7Js2TJ27tyJRqPMapo0aYK9vT0DBw7kjTceTfNs+nJiqFpnb6tBIwRCQFDrmupsykqsGRCsuZJSCNEbxSPI4Ej+hD0NVEwpMukqCpJixZTYhTNnYMUK4uv6UW7b7/QZ0ZvX338VjTaD89fv8+fZOGMqixEd6jF40FusnfspaWmpAGi1WtLS0gr5ZAqG06dPm7kfA3nugmyqttMIJQW5RojCi5R/wrHmar6Lou9fJKW8IITwQqmRoPKEUuTSVRQktrakv/QyP8xayaoZP3K+8f+oFHmItPQ0lnTzp0O9Cjg5OVGpQnlaNWtMqZIlqFChAvXq1aNcuXLG2UDv3r2NUdyZ1UWhoaFFy+j+EFhyP4acXZAzc/36dfr160eNGjVo2LAhTZs2Zc2aNYSHh9O5c2dAuVYdfSviknjFmGDwi0GdKSPin3oPrPwg11c9KeVeoLvJ5zMog4TKE4zhx2bqZeTr4fZU/AgTUzNITtcR49+EGP8mlI4+w5Ffv6fO9vW4AfHFiqN5913FGF26tPG4a9eusXDhQhwcHLh+/Tq7dimlx11cXArpTB6OdK0uV08nS+7HoDgZrFmzxswF2dLgJ6Wke/fu9O/fn+XLlwNw8eJF1q1bR8mSJc329fDwIHzlApYuV1J/L3Wx53+1yjzRaToKi2wHBCHECCnlV0KILyxtl1K+n39iqRQ2T7MbZGZ301vVvNg6Zhp7B4xE99pziPQ0mDQJpk/HJjkZn1q1OHH2LN7e3syZM4dRo0ZRrlw5jhw5wmeffWasd5yUlMS5c+eYO3eusa/u3bsTGxvLP//8U0hn+x+WXI3rVHAls50xPDw8i/sxPJgL8o4dO7C3tzfLWlq1alXeffddwsPDzfbt3Lkzf/75J+fPnqF27drZ+hypPDo5/cINqRsTs1lUngKeRjdIOxsNlZwymDG4GyGDuzG597ME921B8MRBJGdk4JSSQkDFigSkp2MnJc+eOkU5Ozt2zJzJzJkzKVGiBKmpqUydOpVt27YhpSQiIoJGjRoxYcIEY4qFu3fvcvDgQe7du8f58+cL+ayzuhpLCX8dPcuSZSuZNm0aDRs2pG3btrz66qucPXvW6NIJEBcXh7+/Pw0aNCA6OpoBAwYghODAgQN4eXllsSdERkYa8wrlhkaj4f333+ezzz7L0/NVyUpONZUX6P9+nHkB1IxQKk80PZvW4Yd1O3lvwW806tiLlj36M3fVVpycnKhZsyYRV64QkZiIrYODojZKSaFshw40TUsjMSGBxYsWmbW3dOlSzp49a6xWJoRg9erVdOnShT59+pjl789M69at2bJli9m62bNn07FjR4QQfPPNN8b1w4YNIzQ0FFAKFTk7O3P//n3j9pEjRyKEMCulCZYr40kp+fHTd3EuXZ73x43n4MGD1KlTh1u3bmFvb8/MmTONaUnGjRuHnZ0dNWrUMB4/atQogoKC6N69u1l50HStjuS0DGM9B1ACx/z9/WnUqJHFa9CvXz/++ecfLly4kO11Unl0cnztE0JohBClTT7bCyFGAOfyXTIVlULEoDIb8j9P3Is5AHAy9j5p6RlERkZy9Fgkd9Ml0saG7+7e5ZoQ1HZ2JuLWLeylZOaoUWQsXMihf/5BSkn//v1JSEigW7duREZG4u7ubqyE1rdv3xw9cjIXfAHFeDthwgTKli3LV199la1nU82aNfntt98ApTbzjh07qFSpUpb9LFXGOxvxDzY2trw2eR7vDBsBwNdff01UVBReXl706NGDRYsWERERwblz55g8eXIW/f+XX37JjBkzcHZ25tat22ZZSzfu3MuWyFh0OsncuXPZvn07N2/etHgetra2jBkzxmxWopL35JTcrgtwF7ghhNgqhGgCnEbxOOpbMOKpqBQu4advcicxDaF3edRptdjaO/Lm5K+YF36OdK2OYq4laNmyJU1ffpnE0qXRaTT4pafz1jvvMLBLFwTQv08fIiMjGT16NOvWrcPf358zZ87QvHlzatWqhZ2dHcePH7cow0svvcTGjRuND31DRHjlypUpU6YMbdq0YcmSJRaP7dOnj7GWcHh4OM8++yy2tllNh5Zcja9Fn8HDyydbV2PDzOby5cvY2NhQsWJFizIYyoMeuSXNspamp6WyYP58tkYpnkimqSYsMWDAALZt25btoKHy6OQ0Q/gU6AW4AD8AO4BvpJTPSCl353CcispjjSGhX1JaBlHX4hH6kpSpyYlIqUOn1XJofSjfvPsi2ox0MnQSzwbNKVu2LLfj40nT6Thfpw5r7O2potHgCNxYuZKAMmUYPWwYgwYN4vLly9y5c4fq1atTrVo1oqOjs50llCpVisaNGxvdWFesWEGvXr2Mco0bN46QkBC0FhL11apVi5s3b3Lnzh1jEkJLmLoa66QkJV2LRCIhW1fjDh06sHXrVlasWJHF2whg1qxZ+Pj40KRJE8aNn2CmkhJCMDB4LheO7aff84E0amSeasIS9vb2DB8+nBs3bmS7j8qjkaPKSEq5WUqZIqVcDlyXUs4sILlUVAoc00Is88LP8fX2M5y4Gg8ouu5je7dTtUkHHIq5UaJybV4c/ik2tnZInZZk2+IMGRqEEIK2bdty5MQJHEuU4JizM+4lSrChQgUi4uI4cvo0Y/bvJ2zRIjZv3mzMAXXw4MEsaiHTTLOmaiNDSUgDNWrUoEmTJkb3zcy8+OKLrFixgn379tGiRYtsz79N7bJk6HT8G32Lf87f4q5DOS6fPk6b2mUt7m9vb0/Dhg2ZOXOmxfoLo0aNIjIykl9//ZVBb79NfIL5DMC1dFlenziLsYu3snXXHnbu3Env3r3Ncj1lrvkwfPhwYmNj+eCDD3KMXzAwYMAAfvnlF+PnuLg47OzsmD9/frbX4WkmNxuCkxDCWQjhDMRn+qyi8kSR2cvGzkbDnaQ0biemkaHVcXDHBjwatcfRzZ2SVWpzbPfvpKemknT/Hj/P+oianjVwcXHhrbfeAiAjI4PExERiExMJKFUKF3t7HKWk/8qVXDx6lGc++AB+/x2kpHr16ri5ubFv374sA9PcnWdxqfUM27dv59ChQyQlJWUp9vLBBx/w+eefZ3ERBejduzcfffQRbdu2NUYUW2L7qRvYajQ0qlqaZ2qU5oX27dBmpDNmyn/vgUePHuXy5cvGzwa9fqlSpSw1CUDXrl1pGNiQ4+HrLG5/kOh3Q/xCy5YtOX/+vHEgjYmJser4VatW8cwzzxRuIaciTE4Dgh+QYLKYfr6fw3EqKo8dlrxsBEr+nNh7KVyIvc3ZiH84uGw6yffiuPDXOvas/YmyVWrgVrocIZuOE3P1Glqt1ph3f+jQodSuXZuKFSsSceQIiampLF25EtG2LVfatkWEh8MLL4CvL/zwA4f+/psmTZpYdP88ezuDOg2aMnDgQLPZgQFvb2/q1q1rLK1pStWqVZk6dSpDh2Zf18r0/G00ir3E1kbDwOA57Nq5A09PT3x8fJgwYQLlyxtLrOPj40P//v1zvb7Bkyez97cfSc8wV2s9aPR7TvEL1hAWFsbMmTO5cuWK1YPI00RObqcaKaWN/q8m02fLiW5UVB5TLHnZnI9LRIPAq+MbOLmVpnLj9jw7cQUdP1vDjI3HqFa3Pi8P/xjHYq7UqeBKhXJlCQ0N5c8//wQUV8qkpCQzVU5ScjJUrAh//AGHD8Mrryi1oAcOhBo10E6bzvkzMbzfqS4hg7vx+dsvMGNwV3av/oEaz7TlyJEj9O3blz179tCtWzfOnDmDt7c3CxcuZOLEicTExLB27VqEEMTHxxv7TU5OpmbNmmSXGNLS+QO4lS5H73FfcvBYFJGRkWzcuBEvLy+LBnBT9U5wcLBZMryGDRty6dwZ/KuUNKagEIIHjn7PLX7BtEJfQEAA69b9Nyu5fPky165do3HjxvTq1ctobFf5j6cn2khFJQcye9nopOTm/VQ0GkGZ4g7ci9xFq3Yv0N6nPBVLOCGR1G32PNtXLsDRTmN8qHXv3p2kpCR2795N+fLlWblyJRMmTKBmzZo0a9aMX3755b/I3YAAWLoUzp9X0m3Hx2PzwQQG9vsfjkLDx5/MY9x3Gxk87Qei9v/JpdMnuJOYSokSJejXrx8//PADqamp7NmzhwULFhATE4NOpyMgIABfX18aNGhg1O2vWrXKWB0sOjoad3f3HM/flLxKaGhw5Q1qXZMhrTwJal2T9j7lrU5BYU38gmmFvoiICLp27Wrcd+XKlfTq1QtQPKRUtVFWnsC0lSoqD47By+ZYzD1sNIK0DB0ZWh0ajaBscQeGhfxk3LdSCWdeb1qNN5sHZ0npIYTgyJEjxs/PPPOMMadRtlSpAjNnwkcfoZ03n/SZs9Ak32Dg689zqlUnDr78Jr1GfsqsYS/hbP8NM+fOZcCAAcY3ZXd3d7744guCg4N54YUXAGVg+u233/jwww85d+4cbm5u2NnZWX3+BrQ6ia+HW55Gqhui363FNKXG6bQS7Nq5gk6RsbStU465c+cSFxeXYzU6A2FhYcTGxrJs2TIArl69ypkzZ/Dy8nroc3nSUGcIKip6TOsaa6UOW1tlMPAsU8xsPyd7G0o42+d9So8SJbCZMJ692/ajtbPnduUa1NmxnleHdOet2ZPQZKRxJ+4mkZGRWYzKgYGBREZGGj+7urpSuXJljh8/nq1baGZMz/9hVTr5galN5WHiF0BJ1Z2QkMCVK1eMnl0TJkxQZwmZUAcElccaGxsbAgICqFevHi+//LLx4ZCRkUGZMmUYP3682f6tWrUy6tEvXLiAl5eXMS2EqUpjWGsv+jetTnX3YkZ///Fd62cxghrqJwOkpKTQtm1bY3oKUFJMODo6cu+edYUB07U6Ar3Kga0tS79bz8qP53M5oAnVDu8lITkJnnsOLl3iwL591KpVi4sXLxIcHEydOnWIj48nICCA+fPnk5ycjK+vL76+voSGhtKjRw9AcQXNnDzOwKOqdPKDzMb+h4lfAGV2YLgGBnr27KkOCJmwSmUkhGgDeJruL6X8Nr+EUlGxFicnJyIiIgB45ZVXmD9/PqNHj2br1q3UqlWLVatWMW3aNOND3UBMTAwdOnRg5syZtG/f3mybQaXRwac8NhphzP4J2b8xp6Wl0bNnTxo2bGg2IISFhdGoUSNWr16dY/WwzJlGdRJqlXflmdEDKD7xbc5vWA/du1P21CmcdDrmHjrE/vffp6o+VUTXrl2JjIxk9+7dBAcH4+TkRL169bCxseHWrVu4urpafU0fVKWTnxiM3abV+wzxC5aq97Vq1crseENeJ0txEn5+fkRFReWL3I8ruc4QhBChKHWUmwON9EvuCjsVlQKmRYsWnD17FlAexCNGjKBKlSr8/fffZvtdu3aNdu3aMXXqVDOjY2YyvzHb2QiLb8wZGRn07t0bLy8vpk//L+/jg1QPy+xqCnDiajz7L97h7u1bDJ47Fzs7O3YvX86e4sVxkJL7n38OlSuzfdEiFi9ezOXLlwkICCAmJoZ58+YxaNAg7O3tSU1NNQakXbx4kQ4dOphF+xYrVsyiTEWBp7J6XyFijcqoGVBfSvm6lPIN/TIwvwVTUXkQMjIy2LRpE76+vqSkpLBt2za6dOliMXFc//79GTZsmMW3Rkvk9sb8xRdfYG9vz+zZs83WW1s9zFIMRHpaCrOGdqf/Cy1o8/zztGvXTgnKGjKE9Xv2sHLDBt6uVImqCQn8feUKzhkZlIiPJyMhwagCmzhxIs899xybNm0yC0hzc3Nj5szHI+nAU129rxCw5mpezn0XFZWCxZDWwVB8JjAwkCpVqvDmm2+yYcMGWrdujZOTEz179mTt2rVmeX6ef/55li5dapUx0rSv7GjevDl79+7l9OnTZusNuYNMq4dZwlIMwMzNUYyd/xsjvl3H9j37eGvocOzs7GjWrBmLFy+m5QsvsD8mhstSYmNri5ONDV537vDvuXNc9vHBo1gxxowZw8mTJ2nZsiUAU6ZM4ebNm2RkZLBy5Upu377NgQMHSE5Otuo6FBZF1dj9JGLNfOs0sF0IsZb/iuaoNgSVAuX69euMGjWKf/75BxunYqRLG5q/OBCNjS1V6jZg7dLFRlXOkCFD0Ol0FCtWDAcHB5KTk9mxYweurq4cOHCAbdu2ERYWxssvv8xvv/1mMfsnZNXrp2uVlBJt65QzUxu1bNmS/v3707FjR/bs2UOFChUeqHpY5gptBqSUXLmbxI9/R5OWIdFKwdvBXzP93Vf4dMpU3h39Hk5OThQvXhytVssRe3vek5K5W7ZwB0gXgvPnzxPg749XrVpMmTLFKEuLFi346quv6NKlSx58O/nL01y9r6Cx5qo6otQ/8CWPbAhCiBJCiF+EECeFEFFCiKaP0p7Kk41p/pp56/fyzpereO2DL0m6o+jBbyWkGl0Q4+PjuXfvHuvXrzf630+fPp3ly5czbNgwvLy8sLW1Zfbs2bi6uvLmm29azP8DlvX6x2LuGfsypWfPnowdO5YOHTpw9+5dwsLCCA4ONro4Xr16latXr3Lx4sUsx2anFjl7IwGdBI3QGPs/HZdG21FfMndRKG9O+IIMHTzXvR89e/bEKyCAH+7eJfSjj7jh4MDrUlJSSr69fp1VzZpxJy6O9PR0Y+rtJUuWkJj4+BQ/fBqr9xU0uV5ZE7vBG3loQ/gK2Cyl9Ab8AdXUr5Ithvw1b749yKhrL1WuEi26vwYorohR1+JJ1+pYs2YNJUqUwN7ennLlyjF27Fj++ecffvnlF3x8fHBzczMes2TJEq5du8b772ctD25Rr5+azJRX/0fvVgF4eHjw5Zdfmh0zZMgQevToQdeuXVmxYkUWN8cePXpkWxkts1pEJyUajaBmphiI6FuJRN2WvD3lO8JXzken0xKXkMaq1WvZs2cPTk5ODJ89m9s6HX4ffUSlUqUIi4uD0aOZGBWFVqulnJsbzs7OPPvss9mqsQqKqVOn4uPjg5+fHwEBAezbt8/MNdhAeHg4bm5uZmkptm3bVkhSP7mI7N6OjDso/nqDgOf1q/4AFsncDsy+PTcgAqjxMG0EBgbK7PKxqDyZfP3111y4cIHJUz9nXvg5MxfEf/9YzeXTx3nhnYlGF8RWrVoREhJCYGAgOp2Opk2bcuPGDQ4cOEDp0qVz6Ok/7ialZenLgCV3x7wiXasjMTWDdK2ORbsvmPWv1Un2R98iQytpULUkGiH4oLMPZT2qY2tnT/eObXAtXpzDhw9z+PBhNm7cSPXq1Qnw8yN60CBs5s6l2s2bdBOCbu3aMfnmTc7HxnLt2jV0Ol2en0tu/P3334wePZrw8HAcHByIi4sjLS2Nfv36Gb8/A+Hh4YSEhBjTYqs8OEKIg1LKHLU71sy9vgBeBtbql5eBR6ljVx24CfwghDgshFgkhHDJ6QAhRLAQQgoh5NWrVx+ha5XHDdP8NQZd++h2tfngxUBmDeuJQCB1Oqb2e5ZXXlbeyGNjY2nUqBHbtm1Do9Hwzjvv4O3tjbu7uzE3fqtWrahdu7YxD44hlsFAdu6OYzrUYc67L/Js4wZmgXCGADnDYnA/za2fzBjUIiWc7bP0n67Vka6VxKekc+TyXQ5E38bG1h6NYzHS09MID9/FqVOnsLW15c6dOzRt2pQKFSpw5do1Gv32G1y8CKVKgZsbz23ZgvbQIWqnpCgqs4d7v3skrl27hru7Ow4OSolSd3f3bKuuqRQM1gwI7YEOUsplUsplwAtAh0fo0xZoAMyTUtYHEoHxOR0gpQyWUgoppVBvmCcPSw/T7t174Oldj4pVqvPRhxOZv3Ah3j5+6K6dRGg0uBQvwf27t3F2LcHlMydIS04yBp8lJibi7e1tVM9oNBrOnj2Lv7+/Wb/Lli3jyJEjDB06lPfee89sW3Z6fTt7R75ft5PI48ext7c3FloxBMgZFtMI6Zz6yQ5L/dvZaIhPTkeiqLxsbTQIjYbmo+fRa/oqIo4cITQ0lH///Zf79++Tnp5OWloaixcvJiAgAJycqObry2tbtsDatXxYpw5nb9/mfwCNGsHKlZCRYZV8eUG7du24fPkytWrVYujQobnmfMqcyfTcObW0e15jjZeRwFAySkHq1z0sMUCMlHKf/vMv5DIgqDzZmEYbG9gSGcuxmHtcOPYvO1ctJvHeHfzb9qBO/UYIwKNmXc4c+Zeq3r5cv3gaN9fiABw4cID09HRat27Nnj17SE9PJyUlhXv37vHss89a7L9p06bMmDEjy3qDW6PBy8jJ3gaN+G99ixYtOHr0qNXnmV0/2ZG5f3tbDW7Odtn8+JSf6Jo1a3juueeMb90A3bp14/333yc1NVVZodFAt2506taNMt7ecOsWHDoEffpAtWowapSSjjufAtYMajEXJ2cOHjzI7t27jdXSTAP7MtOiRQtVZZTPWDMgbAE26SOWAfoDmx+2QyllrBDishCitpTyFNAGOPGw7ak8eWSXv2bN/Gm88vwidDodqTcvUqdWTUZ2CmChe2muX7/On3/+yciRIxk8eDBxcXE8//zzbNmyhcOHD1OtWrVs+9u8eTPdu3fPst7U3dEQh/ChRqDRCGMgXIcOymTZEA9hYMKECVkSymXXT3ZkdrdM1+pITddy5W4ycQmpZGglL361DfdiDlR0cyIxNYP+/ftnKVhTqlQpY2H6zHmMDp48qfxz5gzMmgU//AAjRkBwMAwZAu++CyYFcR6FzG68TvY21KngStuW/6NVq1b4+vqyZMmSPOlL5eGwZkB4H3gHeFH/eQ2w8BH7fRdYJoSwB84D2Sd5UXniyfwwHT76PZLdGxgNqlJKbmidqN13IlXTdawf8T/up6Rz8eJxfOp4k5ycTKVKlfDx8WHDhg2EhoYSFxdHnz59+Prrr7l37x5Lly7ls88+M+v3lVdeIS0tjYSEhGx1+zqdZMfJG8aHWGJiIg6OTmiEkvJh9OjRgDLLmTNnDqNHjyY+Pp7JkycbE9q98sor3Lhxg/T0dE6dOvXA18dgV0jX6nB2sKVm2eJUdy9GulaHnY0GG41ACB4tjYOXF3z7LXz8sfJ3zhz47DMICYHXXoMxY6BOnYdvn//ceA0V2a5fOs+1SwLwp71PeSIiIqhatarF4jsqBYM1bqc6KeU8KeVL+mW+lPKRXBKklBFSykAppZ+UsruU8s6jtKfyeGKIAM6sf3/tlb5mBtWkNC037qcCwjhI9Bw9nZcGDiMxMZGlS5dabL9x48YcO3aMuLg4atWqlWX7smXLOH/+PP3798+2BKOlWISpv0Ww9sAFli5dyqRJkwBl0OrXrx/z58/n5MmTxqI1t27dYtmyZQwfPpyGDRtaXerREqZ2BcND1UYj8jaNQ5kyMHmyYoCeN0+p1bB4MdStC127wp9/PpQB2pIbb2pyEj/PnED/Ts3x9fPjxIkTxsSAL7zwAh4eHnh4ePDyyy8DWW0IBgcBlbwj21cKIcQIKeVXQogZmNsQAJBSZnXeVlGxgtwigE2Lteh0kjStDo0Q6KSknKsjADYaQY1nO/NhzUp4e3tn29f06dNxdHTMdrsQgk8//RRPT09Onjxp1palh5ih76hr8ZS7e5eS+myj6enpFovWvPjii8Z+OnTowDfffJOlnwfBkl0jX9I4ODvD4MHw9tvw228wYwasX68sjRvD2LHw4otgY101XUtZSyvXqsfw2SuyuPFml57b2hTiKg9PTq8UhjQVCSieQJkXFZWHwpoIYEOgllZKY8nEcq6O1HD/z0PZwa0MbwwakmNfHTt2pHXr1jnu4+TkxJgxY7IYfE1zDOmkJCVd+T9kcDemvtGBdwYN4qOPPgKUAeGnn34yvr2OHz+ewMBAs3xJdnZ2Fvt5EAq8ZoGNjfLg//tv2LMHunWD/fuhVy+oVQvmzgUrop3VrKUKQgjGjBlj/BwSEmKWLn3p0qX4+fnh4+ODv78/b731Fnfv3jVuj4uLw87Ozujdltdk+y1IKRfo/10ppTxpuk0I8XCvNypPPTll9gRwL+5Axw4dmD59Ou19yiOvevDnr/Y0rl4Kjd6ttEY9JbbG8CC5C7i4uBg9UAYMGMCAAQOy9G3IjQ9Z30JNf6QGXBxscbTVcO5mAjfvp5Kh1WHr4ESPT5fh6e5CoHMcr7/+OsePH6dHjx7079+fbt26GY+/d+8exYsXJzAw0CibpX4ehkKpWfDss8py6hR8+SUsWQLDhsGkSRAUpPxftmy28hZUic6ijIODA6tXr2bChAlZ6lpv3ryZWbNmsWnTJipVqoRWq2XJkiVcv36dEiVKAEpt7GeeeYawsDAGDx6c5/JZ8y0st3Kdikqu5JTZc9g3q9m1918z18MObdvw1Q8rzNTWQSE/UbFmPaPevFq1avliiLSz0ZCq1XHtnpIN1Fb/0Lp2L5lUrY4WzZ8lLi6OmzdvUrduXQ4ePGh2/MGDB42F7Z8oateGBQvg0iXQz5D49FPF3vDOO5Ap66sBNWsp2NraMmjQIGbNmpVl29SpUwkJCaFSpUqAEp8zcOBAateubdwnLCyMmTNncuXKFWJiYvJcvmwHBCGEuxCiLuAohKgjhKirX5oCOUYWq6hkx8OoDgrrQZKu1eFgq6G8qyMgydAqvhTlXR1xsNVwLPIEWq2W0qVLExQURGhoqNFb6datW4wbN85inqQnhrJl4ZNPlIFhzhyoVAkWLgRvb+jeHf76y2z3oliis6AwTaEeFBTEsmXLsthEIiMjjTYoS1y+fJlr167RuHFjevXqxcqVK/NczpwUd68AI4GKwO8m6++hpLNQUXlgHkZ1UFjpjxNTM0hJ15m5ea5JT2XtpFeRUrLYxZ4lS5ZgY2NDhQoVWLp0KW+//Tb3799HSsnIkSPN0ktPmTLFrIhOfrzhFQouLorKaPBgWLNGMUD/9puyNG0K772neCjpDdBFqURnfmPJgeLvy0m89tprfP311zg5OVk87tixY7z22mvcv3+fzz77jN69e7Ny5Up69eoFQJ8+fRg4cGCeqSANWJPc7gMp5Wc57lSAqMntHn+yDVDKVGegsEnX6pi786xFL0shIKh1zSde97127Vp69OhBVFSU0TPq33//5f333+fKlSsUL16cChUqMH36dHx9fQFYuGABX372Gdy6hWtiIl8Czb28YPRo6N8fsnkIPokYIu4NLz/ju9Zn6tpDVCuuY2y/jrzxxhtIKQkODqZFixZ88sknZk4Qw4YNIzAwkAEDBtCwYUNiY2Oxs7MD4OrVq0RGRuLl5WWVLHmS3M4wGAghygohqhgWqyRQUbHA46I6sKZ8Y+Z6xKGhocYiOMHBwTg7O2dbv9iQw8ngUTJz5sxCyTqaE2FhYTRv3txYhvT69ev06tWLzz77jDNnznDo0CEmTJhgzCu0YcMGFixcyJ6DBzmZkMD8X36hn4sLsdHRSuRzlSpK8FtcXCGeVcGQk9tyTJINPV96mcWLFxvXT5gwgbFjx5rNHA3V7E6fPk1CQgJXrlwx1tiYMGFCrrW6H5RcBwQhRGshxGWUHESngAuA+oqu8sg8DgVPHtV+4e7unm39YkNAXmRkJFu3bmXTpk18/PHHeSn+I5GQkMCePXtYvHixMVHgnDlz6N+/P82aNTPu17x5c2NKjs8//5wZM2YYPWga9OxJ/1GjmBsUBBMmKMnzgoOVgSEoCJ7gBHWWHCgMJKdpeSdoOHEmA2OnTp0YPnw4HTt2pG7dujRr1gwbGxvat29PWFhYlvoaPXv2zPMBwRqV0UGgL7ASJUvpm0A1KeWHeSqJlagqI5XCwJiQLZP9olixYiQkJBg/h4aGcuDAAebMmWP0Lw8NDeXQoUOUKlXKbP/Mx54/f55GjRoRFxdnzNxamCxbtowdO3awePFimjVrxjfffMPUqVOzuNeaUqpUKS5cuGAsRATw22+/sWTJElavXg0JCUrk86xZSjS0EEqcw3vvQZMmBXVqBUJRUznmVT0EpJSnATupsIhHS3+tovLYkd1sxpCHybAYUlkYKFasGAMHDuSrr77KtY8aNWqg1WrNVEy5kVOgU3BwMCEhIWbbvL29CQgIoFGjRvz4448W2zR4xCxfvpw+ffoAihHT0ttokyZNqFOnDiNGjLBO4GLFlOR5Z89CWBjUrw+//grPPAMtWyqR0EVMbfawWKNyLGpYI1G6/u8VIUQXIYQvUCofZVJRKfJkl4fpk08+ybLv8OHDWbJkCffv389zOQyBTnG56OTnz5/P1q1b+ffff4mIiGD79u1ZaknrdEoKkbk7zzJz/UG2bt/Bq/3foFq1asyYMYOff/4ZHx8fDh06ZDxm3759fPrpp0YXSqvjMWxtlXTbBw7Ajh3QsSPs3q14I/n4wKJFkJLC487jFnthzYDwlRCiJPAhMAvYAUzK+RAVlScT04fmvPBzxjxMOl32qtcSJUrQr18/5s6dm2Pb58+fx8bGhrLZRPtaIqdAJ1M+++wz5s2bh6urKwCurq5Z0mSbphQ59fc2Grbpxns/bGfBxn+4fPky1atXp23btoSGhrJ3717jcabpOd5//33GjRvHrVu3AIiIiCA0NJShQ4daFkwIaN0afv8djh2DAQMUu8Lbbyu1GaZOhdu3rb4eRY3HxYHCQK4JRKSUhnnifqBm/oqjolK0sbOzxd2jBvG3rpOemoJWm8GglzpAejIOGsnFixdJS0vj+++/R6vV0q5dO44fP27MyKrT6fDy8kJKSWJiInFxcbi7u3Pz5k0GDx7MsGHDHth+EBQUhJ+fX7ZBcPHx8dy/f58aNWpk20Zmj5hD4Rt4rtfbxkR+z3mXNRoxV65cybhx47hy5Qply5bF3d3dqCrr2rUrV65coVmzZgghKF68OEuXLqVChQq5n0i9eko9hqlT4euvYf58+PBDJQ33m28qhXuqV3+ga1NUeFxiL7I1KgshnpVS/iWE6GRpu5Tyd0vr8xvVqKxSWKRrdbi4FKNijdo0atudUuU9+O7DQUxcso1DO9YRuWMNnTt3JiwsjFu3bjFw4ECjG+mmTZsoWbIkp06dQqvVUr16dS5duoS3tzcajQZbW1tee+01Ro8ejUaT+8TdYOT2KFuKhIQEJk2ahJ2dHU5OTiQkJBAcHExwcDDFihVj0KBBVK1alTt3ss8yfzcpjXnh58yykRrInI20wLh/X1EdzZoFly8rld5eekkxQAfmaBtVscCjGpUH6P++Z2EZmxcCqqg8TiSmZiB1Wmxs7WjWuS8pSQn4NGlNqXKVSE5OoXe/Vxk7dqzxbfj77783JtSrV68eGo0GKSXh4eF06tQJFxcXXn31VSIjIzly5Ahjx47NdTDITmU1fPgIFi9eTKKFzKOurq4UK1aM8+fPZ9tukcxGWry4Mis4dw6WLQM/P/j5Z6X+c+vWsHHjE2OALipke/dJKd8WQmiAkVLK1pmW5wpQRhWVIoGLgy0ZGelcv3iW6QM78POXH9L2FUU3Hnf5HI0bNcz22NKlS5OYmMidO3cICwujT58+2Nvbc+bMmQeSIbvU4Qevp9OrVy+zQCdTJkyYQFBQEPHx8YASY2DqZVSkPWLs7KBfP6Xu89at0L49hIdD587g66uomQz1olUeiRy/ZX1ltJ8KSBYVlSJLUloGV+8mY2trS4M2XRn//WYGfbaI+eMHMqn3sxz/eztjRo2iU6dOJCUl4erqSrVq1RgxYgQJCQls2LABHx8fVqxYwb59+2jRogXx8fHGWsfWkFvBnuEjR2XrbTRkyBBat25No0aNqFevHi1atMgyG8lrj5jcXGIrVapk5rJ79+5dwsPDcXNzIyAgAG9vb8aOHWvaIDz/PGzeDEeOKKU9T5+GgQMV28L06ZCDWkwld6yZB54VQlSTUkbntzAqKkWNjAwdIVtPcfjyHVLSdOgQnD1xFKnTUd7LD5CULVuOOG2aMWCrb9++NGjQgI4dOxIVFUWxYsVwdnamVKlSfPTRR/Tv35/Tp0+j0+l49tlnrZbFUtWx6esOA0rka7ESpc08fkwLrwgheP/993PMvprXSQRzyv0PMGrUKPMHvp4WLVqwYcMGkpOTqV+/Pj169Mh6nfz84McfFYPzV18p6bgnTFAM0m+9BSNHQtWqDy3704o133Zx4KgQ4nchxM+GJb8FU1EpCoRsPcX+6NsglRrGQmhISU0hfN1y2lRIxwYd1SqVw9XVldDQUC5cuGBMPpZi4kdfuXJlzp07x9SpUxk6dCiTJk1CSmkssWkNBaXnz6uUIta6xGaHk5MTAQEBXLlyJfudPDyU7KqXLyt/3dxg9mzw9FTUTIcPP5zwTynWfONLgXdRUldsNFlUVJ5oktIyOHz5DrYmqhVdehra1BQ2/zCL5k0CKVeuHPv27SMjI4PAwEA+/PBDTp8+zeHDh9m+fbsx0V3lypX5999/+fbbb+nQoQObNm2iTJkyFt+cs6NI6/mzIbvc/wCzZs0yqosslTm9c+cOZ86coWXLlrl35Oam1Hk+f16ZOdStq0RCN2jwn5oplzQ9KtbFISzJ606FENHAfUALZOTmCqWiUhjcvJ9KSprOTEXz8rw9ACSnZTC1hy/eFVzRINm9ezc7d+7kr7/+YuHChYSGhhISEkJgYKDR06ht27aMHTuWcuXK8dprr3Hs2LEHlsmgzzdNHV7UIl8NLrGgeDi9/vrrFnP/Z6cy2r17N/7+/pw5c4aRI0dSvnx56zu3t1dsC6++Cn/8ASEhsG0bbN+uxDmMHQt9+yr7qWQh1wFBCFEB+BowDOE7gBFSymuP2HdrKeWTnwNX5bGlTHEHHO01YPJiKaXkTlIaSela1hyOwfWkvVLLoeX/aNWqFb6+vixZYvkdqm/fvqxYsYJy5crRt2/fh5KpsIoFWYOlYjAGl9jAwIa88cYbVrVjsCFcuHCBZ555hl69ehEQEPBgwgiheCO1b6+ojUJCYOVKJRJ64kQln9KgQcrMQsWINXfST8AxwE+/HEX1PFJ5CnC2t6V+5ZJkmPi630lK435KOuVdHSnuaM/1S+fZse8IW6OuA0qqhqrZGDNffPFFfv/9d1auXGlMGvewFMXU4Q/rEpsd1atXZ/z48Xz++eePJlj9+kocw/nzSlzDvXvw/vtQubIyY7h8+dHaf4Kw5m6qIKX8REp5Vb9MAayIQ88RCfwhhDgohBj0iG2pqOQbY9vWplG1UiAkyWkZJKVrqVTSmWY1FN1/anISP8+cQP9OzfH18+PEiRNG754XXngBDw8PPDw8ePnllylRogRNmzalXLlyOaaReBx5GJdYUxtCQEAA0dHRWdodPHgwf/75p8VtD0yVKvDll8oAMH26knl15kyoUUNRMx058uh9POZYUw/hN2CMlPKs/rMnECKl7JHjgTm3WUlKeUUIURbYCrwrpfwzh/2DgckAFSpU4OrVqw/btYrKQ5GUlsG5GwmsORxDcces+udCS+9QRCiSqS9yIzVVMTyHhEBkpLKubVslNcbzzytqpyeIvKqH4AQcEUJsEUJsAY4Azo/ifiqlvKL/ewNYAzTOZf9gKaWQUoqKFSs+TJcqKo+Es70t3hVccXWy/FArtPQORYQimfoiNxwcFJvCsWNKGozWrZVI6HbtFDXT0qWQnp5rM08S1gwIy4ChwHL9EqT/+1Dup0IIFyFEccP/QDvg+IO2o/L0MXXqVHx8fPDz8yMgIIB9+/ZRrVo1M1VEeHg4nTt3BpRKZWXKlDFGvT6sP7yBx9Ht0xpiY2Pp06cPnp6eNGzYkE6dOnH69GkiIyN57rnnqF27Nl5eXnz66admNRQ2b95M48aN8fb2plHDBqwJeY+4WPOYgcfi2ggBnTopdRn274fevZVB4rXXFHXSzJmgT/nxpFMYbqflgDX6FL+2wHIp5eY87kPlCePvv/9mw4YNHDp0CAcHB+Li4khLS8v1uN69ezNnzhxu3bpF7dq1eemll6hcufJDy/E4uH0+CFJKevToQf/+/Y11k48cOcL169cZMGAA8+bNo127diQlJdGzZ0++/fZbgoKCOH78OO+++y7r1q2jTp06AKxd+xtn791FiEqP77UJDIQVK2DaNCXAbdEixfD86afwzjswfDhUqlTYUuYfUsocF8AdWAHc1C/LgTK5HZdfS8OGDaXK08evv/4qO3funGV91apV5c2bN42fd+7cKV944QUppZQ//PCDDAoKMm5r0qSJ3LdvX57Ik5ahlXcSU2VahjZP2isstm/fLlu0aJFl/aJFi+Rrr71mtu7s2bPSw8NDSinlq6++Kr///nuLbT4p10ZKKeWtW1JOnSpluXJSgpR2dlL27y/lsWOFLdkDAxyQuTxfrZnHLQBOAwFAfeCMfp2KSr5jKFXZus3zXL58mVq1ajF06FB27dr1QO1cunSJlJQU/Pz88kSuouj2+aCka3XsPxRB/foNsmyLjIykYUPz7K2enp4kJCQQHx9PZGQkDRpkPQ6ejGtjpFQp+OADiI6G775TUmIsWaJkWe3YUVEzPUER0NZ8Y55SyklSyitSyhgp5WTgyfKZUylyZM77v2R/LJ/9uJH58xdQpkwZevfuTWhoqMXqYqbrVq5ciZ+fHzVr1mTo0KE4OjoW5GkUCNevX6dfv37UqFGDhg0b0rRpU9asWUN4eDhCCNavX2/ct3PnzuzYsZMtkbHUadCUT6dMZfGSH6lSw4v58x/uPe/WrVsEBARQq1YtQkJC8uq0ihaOjkrSvMhIWLcOWrZU0mG0aaOomcLCICOjsKV8ZKwZEDR691AA9P8/AUO/SlEmc5CTlHDiWgLp5erw8ccfM2fOHH799VdKly5tVgns9u3bZvmBevfuzdGjR9m7dy/jx48nNja2ME4n35BS0r17d1q2bMn58+c5ePAgK1asICYmBgAPDw+mTp1qdsyBi3c4FqPkFmr/+ruUq+zJW18s5b1x44x2mbp163Lw4EGz486fP0+xYsVwdXXFx8eHQ4cOAUqth4iICAYNGkRCQkJ+n3LhotFAly6waxf8849SwS0iQkmkV7Omknn1Mb4G1jzYQ4DDQoiFQoiFwCHgi/wVS+VpxlKQ043L57l97SJR1+JJ1+qMEcGtWrXip5+UwHmtVsvSpUstJkoLDAzktdde46uvviqw8ygIduzYgb29PYMHDzauq1q1Ku+++y4A/v7+uLm5sXXrVkAZQC7dTjRe2yq16pGRnsa/v/+Mjb0jOgRHjx6ldu3a7Nmzh23btgGQnJzM8OHDjemz33//faZOnUpUVJSxX9PU208FTZrAqlVKTYagILhxQ0m7Xbmyoma69qjZfQoeawaE31BcQ4/ql/ZSyqX5KpXKU40h778pqclJhM0Yz2cDOxLg72+MCP7oo484e/Ys/v7+1K9fn5o1a/Lqq69abHfcuHH88MMP3L9/vyBOo0DISZdvYOLEiUyZMgWADJ0kNf2/VBzLPn+PlKQEfl8ym/S0NOoH+DNhwgTKly/PhQsXaNeuHQ4ODpQoUYLy5csTFxfHBx98QOPGjblz5w4+Pj5oNBpcXFxYvnw5586do1KlSuh0yqAthGDzZuucCIUQZt9dRkYGZcqUyeJGXL9+fby8vGjfvj179+41a+PLL7/E29sbX19f/P39GT16NOn5HUvg6Qlz5sClS/Dxx0qFt2nToFo1ePNNMBk0izo5DghCUcb+LaWMlFLO0S+RBSSbylOKpSCnyrXqMXz2Cj74fhMRR46wevVq3N3dcXNzY/ny5Rw5coSjR4/yxRdfGCuBDRgwgDlz5hjbqFixIrGxsRQvXrxAzyc/MBjbM8dEBAUF4e/vT6NGjYzrDOmj9+zZg61G4GD338/+1fEhjF+8iY9X/IVzMVfWb9jIxo0b8fLywsnJibp165Kamsq6des4c+aM0T5TsmRJGjVqxO3bt9HpdNy5c4f+/fsTHh5O5cqV2bVrF2FhYTRv3pywsDCrzsnFxYXjx4+TnJwMwNatW6mUycWzd+/eHD58mDNnzjB+/HhefPFF4yxl/vz5/PHHH/zzzz8cO3aM/fv3U7ZsWWN7+Y67O0yaBBcvwvz5SoGe779XUnEb1ExF3ACdWwlNCVwWQpQsIHlUVJ7YALC8ILOx/UxaCbbv+Qed/lrNnTuX7du3ZynNaZglCCGoUsoly7V1ci1JXV9/Dh3Yb7Hf+Ph4Spb87zFw8+ZN5s2bR4kSJQCwt7enSZMm+Pr6MmTIEJYvX86qVasIDQ1l69atZsWCcqJTp05s3KjEu4aFheWYFbZ169YMGjSIhQsXAkrgYmaZxo8fj6urq1V95xlOTkrMwsmTsGYNNGsGGzZAq1aKmunnn4usAdqaX9Y9FBvCt0KILwxLfgum8nST1/V9nxQyG9trBjzD3ftJjJj830/Ski6/Xbt23Llzh6NHjxJYtSS+Hkra59QMHUKAV2k7rp2PwtPTk6S0DC7eSiQlJYVz587h7e3NW2+9xUcffQQoleAyMjLo0aMHAQEBBAUFAf89wHv06MGaNWuoVq0anp6etGrVyviQt4RhtgPQp08fVqxYQUpKCkePHqVJkyY5Xo8GDRpw8uRJ4uPjSUhIoHr16g92QfMTjQa6d4e//lKWHj3gwAElErpWLUXNlJhY2FKaYU2CkUj9oqJSYBTlvP+FhSVjuxCCN4Pn8tuCaVSvPp8yZcrg4uJiMWX0xIkT6datm/HaVirhxB9zJ7LH2Zm01FRef70/W28480Xov6Sk6dDY2eNWzoPjx0+wf/8+Xn/9dV566SUlgEmjISIiAoAtW7bg7+9PZGQkvXv3xtXVFTc3N3x8fADlIf/jjz/Ss2dPM3ks1U+4ZlOW6OhowsLC6NSpU67XRGajgtmyZQvjxo3j7t27LF++nGbNmll7mfOHZs1g9Wo4c0bJuBoaCu++C5Mnw9ChMGwYlCv8lx1rUld8XBCCqKhYwhDkpPKfsT1zRlHX0mV5+b0QixlFW7VqZfy/a9euZg/QXbvCzfadvimK/dG3sdVo9H0I7t29S8jWU4zv2JS4uDiuXLlCnTp10Gg0XLhwgerVq9O+fXvS0tLo3r07r776Kg4ODsTExPD999+zfv16pJTcunWL+/fvm9lvDLMdG40wq59Q95nnGDt2LOHh4dy6dSvHa3L48GHq1KmDq6srxYoVM5Opffv2dO7c2aoUJwWGlxfMm6cYn7/9VpklTJmi1IN+/XUYMwZq1y408XJ95RJCFNeriQ7ol88NyelUVFQKjvzMKGqpfrQQAic3dzZv3cbho8dJT0/nzz//JDAwkDJlyjBkyBDu3r0LwPLly6lTpw4rVqzgu+++o3Xr1ri4uHDixAkuXrxIz549WbNmjbHtnOon1Hi2Mx9+NAlfX98cZd61axcLFy7k7bffBmDChAlmMkkprbZdFDhly0JwsOKZ9O234OGhREJ7e0O3brBnT6EYoK2Zg38PlAaG65dSwA/5KZSKikpW8tPYbqgfbdZuWioZqUnsmTuOpk0aUaJECSZPnkzVqlUpWbIkbdq0oUmTJtSrV49Vq1bx3HPPUb9+fcLCwnjppZdo3ry5MUq6Z8+eZt5GllyLDTi4leGNQUMsblu5cqUxKvqzzz7j119/NSbXGzJkiFEmPz8/nn32WerXr0/9+vUf+rrkO87OMGQInDoFv/yiGJ3XrYMWLRQ106+/gtbydcoPrCmQEyWlrJPbuoIiMDBQHjhwoDC6VlEpdDLr3Z3sbZSaznXKodFkTeNhLUlpGbwR+i9IC20IyQ8DGuNsn3c1DdK1OubuPGvxJVgICGpd8+m0GUmpGKBnzFAGBlDiHEaPVmo3ODs/dNN5VSDnqhDCmAtACFEauJLD/ioqKvmEwSAc1LomQ1p5EtS6Ju19yj/SYACW60cDZOh01K9cMk8HA1Bdi7NFCGjeHH77TQloe/ttiIlRIqGrVFHUTJm+o7zEmqseh1IxbYEQYgFKxbSbqvupikrhkR8ZRU3rR6eka0FIGlUrxdi2+WPkVF2Lc8HbGxYuVALdPvxQGQj27lXcWfMJa1RGk3PaXtBeSKrKSEUlf0lKy+Dm/VTKFHfI85mBJdK1OtW12BoSEyEuTomAfgisURk9lNupEEIjpcy/eYuKikqh4WxvS9XSBVcDWXUtthIXF2XJR6xxO12ur31s+FwR2JmvUqmoqKioFDjWzM9OAQeEEAFCiI7A30BovkqlovKYM3XqVHx8fPDz8yMgIIB9+/bRqlUrateujZ+fH97e3gwbNszoMw+K3/+YMWOMn0NCQggODi544VWKNNZkhR02bJilQ31NHYQskeuAoFcZDQf+Ar4DOkgp1TgEFZVs+Pvvv9mwYQOHDh3i6NGjbNu2jcqVKwOwbNkyjh49ytGjR3FwcKBbt27G4xwcHFi9ejVxcXGFJbrKY4A1WWEfFmtURpWBT4CfgavAu0IIVeGnopIN165dw93dHQcHBwDc3d2pWLGi2T729vZ88cUXXLp0iSNHjgBga2vLoEGDmDVrVoHLrPJ48SBZYR8Ea1RGfwFzpJRvAM8C6cC+POldReUJI12ro3HzVly6fJlatWoxdOhQdu3aZXFfGxsb/P39OXnypHFdUFAQy5Yt4969ewUlsspjwqNkhbUWa1wJ2kopTwFIKdOBEUKIbrkckytCCBvgAHBFStn5UdtTUSlMMkcQ9/88jIyrUSRciKB3795Mnz7d4nGZ3b5dXV15/fXX+frrr3FycioI0VWKOHmRFdZash0QhBBVpJSXDINBJi7nQd8jgCiggKtXqKjkPZYyd2oq+vBs42b4+fmxZMmSLMdotVqOHTtmzMVjYOTIkTRo0IA33nijQGRXKdrkRVZYa8lJZbTW8I8Q4t9M2xY9SqdCCA/ghUdtR0WlKJA5c+eNy+e5eSUaG40g6lo8Bw8dpmqmYKL09HQmTJhA5cqV8fPzM9tWqlQpevXqxeLFiwvsHFSKJnmRFfZByGlAMJXALodtD8Ns4H3AquA2IUSwEEIKIeTVq1cfsWsVlbwlc+bO1OQkwmaM5/O3OvH52104HhlpdB995ZVX8PPzo169eiQmJvLbb79ZbHPMmDGqt5HKQ2eFDQ0NxcPDw7jExMRY1V+2qSuEEIeklA0y/2/p84MghOgMdJJSDhVCtALGPogNQU1doVLUyI/MncWKFSMhIcH4OTQ0lAMHDlChQgVWrVoFwLFjx4xvhwMHDuT27dsUK1aMsWPHPvzJqBQp8vLeetRsp45CiDpCiLqm/xs+WyWBZZ4FugohooEVwHNCiKWP0J6KSqFSkJk7J06cSEREBBERETg5ORn/Hz58eJ71oVJ0KOissDl5GTkDv5t8Nv3/oUv5SCknABMATGYIr+Z0jIpKUceQodO0ToGauVMlLyjIeyvbAUFKWS3Pe1NReUIx1Cl4zrtsnmTuTE5OJiAgwPj59u3bdO3aNQ8kLTpkVosBBAcH891331GmTBlSUlJo3bo1c+fORaPRMGDAAHbt2oWbmxtSSr788kvatGlTSNIXHHl9b+XYV760aiVSynA1BkHlSeJR6xQYgo9M1UERERF88skneSxp0WXUqFFERERw4sQJjh07ZhbYN2PGDCIiIpg9ezaDBw8uRCkLnvyogZGZgstxq6Kiki2Wgo+2RMY+cmnMx5m0tDRSUlIoWbJklm1NmzblyhW1cGNeo1ajUFEpAhiCj6TELPhoa9T1QpYs7zFNwWCJWbNmERAQQIUKFahVq5aZ6szA5s2b6d69e/4J+ZSiDgj/b+/eg6MqzziOf39cAoiKvVhFKQYYuXSijVasrYZRgYpKRWqLGKugrR2lWKmtijOOorVMUXCsxVqpVpAyCYrEsVhEUaKIdwmCN0gpQREqXrl5AZOnf5x3YRN2Nwlscjbk+czsuOeS9zy7wTznfc97nuNczDLdfPTWhs3sqG66Z1GVlZVRWFhY69WmTRvmz58PwO23307Hjh1r1VYqLy+nS5cuFBYW0rdv3wZPc62piXo9dy76D3eVr97ZC6qpM4MmMWS0ceNGtm3bRmlp6c5tV111Fb1796a4uJhrrrkmC9+AS+YJwbmYpbr56E+PVADw+fZqtn35FaNHj2bq1Km19kl1Qbax9yAMHz681rWKMWPGUFRUxGmnnQZElTT79+/P3Llza/1cUVERy5Yto6Kignnz5rFkyZJ6j9XYXlD79u0ZMmQIzzzzzM51t956K6tWrWLSpElcfPHFjfqsrn6eEJyLWecO7eiU1zbltk55bencoXku9a1atYqbbrqJmTNn0qZNG1avXs3WrVu5+eabKSkpSR1fp04UFhbWO56/J70gM2PJkiX06tVrt21jx46lpqaGBQsWNOITuvp4QnAuZs1981FCYix/R3UNO3bsoLi4mClTptC9e3cASktLGTlyJEVFRaxcuZL339/9TP6TTz6hsrKSAQMGZDxWql7Qji8/58biAdwwcgD53btz2223AbuuIRQUFFBdXc2YMWN2a08S1113HbfccsuefnyXQtrSFbnKS1e4fVHdWUad8trSr+uBTTLLKNWxni/9Cx22b+b++3dVZS0oKKCsrIwjjzySK6+8kp49ezJ27FjKy8sZNmwY+fn5VFZWMm7cOCZOnJjxmE1R3sM1TkNKV/i0U+dyQHPefFS3nHLlshd5av6/mFa2cOc+K1asoLKyksGDBwPRFNAePXrsfFZvUVER8+bNY82aNZxwwgmMGDEi5WyghEQvKHHchOoa46huXTwZ5Aj/LTiXQ5r65qO6Y/mfbdlE6eRrKb56Ems31+wcyy8pKWHChAlUVVVRVVXF+vXrWb9+PWvXrq3VXo8ePRg/fjyTJk2q99iD+x3CUd26IMEXO6qR8PIeOcZ7CM61Iomx/MQsn+fmlbL104+Zc8cEzIx79+9A2zZi06ZNO6eeJgwfPpzS0tLdHtd46aWXMnnyZKqqqsjPz0977ObsBbk949cQnGtFfCy/9drb8tfOuX1MXDOaXMvgv33nWhkfy3fp+DUE51oZH8t36XhCcK6VSsxoci7BTwucc84BnhCcc84FnhCcc84BnhCcc84FLe7GNEkfAGvr3bF+hwHrs9BOc2hJsULLitdjbRoea9PZ03iPMLODM+3Q4hJCtkgyM2sRD6ttSbFCy4rXY20aHmvTacp4fcjIOecc4AnBOedc0JoTwo1xB9AILSlWaFnxeqxNw2NtOk0Wb6u9huCcc6621txDcM45l8QTgnPOOcATgnPOucATgnPOOcATgnPOuaDVJQRJ/5C0UdLrccdSH0nflrRI0puS3pB0RdwxpSOpo6SXJL0WYs35qXyS2kqqkDQv7ljqI6lK0gpJyyTl9EPFJR0kaY6ktyW9JekHcceUiqQ+4ftMvDZLGhd3XOlI+m34f+t1SSWSOmb9GK1t2qmkAcBW4H4zK4g7nkwkdQW6mtlSSQcArwJnm9mbMYe2G0kCOpvZVkntgWeBK8zshZhDS0vSlcBxwIFmNjTueDKRVAUcZ2Yfxh1LfSTNABab2T2S8oD9zOzTmMPKSFJb4D3g+2aWjVppWSXpcKL/p75jZp9LegD4t5lNz+ZxWl0PwcyeAT6OO46GMLMNZrY0vN8CvAUcHm9UqVlka1hsH145e7YhqRtwJnBP3LHsSyR1AQYA9wKY2fZcTwbBQGB1LiaDJO2ATpLaAfvRBAX5Wl1CaKkk5QPHAC/GHEpaYQhmGbAReMLMcjZW4HbgaqAm5jgayoDHJb0q6VdxB5NBD+AD4L4wHHePpM5xB9UAI4GSuINIx8zeAyYD7wAbgE1m9ni2j+MJoQWQtD/wEDDOzDbHHU86ZlZtZoVAN+B4STk5JCdpKLDRzF6NO5ZGOMnMjgVOB34dhj5zUTvgWOAuMzsG2AaMjzekzMKw1lnAg3HHko6krwHDiBLuYUBnST/P9nE8IeS4MB7/EDDLzObGHU9DhCGCRcCQmENJ50TgrDAuXwqcKumf8YaUWThDxMw2AmXA8fFGlNY6YF1S73AOUYLIZacDS83s/bgDyWAQsMbMPjCzHcBc4IfZPognhBwWLtTeC7xlZrfFHU8mkg6WdFB43wkYDLwda1BpmNm1ZtbNzPKJhgqeMrOsn21li6TOYVIBYfjlR0BOzpIzs/8B70rqE1YNBHJuEkQd55HDw0XBO8AJkvYLfxcGEl1TzKpWlxAklQDPA30krZP0i7hjyuBE4AKiM9jE1Lgz4g4qja7AIknLgZeJriHk/HTOFuIQ4FlJrwEvAY+a2WMxx5TJ5cCs8G+hEJgYbzjphQQ7mOiMO2eFHtccYCmwguhv97RsH6fVTTt1zjmXWqvrITjnnEvNE4JzzjnAE4JzzrnAE4JzzjnAE4JzzrnAE4LLulCZs9nvUpZ0WaiwWZGYt5+07QhJj0paHqqGvtKQGCWNltQ7ablQ0og6+1i4m3xv4y+X9KWkryetOzm0P3kP25wQ7sRtyL47f2+h5ERReH+2pFy9Ec5lkScEty/5DXCBmR0TigEm+ysw38yONrOjiMoAbGxAm6OB3knLhcCIlHs2QihQlsrrRDfLJVxENPd8T9u/AWhQQkhmZr80s8Vh8Wxy985ol0WeEFyzkXRhODtfLqlM0rfC+jxJ0yStkvSspKmS5qRpo7+k50Mbz0vqH9bPBnoBMyXNSvGj3YjKGwNRKYhQBgJJB4Yz4pdCu38OhfouIiqPfUe4KfBc4CZgUFi+I0V8fSTNl/SyomdDXJS0zcIZ+8tEf6hTmQFcGPbfHzgJmJ/URltJkxXVxH89vG8btk0Pn2Mx8IqkO8OPPRfiPUhSsaQXQy+qQtLANN9zuaShkk4jqvMzPrRxYehp/Sxp359IynqhNRcDM/OXv7L6AqqAgjrrCojK9XYNy38AZof3lwOPERVG6wi8AMxJ0W4e0S38A8PyoLCcl+64ST97PrAFeBqYBPRP2nYPUc8CopOkEuCSsFwODE3ad3Td2Igqke4f4n8V6BvWHwCsTFo24JoM31s5MDT8tx9R72AKMAGYHPa5DFgYvos84EngsrBtOvAK0XMpasWWtPwNdt2Q2oeo7tBuv7fkzx3aHZu03xBgUdLyk8CwuP/d+WvvX95DcM3lFKIHemwIy3cT/UFPbJtpZl+Z2RekryvTB9huZk8CmNlCYHtYn5GZzQLygTuBzkRlNs4Lm88CrlJUunsp8D1qDxM1VG+iP+Sloa3FQIewLmFGA9qZAYwKr+l1tg0Cplv0nIHtwH3s+h4hSlbbMrTdC1gg6Q1gNnCopEMbEFOyBUBXSf0k9QttepmSfUC6cUzn9jlm9hHwAPCApHfZVdRMRE+i++9eHkLAhxaVAE9na4ZtCQ8CbxCV6F4h6ZxGxFBf+yXA78zsYUltgM+IemUNZmYmaSowJqy628yqG9OGy03eQ3DNZRFwRtLZ6CXAE+F9OXC+pHaKnhN7bpo2VgJ5kk4BkHQq0ZPZVtZ3cElnhrYTj0s8GlgTNj9CNEaeGIv/pqQeYdtmoEtSU3WX68b3maQLko7bV9KB9cWXzKInz10N/D7F5oXAKEntFZVGH8Wu7zGVLXXiPYhdn/tioh5MfVJ95hlEF5vPxZ86t8/whOCaykJF1WTXSVpHdEF3PPCEoiqY3wWuCPv+jegpUG8CTxGV9d1Ut8EwRHIOMDG08Ufgp2F9fU4GKsLPrQCqgevDtnFh+TVJK4iuZyQeVToNuD5cUB1ENF7eOVwwrnVR2cy+An4MjAwXp98gmt20J7N8ZpvZ0yk2TQOWAxXhtRz4e4ampgBPJS4qh8/6sKSlQE/gowaEMxMoTlxUDvFtIfqeHjezDxr2qVyu82qnLidIOsDMtkjqQHTG/qCZ+ZlnjgrTWpcDo8zs5bjjcdnhPQSXKxaGC7GvAZXsfjHV5QhJZwGriXoHngz2Id5DcM45B3gPwTnnXOAJwTnnHOAJwTnnXOAJwTnnHOAJwTnnXOAJwTnnHAD/B6FrVZ4/poMDAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ], "metadata": { "collapsed": false } }, { "cell_type": "markdown", "source": [ "Let's estimate the first stage relationship for the univariate regression model.\n", "To do this and get results similar to AJR we will need to work on a subset of the data where `baseco == 1`, essentially where there are complete data on observations." ], "metadata": {} }, { "cell_type": "code", "execution_count": 31, "source": [ "data = data.query('baseco ==1')\n", "\n", "first_state_reg = smf.ols('avexpr ~ logem4', \n", " data=data)" ], "outputs": [], "metadata": { "collapsed": false } }, { "cell_type": "code", "execution_count": 32, "source": [ "print(first_state_reg.fit(cov_type='HC1').summary())" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: avexpr R-squared: 0.270\n", "Model: OLS Adj. R-squared: 0.258\n", "Method: Least Squares F-statistic: 16.32\n", "Date: Mon, 13 Sep 2021 Prob (F-statistic): 0.000150\n", "Time: 12:26:10 Log-Likelihood: -104.83\n", "No. Observations: 64 AIC: 213.7\n", "Df Residuals: 62 BIC: 218.0\n", "Df Model: 1 \n", "Covariance Type: HC1 \n", "==============================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept 9.3414 0.704 13.264 0.000 7.961 10.722\n", "logem4 -0.6068 0.150 -4.040 0.000 -0.901 -0.312\n", "==============================================================================\n", "Omnibus: 0.035 Durbin-Watson: 2.003\n", "Prob(Omnibus): 0.983 Jarque-Bera (JB): 0.172\n", "Skew: 0.045 Prob(JB): 0.918\n", "Kurtosis: 2.763 Cond. No. 19.4\n", "==============================================================================\n", "\n", "Notes:\n", "[1] Standard Errors are heteroscedasticity robust (HC1)\n" ] } ], "metadata": { "collapsed": false } }, { "cell_type": "markdown", "source": [ "So the first stage seems OK with a t-state > 3 (approx F > 10). \n", "\n", "Let's do an IV2SLS estimation.\n", "`statsmodels` does not have the functionality to do IV regression, so we'll need a new package -- `linearmodels`.\n", "We'll only borrow the `IV2SLS` function from the `iv` module:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 33, "source": [ "from linearmodels.iv import IV2SLS" ], "outputs": [], "metadata": { "collapsed": false } }, { "cell_type": "code", "execution_count": 34, "source": [ "iv_regression = IV2SLS.from_formula('logpgp95 ~ 1 + [avexpr ~ logem4]', \n", " data=data).fit()" ], "outputs": [ { "output_type": "stream", "name": "stderr", "text": [ "/home/lachlan/anaconda3/lib/python3.7/site-packages/linearmodels/iv/data.py:25: FutureWarning: is_categorical is deprecated and will be removed in a future version. Use is_categorical_dtype instead\n", " if is_categorical(s):\n" ] } ], "metadata": { "collapsed": true } }, { "cell_type": "markdown", "source": [ "Somewhat differently from `statsmodels`, the `summary` which prints out the STATA-esque regression output is an **attribute** of the fitted model:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 35, "source": [ "iv_regression.summary" ], "outputs": [ { "output_type": "execute_result", "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
IV-2SLS Estimation Summary
Dep. Variable: logpgp95 R-squared: 0.1870
Estimator: IV-2SLS Adj. R-squared: 0.1739
No. Observations: 64 F-statistic: 28.754
Date: Mon, Sep 13 2021 P-value (F-stat) 0.0000
Time: 12:26:20 Distribution: chi2(1)
Cov. Estimator: robust
\n", "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Parameter Estimates
Parameter Std. Err. T-stat P-value Lower CI Upper CI
Intercept 1.9097 1.1740 1.6267 0.1038 -0.3912 4.2106
avexpr 0.9443 0.1761 5.3623 0.0000 0.5991 1.2894


Endogenous: avexpr
Instruments: logem4
Robust Covariance (Heteroskedastic)
Debiased: False" ], "text/plain": [ "\n", "\"\"\"\n", " IV-2SLS Estimation Summary \n", "==============================================================================\n", "Dep. Variable: logpgp95 R-squared: 0.1870\n", "Estimator: IV-2SLS Adj. R-squared: 0.1739\n", "No. Observations: 64 F-statistic: 28.754\n", "Date: Mon, Sep 13 2021 P-value (F-stat) 0.0000\n", "Time: 12:26:20 Distribution: chi2(1)\n", "Cov. Estimator: robust \n", " \n", " Parameter Estimates \n", "==============================================================================\n", " Parameter Std. Err. T-stat P-value Lower CI Upper CI\n", "------------------------------------------------------------------------------\n", "Intercept 1.9097 1.1740 1.6267 0.1038 -0.3912 4.2106\n", "avexpr 0.9443 0.1761 5.3623 0.0000 0.5991 1.2894\n", "==============================================================================\n", "\n", "Endogenous: avexpr\n", "Instruments: logem4\n", "Robust Covariance (Heteroskedastic)\n", "Debiased: False\n", "\"\"\"" ] }, "metadata": {}, "execution_count": 35 } ], "metadata": { "collapsed": false } }, { "cell_type": "markdown", "source": [ "Unfortunately there is no way to currently get the outputs from a regression estimated via `linearmodels` into a regression table!" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "### Final Challenge\n", "\n", "1. Estimate the IV model for the regression that includes 'asia', 'africa' and 'other' as exogenous explanatory variables\n", "2. Use a robust Hausman test to conduct a formal test for endogeneity. Was the IV strategy appropriate based on this test result?" ], "metadata": {} }, { "cell_type": "code", "execution_count": null, "source": [], "outputs": [], "metadata": { "collapsed": true } } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "name": "python3", "display_name": "Python 3.7.6 64-bit ('base': conda)" }, "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.7.6" }, "interpreter": { "hash": "bf5f67538785b35c4702aeb67651fc986043af768008a1d9e9d94d6c2205c8f1" } }, "nbformat": 4, "nbformat_minor": 2 }