{ "cells": [ { "cell_type": "markdown", "id": "aa6ee17c", "metadata": {}, "source": [ "# Solution sketch, lab 7" ] }, { "cell_type": "code", "execution_count": 7, "id": "7a9a2d17", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.stats as spt\n", "\n", "#new packages\n", "import statsmodels as sms\n", "import statsmodels.formula.api as smf\n", "import seaborn as sns\n", "\n", "\n" ] }, { "cell_type": "markdown", "id": "aa7efac9", "metadata": {}, "source": [ "### 1. Running a simple regression\n", " \n", "In Penn World tables data (\"http://jmaurit.github.io/anv_statistikk/data/pwt100.csv\") consider the relationship between average hours worked for employed persons (avh) and real per capita GDP in 2019 (rgdpe/pop). \n", "\n", "- Run a simple regression of average hours worked on rgdpe with statsmodels\n", "- plot the data with the regression line overlayed\n", "- What is the 95% confidence interval of the relationship? \n", "- What is the estimated $\\sigma$ of the model?\n", "- Briefly discuss whether you think the relationship is causal? " ] }, { "cell_type": "code", "execution_count": 8, "id": "b7de0571", "metadata": {}, "outputs": [], "source": [ "pwt = pd.read_csv(\"http://jmaurit.github.io/anv_statistikk/data/pwt100.csv\", sep=\";\", decimal=\",\")\n", "pwt[\"rgdpe_per_capita\"] = pwt[\"rgdpe\"]/pwt[\"pop\"]\n", "pwt2019 = pwt.loc[pwt.year==2019, :]" ] }, { "cell_type": "code", "execution_count": 9, "id": "e57fc223", "metadata": {}, "outputs": [ { "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", "
OLS Regression Results
Dep. Variable: rgdpe_per_capita R-squared: 0.254
Model: OLS Adj. R-squared: 0.242
Method: Least Squares F-statistic: 21.78
Date: Mon, 10 Oct 2022 Prob (F-statistic): 1.61e-05
Time: 11:22:41 Log-Likelihood: -745.28
No. Observations: 66 AIC: 1495.
Df Residuals: 64 BIC: 1499.
Df Model: 1
Covariance Type: nonrobust
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [0.025 0.975]
Intercept 1.142e+05 1.7e+04 6.732 0.000 8.03e+04 1.48e+05
avh -42.3617 9.078 -4.666 0.000 -60.497 -24.227
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 34.079 Durbin-Watson: 2.164
Prob(Omnibus): 0.000 Jarque-Bera (JB): 79.123
Skew: 1.675 Prob(JB): 6.59e-18
Kurtosis: 7.189 Cond. No. 1.31e+04


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.31e+04. This might indicate that there are
strong multicollinearity or other numerical problems." ], "text/plain": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: rgdpe_per_capita R-squared: 0.254\n", "Model: OLS Adj. R-squared: 0.242\n", "Method: Least Squares F-statistic: 21.78\n", "Date: Mon, 10 Oct 2022 Prob (F-statistic): 1.61e-05\n", "Time: 11:22:41 Log-Likelihood: -745.28\n", "No. Observations: 66 AIC: 1495.\n", "Df Residuals: 64 BIC: 1499.\n", "Df Model: 1 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept 1.142e+05 1.7e+04 6.732 0.000 8.03e+04 1.48e+05\n", "avh -42.3617 9.078 -4.666 0.000 -60.497 -24.227\n", "==============================================================================\n", "Omnibus: 34.079 Durbin-Watson: 2.164\n", "Prob(Omnibus): 0.000 Jarque-Bera (JB): 79.123\n", "Skew: 1.675 Prob(JB): 6.59e-18\n", "Kurtosis: 7.189 Cond. No. 1.31e+04\n", "==============================================================================\n", "\n", "Notes:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "[2] The condition number is large, 1.31e+04. This might indicate that there are\n", "strong multicollinearity or other numerical problems.\n", "\"\"\"" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ex_mod1 = smf.ols(formula='rgdpe_per_capita~avh', data=pwt2019).fit()\n", "ex_mod1.summary()" ] }, { "cell_type": "markdown", "id": "6f0d5112", "metadata": {}, "source": [ "We have a negative and statistically significant relationship\n", "\n", "I'll use seaborn lmplot to give a quick visualisation of the regression line: " ] }, { "cell_type": "code", "execution_count": 10, "id": "b13ff880", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABE/0lEQVR4nO3deXzU53Xo/8+ZTTMS2thBwsbY2NjEG2DHcR1C7CR2Gsd2DEnwTW+cxjem+eWmbnN7G6ftdVO3/TVOe9M6XVKTps3uJeDEZHEcJ4SQxRgDxgsGBww2khAIhNA2+8y5f3y/I0YgJI00o++MOO/XSxnp0SxHxDp65nyf5zyiqhhjjJl4Pq8DMMaYs5UlYGOM8YglYGOM8YglYGOM8YglYGOM8UjA6wDKyU033aQ//vGPvQ7DGDP5yFCDNgPOc+zYMa9DMMacRSwBG2OMRywBG2OMRywBG2OMRywBG2OMRywBG2OMRywBG2OMRywBG2OMRywBG2OMRywBG2OMR2wrcgXZtKeDhzbvp6UryrzGatYsX8CKRTO9DssYM0Y2A64Qm/Z0cN+GXXT0xmmIBOnojXPfhl1s2tPhdWjGmDGyBFwhHtq8n6BfqA4FEHFug37hoc37vQ7NGDNGloArREtXlEjQP2gsEvTT2hX1KCJjzHhZAq4Q8xqriaUyg8ZiqQzNjdUeRWSMGS9LwBVizfIFpDJKNJlG1blNZZQ1yxd4HZoxZowsAVeIFYtmcv8ti5lZG6Y7lmJmbZj7b1lsqyCMqWC2DK2CrFg00xKuMZOIzYCNMcYjloCNMcYjloCNMcYjloCNMcYjloCNMcYjloCNMcYjloCNMcYjloCNMcYjloCNMcYjloCNMcYjloCNMcYjloCNMcYjloCNMcYjloCNMcYjloCNMcYjloCNMcYjloCNMcYjloCNMcYjdiSRKWub9nTw0Ob9tHRFmddYzZrlC+xYJjNp2AzYlK1Nezq4b8MuOnrjNESCdPTGuW/DLjbt6fA6NGOKwhKwKVsPbd5P0C9UhwKIOLdBv/DQ5v1eh2ZMUVgCNmWrpStKJOgfNBYJ+mntinoUkTHFZQnYlK15jdXEUplBY7FUhubGao8iMqa4SpqAReQ/RaRDRF7OG5sqIk+LyF73tjHve58RkX0i8qqI3Jg3vlREXnK/90UREXe8SkQedcefFZH5eY+5032NvSJyZyl/TlMaa5YvIJVRosk0qs5tKqOsWb7A69CMKYpSz4C/Ctx0yti9wM9UdSHwM/drROQSYDWw2H3Mv4lI7v3nl4C7gYXuR+457wK6VPUC4B+BB9znmgr8JfBm4GrgL/MTvakMKxbN5P5bFjOzNkx3LMXM2jD337LYVkGYSaOky9BUdXP+rNR1K7DC/fxrwCbg0+74I6qaAA6IyD7gahF5HahT1WcAROTrwG3Ak+5jPus+1zrgX9zZ8Y3A06p63H3M0zhJ++Fi/4ymtFYsmmkJ10xaXtSAZ6lqO4B7m/vtagJa8u7X6o41uZ+fOj7oMaqaBrqBacM812lE5G4R2SYi244ePTqOH8sYYwpTThfhZIgxHWZ8rI8ZPKi6VlWXqeqyGTNmjCpQY4wpBi8S8BERmQPg3uZW1bcC8/Lu1wwccsebhxgf9BgRCQD1wPFhnssYY8qGFwl4A5BblXAn8ETe+Gp3ZcN5OBfbtrplil4Rucat7374lMfknmsVsFFVFXgKeJeINLoX397ljhljTNko6UU4EXkY54LbdBFpxVmZ8DngMRG5CzgIvB9AVXeJyGPAK0Aa+ISq5haBfhxnRUUE5+Lbk+74V4BvuBfsjuOsokBVj4vIXwPPufe7P3dBzhhjyoU4E0YDsGzZMt22bZvXYRhjJp+hrkuV1UU4Y4w5q1gCNsYYj1gCNsYYj1gCNsYYj1gCNsYYj1gCNsYYj1gCNsYYj1gCNsYYj1gCNsYYj1gCNsYYj1gCNsYYj1gCNsYYj1gCNsYYj1gCNsYYj1gCNsYYj1gCNsYYj1gCNsYYj1gCNsYYj1gCNsYYj5T0UE5jjBmvTXs6eGjzflq6osxrrGbN8gWsWDTT67CKwmbAxpiytWlPB/dt2EVHb5yGSJCO3jj3bdjFpj0dXodWFJaAjTFl66HN+wn6hepQABHnNugXHtq83+vQisISsDGmbLV0RYkE/YPGIkE/rV1RjyIqLkvAxpiyNa+xmlgqM2gslsrQ3FjtUUTFZQnYGFO21ixfQCqjRJNpVJ3bVEZZs3yB16EVhSVgY0zZWrFoJvffspiZtWG6Yylm1oa5/5bFk2YVhC1DM8aUtRWLZk6ahHsqmwEbY4xHLAEbY4xHLAEbY4xHLAEbY4xHLAEbY4xHLAEbY4xHLAEbY4xHLAEbY4xHLAEbY4xHLAEbY4xHLAEbY4xHLAEbY4xHPEvAIvLHIrJLRF4WkYdFJCwiU0XkaRHZ69425t3/MyKyT0ReFZEb88aXishL7ve+KCLijleJyKPu+LMiMt+DH9MYY87IkwQsIk3AHwLLVPVNgB9YDdwL/ExVFwI/c79GRC5xv78YuAn4NxHJtcn/EnA3sND9uMkdvwvoUtULgH8EHpiAH80YY0bNyxJEAIiISACoBg4BtwJfc7//NeA29/NbgUdUNaGqB4B9wNUiMgeoU9VnVFWBr5/ymNxzrQNuyM2OjTGmHHiSgFW1DfgH4CDQDnSr6k+AWara7t6nHcg1AW0CWvKeotUda3I/P3V80GNUNQ10A9NK8fMYY8xYeFWCaMSZoZ4HzAVqROT3hnvIEGM6zPhwjzk1lrtFZJuIbDt69OjwgRtjTBF5VYJ4B3BAVY+qagp4HLgWOOKWFXBvO9z7twLz8h7fjFOyaHU/P3V80GPcMkc9cPzUQFR1raouU9VlM2bMKNKPZ4wxI/MqAR8ErhGRarcuewOwG9gA3One507gCffzDcBqd2XDeTgX27a6ZYpeEbnGfZ4Pn/KY3HOtAja6dWJjjCkLnpwJp6rPisg6YAeQBp4H1gJTgMdE5C6cJP1+9/67ROQx4BX3/p9Q1dxZ1R8HvgpEgCfdD4CvAN8QkX04M9/VE/CjGWPMqIlNCk9atmyZbtu2zeswjDGTz5ArsGwnnDHGeMQSsDHGeMQSsDHGeMQSsDHGeMQSsDHGeMQSsDHGeMQSsDHGeMQSsDHGeMQSsDHGeMQSsDHGeMQSsDHGeMSTZjxmYmza08FDm/fT0hVlXmM1a5YvYMWimSM/0BgzIWwGPElt2tPBfRt20dEbpyESpKM3zn0bdrFpT8fIDzbGTAhLwJPUQ5v3E/QL1aEAIs5t0C88tHm/16EZY1yWgCeplq4okaB/0Fgk6Ke1K+pRRMaYUxVUAxaR9+AcDR/Ojanq/cUOyozfvMZqOnrjVIdO/l8cS2Vobqz2MCpjTL5Rz4BF5N+BDwKfxGku/H7g3BLFZcZpzfIFpDJKNJlG1blNZZQ1yxd4HZoxxlVICeJaVf0w0KWqfwW8hcEHZZoysmLRTO6/ZTEza8N0x1LMrA1z/y2LbRWEMWWkkBJEzL2NishcoBPnWHlTplYsmmkJ15gyVkgC/oGINAB/j3OYpgL/UYqgjDHmbFBIAv68qiaA9SLyA5wLcfHShGWMMZNfITXgZ3KfqGpCVbvzx4wxxhRmxBmwiMwGmoCIiFzJyeOV6wBb02SMMWM0mhLEjcBHgGbgC3njvcCflSAmY4w5K4yYgFX1a8DXRGSlqq6fgJiMMeasMJoSxO+p6jeB+SLyqVO/r6pfGOJhxhhjRjCaEkSNezullIEYY8zZZjQliIfc278qfTiVpy+RZkqVtVU2xhSukF4QC0Tk+yJyVEQ6ROQJETnrGwt09MQ53B0nncl6HYoxpsIUsg7428BjwBxgLvAd4OFSBFVposk0rV0xeuIpr0MxxlSQQhKwqOo3VDXtfnwTZzuyAbKqHOtN0N4dI5m22bAxZmSFFC9/LiL3Ao/gJN4PAj8UkakAqnq8BPFVnFgyQ1sqRm04QGN1CL9PRn6QGZKdaWcmO1Ed3SRWRA4M821V1YqvBy9btky3bdtW0GP2H+074/f8PqGxJkRdODje0M46uTPtgn4hEvQTS2VIZdRaappKNeRMbNQzYFW11pMFymSdskRfPM20KSGqAv6RH2SAwWfaAVSHAkSTaR7avN8SsJk0Cj2S6E3AJQw+kujrxQ6qUmza08E//XQv7T0x5tRFWH3VPK5eMPW0+8VTGdq6YtRFglaWGKWWrigNkcHvHOxMOzPZFLIM7S+Bf3Y/3g58HrilRHGVvdxb5M7+BHXhAJ39CR7cuJet+89cCu+JpWjtitpqiVGY11hNLJUZNGZn2pnJppBVEKuAG4DDqvr7wOVAVUmiqgC5t8iRoB/BuQ34hEeeaxn2cbmyRNuJGPFTEow5yc60M2eDQhJwTFWzQFpE6oAO4Kz9bRjq2Pdw0MfhntgZHjFYIpXh0IkYHT22iWModqadORsUUgPe5h5J9GVgO9AHbC1FUJUgd+y7X07Wc+OpLLPrIgU9T18iTTSZoaE6SH0kiIjVh3PsTDsz2Y16Bqyq/5+qnlDVfwfeCdzpliLGREQaRGSdiOwRkd0i8hYRmSoiT4vIXve2Me/+nxGRfSLyqojcmDe+VERecr/3RXEzmIhUicij7vizIjJ/rLEOJfcWOZbKoDi36ayy+qrCD4rOqnK8P0lrV4xY0soSxpwtCrkI9z4RqQdQ1deBgyJy2zhe+0Hgx6q6CKeevBu4F/iZqi4EfuZ+jYhcAqwGFgM3Af8mIrn3/18C7gYWuh83ueN3AV2qegHwj8AD44j1NLm3yNNqquiNp5lWU8U91y8cchXEaKUyWdq7YxztTZDN2iZDYya7QjZi7FTVK04Ze15Vryz4RZ0a8gvAAs0LQEReBVaoaruIzAE2qepFIvIZAFX9O/d+TwGfBV4Hfu4mcUTkDvfxa3L3UdVnRCQAHAZm6DA/cLE3YoxH0O9j+pQqIiFbO2xKw3YaTqgha4uFXIQb6r5j7cO4ADgK/JeIPC8i/yEiNcAsVW0HcG9z/zU0AfnLC1rdsSb381PHBz1GVdNANzBtjPFOuNxs2C7SmVLILaPs6I3TEAnS0Rvnvg272LSnw+vQziqFXoT7AvCvOL0gPolzMW6sr7sE+KSqPisiD+KWG85gqL8eOsz4cI8Z/MQid+OUMDjnnHOGi7notu4/ziPPtQy7kSN3ka6xOkR9tW1pNsVhOw3LQyEz4E8CSeBRnLaUMeATY3zdVqBVVZ91v16Hk5CPuKUH3NuOvPvnX91qBg65481DjA96jFuCqAdO2yWhqmtVdZmqLpsxY8YYf5zCbd1/nAc37h3VRo6sKp39CVq7orZ22BTFUMsobafhxCtkFUS/qt6bS1aq+meq2p/7voj8cwHPdRhoEZGL3KEbgFeADcCd7tidwBPu5xuA1e7KhvNwLrZtdcsUvSJyjbv64cOnPCb3XKuAjcPVfyfaI8+1EPAVtpEjmc46a4d742TsIp0ZB9tpWB6KeZbO7xR4/08C3xKRELAf+H2cPwiPichdwEHg/QCquktEHsNJ0mngE6qa+6/n48BXgQjwpPsB8BXgGyKyD2fmu3qMP1dJtPfEqAsP/ucf7UaOvniaWDJDQ3WI+oiVJUzh1ixfwH0bdhFNpgd1m7OdhhPLs8PMVHUnsGyIb91whvv/LfC3Q4xvA940xHgcN4GXozl1ETr7E4PeBhaykSOTVTr7EvTGU0yfUkU4aKslzOitWDST+3Fqwa1dUZptFYQn7DRJj6y+ah4PbtxLLJUhHPQRT2XHtJEjV5aYEg4wtTpEwF9IWd+czWynofeKmYDPqj20uTWU+4/1DduK8lT5Kx9qQgFQpTeeZnYBzzGUvniaaMJZLVEXCRS8pdnWhBoz8UaVgN1dZ59T1f89zN0eLE5I5S//tIb8FQz3MPxOuNzKh4DPeVxu1vtHN1w4rh10ObnVEj1uWWK0mzhyP08qk6E7mqK9O8aOg118YsX5/OE7Lhx3XMaYoY3q/ap7wWupDDOtUtWvFiuocpdbQ5lIZ0mlddStKMey8mEscps4jvWNbkvzQ5v3k8pk6OxLkVFnF15WlX/d9JotzDemhAopQTwPPCEi3wEGlp+p6uNFj6rMtXRFqQn52d8VQ4GakJ+G6iDt3cOvoRzPyoex6ImliCUzTJsSGlhwP5SWrijd0RQi4HP/xvoF0lm1hfnGlFAhV2ymAp3A9cB73Y+bSxFUuZvXWE1fIk3A7ySr/mSGthNxumNpfvRS+xmPpZ9TFyGeGvy9sbSwLEQqk+Vwd5wjw2xpntdYTSKTJf/9jSpUBXy2MN+YEipkI8bvD/Hx0VIGV67WLF+A3+djXmOEOfVVVAWcf8Z4Oss//OS3rF67hf/69QGO9ycHPW71VfNIZ4vTwrJQ/Yk0rV0xTkSTnLofZc3yBQR8PjJZRVXJZpUsSm04YAvzjSmhQtpRXigiPxORl92vLxORvyhdaOUr14pyVl0EVbh4dh2ffPsFvPOSWQR8wolYim9sOcgdX97CAz/ew2sdTse0qxdM5Z7rFxa1hWUh8vsO9yfSg36eT6w4H58I6awS8AvTakKEAn5bmG9MCRXSjvIXwP8GHsq1oBSRl1X1tE0QlaoY7SiP9SV4Yuchvv/CIXriJ5PcFfMaWLmkibecP22gzuq1cNDP1JrQwCaO3FI0W5hvTNEN+UtfSAJ+TlWvyu8BPFSP4EpWzH7AiVSGp3d3sH5HK290nqyjNjVEuH1JEzctnl02vX6nhANMq6nC7yuPPwzGTEJD/nIVsgrimIicj9vSUURWAe1FCGxSqgr6ufmyObzn0tlse6OLddtbee71LtpOxPjnjfv4z18f4D2XzuF9VzYxqy7saaz5mzis5aUxE6eQGfACYC1wLdAFHAA+pKpvlC68iVXqEzFe7+zn8R1t/OSVIwMrJXwCyxfOYNXSZi6ZW1fQa5dCVdDPtLyyhDGmKMZXghh4gHNyhU9Ve4sRVTkpNAFv2tPBP/1077AN1YfSHU3x/RcP8cTOQ3TmrZS4ZE4tK5c0s/zCGZ6XA+ojQabWhOyUZmOKY9w14GnAXwLX4ZQhfgXcr6qdxYrQa4Uk4Nz2XVUd1EynkFUNqUyWn796lHXbW9nXcXImPbO2ituubOI9l86mNlyaksBoTuMI+n3MqLVOa8YUwbgT8NPAZuCb7tCHcA7AfEdRwisDhSTgO9ZuoaM3jj9vhhhLZZhWU8UXPnh5Qa+rqrzY1s267a38Zl/nwLlJ4aCPGxfPZuWSpqKux83vSTGaPx5TqgJMrSlup7VKa/5TafGasjPuBLxdVZeeMrZNVYfq6VuRCknA1z2wkYZIcNCuN8XpbPbtj10z5hgOnYjx+PNtPPnS4YETCwS4ZsE0Vi5t4sp5DeMuC3zq0RdO60U80h8Pn0jRLtLlNzPKbwZ+/y2LyzKpVVq8piyN+1Tkn4vIahHxuR8fAH5YnNgqz1BHuhRjW/Hchgj/8+0X8Oiaa/j42xYwq64KBZ7Z38mffOdF7v7Gdn788uEzbncejfaeGOHg4P/rR+pJUcxz6fIPhBRxboN+4aHN+8f1vKVSafGaylFIAl4DfBtIuB+PAJ8SkV4R6SlFcOVszfIFpDKl21Y8pSrA+5fN45t3vZnPvvcS3uSukHjtaD+ff+pV7vjyFr72m9fpiiZHeKbTjacnRf65dGfqLbFpTwd3rN3CdQ9s5I61W07rqFZpB0JWWrymchTSC6JWVX2qGnQ/fO5YrarWicjiUgZabnLbkUu9rdjvE5ZfOIMv3nEl//ahK7lh0Uz8PqErmuJrz7zB6rVb+PunXi1oOVwxelL0xYfuLZF7u97RG6chEqSjN859G3YNSsKVdiBkucc70h88U74KXoZ2xicS2aGqS4ryZB4p9TrgYjnam+B7O9v4wYvt9OZtd156TgMrlzZz9XlTR9zunFsFcbgnNu7TOIJ+30DLy9zFyfz2l9Fkmpm1YR6+26mNV1pNtZzjLefYzCDFWQd8xmfP26JcqSolAefEUhmefuUI67e30tJ1sn47rzHC7UuaedfiWae9dS6lSMjPqi/9hsbqweuHVZXuWIpffvr6gbFK6ztRrvGO5g+eKQvj3oo8kuJkcjNqkaCfWy6fy82XzWHrgeOs397K9oMnaOmK8eDP9vKfvz7AzZfN4bYrmphRW1XyeGLJDDOmhOmKJamtOnku3VBv1yvtQMhyjbelK0pDZPDKFKtPVw47FXkS8IlwzYJpXLNgGgeO9bN+eytP7z5CbzzNw1tbeGxbKysunMHKpU0smn1yu/NoNmMUKnfac4+mqAkFSGaypDJqbS1LZF5j9Wkz4HKqT5vhFbMEsUVVK/o9T6WVIIbTFU3y/Rec7c5d0dTA+OK5daxa2kzI5+NfNu0b9WaMQuTXl+fUR7j7rQt492VzxvsjmSFYDbhijHsjhuDsflugqveLyDnAbFXdWrwYvTWZEnBOMp1l454O1u1oZf/RgaP8CPl9VId8g9pQjnUn32iEAs5rlUsLzsmkXOvTZpBxJ+AvAVngelW9WEQagZ+o6lXFi9FbkzEB56gqO1tOsG57G1v2n9zuLAL14SANkSDBgBS8k6/QMkYk5DSBrwpYIjZnlXFfhHuzqi4RkecBVLVLREJFCc2UnIhw5TmNXHlOI61dUXc7chJVOBFLcSKWIhL00dRQjaqOartzfk+JunCAzv4ED27cyz2cuYwRS2ZoS8aYEg4wtbq4/SWMqTSF/NefEhE/Jxuyz8CZEZsK09xYzZ+86yJm1lVRHwkQGChBZNl3tI8139zBT145QuoMO91yHnmuhYDPqT0Kzm3AJzzyXMuIMfTF07R0xejqP/2QUGPOFoXMgL8IfBeYJSJ/C6wCzspDOb1SzFULVy+Yyh9zofN83VHCwQCqysGuGPs6+vjck3tYu3k/t14xl1sumztkE572nhh14cH/CY3UUyKfqtIVTdIbT1MfCVIbDuCzY5HMWaSgVRAisgi4wf1yo6ruLklUHinnGnChLSTH6pVDPazf0covfnuUrPufRijg450Xz2Ll0ibmT6sZuO9YuqoNxydCbThAXSRI0EoTZnIZdzc0gGrA7z5ufG2/JomJOjFiPG/3C3HJ3Dr+z82X8K3/8WY+uKyZmio/yXSWH77Uzke/uo1Pr3+RrQeOo6pF6SmRL+vumGvtitHREx9XxzdjKkEhqyDuA94PrMfJ5rcB31HVvylZdBNsLDPgrJuAoskM8VRmxLrpWN3x5S3UhQNI3h/SYvQfHkksmeGpXYd5/Pk2WvO2O587tZqVS5tojIRYv6OtKD0lhlIbDtJYHbSLdabSjXsZ2m7gSlWNu19HgB2qenHRQvTYWBLwqZLpLLFUhkQqQyyVIZMtzgWmYr/dL1RWlS37O1m3vY2dLScGxuvCAd57+VxuvWIu06eUZruzuKWJeitNmMo17mVorwNhIO5+XQW8Nr6YJp9QwEco4AN3f34ynSWedmbHybSzLXcsV/1zW3xjqcygGnCx+g+PxCfCtedP59rzp/NaRx/rdrSycU8HPfE033r2II8+18KKi5zTnS+cVVvU11ZVemIpeuNpakJ+6iLBkp9TZ0cQmYlQyAz4e8BVwNM4S9HeiXMwZweAqv5haUKcOMWYAY9EVUmks8RTGRLpLIlUlnR2dGWLYraQLIbj/Uk27DzEhhcOcSJ2crvzZc31rFzSzLXnTyvZ6c6hgI+G6hBTqorfzsS295oSGHcJ4s7hvq+qXxtDUGVlIhLwUDJZJZHOkEg5s+VEKku2gtbGJtNZfrr7COt3tHHg2MntznPqw9y+pIl3v2n2oGYxY3GmJXihgI+6SJApoeItYbMWj6YEStsPeDLwKgEPJZnOOknZnS2PtXQxkVSVHQdP8B+/PMCrR3oHxqsCPt57+Rxuv7KZg53Rgtcyj2YJnk+E6io/deHxlydyB66O1NPYmAKMrQYsIi8xTK9fVb1sHEGZM8jVknPV1PzSRW7FRbkRETIZpSeeYnZdFbFUhp5YmkQ6y7rtbazf3kZVwEddJDDqrcsweAkeMFAWeOS5loHHZVXpi6fpi6epCvqpjwSpCfnHtEzQWjyaiTKa94U3u7efcG+/4d5+CLCuzxNERAgH/YSDfhqqnbKFs/wtTSxZvNUW43UyWQaoCweZXqMc608QTTqz+Hg6S7w3SXcs7Swvc9cyD5eA83fc9SfTHO9PkkxnOdwTZ+v+46c9NpHK0JHKEPDlkn2woPLEmuULuG/DLqLJ9KAa8Fh6GtvFPDOcEdf0qOobqvoG8Duq+qeq+pL7cS9w43heXET8IvK8iPzA/XqqiDwtInvd28a8+35GRPaJyKsicmPe+FIRecn93hfdtpmISJWIPOqOPysi88cTa7nx+4QpVQFm1oY5d1oNcxsiNFaHnBUYHjr1yHu/T5hZW8XU6iB1kQBVbnyJdJbDPQkO9yTYd7SXnryLeKfKneLcn0zT0ZMgnVV84nRye3DjXrbuPz7k49LZLMf7kxw8HuV4f3LUf6RyB67OrA3THUsxszY8pgtwozmg1JzdCrkyUiMi16nqrwBE5FqgZoTHjOQeYDeQO6bhXuBnqvo5EbnX/frTInIJsBpYDMwFfioiF6pqBvgScDewBfgRcBPwJHAX0KWqF4jIauAB4IPjjLds5WbHjTUh0pks0VSGaMJZi3ym2nEpTsSYUxc5bb1yPJVlTr3z9v1YXxxB6Iol6Us4O+f6Ehk+uHYL71o8i5VXNnPOtMFv9XNL8Dr7EoCCCgrMcHsZjzSDzqpyIpqkO5aiJuSnNhwcsS9xMY4gemjzfoJ+GShlVIcCRJNpHtq832bBBihsK/JdwL+KyOsi8jrwb8BHx/rCItIMvAf4j7zhW4Hcaoqv4ey2y40/oqoJVT0A7AOuFpE5QJ2qPqNOlvn6KY/JPdc64Ibc7HiyC/h91IWDzK4PM39aNbPrw9RFggR8J//vzl3Y6uxPDKrHnmk2OVq57cnH+xO0dEV57Wgfh3viXDmvntVXzSOjgDgrJObUVVET8lMV8JFIZ/n+C+185KvPce/jL7Ht9eMDfziuXjCVe65fiOJcjAj4hZm1YaZUBQpu/tOXSNPeHaPleJTuaKqkpZuWruhph6LaeW0m36hnwKq6HbhcROpwVk90j/O1/wn4UyB/1f4sVW13X69dRHLThCacGW5OqzuWcj8/dTz3mBb3udIi0g1MA46NM+6KIuLMwKpDAZgCiXSGWDLDY9tHvrA1FlcvmMpNh2fxza0HyWSVkN/ZxfbjV45w0ew67rl+4aC1zPfcMI83Ndfx45cPs35HG+3dcbYeOM7WA8c5b3oNK5c08Y6LZ3H1gqksnlM/5Ox6dl3hbUlSmSyd/QmOR5NUh/xMqQpQPcaLdmdiF/PMSEadgEXkU6d8DdANbFfVnYW8qIjcDHSo6nYRWTGahwwxpsOMD/eYU2O5G6eEwTnnnDOKUCpbVcBPVcBPR2+C+nAARciqklUtaDY5nOdbupldFz5t2/Qjz7XwhQ9ePmSCv31JM7de0cQzr3WybkcrL7Z2c+BYP//wk9/yH788wC2Xz+U9l87mP3/z+qDdgH2JNEGfcMeXt4ypjKKq9CfS9CfSA3X12nCwKLX0Yl7MM5NTITXgZe7H992v3wM8B/yBiHxHVT9fwHP9DnCLiPwuzvbmOhH5JnBEROa4s985uLvscGa2+Xtum4FD7njzEOP5j2kVkQBQD5z2/lpV1wJrwVkHXMDPUNHyZ2d+xElEyTTNjdWEAr5xrTsea59gv0+4buF0rls4nd8e6WX9jjZ+vqeDE7EUX9/yBkG/cGlTPf3xDN3x5MDMMpXVgpa1nUkm66z17Y6lqA45vSfGc4bdikUzuR/svDZzRoXshHsKWKmqfe7XU3Bqq+/DmQVfMqYAnBnwn6jqzSLy90Bn3kW4qar6pyKyGPg2cDXORbifAQtVNSMizwGfBJ7FuQj3z6r6IxH5BHCpqv6BexHudlX9wHCxlNNGjFIbzXbbbFZJZ5VUJks6o0RTaeKp7IiJuZiNg471JXhi5yG+/8IheuLpgfEr5jXQHU2RymROe4tfzAZFuZ12tVWBCWs9aialcTfjOQdI5n2dAs5V1ZiIJMYTWZ7PAY+JyF3AQZz2l6jqLhF5DHgFSAOfcFdAAHwc+CpOf+In3Q+ArwDfEJF9ODPf1UWKcVIYzezM5xNCPhl4O15PkEzWmSn3xtMkzrAZpJiNg6ZPqeKu687jQ28+x9nuvL2NN45HBzqyBf1CQyRLvbvWdzxllDOtCjnWm6CrP1nU8oQxMMoZsLt64P/grDB4wh1+L7AB+L/AWlX9UIlinDBn0wy4GJLpLLFkhkQ64ybZk02FStU4SFXZ9kYXX/nlAX7bcfI0EgEaqoOEgz5m1UYKngEXcuJIVdBPbThQ1P4TZtIb+wxYVVVEbgU+BlznPtkfqGouW1V88jWFG2i9idN6M5VxtkrHUhnecsG0knRqExE0C72JNI3VAU5E0wPL07qizmaOGVPC7DrUzeK59aN+3tFsd85JuP2eOyVJJOgnEvJTE/Jb03hTsEJKEFsAn6o+WKpgTGUL+n0E/T5qw05Cjqcy9MRT9CfOvBlkLHLJcsaUMNWhNJ19CRJp5/kVeLGtm08+vJOL59Syakkzyy+cMWJbzLFcOFRVosk00WSaTpyNFrXhAJGg32bGZlQKScBvB9aIyBtAP84sWK0ZjzmT3O68TI2zASLmNhEab6vN/GRZEwpQMzWA4jRt/+h1C1i3vZW9HX3sbu/lr3+4m5mb93PblU3cfOkcpoSH/k/+TDv4ClljnEvG4Pwxirjri0vdPN5UrkIS8LtLFoWZ1Pw+oT4SpD4SdJe7ZeiNp4glx9bRbbjtzu+8ZBbvuHgmL7Z1s257K7/Z10lHb4K1m/fz9Wde58bFs1m5pOm0zRDFPnEklcmSimXpiaUIBXxUh5xde1UBf8ma1JvKY/2A89hFuImVymQHjhoqZFZcyAWzthMxvrujjSdfPkzMXbUhwDULprFqaRNXzGsYWF42ESeOiDh15tpw8XfembJmDdlHYgnYG1l3aVuuTDEahSbLvkSaJ19q5/Hn2zjSc3LVZMAnzKkPs+atC7h24fRx/yyF8PucbeJTqgLj2vBhKoIl4JFYAvZeOpOlP5GhL3nmdcbjkckqX/316zy2vYVU5uR/+z6B6y+aycfffj6N1aGiv+5IAj4fNVV+aqxmPFmNeyOGMSUX8Puor/ZRXx0kmc7Sn3BmxqnM6A4uHYnfJ+w61MPsujCCs3StN5Emq/DTPR38Yu9R3nHxLFYuaWLBjClFec3RSGezdMeydMdS+N3lcFVB/0Dd2ExOloCNJ0ZzUoSzzjhEY01o4Cimvnh61KdIn0luFYUgzKn3Mz2T5UQsxYloilRGefLlwzz58mGWntPAyqXNXH3eVHwlrtUO15vZ7zt5GkpVwEdVwGe140nCShB5rAQxMcZz7Huup293LEUyPbZEfKZeFQ2REDdcMpP121tp6Tq5/ndeY4TblzTzrsWzTuvvWwyFXFQE50Je7t8uEvITDti64wpgNeCRWAIurjPNcot17HssmaE7lhpYeztaIyW8rCpbDxxn/fZWth88MfC42nCAmy+bw21XNDGjtqqg1xzOeJsXiQhVAd/JhGw15HJkNWAzcfJnufnnod2Pc1JEQyQ46P5jOSkiEnISTjLt1E77EulR7bi7esFU7mHhGVdR+ES4ZsE0rlkwjQPH+lm/vZWndx+hN57m4a0tPLatlbddOINVS5tYNLtuhFcb2Vjbd+aoKvGUs8mlK+rEn0vEkaDfmgeVMUvApiSGOw+t2CdFhAI+ZtRWMbUmRF88TU88NeqLdiOl6/Om1/AnN17EXW89j++/cIgndh6iK5pi454ONu7p4E1z61i1tJnfuWD6mDdYFGMXXr5sXpN5cFZYhEM+wkE/Ib+PkN9nJYsyYSWIPFaCKJ7rHthIQyQ46GKRqtPw/K9vfdOYa8CjNVwfikJrrvmS6Sw/f7WDddtbee1o/8D47Low77tyLu++dA5Tqgqb14wnnrEK+HwEA0LI7xu4wGc79ErKasAjsQRcPCPVeXP14VKfFJHJ6mmz4mI0jFdVdracYN32Nrbs7xyYSUeCft79ptm8b0kTTQ2jn8FOxC68keS2TFe7B6XaSouisgQ8EkvAxTOelQ6lEk2m6YmlufVffzWwDC1HUXrjab79sdFfBMxp7YqyfkcbT718mLi7MkOAay+YxqqlzVzWVF9xyUxEBtYg55a+WbvNcbEEPBJLwMU1UbPcQn3woWc43BOnKuAbKAIX4yij3niKH750mO8930ZH78ntzgtnTmHl0mbeftEMghWcxPw+OZmQrbFQoSwBj8QS8NkhNzsP+CDk9xF1Z+fFqrlmssov9x5l3fZWXmnvHRifVhPi1ivm8t7L5lJfHRzmGSpHVdBPtbv8zcoWw7IEPBJLwGeP/Nl5U0OED7/lXC5tbijaluecVw71sG57K5v3HiXr/qqFAj7eefEsVi5tYv60miEfN9zOuHIWCjirLIJ+57QUp0m/WGK2BDwyS8CVYTTbmMeq391lFx9HI6Chkue506v57vNt/PCldvoTJ5/7qvmNrFzSzFXzGwe1xZzoVRGlJCIEfM5mkVxiziXns4gl4JFYAi5/E3VxL5bM0BVNFpyIR0qesWSGp3YdZv2ONtpOnNxoce7UalYubeKdF8/iM4+/PO5VGpXAJzJQS64K+Cb7UjhLwCOxBFz+irWNebTiqQy9cWdTw2iaxo92iVtWlS37O1m3vY2dLScGxuvCATJZZWZtiKD/5HOMZ5VGsUxEWSTorkt2elxMqpUXthXZVL7RbmMuVpkit0lhWk2IE7EU3bHUsNudR7ut2CfCtedP59rzp/NaRx/rdrSycU8HPXFn99qBzhi1VQEaq4OEg/5x7YwrhvyZfV04QGd/ggc37uUeilsWSWWypDJZeuPOCde5hFwV9BEOTL5t1ZPrpzGT3rzG6oGjhXJO3cacK1N09MYH9aHYtKdjzK/r8wlTa0LMa4xQGw6e8aLSnLoI8dTgC3kjJc/zZ07h0zct4uGPXcOHrzmXGvd0jN5EmoNdMd44HqUvkeYDy5rHHP945U6ijgT9CM5twCc88lxLSV83l4yP9SZo7Ypy4Fg/bSdiHO1NDNTqs9nKfRdvCdhUlDXLF5DKOMfB546FT2WUNcsXDNznoc37SaYzHO6O8+qRXg53x0mmMzy0ef+4Xz/gd/pOzGuMUB8JntYnePVV80hnlVgqg+LcjvZwz6k1IT7yO/NZ//FrWXllE2F3tpdIZ+mKpvjnn+9j/Y7Wgru/FUN7T4xwcHC6KKRhULGoKomUc6hrZ1+CQydivN7ZT8vxKB29cU5Ek0STadJFXs1SKlYDzmM14Mow0gaPpX/9E3riaXwIIqAKWZT6cIBt/+ddRY0lm1V6E2l6Yie3OhdrW7GqsuPgCdZtb+XZA8cHxmtCft596Wxuv7KZ2fXhov0swynG9u2J5vfJwLK43MqLkN+ztcp2EW4kloAnh8s++xSxVIaA7+SMLZ3NEgn6efGzN5bsdfsSabr6k0VfSwxwsDPK+udb+cmuIyTc7c4+gesumM6qpc0snltX0sQyWZbG5ZrZhwI+qvxubXli+idbAh6JJeDJYdnfPE13NIXPlzcDzir11UG2/cU7S/76vfEUXf2pcR+dNJSeWIofvNjO93a2cawvOTB+0axaVi1t4m0XzijZyoFyaBhUCgNn8AWchFyidp2WgEdiCXhyuGPtFl7v7KMnliaZyRLy+6iLBJg/bUpJlqoNRVXpiaU5EUuSKcFFonQmyy9+e5R129t49UjeducpId53RRM3XzaHusjk2O7shVy7zoDvZAkj6Jfx/HGzBDwSS8CTQzl1Yssl4u5YaWbEqsoud7vzr/YdG9juHA74eNfi2dy+pIlzpo6t0b05XcDn45xpY/r3tAQ8EkvAk0e5dWLLHSbqnLxcmiv07d0xvvf8IX70Ujv9yZNL9d583lRWLmli6bmN1pNhnHwizJ8+dP+OEVgCHoklYDOUYvee6I2nSpqI+xNpfrzrMI/vaKO9Oz4wft70GlYuaeIdF8+adBsaJool4BKyBGxOVapyRu54phPR1Ki2OI9FJqv85rVO1u9o5cXW7oHxhkiQ914+h1uvaGJqTagkrz1ZWQIuIUvAwytlF7JyVereE+lMluP9SfoSpd1c8dsjvazf0cbP93SQdgvFQb9w/aKZrFrSzPkzp5T09ScLS8AlZAn4zMrpwtZEGu5w0V9++vqivU4smeFYX6JkZYmcY30Jnth5iO+/cGig7wTAFfMaWLW0iWsWTDttd585qdgJ2JrxmAHDzXCHO2Z+MifgeY3Vp82AT+09UQyRkJ/mxgg9sTRd0WTJyhLTp1Rx13Xn8aE3n8NPdx9h/fY23jgeZWfLCXa2nKCpIcLtS5q4afFsIqEJ2aBwVrNKvAFGbmDT0hUdtA0Vhu5CNtmMpvdEsYgI9dVB5k2tpjZc2jW84aCfmy+by39+ZBkPrLyUq+Y3AtB2IsY/b9zHB9Y+w7//4jWO9MRHeCYzHjYDngSKUZsdaYY7UTPBcrNi0Uzuhwld0ub3CTNqq6iPBOmKJukvYX1YRLhq/lSumj+V1zv7Wb+9jad3H6E/keGxba2s297K8oUzWLm0icVz60sWx9nKasB5KrEGXKza7Ei1zrO1BlwO4qkMnf1JEuM4JqkQ3dEU33/xEE/sPERn/8ntzhfPqWXVkmaWXzhjMp9cMaxi14CtBFHh8meuIs5t0C8Ft14cqc/uikUzuf+WxcysDdMdSzGzNmzJd4KEg36aGiLMqK0a1GCoVOqrg/zeNefy7Y+9mc+8exEL3RUSu9t7+esf7uZD//Esj2w9ONA03YydJzNgEZkHfB2YDWSBtar6oIhMBR4F5gOvAx9Q1S73MZ8B7gIywB+q6lPu+FLgq0AE+BFwj6qqiFS5r7EU6AQ+qKqvDxdXJc6Ai3WV3ma4lUFV6Ymn6Y6WZmvzmV7zxbZu1m9v49f7jpHLGOGgjxsXz2blkqZJX4rKmSyrINLA/1LVHSJSC2wXkaeBjwA/U9XPici9wL3Ap0XkEmA1sBiYC/xURC5U1QzwJeBuYAtOAr4JeBInWXep6gUishp4APjghP6UE6BYtVkvap2mcCJCfSRIbVVgVEckDaXQs91EhMubG7i8uYFDJ2I8/nwbT750mFgqwxM7D7Fh5yGuWTCNVUubuGJeg213LkBZ1IBF5AngX9yPFaraLiJzgE2qepE7+0VV/869/1PAZ3FmyT9X1UXu+B3u49fk7qOqz4hIADgMzNBhfuBKnAHbzPXsk3/RtanBSaCXz2sY1WOL1de3L5HmyZfaefz5No70JAbGz59Rw8olzVy/aOak3O486WrAIjIfuBJ4Fpilqu0A7m0ugzQB+YdPtbpjTe7np44PeoyqpoFuYNoQr3+3iGwTkW1Hjx4t0k81caw2e3Y5dbngsb4E//jTvew70kfVKBqLF+tstylVAd6/bB7fvOvN/OV7L+FNc+sAeO1oP59/6lXu+PIWvv7M63RFkyM809nN02VoIjIFWA/8kar2DPPWZahv6DDjwz1m8IDqWmAtODPgkWIuRysWzbSEe5Y403LB//rN6zx89zUjNoMf7anNo+X3CW+7cAZvu3AGu9t7WL+jjV/89ihd0RRf/c0bfOvZg7zz4lmsXNrMeWObOU5qniVgEQniJN9vqerj7vAREZmTV4LIHWPbCuSfatgMHHLHm4cYz39Mq1uCqAeOY0wFa+mK0nBKo/X8DTG14SBTqgJ0RYeuD8+pi9DZnyCrynH3+CSfSFEuol08p46/eE8da5Yv4LvPt/HDl9rpjaf50cuH+dHLh1l6TgMrlzZz9XlTbbuzy5MShDhT3a8Au1X1C3nf2gDc6X5+J/BE3vhqEakSkfOAhcBWt0zRKyLXuM/54VMek3uuVcDG4eq/xlSC4ZYLbtrTwR1rt/DWz/+cT3xrB/uO9J5Wllh91Tz6EmkOd8dJZ7IITte0rmiSrfuLMz+ZUVvF3csX8Mjd13DPDQtpbowAsP3gCf7suy/z+//1HE/sPHTaz3E28moZ2nXAL4GXcJahAfwZTh34MeAc4CDwflU97j7mz4GP4qyg+CNVfdIdX8bJZWhPAp90l6GFgW/g1JePA6tVddjFsZV4Ec6cXc500XXVkibW7Wgb8mLslec0cjyaHJgN/4+vPkdrd4xsVgn6fTRWh/D7pGQnHGdV2XrgOOu3t7L94ImB8dpwgJsvm8NtVzQxo7aq6K9bCtYNrYQsAZtKMNRpHw9t3j9s28xUJsuxvgSxZIY7vryFunAAycsJitIbT/Ptj5X2zLz9R/t4fIez3TmVcXJPro68amkTi2bXlfT1x2uyrAM2xozRUBdd/+KJl4etDQf9PubUR+hLpJlbH+FYX2JQc6V4KsvsukjJY18wYwp/cuNF3PXW8/j+C852565oio17Oti4p4PFc+tYtbSZ6y6YflZsd/Z8GZoxZvxG2kqeM6UqwCevv4CsOt9XlFgqQzqrrL5qHhOlsTrEh98yn4c/dg2fvukizp/hzCp3Herhr77/Cr/3lWd5bFtLyRvVe81KEHmsBGEqVaEbcjbt6eBLv3iNg8ejzKoNj7gbrtRUlZ0tJ1i/o41nXuscWC8aCfp595tm874lTTQ1lH6GPhKrAZeQJWBTycZ6EnR/Ik1nX3LCekuMpLUryuM72vjxrsPEU05MAlx7wTRWLWnmsuZ6z7Y7WwIuIUvApXM2nidXSbJZpbM/WVYdzvriaX74Ujvffb6Njt6T250vmDmFVUubeftFMwj6J7aKagm4hCwBl4b1q6gcT73czr//Yj+HukfXqGciZLLKL/ce5TvbW9nd3jswPq0mxK1XzOW9l82lvrq0J4jkWAIuIUvApVHqk4XLUSXO+HN/KAM+CPl99CczY2rUU0qvHOph3fZWNu89inu4M6GAj3ddMovblzQxf1pptzvbMjRTcUbaPjvZ5M/488/Xux/KOgmf2meizuejN5HikedayiYBXzK3jvvmXsKRnjjfe76NH7zUTn8iww9ebOcHL7Zz1fxGVi1tZtm5jRXRFtMSsCm5yXKe3GhntZV6gvSpfyh9PqEuHOTg8T7+12MvlFVZYlZdmDVvO58Pv2U+P951mMd3tNF2IsZzr3fx3OtdnDu1mpVLm3jnxbNG1SXOK7YO2JTcRJ4sXCojnRqdr9xPkM71jLjugY3csXbLwM8w1FriY30J+pNZumNJplaHOB5N8ODGvUXrGzFekZCf913ZxNc+ehV/c9tirnD7Ir9xPMoXnt7LB9du4Su/OsCxvsTwT+QRS8Cm5CZDz+JCzt4b7aYILwz3h2SoP5Rd0RRTa4LUVAUJ+H3UhYNUBXw8sq2w/sGl5hPh2vOn84UPXM6X//tSblw8i6Bf6Imn+dazB/lvX36W//9Hu/ntkd6Rn2wCWQnCTIhK71lcSB17zfIF3LdhF9FketCqj+Fm/BN10W648sjDd19z2rFU3bEU02pONsoREaZUBTjWGycS8hNLll9Hs/NnTuHTNy3iY29dwIadh9jwwiFOxFL8dHcHP93dwWXN9axc0sy150/zfLuzJeBJqhKvwpezQurYhZ6vN5EX7Ub6Q3LqH8qhVrDEUhnmTa1hTn2E3niK4/1JMtnyW001tSbER35nPv/tzefw091HWL+jjQPH+nmxtZsXW7uZUx/m9iVN3LR4NjVV3qRCW4aWZ7IsQ7N1t8VXyn/TiVymV+hrjebnzmSd5u7ltIljKKrKjoMnWLe9lWcPnKxh14T8vPvS2bzvyibm1A+/3XnSnQlniq+QeqUZnVLWsSfyol2hF0RH83P7fcKM2irmNkTK+iBOEWHpuY383e2X8tXfv4pbLp9LVcBZ77xuexv//Stb+eyGXbzc1l3wSdNjjslmwCdNlhnwdQ9spCESHLQOUlXpjqX45aev9zAyM5SJ3qgy1p4Ro9UdTdEVTZKtgNzSE0vxgxfb+e7ONjr7Th4getHsWlYtaeZtF04nkLfd2XbCldBkScBn486zSjYZS0bpTJbj/cmKaSeZzmT5xW+Psm57G6/mrZSYPiXEbVc0cfNlc6iLBC0Bl9JkScCT8Rd6siv1rNQr0aTTaS2VKY9OayNRVXa5251/te/YwHbnqoCPdy2exfuXzmP5hTPG8tSWgEcyWRIwTN5faFN5VJUT0RQnhjiluZwd7o7z3efb+NFL7fTnLbd7+0Uz+IubL+H8GVMKeTpLwCOZTAnYmLEo5fLFZDpLZ3+iLNcOD6c/keapXYdZv6ON9u44AZ/wq09fz+z6cCFPYwl4JJaAzWQ3XIKdqNJVOa8dHk4mq2zZ30lPPM0n3n5BoQ+3ZWjGnM1G6mcxUcsXa8NB5jVWUxuemB6+xeL3CW9dOGMsyfeMLAEbc5YYKcFO5Hpkn7t2eE59ZMJPtSgnZ+9PbsxZZNOeDnYc7OKNzn72H+0b2LWWn2C9aCIUCflpbozQWB2qiP69xWYJ2JhJLld6EHE2EqSzyqETcXrjqUEJ1qu2oSJCY02I5sYIkVD59u4tBUvAxkxyudLDrNqwc9y7+z+Hu+ODEqzXbUODfh9z6iPMqK3yvEvZRLFuaMZMcrkOaBJyktqxvgTJjKJwWoIth7ahteEg1aEAXdEkPbHybvAzXpaAjSmBcmoHmt9Ksy4SpC4SHNia7nWyPRO/T5g+pYracIDOviTx1MStHd66/ziPPNdCe0/pj2CyEoQxRVbI8UUToZKPhKoK+JnbUPyyxNb9x/nUoy9wx5e38KlHXxg4Ymnr/uM8uHEvnf0J6sIBOvtLewSTzYCNKbJyO5Sz0Abx5ag2HKQmFKCzCH2Hc0k24JNBSfYeFvLIcy0EfDKwHC+3IeWR55wjmB7Z1sKxvkTR3tXYTrg8thPOFIO1Ay2tWDLDsb7EmBv8fOrRF+jsTwxa8xxLZZhWU0V7T4y6cADJ27imKMf6koSDfoJ+56ToMewStJ1wxkyEcj6UczLIrR2uj4xtJ117T4xwcHDqCwd9HHZrvvHU4MQeT2VJprMEfEI46C/qLkFLwMYUWSXXXCuFiDBtythO4ThTkp3tXnBLZ5VYKoPi3KazStAvpyXtYuwStARsTJF5vZ72bBIO+mlurGZaTRW+Ue6kO1OSza12uOf6hUyrqaI3nmZaTRX3XL+Q+dOmnJa0i/GuxmrAeawGbEzlSmeydPYn6R/FKRy5pWaHe2IDM9/hlprlLtwVuwZsCTiPJWBjKl+pTuHYuv84j2xrobMvMZaVJEMmYFuGZoyZVKpDAcINfmcnXTxdtFM4rl4wlWvOnzbWM+GGZDVgY8yk4/PlLtKFCQfLt8GPJWBjzKSV20k3vUwb/Ez6BCwiN4nIqyKyT0Tu9ToeY8zEqwsHaW6sZkq4vKqukzoBi4gf+Ffg3cAlwB0icom3URljvOD3CTNrw8ypjxDwlUfqK48oSudqYJ+q7lfVJPAIcKvHMRljPJTbSVcOZ9JN9gTcBLTkfd3qjg0QkbtFZJuIbDt69OiEBmeM8Ua5nEk32RPwUFX3QWtSVHWtqi5T1WUzZsyYoLCMMeUgNxsuZCddMU32BNwKzMv7uhk45FEsxpgyJCLUVweZN3XiL9JN9gT8HLBQRM4TkRCwGtjgcUzGmDKUu0g3lgY/YzWpE7CqpoH/CTwF7AYeU9Vd3kZljCln4aCfpoaJKUuU16K4ElDVHwE/8joOY0zlyJUlaqr8o27wMxaTegZsjDHjEfD7mFUXZlZduCRrhyf9DNgYY8arpipAJOinOza+8+hOZTNgY4wZBZ9PaKwJFfc5i/psxhhjRs0SsDHGeMQSsDHGeMQSsDHGeMQSsDHGeMQSsDHGeMQSsDHGeMQSsDHGeMQSsDHGeMQSsDHGeMQSsDHGeMQSsDHGeMQSsDHGeERUdeR7nSVE5CjwxjifZjpwrAjheKFSY7e4J1alxg3exX5MVW86ddAScJGJyDZVXeZ1HGNRqbFb3BOrUuOG8ovdShDGGOMRS8DGGOMRS8DFt9brAMahUmO3uCdWpcYNZRa71YCNMcYjNgM2xhiPWAI2xhiPWAIeBRH5TxHpEJGXh/jen4iIisj0vLHPiMg+EXlVRG7MG18qIi+53/uiiIgXcYvIJ93YdonI5yshbhG5QkS2iMhOEdkmIleXYdzzROTnIrLb/be9xx2fKiJPi8he97axnGIfJu6/F5E9IvKiiHxXRBrKKe7hYs/7ftn+fgKgqvYxwgewHFgCvHzK+DzgKZzNG9PdsUuAF4Aq4DzgNcDvfm8r8BZAgCeBd0903MDbgZ8CVe7XMysk7p/kXhf4XWBTGcY9B1jifl4L/NaN7/PAve74vcAD5RT7MHG/Cwi44w+UW9zDxe5+Xda/n6pqM+DRUNXNwPEhvvWPwJ8C+VcybwUeUdWEqh4A9gFXi8gcoE5Vn1Hn/+2vA7d5EPfHgc+pasK9T0eFxK1Anft5PXCoDONuV9Ud7ue9wG6gyY3xa+7dvpYXR1nEfqa4VfUnqpp277YFaC6nuIeL3f12Wf9+gpUgxkxEbgHaVPWFU77VBLTkfd3qjjW5n586PtEuBN4qIs+KyC9E5Cp3vNzj/iPg70WkBfgH4DPueFnGLSLzgSuBZ4FZqtoOTsIAZrp3K7vYT4k730dxZoVQhnHD4Ngr5fczUOoXmIxEpBr4c5y3aKd9e4gxHWZ8ogWARuAa4CrgMRFZQPnH/XHgj1V1vYh8APgK8A7KMG4RmQKsB/5IVXuGKSWWVeynxp03/udAGvhWbugM8ZXFvzlOrBXx+2kz4LE5H6d+9IKIvI7z1myHiMzG+cs5L+++zThvl1s5+RYuf3yitQKPq2MrkMVpUFLucd8JPO5+/h0gdxGurOIWkSBOIviWqubiPeK+xcW9zZV9yib2M8SNiNwJ3Ax8yH1rXlZxuzGeGnvl/H6Wusg8WT6A+ZxyES7ve69zssi/mMFF/v2cLPI/hzPzzBX5f3ei4wb+ALjf/fxCnLdjUgFx7wZWuJ/fAGwvt39v93W+DvzTKeN/z+CLcJ8vp9iHifsm4BVgxinjZRH3cLGfcp/y/f0s9QtMhg/gYaAdSOH8pbzrTP8Hu1//Oc7V1VfJu5IKLANedr/3L7g7EScybiAEfNONYwdwfYXEfR2w3f3leRZYWoZxX4fztvVFYKf78bvANOBnwF73dmo5xT5M3Ptw/kDnxv69nOIeLvZT7lOWv5+qaluRjTHGK1YDNsYYj1gCNsYYj1gCNsYYj1gCNsYYj1gCNsYYj1gCNmacRGT+qR3njBkNS8DGGOMR6wVhzBmIyPdwtq2GgQcBP3Ceqv6p+/2PAEuB/wv4ReTLwLVAG3CrqsY8CNtUENuIYcwZiMhUVT0uIhGcbao3AL9W1Qvc7z8J/C3Obr19wDJV3SkijwEbVPWbXsVuKoOVIIw5sz8UkRdweuHOw+0dICLXiMg04CLg1+59D6jqTvfz7Ti9LIwZlpUgjBmCiKzAaXf5FlWNisgmnFLEo8AHgD3Ad1VV3XaTibyHZ4DIRMZrKpPNgI0ZWj3Q5SbfRThdssBpiXkbcAdOMjZmzCwBGzO0HwMBEXkR+GucMgSq2oXTovFcdfopGzNmdhHOGGM8YjNgY4zxiCVgY4zxiCVgY4zxiCVgY4zxiCVgY4zxiCVgY4zxiCVgY4zxyP8DXewYq/vZgWUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sns.lmplot(x=\"avh\", y=\"rgdpe_per_capita\", data=pwt2019)" ] }, { "cell_type": "code", "execution_count": null, "id": "cff81d57", "metadata": {}, "outputs": [], "source": [ "pwt2019" ] }, { "cell_type": "markdown", "id": "02bad4a8", "metadata": {}, "source": [ "The 95 % confidence interval is displayed in the summary output, but we can also quickly make the calculation ourselves:\n", "\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "dec0de50", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[-60.5175164514408, -24.205978222312357]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ex1_cov = ex_mod1.cov_params()\n", "SE_b = np.sqrt(ex1_cov.iloc[1,1])\n", "\n", "\n", "CI = [ex_mod1.params[1] -2*SE_b, ex_mod1.params[1] +2*SE_b]\n", "CI" ] }, { "cell_type": "markdown", "id": "bb731c5f", "metadata": {}, "source": [ "We can see that there is quite a bit of variation around the regression line. We can get the estimate for the $\\sigma$" ] }, { "cell_type": "code", "execution_count": 10, "id": "ca8b1905", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "19705.147498247465" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.sqrt(ex_mod1.mse_resid)" ] }, { "cell_type": "markdown", "id": "8fa9f5e9", "metadata": {}, "source": [ "Initially, we might postulate that higher average hours would causally increase per-capita GDP. Instead, the results of this regression suggests a reverse causality: as GDP grows, people have a preference for more free time and the average hours worked goes down. But again, notice the very high variance around the regression line ($\\sigma$)" ] }, { "cell_type": "markdown", "id": "d944104c", "metadata": {}, "source": [ "\n", "\n", "## 2.) Fake data simulation with a lurking variable. \n", " - Generate a model based on a lurking variable - generate X, Y and Z variables where both X and Y depend on on Z (but not on one-another). \n", " - Run a regression of Y on X, does the regression suggest a relationship between the variables. \n", " - Now add Z as an exogenous regressor in the regression. How does this change the suggested relationships." ] }, { "cell_type": "code", "execution_count": 12, "id": "d83e4dd1", "metadata": {}, "outputs": [], "source": [ "n=200\n", "Z = spt.norm.rvs(0,10,n)\n", "X = 3 + .3*Z + spt.norm.rvs(0,3,n)\n", "Y = 5 + .5* Z + spt.norm.rvs(0,3,n) \n", "\n", "fakeData1 = pd.DataFrame({\"X\":X, \"Y\":Y, \"Z\":Z})" ] }, { "cell_type": "code", "execution_count": 13, "id": "deee559b", "metadata": {}, "outputs": [], "source": [ "as_mod2 = smf.ols(formula='Y~X', data=fakeData1).fit()" ] }, { "cell_type": "code", "execution_count": 14, "id": "295854b0", "metadata": {}, "outputs": [ { "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", "
OLS Regression Results
Dep. Variable: Y R-squared: 0.379
Model: OLS Adj. R-squared: 0.376
Method: Least Squares F-statistic: 121.0
Date: Thu, 18 Aug 2022 Prob (F-statistic): 2.89e-22
Time: 13:42:07 Log-Likelihood: -568.23
No. Observations: 200 AIC: 1140.
Df Residuals: 198 BIC: 1147.
Df Model: 1
Covariance Type: nonrobust
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [0.025 0.975]
Intercept 3.4916 0.355 9.838 0.000 2.792 4.191
X 0.7660 0.070 10.998 0.000 0.629 0.903
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 0.901 Durbin-Watson: 2.237
Prob(Omnibus): 0.637 Jarque-Bera (JB): 0.932
Skew: -0.159 Prob(JB): 0.628
Kurtosis: 2.898 Cond. No. 6.21


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified." ], "text/plain": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: Y R-squared: 0.379\n", "Model: OLS Adj. R-squared: 0.376\n", "Method: Least Squares F-statistic: 121.0\n", "Date: Thu, 18 Aug 2022 Prob (F-statistic): 2.89e-22\n", "Time: 13:42:07 Log-Likelihood: -568.23\n", "No. Observations: 200 AIC: 1140.\n", "Df Residuals: 198 BIC: 1147.\n", "Df Model: 1 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept 3.4916 0.355 9.838 0.000 2.792 4.191\n", "X 0.7660 0.070 10.998 0.000 0.629 0.903\n", "==============================================================================\n", "Omnibus: 0.901 Durbin-Watson: 2.237\n", "Prob(Omnibus): 0.637 Jarque-Bera (JB): 0.932\n", "Skew: -0.159 Prob(JB): 0.628\n", "Kurtosis: 2.898 Cond. No. 6.21\n", "==============================================================================\n", "\n", "Notes:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "\"\"\"" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "as_mod2.summary()" ] }, { "cell_type": "code", "execution_count": 15, "id": "8d9b5885", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABLeklEQVR4nO29eXRb93nn/fld7CDBVaRWSiJtOfISL7GsyrasKE7aOnFrN81mt3GzOJbfafLG7UznJDNtOhl3OhO3TTvpTNvXcuw4W+2kbnPiZnEWO4ok24rkJY5lW95ISaREifuGHbi/948LQCAJkAAJ4ALg8zlHh+LlxcUPy/3e5z6r0lojCIIgVB7D7gUIgiCsVESABUEQbEIEWBAEwSZEgAVBEGxCBFgQBMEmnHYvoBTccMMN+rHHHrN7GYIgCPlQuTbWhQU8MjJi9xIEQRCKpi4EWBAEoRYRARYEQbAJEWBBEASbEAEWBEGwCRFgQRAEmxABFgRBsAkRYEEQBJsQARYEQbAJEWBBEASbqItSZEEQhHKy79gQ9+7vpX88RFernzt39bB7a+eyjysWsCAIwgLsOzbEnz/6EkPTEVp8LoamI/z5oy+x79jQso8tAiwIgrAA9+7vxeVQ+N1OlLJ+uhyKe/f3LvvYIsCCIAgL0D8ewudyzNrmczkYGA8t+9giwIIgCAvQ1eonHE/O2haOJ9nQ6l/2sW0TYKVUl1LqZ0qpV5RSLyml7kptb1NK/UQp9XrqZ6tdaxQEQbhzVw/xpCYUS6C19TOe1Ny5q2fZx7bTAk4A/0lrfSGwA/ikUuoi4LPA41rrLcDjqd8FQRBsYffWTu6+6WI6A14mw3E6A17uvunikmRB2JaGprUeBAZT/59WSr0CrAduBnandvsqsA/4jA1LFARBACwRLoXgzqUqfMBKqc3AFcAvgNUpcU6LdM5XrZTao5R6Rin1zPDwcMXWKgiCUCpsF2ClVCPwr8Afaa2nCn2c1nqv1nqb1npbR0dH+RYoCIJQJmwVYKWUC0t8v6m1/rfU5rNKqbWpv68Flp/tLAiCUIXYmQWhgPuBV7TWf5v1p0eBj6T+/xHgu5VemyAIQiWwsxfEtcBtwItKqV+mtv1X4AvAt5VStwMngQ/YszxBEITyYmcWxEFA5fnzOyu5FkEQBDuwPQgnCIKwUhEBFgRBsAkRYEEQBJsQARYEQbAJEWBBEASbEAEWBEGwCRFgQRAEmxABFgRBsAkRYEEQBJsQARYEQbAJEWBBEASbEAEWBEGwCRFgQRAEmxABFgRBsAkRYEEQBJsQARYEQbAJEWBBEASbEAEWBEGwCTtnwgmCsELZd2yIe/f30j8eoqvVz527eti9tdPuZVUcsYAFQago+44N8eePvsTQdIQWn4uh6Qh//uhL7Ds2ZPfSKo4IsCAIFeXe/b24HAq/24lS1k+XQ3Hv/l67l1ZxRIAFQago/eMhfC7HrG0+l4OB8ZBNK7IPEWBBECpKV6ufcDw5a1s4nmRDq9+mFdmHCLAgCBXlzl09xJOaUCyB1tbPeFJz564eu5dWcUSABUGoKLu3dnL3TRfTGfAyGY7TGfBy900Xr8gsCElDEwSh4uze2rkiBXcuIsCCIMxCcnQrh7ggBEHIIDm6lUUEWBCEDJKjW1nEBSEIQob+8RAtPtesbSs1R7cSrhixgAVByCA5uhaVcsWIAAuCkEFydC3x/fTDz3N6IsyZyQjTkUTZXDEiwIIgZFjpObppyzcYS+AwIJHUnJ4MMxWOl8UVIz5gQRBmsRJydPP5d9NBSK/TQcLUGIYCE0ZmojgdquSuGLGABUFYUSzk3003CuoIeNAaTK1BaaIJsyyuGBFgQRBWFAul2qWDkAGvi3UtXpyGImFq/G5HWVwxtgqwUuoBpdSQUupo1rbPK6VOKaV+mfr3HjvXKAh2s+/YELfuPcTOe57g1r2HpChimSzUDjM7CNnocbKm2cv6Fj9/f8sVZXHL2G0BPwjckGP732mtL0/9+0GF1yQIVYNUppWehVLtKh2EtDUIp7Xer5TabOcaBKGayb5dBvC7nYRiCe7d31v3gbJyceeuHv780ZcIxRL4XA7C8eQs/24lg5B2W8D5+JRS6lcpF0Vrrh2UUnuUUs8opZ4ZHh6u9PoEoSLI9IjSU02pdtWYhvZPwF8AOvXzi8DH5+6ktd4L7AXYtm2bruQCBaFSdLX6GZqOZCxgWJmVaaWmWlLtqs4C1lqf1VontdYmcB+w3e41CYJdSGVafVN1FrBSaq3WejD163uBowvtLwj1zO6tndyN5QseGA+xoQxNYaT/r30ore27e1dKPQTsBlYBZ4H/lvr9ciwXxHHgzixBzsm2bdv0M888U8aVCkJ9ks6ycDnUrIDUSio/rhAq10a7syBuzbH5/oovRBBWKJJlYS9V54IQBKFy2NH/V1we56i6IJwgCJWj0v1/a7WwJGlqhqejxJNmSY8rAiwIK5hKZ1nU4sijmWiCgfEQ05F4yY8tLghBqBDVeOtdiSyLbKph5FGhn0M8aTI6EyMUS5RtLSLAgrBEihHU7GyD7Fvvu6EqRDi9hvRr+rPvHi3LRcLuwpJCP4fJUJzxUMxqR1lGRICFuqTc1mYhJ3L2GqbCcfxuB80+L1Cd2QaVuEgs1oeh3CyW9RFNJBmZiRGd4xcvF+IDFuqOSgR6FvNlzl1DMJZgNBhjKnzOj1htPR0q4Z+1uw9Dvt4a/WNBRmeinJ6IVEx8QSxgoQ6pRG7rYr7MuWvwOh3EkiYjM1GaUo+rtp4OlfLPlrsPw0J3P7lcIMFYgo7UBaHSiAUs1B2V6CC2WPrW3DV0BDygIZowq7anQz2MpF/s7ic768M0TabCMSJxkw9t67JlvSLAQt1RCSFZLH1r7hoCXherAm78boftLRDzUQ+NfxZzo6RdIG0NbsaCMVr9Hu66fgvbe9psWa+4IIS6oxKBnsXSt3KtweVw8Pe3XFpVoptNpVPSysFibpRYwuQtawPc875L7VjePESAhbqjUkKykC+zVsWskn1yy5Gpki/NbX2Lj4lQjPFQHDsbkM3F1m5opUK6oQm1TjUWaZSTcnVhy3XcaMLkruu3cOXmnMN1CsbUmmA0yWVdLUt5ePV1QxOESlKtIlfNRRrlothMlUI/u+w7j/6xIKubfHzgyg3LEl+tNUeOj3PfgV4GxsPs+8+7WdvsW/LxshEBFlYE1SxyK7ElZDEpb8V+dru3dnJVdxujMzES5vKa57wyOMV9B/r4Zf9EZtvjrwzx4R2blnXcNCLAwoqgmkWuGvojVJpiSpKL+ewSSZPRYIxgdHn9G/rHQtz/ZB/7XxvJbLtiYwt/duOFXLmpdBkTIsDCiqCaRc7u/gh2UEymSqGf3VQkztjM8vo3jM5E+drTJ/j+i4OYqcOc39HIHbu62baplY3tDUs+di5EgIUVQTWLnN39EeygmCyRxT67WMKqMIwso4R4JprgW0f6+ddnB4gkLLfF2mYvH7+2m3ds7cBQOWNoy0YEWFgRVLPIFSpGpQwiVkNAstCUt3yf3Z7rupkIxXjs6BkePtzP4FSYtU0+brmqq+DCiljC5LsvnOabh04wFbHcFi0+Fx/esYnfvmwtLkd5a9UkDU1YMaRFp5byctOUMm2rGgdxLnZBmPvZfeyazVy0vomDr43wpSdex2kovC6DSNwkYepFq9uSpuanr5zlK08eZ2g6CoDXZfDBK7v44FUbZlnb2XS1+ZcqyjlNaBFgQbCJYqzQW/cemncbHool6Ax4eWjPjqKet5THKgXFXBBMUzMeimUa5/zHb73AaDA6q+9GOJ6kvcHD337osnnPpbXmF31j3Hegj76RIAAOQ/Hbl67lwzs20dbgXnCtpRZgcUEIgg0Um1pVyiDiUo5VTpdFoVkOwWhiXmrZ4FSYJu9sGfO6DM5Mhec9z8unp9h7oJdfDUxmtl2/tZOPXbuZ9S3z83oP947x8JHZro2uttLGDESABcEGik2LK2UQsdhjlTuHerELwkKjgdY2+eZZwJG4yZqmc4J6cjTElw/2cfCNcyllV25q5Y7rurlgdSDnmg73jmVcG01eJ6PBKF964nVWBdy866I1y3q92YgAC8IClCvwNTwdZU2TZ9bfF7JCSxlELPZY5c6hXuiCMBmKMxaK5e3fcMtVXXzpidcJx5OzfMC3XNXF8HSUrz59nMeOnsmklG3pbGTPrh6u3LRwZdzDR/pxGioj7On36b4DfSUVYGlHKQh5KOVkjbnHUgpOTURmTdpdyAot5SSJYo9V7v7KudpgxhIm73/bekaD0QWb52zvaeOu67fQ3uBhOpKgvcHDnp09vHBqgtseOMwPXrTEd12Ll8/deCH/9OG3LSq+YLk2vK7Z8uh1GZwqcd64WMCCkIdSWn5zj7U64OXURJgzkxEaPc6CLNpSdior5ljlzqGe279hTbOP91+5gSsKEEqwRHh7TxvReJLv/PI0f/v4a0ynUspa/S5u27GJGy8tLqUsn2tjfYnzxkWABSEPyw18LeRysMYSac5MRZkMx6s6La5S/ZV3nNfO8HSUeLK4/g1JU/Pjl8/y4JPHGZ6xUsp8Lge3XNXF+6/cgM/tWOQI88nn2rjjuu6ij7UQIsCCkIflWH5zA1cjM1FOTURQqYqq4ekokUSSBreTv7j5kqoS3lx+77tvurhsOdRJUzMajDITKa5/g9aap94c5csH+zgxal0UnYbipsvX8eFf20iLf+GUsoXY3tPGXWzh4SP9nJkKsyaVBfH2t5T2cxIBFoQUc4Xn6p42Hnnu1JIsv3wuh9MTYbQGFBhK4Xc7qqYrG5y7cMSTSSZDcQYnwzx3cpxP7j6vLDnCk+E4E6EYSbO4eoSjpya5d38vL52eAqwk23deaKWUlapVZNq1YShFo9dJk9dV8so4EWBBIHeq1SPPneL9b1vP071jRVt+c90XaZfDyfEwDqXwOAxWNXpo8rkq0pWt0GyOe/f3Ek8mGZ2JoxS4HAZJU/MP+97k0g0tJVtjNJFkZCZW9Aj4vpEg9x/s46k3RzPbrtrcyh3X9XB+Z2NJ1pbGaRg0+SzhNQzpBSEIQHmKAvIF3J7uHVuS5ZfLfeF0GLgdBls6GzOuCCh/V7Zi8nj7x0NMhizxTTegcShImLokFwnT1EyE40yGixsNNDQV4cGnTvDjl8+llL1lTYA913VzxcblTbqYi8th0Ox3EfA4Z31O5UAEWKgpylUUUOp2lfkCV93t1rTkYv3Ky7noFJPN0dXqZ3AyPOtWW2vwOI2i34u5a/6Dqzdx8brmgpukH+4d4xu/OEHvSJBwLElarje0+rh9Zze7tqwqmUAe7h3j28/2c2Yqwqa2hooFRCUPWKgpFhs7vlTyjbJvcDu4de8hdt7zBLfuPVRwDnC+XNvPvvvCoke/LzcfuZg83jt39eA0LLeD1hrT1JhoAl5nUWln2Wtu9jo5PRHi7u+9zFNZ1WgLcfC1Ef7yB6/w0uAUoZT4Ggree/k6HvjINt5+QUfJxPeFkxP8331vMBmO0+Z3Lyvfu1jEAhaA6mhPWAjlaqyey2KdDMdRQNzUC1rb+d67XLm2+44N0eB20JtqBNPd7udzN1604Hu93HzkYrI5dm/t5JO7z+Mf9r1JLGFmWsiMh+L8XoEtHtNrdhrgdhrETY3H6cDUSR4+0r9ol7LHjp7h7x9/nXjK12AoaPW78boM+kZCOEsQCFNK0eBx0Oxz8affOYrHadgyLUUsYKGkFV/lJp+lutyigFwWa0cqSLaQtV3Me5feN5Y02dLZyIZWH6H44rfjy61Ey1VptpDV/el3XcAnd5+Hy2mgFHidBm0NLh557lTB34mTY0HLkk5q0r6DfE1ywEopO/D6CLd/9Rm++JPXiJsaBbT6XHS3N9De4MbvduR9fKEYStHsc9HV6qMz4MXjdJS90m8hxAIWqnpe2lzKWRQw12Ldec8Ti1rbxbx3S32fl1uJVsz0iTRP946xodU3r2XlYmtNmpqxYIyORu+iTXLSvDAwwX37e3l5cBqwjO4Wvwuv0yDgdS36+EJwGIomr4smnwvHnIwGO6el2CrASqkHgN8ChrTWl6S2tQHfAjYDx4EPaq3H7VrjSqCa56XNZSlislQKOTFzvXeJpMlzJ8fZec8Ts1wSS32f5150RoNRxoJW/uytew8V9PqLLWNeylqnI3HGglZO70JNctL0Ds/w5YN9HOody2zb0dPGJ3Z2MzIdW/TxheByGDT5XDR582c02DktxW4L+EHg/wJfy9r2WeBxrfUXlFKfTf3+GRvWtmKo5nlpuShlT4SFKOTEnPveTYXjnJqI4DTmZ2lk7zsdic+qhtt3bCjva8q+6Lw+NM10JEGr38WqRk/JW0Pme12Q/zuRayZbvkqy7T1tnJmK8OCTx/nJy2czmQ0XrW1iz65uLt3QAkBPB3kfXwhup0GL302jZ3GJq+RFfS62T8RQSm0GvpdlAb8K7NZaDyql1gL7tNZvWegYMhFjeVTjiJpqYbExRnPfuzeGZ0gkNRtafZnb54HxIFORJFprNNDgMogmdCbA1d7gxu10FPR+V2qaRSHfCa01k+E446HCcnonQ3G+efgE3/3laeJJa/+NbX5u39nNzvPbS5LV4Hc7afa5ltT/oczUzESM1VrrQYCUCK9sBagAdloA1U4+azs78yHgcWbESGtY3+LNiO/ZyTDjIavHgcepSCQ1MzEr8GYoKxo/HUkQ8FKQz71S7qLFvhOReJKRmaiVKbEI4XiSf312gG8d6ScYs6zkVY1uPnL1Zm64ZM08n2yxOAxFo8dJwOvC7aytvIJqFOCCUErtAfYAbNy40ebV1D6Vuq2vB+YWg6Stw7+4+RLu3d/L0HQks+9IMAZYYmsoA7eTTBaHy7CyDBJJzWgwRiI5tehzV9JdlOs7MXcm20IkkiY/PHqGrz59grHU+9DocXLr9i7ee8V6vK7lWak+t4OA10WD21H2irVyUY0CfFYptTbLBZEz70VrvRfYC5YLopILFKqPSuYxL5TNMNdvnC6bdRrzLbN0fwGlLGGLJRf/GtsZMJqKxBkPLt44R2vNz18b4YEn+xgYt9LGXA7Fe69Yz+9t35jqi7E0HIYi4HUR8DrLPjK+ElSjAD8KfAT4Qurnd+1djlDtlHtm2VwWcgPMvXV3GAq0znmbbWqNsv4MmoJun+1wF0XiSUaDhTXOef7kOHsP9PHqGSulzFDwmxev4SNXb6KzybvkNfjdTgJeJ/4atnZzYXca2kPAbmCVUmoA+G9YwvttpdTtwEngA/atUKgFypXHnM+qXswNkH3r/vc/fY0vPfEGCdPEUGQs4kaPA60hljRxOwyaGlxsbi+sm1el3EXpnN7ssUn5eGNohvsO9HLk+LmM0WvPa+f267rZ3N6wpOd3GgaNXmfdWLu5sFWAtda35vnTOyu6EKGmKUdgaiGruhg3wKffdQEAXz7YRzBm9Zb4rbd28OzJyXkZBpVwIxTqqinU3XB6IsxXnjzO41kVcpesa2LPrh4uWd+8pDX63U6afM5ZF7h6pf5foVD3lCMwtZBV/dCeHUW5AT79rgsyQpxmsfS2clCIqyaaSDI6E5uV05uLiVCMbxw6yaMvnCaREulN7X7uuK6bq3uKTylTyspkaPYtnslQK31LCkEEWKh5yhGYWsyqXq4bwI6sk4UuKrsu6Fg0u+Fw7xjf/MVJekdnCMeSGXdKZ8DDR6/ZzK9ftLrolLJ00/OAd36JcC4q7e8vNyLAQs1TqsBUtmU1FY6TNE1WNZ4LHM21qmvNEst3UTk5FmRgPLxgn96n3hjhr370KjPRREZ4lYJ3X7yGT79zS9H5tx6X1Yms2BSyWupbUggiwEJdsFyLct+xIf7kkReYiSYsv6fWzEStAor2Bs88q7oWLbG5rhqtNdOROB2N3rzia2rNvleH+avHXiWWmlacbpbjdzk4PREpSnz9bie/6p/gK08dX9KFqxh/fy1cIEWAhRXPvmND/OE3nyMUT6IAp2EVTSRNk8lQAqdhzLOqK2GJlVpA0q6aYDSO22EQjCUXbHDz7Ilx9u7v5fWhmcy2Zq+TtgY3LoeBRhfUHjLde7fF5+apN0b4Hz94ZckXrkL9/bVygRQBFlY06RM1lAo6aSBugsvQOB2KpNYc+Mz18x5X7pLgcgjI7q2d/Jd4kr37exmczN/g5rWz09y3v5dnT05ktjX7XPhcBk1FtIdUShHwOmnxuTJN1Jd74SrU318rrgoRYGFFkz5RM17IVGFEUmucC/gmy10SXGoBSWc3vGVNgC9+8LKc+5waD/PAk3387NXhzLZLNzSz57oeZiKJgttDGkrR5HPRnKP37nIvXIX6+2ulxaoIsFAzlMOnlz5RPU6DSMJEpQJMprb+benIXURQ7pLgUgmIaWrGQjGmFshuGAvG+PrTJ/jei4OZvN+eVQ184rpufq27LRMky9ce8nDvWGb7hlY/f7j7PK6/cHXO5yrFhasQf3+ttFgVARZqgnL59NIn6ppmLwPj4YwAWXPIXHzmhq05H1fukuBSCMh0JM54MD4rwJYWy8GpMJ2NXjoCHp7qHSGSGo3UGfDw8Ws3884L56eUbe9pm+euONw7xt8/8Toep0F7g5vxUIzP//vLGErlfC8q1cvCzp4ZxWB7P+BSIP2A659y9cHNFvZE0uTsVJS4abKlo5HPvvvCsvoLF7Lo5/bjHZmJMh6KE/A62dIZWFDs8xVTHO4d40tPvI5DQTRhMhqMZVLKmrxOfn/HJm6+bF1BWQ3pwolP/fNzjMxEi/pcKlWEYkexywLUTD9gQZhHuXx6cy3ZKza2VkVV2qwpGGenmI4maWtw0d6QfwqG1lbvhqlIImeD9IcOnyQaTzIdSWQmDissq/e+j2wraHrE3MDaqYlw0Z9LpYpQaqHFqgiwUBOU06eXnVrWPx7KTD0u58lbSJAtLSBzrX+/28nITIRPP/w8TT4XXa1+Pnr1Ji5a30w8OT+fV2vNkePjvDQ4lSkbBmj2OWnzuwjHzUXFN90Gcm5grVZ8rdVKfbYYEuqOYkerF8q+Y0Pc8Hc/5/avHeH5k+M4FAuOli/22LfuPcTOe57g1r2HZh2vmFHoc/edjsQZmY4RiiVp9jo5PRHi8997mSdfH5n32GNnpvhP//IrPvtvL2bEt9HjZHObn9UBLwmTBVPJ3E6DVQEPG9v8tDW45/mFy/W5rBREgIWaYPfWTu6+6WI6A14mw3E6A95lz6xLuwGOj4VwKIUGBiejJE2Ny6EylvByjj00HZnlYkiLcFerPzMZI00+y3HuvsPTUTQat0MRNzUepwOnoXj4SH9mn/6xEJ//95f4w28+zy/7JwArs6G9wU2r34XLqQjH8xdieF0O1jR72dDqp8nrylsuXI7PZSUhLgihZii1Ty/tBkiaGodSKKUw0QxPR+le1bAs//JiLoZiovS59jUUtPjdpMcKe10GZ6bCjM5E+drTJ/j+i4OZANvaZi8+l4OZaJxmrxNSc+hyFWI0epw0+VxFjQuqBV9rtSICLKxY0oE9t8MgkbSmUyhlNUlfrh+zkG5qhaaxpff9p5+/yckxyx3hcxmz/LahWBKt4bb7DxNJDcpc2+zl7Vs62PfaEPGkSbPPlSqgMPmjd16QEV4jFVhrzqpYEyqDCLBQ9yw22WJVo4fTk2EwQWNZw8v1YxYSnCrUcjRNzaVdLXzhfZeitc6kk4XjSdxOxehMnMmINZEZoMXn4sM7NvLbl63jM4+8iMthZHzIaQv64SP9XHP+Kpp8Tpq8rsx8OqGyyOVOqGsW8sWmA0hOh2JdsxdlQFJD96qGZfsxiw1O5QvYTUfiDIyHmQjFMqll23va+H/fcT4KRd9IiImwJb5el8Ef7NjENz6xnd992wZcDoPBqTBe1+zT3OdyMDwToavNR4vfLeJrI2IBC7ZTzraBxUy2uKKrdDnAxbgY5uYEHx+dYc/Xn8HvdrC5vXGWn1ZrzS/6xvjywT5OptwZDkPxW29dy21Xb6KtwT3r2GubfIwGo/hcVt9dp6GIJJJsbGsouA9vLbR1rFVEgAVbKXfbwHJPtliIQo+dfZGYDMcYmY6h0alqtShfeuJ17mILDV4He/f38eKpycxj3/GWDj5+bTfrW61UsuxS4wa3k2AkxnAwjstQrGn24nQYRblXaqWtY60iAizYSqm7fs211gIeJ+F4sqoLBfrHQzR5ncSTJsPTUVBgoIgnTXwuB1OROP/rh8eYzJpOfOWmVu64rpsLVgcy29K+YaehcCg4MRoEoL3BxXQ0ycBEhAs6G/ncjVsLfm9rpa1jrSICLNhKKUuMc1lrk+F4pgi/Uk1ZirllTyRNVge8DE1H8LkcxJMmhqHQJjgNxZmpCFORRGb/C1Y3ct35HTx7YpzP/uuviCc1bqfBprYGpiJxnIYllv1jQSujQUMoZrKlM0AolqDF7y5KOGulrWOtIkE4wVaKKUhYjGxrTSnrZ7PPRXuDu2KFAosVYKTRWjMRijEwHuYDV24gYWrC8SROQ5FMahKmJpLQGfF1Oww+d+OFfPTqzfzg6CAD40FmogmiiSTT4TinJoIcHw0CliDHzdlpdbA04Szl5yPMRyxgwVZK2TYwn7U2GY7z2B8vr2NaoRZtIbfsc9tEbu9p4y628M+HT3J2KkIyq4+OUhDwOPnMb27l6vPb+Y/fegGnoZiIJVEoy1rWmmA0idtpMDwdo8XvsXKbTU3S1Jhac+zMFA5DsbmtOOGslbaOtYoIsFA2ChGuUvbVLUdjmGKDULkuAomkyXMnx7nmC4+zJuDlg9tmV58lTc1oKMaJsRCxLPVVgEPB+65Yz9XntwMwOBXO+IsdyprkoQyrJHlds5eBiTChWIJVjW6rv7G2jqGARFIzGoyx79hQ3myMXJ9XOfser3REgIVlk+vEBQoWrlJlIpTDWis2CDX3IpDO43U6FA1uB8Mz57Iarupu5ak3R/nywT5OjJ5zDbT4XLQ1uHAaBuF4kuf7J7kt9be1TT7GQ1HcDoOkabWHNLXG7TBwOgy2dDTS2uBhYDyEx+UgnkiCUrgcBh0BDw5D5Vx7Ie0xhdIjDdmFZTG3cXha9Pwuyw9Z6gbqhaynlNbaznueoMU3uxmN1prJcDznsM7s98PrNHhjaIaEqVnd5M2UDofjSTxOB16XwdHTU4BloXpcBmuaPLgdVtVaMJZgdCZKLKm5dH0LH71mEw0eJ3/x/VeIJZKMBmOpBcGqgBuXwzHLv13M2t/9v/fTNxIkmRLztFiX+/NaQUhDdqH05LMQ+0ZDbOlsnLXvcqPnhbo0SmmtFevW2L21k89rzT/ue5OB8RAmsLrJkxHfaCLJWDBGMHYusLV9cyufuK6Hf/zZm4wGo+CwxHdoyup65nEaDE6G+My/vUjA66Sj0YPfZZBImsRSWRCb2xvnvR/FjHB/bWgGh7KKOhKm5vREhLXNHsl2KDMiwMKySPs8p8JxRmaixJJm6vbYLGn+rV0FAcW4NbS2shbO62zkr95/KQD/8VsvMBqMEk9aI4CyU8resibAnuu6uWJjKwC3XNWV6fEwOmOJL1gdysZCVg5wKJog5nMRT2r+5gOXL/jaix3hrk1QKJQCE83ZqWhmbXZTr9V4IsDCsuhq9dM3MsNoMIaBwqEUsaSJUorJ1CTeYv2xuU42uwoCCg1CTUXiTIbi8yZS3HTZWr74k9dmWbwOQ3HLti4+vnPzLPdAOhvi28/2c2oijM/lYFWjh5GZKAYKZZBx6xTy2osZ4b464OH0ZARMK/NCa01CL81/XmqxrOdqPPEBC8ti37Eh7vzGs5ha4zAUWoPW0N7oosXnpsXvLsofm8+nHIolWNPkLdgXWylmognGg7F5whuJJ/m3507x0JGTBKOW+BoK1jX7+A9vPy+T1ZCN1+Wg1e/G53bMGkN07MxUpmG801D0dDSW9LWnnyuR1Jm7mHTK2mN//PaijpXv81tO7nW5BrJWGPEBC6Vn99ZOAl4noag16DEdwGn0OJkMx/nhH+0q6nj5LN1YorQujeUSjCYYD8WIJWYLb9LUPHb0DA8+fZzRGStI1uB2cMv2Ln73bRvmjSGC2cKbJtt94HYY1l0Fio6AByjta08/l8uh6F7VkBHNz777wqKPVY47lXquxhMBFpbNls5ATgtlKQKR72RzO1TGErazICAUSzAeihOdUx2mtebgG6Pcf7CPk2OWMLgcipsvX8fvb99Es98171gel4OjAxM8+NSJebfr2e6DyVCMhKlpa3DR6HGWfO5aKXN9yyGW9Tz4UwRYWDalzL/Nd7JtWd2U8QXbURAQjCaYCM8XXoAXBia4b38vLw9OA9a95q9ftJqPXruZNU3eefu7HAatDW6e6RvjL39wbMHc2/TrW2p6XaH+2GKyRxY6ZjnEsp6r8cQHLJSEUuXflsOHuBzyuRoA3hye4csH+vhF31hm246eNj6xs5uejsZ5+zsMRYvPTZPP6lVRbt9mOd7LxY5Zrs+v1PndNpDTBywCXKfUctrO3JPt6p42nu4dK/q1LOc9mIkmmMgjvGemIjz45HF+8vLZ9ExMLlob4I5dPVy2oWXe/g5D0exzzRv9U2yRR7GUQ+ALOWYdiGU5qK0gnFLqODANJIGE1nqbvSuqHWo9bWfurfdSXstSHqe1ZjqayJlOBjAZivPNwyf47i9PE0/1bNjY5uf2nd3sPL993oQJQymafC5afLlnrpXbt7mYP3YpF6hCfLxSulw4edtRKqV+oJTaXMG15OIdWuvLRXyLI1dbRpfD6gFQayz1tRT7uJlogoHxMCPT0XniG44n+cahE3z4/l/wyLOniCc17Y1u/tOvX8D9H9nGdVtWzRJflRLeDa0+2hryz1wrdm5csSzUSrLQtpnFHLMS5JudV6ss1A/4QeDHSqk/VUrND+EKVUv/eGheulOp0nYqfQIs9bUU+rhQLMHAeIihqcg84U0kTf79hdPcdv9hHnjyOMFYkgaPgzuu6+brH9/OjZeuxTFHXBs9Tta3+FjV6Fl0xPvurZ3cfdPFZetVvJDAL/XCVu6LxkIs9aJRzeR1QWitv62U+j7w58AzSqmvA2bW3/+2zGvTWBcADdyrtd5b5uerG8p1a2uHayPgcfLG0My8JjGLvZbF3oNI3OrJEMmR1aC1Zv/rI9x/sI+B8TBgpZT97hXruXX7Rpp88+0Rr8tBW4Mbb44834Uo90y6fOllf/bdo0tKF7OzPWU9jkdazAccB4KABwiQJcAV4Fqt9WmlVCfwE6XUMa31/vQflVJ7gD0AGzdurOCyqp9ype1U+gTYd2yI4ZkoCVNjKIgnTQbGw7T6XXzuxosWfGy+9+Bj12xmcDJMODZfeAGePznOfQf6OHbmXEpZq9+NYcCrZ2Y4Njg9q5evy2HQ3uieJfR2kc+nW0jbTCj8Im2Xj7ceCzLyfmuUUjcAfws8CrxNa13RV6m1Pp36OaSU+g6wHdif9fe9wF6wsiAqubZqp1xWSqVPgHv399Lsc9HgdmZKZJ0ORXvD4nPN5r4Ha5t9fOiqDZy/ujGn+L45NMNf//hVXjs7k9nW1eojHE/ic1mtI7MnFF9z/ipaG1wEvNXhnSv27uTOXT38ySMvcGoiTNK0ysgbPc5FL2zZz1fpLJt6LMhY6LL9p8AHtNYvVWoxaZRSDYChtZ5O/f83gLsrvY5aphxWSqlPgMVO4rTgK7fK3Pan07QK4e1v6WBbd1vedDKAwckwX3nyOI+/MpRJKfO6DJq8Ts5MRWj2OjO+ZJ/LwVgoxv/4wcs0+1xsbGuoiPAUInZLuTtRANp6T9Eqd55UnvVUe2e6WiFvlEBrfZ0d4ptiNXBQKfUCcBj4vtb6MZvWIqQoZQCmkIDKUiPupqmZDFmTKIamIjnFdyIU4/8+8QYfeeAIP02Jr9NQrGv20tXio8XnJmnqc13MFIRiSSZCMcLxJK1+d0WCQIUGnooNVt67v5cmn4stqwNcuLaZLasDNPlcBWXK2JVlU+6gpR3Y77jKgda6F7jM7nUIsymla6MQi61YiyeRNJkMx5mOJDDzFBiFY0n+5dl+vnVkICPuHY0eookknQE3hjpnk7gdiljCGhPvNBSnQmGUUngcRkZ45q651LfmhVq2xd6dLMedVMhjy+WiqLcc46oUYKF6KdUJUGhCfyGCH0+aTITizEQtyzwX8aTJ9381yNcPnWA81dw84HXye9s38juXr+O//NtRRoNRspfU5HMxkSrKcBoOogkTBZmOZHPXXI5b80KFstiL1XLcSYs9ttYLgSqJCLBgC4UKwEKCH0uYTIRjBKPJvMJras3PXx3m/if7OD0RAcDtNHjf29Zzy1VdmSBa9jQKn8tB3DTxupx86h2beLp3jIHxEH63gwaPY1bgLXvNxfhhC7UQi3mfirk7WY4/dbHH1mO6WLkQARZsoVgByBas9S0+br2qi0u7WhZ8judOjLP3QG8ms8FQ8O5L1vIHV2+aZcWCNY3ij9UW/uXZAQYnw3RlBdg+nbWGhdZcqLVajIVYzPtUzN3JctxJiz22HtPFyoU04xFso9CmLWnBchpW3m0oliRhau66fsusnNw0r52d5r4DfTx7Yjyzbef5q/jEzm42tue+xW7yuWj1u+dVthWz5kKb3xTbJKfWmtvUyQSLUlNbzXiE+qdQi+0f972JQuM0HGhtTZ0YDUb53KNHuXhtM7dc1cX2njZOTYR54GAfP3t1OPPYSzc0s+e6Hi5a15Tz2D63VcHmcc6vYCumsAEKt1aLtRBrLfBUj+li5UIEWKhapiNxJkJxTowFafJaX9WZaIKh6YglxBp+dWqCo6cnOL8jwOvDMyRN646uZ1UDn7ium1/rbpvXpQzOVbAd7h3LKbJLCSQVeltfjwUF2dhZrlxriAtCqCpM02oJORU+1xIyPdrd53LQPx4iljBJdYPEYHZ9fIvPxf/z9h7eeeHqnO4EQyla/VZT9J+/Opy3efi9+3vLdhtdbU3nhYqQ0wWxcLsmQagQ8aTJ6EyUk2MhRmdmt4S85aouEqZOCdU58YXZ4utQVhVbi2++L1cpRcDroqvNT7PfaoKeq6Aglkjy6Yef5/DxMc5MRpjKqrorVSCpHgsKhKUhLgjBViLxJFNhK4c3H9t72riLLTx0+CSnJ8J591PKci08fKR/VnDO73bS1uDG7Zxtb8z1xU6F44wGY5ha40lNIj49aT1fk881q5fucosMyuHXreUpKCsVEWDBFkKxBBOheM52kHPRWqMMCMaSs6zfuSRMGJ6OZMTc5TBY1eiZNe49m7m+2JGZKABep4OOgIfTExE0mpGZKM7UVOare9pKWmRQKtGU4ofaRARYqBiLjfzJxbEzU9x3oI/nT05ktm1u8zMwEWJuiweXAQlTE4wmeOX0FO+5dG3OABxYgjUejHJ8NIjLMFjd5JlV6RbwuljXAoMTYYKxJAPjYbrb/fzw6JmSFRnMFc3jozPc+Y1nafQ4uCA1BbrQY0rxQ20iArwIclu3fBJJk6lIgulIPJOlsBgD4yHuP3icn792LqXs8q4W7riumwvXNnG4d4y9+9+kd9TyyToVWHEOTavfxTd+cZIbL1uX89jZwrehxcfZ6SgDExE8ToNmn3NWpVtSW37l8zsaCceTHB8NsqFl9qj5hXzDC31/skVzOhJndCaORhOJm0VbsHb2ZygXtbbepSACvAByW7c8IvEkU5H4gqXCcxmdifK1Qyf4/q8GSWv1eR0N7NnVw7ZNrRmLdntPG9t72vitv99POK5JaFBas6rRTUfAu2CwbK612ORzE4oliMYSDAfjDE1H8TgMEqkFrA54Z3X9OjsVpcnnzhwvXwrZYt+fbNEcno6iFBgoYkmzaAu23voz1Np6l4pkQSxAPQ23rCTBaILTE2FOT4SZieRvkJPNTDTB/Qf7uO3+w/z7C5b4rm328qfv2cq9t13JVZvn5/Me7h0jYYLDAI9D4XYYTITjjMxEF8ypzdW6MZE0GZyO0up34XU6iCU1saSm1e+cNYJodcBD3DQJxRJMhWO8fnaa46NBxoPReS0i831/7nnsGLfuPcTwdJQ3hmaYjsSJJU2UAq3BnZolV0zWxWKtQmvtu1xr610qYgEvgNS0F06u/N1CiCVMvvvCab556ARTESt41uxzcduOjfz2ZetwLTDY8l+eHaCtwcVYMAHKyoLAhPFQnC8sUHWVy1o8Ox3FZRh0BLx0BKxtr5+dZjqaZE3WY50Ogy0djSilOD4ayrgx4qaeZ6Hl+v4kkibHR0NsbvezpsnDqYkIA+NhHIqMe2ZVo+XiKKY4o976M9SjSyUXIsALYEfFUq19qdL+3alwPG8P3lwkTc3jr5zlK08d5+xUKvvAZfDBK7v4wLYNNHjyfzU9LgftDW7OTkdY1ejF60owPG2NLHIZCr/HueB7lq9Udq5vd3WTh4GJyLz9PnfjRdy7v5fN7f55hRrZLoOcQj8VneX+AMXZ6QjJpMYwFK1+FwGvc0nN7hdKbau16rt6c6nkQwR4ASpd027Hl2qpgp/O3w3GCvfvgpUJ8fWnT/DwM/1E4palbCj47cvWcduOTbQ1uPM+1mEoWhvcNKWCZOmTNOA9N5stFEvgdhjcuvfQrNcEzHqd73/b+kybyQ2tftypvN9snA6DCzobafG7lzRVOOf3xzTZ0OLL7NPkswR3MhznL26+pGzlu7XWn2GltLwUAV6ASte02zF1uBjB/9krZ/nHn79J/3iINQFfpglOobx8eoov/vg1+kaDmW1+l9Vj9+ru9rziq5SiyetMTSc+5wfOdZJOheNoIJY0M6/pTx55AYUldultXzt0go5GT2YO3LsvWcMjz53KYe1uXXKf3lzfH5ehiM/JBEk/zq4R9dVIvblU8iECvAiV7ERlx9ThQgQ/kTT54YuD/M8fHsNpKAIe56wJwYuJ8MnREF8+2MfBN0Yy2/wuB6sa3XhTYje3ei2Nz+2gvcEzr4oNFha47Nd0aiIMGtY0+1KvR1sTNCIJzu9sZGg6wiPPnZpnFZeiofnc789iPYXLSa11Vasnl0o+RICriEp/qRYT/OyJE/cfPI7TULMmBC8knGClVn3t6RP88Oi5lDKnoehs8tCY9Rq9LoMzU7NLjB2Gor3RQ+MCvmCYf5LuvOeJea8paepZbpKRmSiGgqTWs2a7Pd07VnCjnaValHZYorUWVyiEWnOp5EMEuIqo9Jcqn+CvbfZxZtIKPqUZnApnWkKmySWcADORBA8dOcm/PncqM5G4vcFNo8fJqYmQVfLbQEZcI3GTNU2WdXq4d4wvH3yTk2PWAMyeVQ185obcboBCX5PDUKDPuS5iSavizZ2VYbGUO42lWpRLfdxShLReglVzqTWXSj4kD7iKqHSXrOzcUdM0mYnECceS/O4V62eJL8DaJl8maJYmWzjBspi/daSf37//Fzx0uJ9YwqTV7+Lmy9bhcihMrekMeEgmNWenIsxE44Tj1nSLW67q4vmT4/z1j49xfNQSda01rw/N8J8feaHg0e+58mEbPc5MZoHWGoehMPXs4ZrVfvta6Hj6udRzPu3urZ08tGcHBz5zPQ/t2VFz4gtiAVcdpfTT5bOYsrc3uAw0irFgjNVN+QNrt1zVxT0/OsbZqQhJ0xKxBo+TT+4+n6Sp+cnLZ3nwqeMMTVspZT6Xg1uu6uL9V27gT79zFJfDSLkvHKhmxfB0lOGZGBevbeb3tnfx7kvXsudrzxKMJXEolQm2Ka2ZjhQeiMxlGX3uxosga9vmNj+jwRgOQ6G1LuhOw+7b+KUGaOslWFWviADXKfluPd8/MMEjz53CaQBa0zcWIpHUbG5vKCyrIV3woCwL9eXBKf6//W9yPN2TwVDcdNk6PrxjIy1+K6thrvuiwe3E3+5gOpLgyx/dRlsqu6F/PGSJe1bFm1JWELAYwch1EUtbihpobfDwnreuLTjgVg238UsV0moIVtl98apmRIDrlHwW030HemlrcBOJa4amIygUDmU1v8mV1XC4d4yHj/Tz0uAkSsGqBiswFo4lOTtt9W0Aqw3OOy/s5GPXbmZts2/WWtY2+TITLdLEEiab2htY1XjODdDV6mdkJkoiqTG1idbWgV2GWpZg5BLQR547VbB7pxpyTpcqpHYHq6rh4lXNiA+4Tsnud6C1zliWwVgSt9NgPBRDYd3qK8PyzzoNxcNH+jPHONw7xpeeeJ3RYBTT1GjT8t2eHAvRPxHOFC5s39zKvbddyX99z4XzxBdmT7QAa/qFRvEf3n7erP3u3NWDy1AkTI2pLWtVa6vF5NVF5BuDdeLfuvcQO+95gk8//DzxZHLJftBcvSMqfRu/WK+HfNg9faOefdClQCzgOqWr1c/ZqTAep4OkttQsHE/iczmIxE3iSRMjdauvtdW8fG5Ww8NH+jOpZy6HIprU1jDMVGaD26HoavXzhfdduuBa0hMtHnlugDOTYbraGnLehu7e2smGVj9vDE1nev16nAYtfhdP947x6QJf+1yr68xkhGA0wXgwTlJr3A6DVY3uggW0Gm7jlxP1tzP/V3zQCyMCXIdEE0k+tM0KmsWTGq/LIBI3SZiaD165gcdePothKLSprUowTcotMTurYXAqTIPbwfBMlEhidvVWe4Mbl0Nxx3WL38p6XQ7ee+V6bvm1jYvuOx1N8JY1TbM6n2mtizph57oMnIYiktCYSROP02ozeWoiwpbOxoKOZ/dtfJpaK6SA6rh4VTPigqgj0m0gT42HubSrmbuu30J7g4fpSIL2Bg93Xb+F267ZzF3Xb2FDs4+kJjUBwo2hVCYdDCAaT2Kg6BsNMR6yBlMayhp86TSsE+uP3nnBgkE7h6HoCHhY1+LD48w9FihN2mWQ3aIxTbEn7FyXQboIQ2tSfo3Z2xcj1238+9+2nnv397Lznie4de+hgtPkVhpLdZ2sFGQsfY2z1DaQcC7AdmYqzJpUCtqVm1v50UtnePCp44zMxABLpNsaXHhdDkwNd11/LlCXPsbgVJi1WWlsjV4n7Q2enKPh55LtMkgkTU5NRABY3+LF6TCKHtl+695Ds6yuY2emLB821kUh7YIwNRz4zPVFvWdz11vqsfL1mDGQfk21XDBRAnKeCCLANUo8aTIVjjMdSRTVBjIfWmsOvjHK/Qf7ODlm3e67HIrtm9uYDMUZCUYzIp0tvvf86BjBaCKTG9zocfK/fvet/OYlawt+7rmCORWOc3Y6gtbwto2tOU/YhYRqrkC+MTRDwtSsb/FlmquHYgk6A96CS48XWu9yj5f9msol7ILt5BRg8QHXGIWMcc9FPksV4IWBCe7b38fLg1OA9U1510Wr+dg1m1nT7M17zL0HepkKxzEMhcOhQMNkOM7//unrRQnw3EBNdovGXIK2WGrT3IBV96oGhlOTjQstvChmvVCawFI1pLsJlUUEuEaYiSaYDMeJFjDGfS7pdDKnYbV1THcy+9BkF4eOj3Kodyyz746eNm7f2c15HYsHqPrHQym/cKpqzbBGUvSOBBd+4ByKDdQUIlS5upCV6ja4XIElyRhYeUgQrooxTc1kKM7J0RBDU5EliS/MTidTKJyGsizVJ17PiO9FawP83Ycu43++960FiS+cu6dSSuUd/14IxQZqlpKXm+4b8Bc3XwLAn3336JKDZ+UKLHW1+jO50mkkY6C+EQGuQuJJk5GZKCfHQowGoyTM4oJrcxmcCuN1GSRNzfB0lOOjIUIx60TvavXx32+6mP9z6xVctqGloOMZymoV2bOqAY2yAlxaY6YKKLrbixOMYosFlipUS21os9z1FopkDKw8JAhXRYRj6THuxfl3F+Ouh37JybEg09FEpi+vQ8G6Zh8PfOyqgjIV0mRnN+w7NsSfPPICM3OCcH/z/styilGhEf7F9ltqsKpcwbNSIhkDdYsE4aoRrXXGv5vunVsqEkmTHx49w/FRS3zByuUNeJ34XA4++Y7zCxZfl8OgI+DBm3Xrv3trJ3/z/ssKEoxCewIUst9Sq8Jqwcdai8UWwtKpWgFWSt0AfAlwAF/WWn/B5iWVFNO02ixORYrP310MrTX7Xx/h/oN9DIxbpcVOQ9Hqd6OUZl2zv6h5bi1+N61+V04/b6GCUWiEv9D9liJUUpVVfuoxj7mcVKUAK6UcwD8Avw4MAEeUUo9qrV+2d2XLJ2lqJsNxpiNxkmbp3T+/7J9g7/5ejp2ZBiyL9zcuWsNHr9lEZ1P+lLJcuJ0GqxpnW71LJW19ToXjjMxYI+TdDoPJUCznftmUykqtlpLiekU6nxVPVQowsB14Q2vdC6CUehi4GahZAV7qGPdCeXNohvsO9HL4+Hhm2zXntXP7zm66VzUUdSylrBOoJY/VuxS6Wv30jcwwGoxhoHAoRSxp9afYd2woc4KW00qtlzE21YrkMRdPtQrweqA/6/cB4Neyd1BK7QH2AGzcuHiTFzsop383zeBkmK88eZzHXxnKjFi/eF0Te67r4a0bmos+ntfloL3RvWjvhmK5c1cPd37jWQCUYfVlUCjaGlyzTtByW6niYy0fteBjrzaqVYBzmV2zzEat9V5gL1hZEJVYVKEkTZ0pE15uClk+JkIxvnHoJI++cJpEypWxqd3PJ3Z2c8157UVbrk7DoLXBRcDrWnznJbB7a6c1ly2aIG5aLSE7AlZz9+wTVKzU2kV87MVTrQI8AHRl/b4BOG3TWgommkgyGY4TjJbHzQBWqtojzw7wrWf6M7m8HY0ePnrtZn7jotVFpZSlafK5MmOBlkohwZctnYGcaWBzT1CxUmsT8bEXT7UK8BFgi1KqGzgF3AL8nr1Lyk8wamUzhGNLq1QrhETS5PsvDvK1p09k2kMGvE5+b/tGfufydXgWCZTl6gWx84JVJQmyFRp8kRO0vpG7l+KpSgHWWieUUp8CfoSVhvaA1volm5c1i+W0gSzqebTm568Oc/+TfZxOtWl0Ow3e97b13HrVRhq9i3+EuXpB/J+fvcHqJg8bLlz+7WExqWNygtY3cvdSHFUpwABa6x8AP7B7HXNJJE2mIomypZFl89yJcfYe6OW1szOAlVJ2wyVr+MjVm+kIeBZ59Dlm9YJQiiafg0g8yd4DfbzjwtXLXmcxwRc5QQXhHFUrwNVGudPIsnnt7DT3Hejj2RPnUsquPb+dT+zsZlN7cSllcG4svMNQOAyrcU4po9MSfBGEpSECvAjBVBpZZImdyIrh1ESYBw728bNXhzPb3rq+mT27url4XfEpZWnWt/iYCMXwus71XiqlQOby7U6G47gdBjvveUIqogQhDyLAOaiUfzfNWDDG1w+d4Hu/Gsy4NXpWNfCJ67r5te62JRdDOAxFa4ObT1+/pey5tdm+3Qa3AwXEkqZURAnCAkg3tCxKPeZnMUKxBN8+MsC3n+0nEreEvjPg4ePXbuadFy4tpSzN3JlsleyyVQtdxwShwkg3tHxE4un83dK2gcxHPGny7y8M8o1DJ5gIWyllTV4nv79jEzdftg63c+ltml0Oq3+Dzz07taySwS+piBKEwlixAlyJMuG5mFrzs2NDPPDkcQYnrZQyj9Pg/Vdu4ENXddHoWfrHUY7+DUtFgnK5kU5hwlxWrABHEybD09GKPJfWmmdOjHPf/j7eGD6XUnbjW9fyB1dvor2x8JSyXPjcDlY1enA5qmPASSEFF4WIUT0JlnQKE3KxYn3AkXiS0xPhMq3oHMfOTHHfgT6ePzmR2bbrglXcfm03XW3LswidhkFbo3tZlnO5WMjnXMhEi3ob0S5+8RWP+IArycB4iPsPHufnr51LKbu8q5k7ruvhwrVNyz5+Kfo3lJNcPue0KD93chwFrGn2opTKWTm33NaG1WY9i19cyIUIcIkZC8b42tMn+P6L51LKzutoYM+uHrZtal22f7aUTdIrSbZFa2qNAk5PRFjXAgGva54YLUewqvF2X/ziQi5EgEvETDTBt47086/PDhBJBfXWNnv52LWbuX5rJ8YyhddQVk5vs6887SLLTbZF63YYJJIaFAxPRwl4XfPEKJdgjQajBKPJRYs7qrExuDQiEnIhArxMYgmTR184zTcOnWAqYqWxNftc3LZjI7992bqSBMYaPU7aGtw4qyTIthSyLdpVjR5OT4ZR2gqG5hq/PlewRoNRhqZjdDS6F7Vqq/F2XxoRCbkQAV4iSVPz+LEhvvJkH2enrGwKr8vgg1d28YFtG2goQWDM5TBob3TPsgKrkUL8rdkWbVNKHM9OR1Ba0RnwznvMXMEKRpN0NLrpCFhz7Rayaqv1dl8aEQlzqe4zuwrRWvOLvjG+fLCP3uEgYJX8/tZb13Lb1Ztoa3CX5HmafS5aqzjIlmapvYCdDkt4F8pqyBasnfc8UbBVK7f7Qq0gAlwErwxOsXd/Ly8MTGa2veMtHXz82m7Wt/pK8hy1FmSrVC/gYqxaud0XagUR4AI4ORri/if7OPD6SGbblRtbuGNXDxesDpTseZp9Ltoa3LZXshVDpXoBF2vVyu2+UAuIAC/AyEyUrz51gh8eHSTde31LZyN3XNfNts1tJXseV2pAZa1YvdlUyt8qVq1Qj4gA52AmkuChIyf5t+dOEU2llK1r8XL7td28/S0dy04pyybgddHeUP2+3nxUyt9abYUVglAKRICziCVMvvP8Kf758EmmUyllrX4Xt+3YxI2Xri1prwWnYbAqUP0ZDotRCcu0GgsrBKEU1PbZXyKSpuYnL5/lwaeOM5Rq0ONzOfjQVRv4wJVd81o7LpdqLyMulnL7W+0urBDrWygXK1qAtdYc6h3jvgO9HB+1gkZOQ3HTZev4/R0bafWXJqUsTa1lOFQLdhZWiPUtlJMVK8DPnRznf/3gFV48NZXZ9q4LO/noNZtZ11KalLI0hlK0+t00+4svIxbry97CCrutb6G+WZEC/A8/e4O//tGrmd+v2tzKJ3Z2s6WEKWVpllNGXO/WV6EXFzsLK6qxrFmoH2q3ucAyeOeFnSgFb1kT4G8+cCn3vO/Skouvy2GwptlLZ5N3yT0csq2vdNtGl0Nx7/7ekq7VDtIXl6HpyKyLy75jQ/P23b21k7tvupjOgJfJcHzRCrpS0tXqJzxnInY1lDUL9cGKtIC3rmniW3t20F6GogelVKqMePmjgerZ+ir21t6uwgopaxbKyYq0gAEu3dBScvH1uhysa/GWrJqtnq2v/vEQvjnByGq8uNhpfQv1z4q0gEuNw1C0NbgJeEvbq7eera9q7ViWCylrFsrFirWAS0XA62JDq7/k4gv1bX3duauHeFITiiXQWufsCSwI9Y5YwEukUjm9y7W+qjWNTXo7CIJMRS76uWppNFC9TRYWhBomZ1BIXBBF0Oh10tXmrwnxhfpOYxOEekBcEAXgcljuhlL3hCg39ZzGJgj1gFjAC6BSJcQbWn01J75Q32lsglAPiADnweNysL7FR2uNTajIRjINBKG6ERfEHJRStC2xcU61IZkGglDdiABn4XM7WNXoKWnjdbuRIoLSU62pfULtUXVKo5T6vFLqlFLql6l/7yn3cxpKsSrgYW2zr67EVyg9xTQREoTFqFYL+O+01n9TiSfyu52salxau8jFEEup/pD+wEIpWbHmnqEUHQEPa5qX3i5yIcRSqk9qpYmQUBtUqwB/Sin1K6XUA0qp1lw7KKX2KKWeUUo9Mzw8XPQTuJ1GWfo3pJEiiPpEUvuEUmKLACulfqqUOprj383APwHnAZcDg8AXcx1Da71Xa71Na72to6OjcosvELGU6hNJ7RNKiS0+YK31uwrZTyl1H/C9Mi+nLNRSu0WhcCS1TyglVReEU0qt1VoPpn59L3DUzvUslXru5bvSkdQ+oVRUnQADf6WUuhzQwHHgTltXs0TEUhIEYTFWbDtKQRCECiLtKAVBEKoJEWBBEASbEAEWBEGwCRFgQRAEmxABFgRBsAkRYEEQBJsQARYEQbAJEWBBEASbEAEWBEGwiWosRRZKjDSGF4TqRCzgOkcawwtC9SICXOdIY3hBqF5EgOscaQwvCNWLCHCdIyN0BKF6EQGuc2SEjiBULyLAdc7urZ3cfdPFdAa8TIbjdAa83H3TxZIFIQhVgKShrQBkhI4gVCdiAQuCINiEWMCCgBSrCPYgFrCw4pFiFcEuRICFFY8Uqwh2IQIsrHikWEWwCxFgYcUjxSqCXYgACyseKVYR7EIEWFjxSLGKYBeShiYISLGKYA9iAQuCINiECLAgCIJNiAtCkCowQbAJsYBXOFIFJgj2IQK8wpEqMEGwDxHgFY5UgQmCfYgAr3CkCkwQ7EMEeIUjVWCCYB8iwCscqQITBPuQNDRBqsAEwSZssYCVUh9QSr2klDKVUtvm/O2/KKXeUEq9qpT6TTvWJwiCUAnssoCPAr8L3Ju9USl1EXALcDGwDvipUuoCrXVy/iEEQRBqG1ssYK31K1rrV3P86WbgYa11VGvdB7wBbK/s6gRBECpDtQXh1gP9Wb8PpLbNQym1Ryn1jFLqmeHh4YosThAEoZSUzQWhlPopsCbHn/5Ua/3dfA/LsU3n2lFrvRfYC7Bt27ac++RDeh8IglANlE2AtdbvWsLDBoCurN83AKdLsyKLdO8Dl0PN6n1wN4gIC4JQUarNBfEocItSyqOU6ga2AIdL+QTS+0AQhGrBrjS09yqlBoCrge8rpX4EoLV+Cfg28DLwGPDJUmdASO8DQRCqBVvS0LTW3wG+k+dvfwn8Zbmeu6vVz9B0BL/73EuX3geCINhBtbkgyo70PhAEoVpYcQIsvQ8EQagWVmQvCOl9IAhCNbDiLGBBEIRqQQRYEATBJkSABUEQbEIEWBAEwSZEgAVBEGxCBFgQBMEmRIAFQRBsQgRYEATBJkSABUEQbEJpXVQv86pEKTUMnLB7HTawChixexE2s9LfA3n9tfH6R7TWN8zdWBcCvFJRSj2jtd62+J71y0p/D+T11/brFxeEIAiCTYgAC4Ig2IQIcG2z1+4FVAEr/T2Q11/DiA9YEATBJsQCFgRBsAkRYEEQBJsQAa5xlFKfV0qdUkr9MvXvPXavqRIopW5QSr2qlHpDKfVZu9djB0qp40qpF1Of+zN2r6fcKKUeUEoNKaWOZm1rU0r9RCn1eupnq51rLBYR4Prg77TWl6f+/cDuxZQbpZQD+Afg3cBFwK1KqYvsXZVtvCP1uddsLmwRPAjMLWb4LPC41noL8Hjq95pBBFioRbYDb2ite7XWMeBh4Gab1ySUGa31fmBszuabga+m/v9V4HcquablIgJcH3xKKfWr1C1aTd2CLZH1QH/W7wOpbSsNDfxYKfWsUmqP3YuxidVa60GA1M+amrYrAlwDKKV+qpQ6muPfzcA/AecBlwODwBftXGuFUDm2rcR8ymu11m/DcsV8Uim1y+4FCcWxIsfS1xpa63cVsp9S6j7ge2VeTjUwAHRl/b4BOG3TWmxDa3069XNIKfUdLNfMfntXVXHOKqXWaq0HlVJrgSG7F1QMYgHXOKkvXZr3Akfz7VtHHAG2KKW6lVJu4BbgUZvXVFGUUg1KqUD6/8BvsDI++7k8Cnwk9f+PAN+1cS1FIxZw7fNXSqnLsW7BjwN32rqaCqC1TiilPgX8CHAAD2itX7J5WZVmNfAdpRRY5/E/a60fs3dJ5UUp9RCwG1illBoA/hvwBeDbSqnbgZPAB+xbYfFIKbIgCIJNiAtCEATBJkSABUEQbEIEWBAEwSZEgAVBEGxCBFgQBMEmRIAFAVBKdSml+pRSbanfW1O/b7J7bUL9IgIsCIDWuh+rrPsLqU1fAPZqrU/Ytyqh3pE8YEFIoZRyAc8CDwB3AFekuq0JQlmQSjhBSKG1jiul/jPwGPAbIr5CuREXhCDM5t1YXeUusXshQv0jAiwIKVI9NX4d2AH88ZxGR4JQckSABQFQVlebfwL+SGt9Evhr4G/sXZVQ74gAC4LFHcBJrfVPUr//I7BVKfV2G9ck1DmSBSEIgmATYgELgiDYhAiwIAiCTYgAC4Ig2IQIsCAIgk2IAAuCINiECLAgCIJNiAALgiDYxP8PpHTvpvTDi+oAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sns.lmplot(x=\"X\", y=\"Y\", data=fakeData1)" ] }, { "cell_type": "markdown", "id": "ab49b965", "metadata": {}, "source": [ "We estimate a strong relationship between X and Y, but we know in fact that X and Y do not depend on each other in reality. Instead it is the \"lurking\" variable Z that causes the relationship. \n", "\n", "But if we *control* for Z in our regression:" ] }, { "cell_type": "code", "execution_count": 17, "id": "170995eb", "metadata": {}, "outputs": [], "source": [ "as_mod3 = smf.ols(formula='Y~X + Z', data=fakeData1).fit()" ] }, { "cell_type": "code", "execution_count": 18, "id": "6589e3b2", "metadata": {}, "outputs": [ { "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", "
OLS Regression Results
Dep. Variable: Y R-squared: 0.725
Model: OLS Adj. R-squared: 0.722
Method: Least Squares F-statistic: 259.4
Date: Thu, 18 Aug 2022 Prob (F-statistic): 6.43e-56
Time: 13:42:37 Log-Likelihood: -486.89
No. Observations: 200 AIC: 979.8
Df Residuals: 197 BIC: 989.7
Df Model: 2
Covariance Type: nonrobust
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [0.025 0.975]
Intercept 5.1078 0.258 19.780 0.000 4.599 5.617
X 0.0191 0.066 0.287 0.774 -0.112 0.150
Z 0.4624 0.029 15.727 0.000 0.404 0.520
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 2.440 Durbin-Watson: 2.126
Prob(Omnibus): 0.295 Jarque-Bera (JB): 2.054
Skew: -0.185 Prob(JB): 0.358
Kurtosis: 3.330 Cond. No. 13.7


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified." ], "text/plain": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: Y R-squared: 0.725\n", "Model: OLS Adj. R-squared: 0.722\n", "Method: Least Squares F-statistic: 259.4\n", "Date: Thu, 18 Aug 2022 Prob (F-statistic): 6.43e-56\n", "Time: 13:42:37 Log-Likelihood: -486.89\n", "No. Observations: 200 AIC: 979.8\n", "Df Residuals: 197 BIC: 989.7\n", "Df Model: 2 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept 5.1078 0.258 19.780 0.000 4.599 5.617\n", "X 0.0191 0.066 0.287 0.774 -0.112 0.150\n", "Z 0.4624 0.029 15.727 0.000 0.404 0.520\n", "==============================================================================\n", "Omnibus: 2.440 Durbin-Watson: 2.126\n", "Prob(Omnibus): 0.295 Jarque-Bera (JB): 2.054\n", "Skew: -0.185 Prob(JB): 0.358\n", "Kurtosis: 3.330 Cond. No. 13.7\n", "==============================================================================\n", "\n", "Notes:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "\"\"\"" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "as_mod3.summary()" ] }, { "cell_type": "markdown", "id": "3a83c772", "metadata": {}, "source": [ "Now the relationship between X and Y is estimated closer to 0 and is not estimated to be significantly different that 0. " ] }, { "cell_type": "markdown", "id": "af6a86c2", "metadata": {}, "source": [ "### 3.) Monte Carlo Experiment\n", " \n", "Let us say we have data on 10 recent graduates from a Norwegian business school who have gone to work at a consulting company. The company has a first-year associate review where associates are graded on a scale from 1 to 10 on a series of evaluations. The company also has the business school GPA (grade point average) from 1 to 5. \n", "\n", "The data for these 10 students: GPA (X) and first-year evaluation (Y) is generated below: " ] }, { "cell_type": "code", "execution_count": 17, "id": "8645d063", "metadata": {}, "outputs": [], "source": [ "#initial \n", "n=10\n", "a= 5\n", "\n", "sigma = 1.5\n", "\n", "b=.5\n", "x = spt.uniform.rvs(1,4,n,random_state=4321)\n", "\n", "epsilon = spt.norm.rvs(0,sigma,n,random_state=1234)\n", "\n", "y=a+b*x + epsilon\n", "\n", "gpaFakeData = pd.DataFrame({\"x\":x, \"y\":y})" ] }, { "cell_type": "code", "execution_count": 18, "id": "554ccbc8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x [1.2832115 4.26025604 4.07161984 2.14541801 1.77237725 4.91564884\n", " 2.62491485 4.03107147 1.35660706 2.23953387]\n", "y [6.3487585 5.34366448 9.18487037 5.60373116 4.80530552 8.78856883\n", " 7.60184005 6.06075048 5.70184809 2.7557395 ]\n" ] } ], "source": [ "print(\"x\", x)\n", "print(\"y\", y)" ] }, { "cell_type": "code", "execution_count": 19, "id": "b3c8865f", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/johannesmauritzen/opt/anaconda3/lib/python3.9/site-packages/scipy/stats/stats.py:1541: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=10\n", " warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n" ] }, { "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", "
OLS Regression Results
Dep. Variable: y R-squared: 0.269
Model: OLS Adj. R-squared: 0.178
Method: Least Squares F-statistic: 2.944
Date: Mon, 10 Oct 2022 Prob (F-statistic): 0.125
Time: 11:33:39 Log-Likelihood: -18.556
No. Observations: 10 AIC: 41.11
Df Residuals: 8 BIC: 41.72
Df Model: 1
Covariance Type: nonrobust
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [0.025 0.975]
Intercept 4.0833 1.360 3.002 0.017 0.947 7.219
x 0.7443 0.434 1.716 0.125 -0.256 1.745
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 0.874 Durbin-Watson: 2.356
Prob(Omnibus): 0.646 Jarque-Bera (JB): 0.725
Skew: -0.507 Prob(JB): 0.696
Kurtosis: 2.156 Cond. No. 8.47


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified." ], "text/plain": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: y R-squared: 0.269\n", "Model: OLS Adj. R-squared: 0.178\n", "Method: Least Squares F-statistic: 2.944\n", "Date: Mon, 10 Oct 2022 Prob (F-statistic): 0.125\n", "Time: 11:33:39 Log-Likelihood: -18.556\n", "No. Observations: 10 AIC: 41.11\n", "Df Residuals: 8 BIC: 41.72\n", "Df Model: 1 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept 4.0833 1.360 3.002 0.017 0.947 7.219\n", "x 0.7443 0.434 1.716 0.125 -0.256 1.745\n", "==============================================================================\n", "Omnibus: 0.874 Durbin-Watson: 2.356\n", "Prob(Omnibus): 0.646 Jarque-Bera (JB): 0.725\n", "Skew: -0.507 Prob(JB): 0.696\n", "Kurtosis: 2.156 Cond. No. 8.47\n", "==============================================================================\n", "\n", "Notes:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "\"\"\"" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "as_mod4 = smf.ols(formula='y~x', data=gpaFakeData).fit()\n", "as_mod4.summary()" ] }, { "cell_type": "code", "execution_count": 21, "id": "3bcb378b", "metadata": {}, "outputs": [], "source": [ "nsim = 1000\n", "\n", "b_hats = []\n", "cover_95 = []\n", "\n", "for s in range(nsim): \n", " epsilon = spt.norm.rvs(0,sigma,n)\n", " # the epsilon part of the fake data is now within the loop, \n", " #so we create a new dataset in each loop. \n", " y=a+b*x + epsilon\n", " fakeDF = pd.DataFrame({\"x\":x, \"y\":y})\n", " as_mod4 = smf.ols(formula = \"y~x\", data=fakeDF).fit()\n", " b_hat = as_mod4.params[1]\n", " b_hats.append(b_hat) # save a list of all our estimated b's\n", " CI_b = as_mod4.conf_int().iloc[1, :] #We take a short cut of just getting the CI from the model object\n", " cover_95.append(int((b>CI_b[0]) & (b" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAASPUlEQVR4nO3dfYxcV3nH8e/jYMvEBjmu8eLy5qJaCQGJQEYkKAKNCaEhRbUrNSikLysUadWqRSC1tG6RqPpXTSUq0YiqWgHqohJMRJPa4q2kbkdRJZoS00AIcWSIQkjj2sqLgbVr7NhP/9ibsrFndu7uztzZM/v9SNbMvXNn7nP2zPx89uy9cyMzkSSVZ82oC5AkLY0BLkmFMsAlqVAGuCQVygCXpEK9qMmdbdmyJbdv397kLjl58iQbNmxodJ9NGvf2wSpo4yOPcO7cOS658spRVzJUY9+PDK+Nhw4deiozX3bh+r4BHhGXA1+Yt+q1wEeBz1brtwOPAe/NzGcXeq3t27dz//331696ADqdDu12u9F9Nmnc2weroI3tNidOnGBTw5+Npo19PzK8NkbED7ut7zuFkpmPZOZVmXkVcDVwCrgb2AMczMwdwMFqWZLUkMXOgV8P/CAzfwjsAmaq9TPA7gHWJUnqY7EBfgvw+er+RGYeBahutw6yMEnSwqLuqfQRsQ54Enh9Zh6LiBOZuWne489m5mVdnjcFTAFMTExcvW/fvoEUXtfs7CwbN25sdJ9NGvf2wfi38aoPfYhz587x4O23j7qUoRr3foThtXHnzp2HMrN14frFHIXybuBbmXmsWj4WEdsy82hEbAOOd3tSZk4D0wCtViub/iPGuP/hZNzbB6ugjZs2ceLEifFuI6ugH2m+jYuZQnkfP58+ATgATFb3J4H9gypKktRfrQCPiEuBG4C75q3eC9wQEUeqx/YOvjxJUi+1plAy8xTwCxese5q5o1IkSSPgqfSSVKhGT6WXxskd9z3edf2t17y64Uq0WjkCl6RCGeCSVCgDXJIKZYBLUqEMcEkqlAEuSYUywCWpUAa4JBXKAJekQhngklQoA1ySCmWAS1KhDHBJKpQBLkmFMsAlqVAGuCQVygCXpEJ5RR6pIT2v4NNwHRofjsAlqVAGuCQVqlaAR8SmiPhiRByOiIcj4q0RsTki7omII9XtZcMuVpL0c3VH4J8AvpaZVwBvBB4G9gAHM3MHcLBaliQ1pG+AR8RLgbcDnwbIzDOZeQLYBcxUm80Au4dToiSpm8jMhTeIuAqYBr7H3Oj7EPBB4L8zc9O87Z7NzIumUSJiCpgCmJiYuHrfvn2Dqr2W2dlZNm7c2Og+mzTu7YPRtvGZk2cW/ZzNG9Yt6rXe8ZE/5ty5czx4++2L3ldJfK8u3c6dOw9lZuvC9XUOI3wR8GbgA5l5X0R8gkVMl2TmNHP/AdBqtbLdbtd96kB0Oh2a3meTxr19MNo29jr0byHta169qNfatGkTJ06csB/HQNNtrDMH/gTwRGbeVy1/kblAPxYR2wCq2+PDKVGS1E3fEXhm/k9E/CgiLs/MR4DrmZtO+R4wCeytbvcPtVJpTB37yc9Ycy4vGqHf2mMkLz2v7pmYHwA+FxHrgEeB9zM3er8zIm4DHgduHk6JkqRuagV4Zj4AXDSBztxoXJI0Ap6JKUmFMsAlqVAGuCQVygCXpEIZ4JJUKANckgplgEtSoQxwSSqUAS5JhTLAJalQBrgkFcoAl6RCGeCSVCgDXJIKZYBLUqEMcEkqlAEuSYUywCWpUAa4JBXKAJekQhngklSoWlelj4jHgJ8C54DnMrMVEZuBLwDbgceA92bms8MpUyrHHfc9PuoStEosZgS+MzOvysxWtbwHOJiZO4CD1bIkqSHLmULZBcxU92eA3cuuRpJUW90AT+DrEXEoIqaqdROZeRSgut06jAIlSd3VmgMHrsvMJyNiK3BPRByuu4Mq8KcAJiYm6HQ6i69yGWZnZxvfZ5PGvX0w2jauP3lm6PtYc/YUkedZ/9QLP1adzqNdt3+mR02bN6wbeG2D5Ht18GoFeGY+Wd0ej4i7gbcAxyJiW2YejYhtwPEez50GpgFarVa22+2BFF5Xp9Oh6X02adzbB6NtYxN/kDy/9lLWnD3F6S1XvGB9+5pXL6qmXtuvFL5XB6/vFEpEbIiIlzx/H3gX8F3gADBZbTYJ7B9WkZKki9UZgU8Ad0fE89vfkZlfi4hvAndGxG3A48DNwytTknShvgGemY8Cb+yy/mng+mEUJUnqzzMxJalQBrgkFcoAl6RCGeCSVCgDXJIKZYBLUqEMcEkqVN3vQpG0wvU6xf7WFX6KvZbOEbgkFcoAl6RCGeCSVCgDXJIKZYBLUqE8CkVaoYZ9MQmPWimfI3BJKpQBLkmFcgpFq8pqnDZo4rqeGg1H4JJUKANckgplgEtSoQxwSSqUAS5Jhaod4BFxSUT8V0R8qVreHBH3RMSR6vay4ZUpSbrQYkbgHwQenre8BziYmTuAg9WyJKkhtQI8Il4J/CrwqXmrdwEz1f0ZYPdAK5MkLSgys/9GEV8E/hJ4CfBHmfmeiDiRmZvmbfNsZl40jRIRU8AUwMTExNX79u0bVO21zM7OsnHjxkb32aRxbx8Mto3PnDzTdf3mDesWtf0gvW3Ph4k8z70f+/jQ97UcvX5GdfleXbqdO3ceyszWhev7nokZEe8BjmfmoYhoL3bHmTkNTAO0Wq1stxf9EsvS6XRoep9NGvf2wWDb2OusxHaPMzGbOIvx/NpLWXP2FKe3XDH0fS1Hr59RXb5XB6/OqfTXAb8WETcB64GXRsQ/AMciYltmHo2IbcDxYRYqSXqhvnPgmfmnmfnKzNwO3AL8a2b+FnAAmKw2mwT2D61KSdJFlnMc+F7ghog4AtxQLUuSGrKobyPMzA7Qqe4/DVw/+JIkSXV4JqYkFcoAl6RCGeCSVCivyKOx5FVomrPQz3qcr3S0EjgCl6RCGeCSVCinUCScclGZHIFLUqEMcEkqlAEuSYUywCWpUAa4JBXKAJekQhngklQoA1ySCmWAS1KhDHBJKpQBLkmFMsAlqVAGuCQVygCXpEL1DfCIWB8R/xkR346IhyLiL6r1myPinog4Ut1eNvxyJUnPqzMC/xnwjsx8I3AVcGNEXAvsAQ5m5g7gYLUsSWpI3wDPObPV4trqXwK7gJlq/QywexgFSpK6i8zsv1HEJcAh4JeBT2bmn0TEiczcNG+bZzPzommUiJgCpgAmJiau3rdv36Bqr2V2dpaNGzc2us8mjXv7YK6NZ2Jd18c2b+i+/pmTZ4ZZ0kC9bc+HiTzPvR/7+KhLWdBSftbzn7Na3qvDaOPOnTsPZWbrwvW1LqmWmeeAqyJiE3B3RLyh7o4zcxqYBmi1Wtlut+s+dSA6nQ5N77NJ494+mGvjky9+bdfH2j2uel7SJdLOr72UNWdPcXrLFaMuZUFL+VnPf85qea822cZFHYWSmSeADnAjcCwitgFUt8cHXZwkqbc6R6G8rBp5ExEvBt4JHAYOAJPVZpPA/iHVKEnqos4UyjZgppoHXwPcmZlfiohvAHdGxG3A48DNQ6xTknSBvgGemd8B3tRl/dPA9cMoSpLUn2diSlKhah2FIq1UJR1tshrN75/1J8/8//KtPY5o0eI4ApekQhngklQoA1ySCmWAS1KhDHBJKpQBLkmFMsAlqVAGuCQVygCXpEIZ4JJUKANckgrld6FIqsXvnVl5HIFLUqEMcEkqlFMokhrXazrGr5ldHEfgklQoA1ySCuUUilaUbr9arz95Bl48gmKkFc4RuCQVygCXpEL1DfCIeFVE/FtEPBwRD0XEB6v1myPinog4Ut1eNvxyJUnPqzMCfw74w8x8HXAt8PsRcSWwBziYmTuAg9WyJKkhfQM8M49m5req+z8FHgZeAewCZqrNZoDdQ6pRktTFoo5CiYjtwJuA+4CJzDwKcyEfEVt7PGcKmAKYmJig0+ksp95Fm52dbXyfTRq39q0/eeaidWueO836pw6PoJpmrDl7isjzY91GqNePd325++ObN6wbRkkD1/TnsXaAR8RG4B+BD2XmTyKi1vMycxqYBmi1Wtlut5dQ5tJ1Oh2a3meTxq19XQ8jfOowp7dcMYJqmnF+7aWsOXtqrNsIy+vHdiFnaDb9eax1FEpErGUuvD+XmXdVq49FxLbq8W3A8eGUKEnqps5RKAF8Gng4M/963kMHgMnq/iSwf/DlSZJ6qTOFch3w28CDEfFAte7PgL3AnRFxG/A4cPNQKpQkddU3wDPz34FeE97XD7YcSVJdnokpSYUywCWpUAa4JBXKAJekQhngklQoA1ySCuUVeTRUXrxWGh5H4JJUKANckgplgEtSoQxwSSqUAS5JhTLAJalQBrgkFcoAl6RCeSKPRqLXCT6S6nMELkmFMsAlqVAGuCQVygCXpEIZ4JJUqL5HoUTEZ4D3AMcz8w3Vus3AF4DtwGPAezPz2eGVqZXOo0o0Cqv964rrjMD/HrjxgnV7gIOZuQM4WC1LkhrUN8Az817gmQtW7wJmqvszwO7BliVJ6mepJ/JMZOZRgMw8GhFbe20YEVPAFMDExASdTmeJu1ya2dnZxvfZpJXSvvUnzwzttdc8d5r1Tx0e2uuP2pqzp4g8P9ZthOX1411f7v689T2273QeXdJ+lqvpz+PQz8TMzGlgGqDVamW73R72Ll+g0+nQ9D6btFLaN8w58PVPHeb0liuG9vqjdn7tpaw5e2qs2wjN9mN7RHPgTX8el3oUyrGI2AZQ3R4fXEmSpDqWGuAHgMnq/iSwfzDlSJLq6hvgEfF54BvA5RHxRETcBuwFboiII8AN1bIkqUF958Az8309Hrp+wLVIkhbBMzElqVAGuCQVygCXpEIZ4JJUKANckgplgEtSoQxwSSqUV6XXovi939LK4QhckgplgEtSoZxCWcVW++WoNL4GNdW30j8LjsAlqVAGuCQVygCXpEIZ4JJUKANckgrlUSgr2KCOEvHkG2k8OQKXpEIZ4JJUKKdQdBGnXKSF9fqMrD95ptET5ByBS1KhDHBJKtSyplAi4kbgE8AlwKcyc+9AqupiHL63Y1BTE/NfZ6Ff2SQtz0r/bC15BB4RlwCfBN4NXAm8LyKuHFRhkqSFLWcK5S3A9zPz0cw8A+wDdg2mLElSP8uZQnkF8KN5y08A11y4UURMAVPV4mxEPLKMfV7kN/tvsgV4apD7XGHGvX2wWtp47WvGv42roR97tLFGVi3kNd1WLifAo8u6vGhF5jQwvYz9LEtE3J+ZrVHtf9jGvX1gG8eFbRy85UyhPAG8at7yK4Enl1eOJKmu5QT4N4EdEfFLEbEOuAU4MJiyJEn9LHkKJTOfi4g/AP6ZucMIP5OZDw2sssEZ2fRNQ8a9fWAbx4VtHLDIvGjaWpJUAM/ElKRCGeCSVKixCvCIuDkiHoqI8xHR81CeiLgxIh6JiO9HxJ4ma1yuiNgcEfdExJHq9rIe2z0WEQ9GxAMRcX/TdS5Fv36JOX9TPf6diHjzKOpcjhptbEfEj6t+eyAiPjqKOpcqIj4TEccj4rs9Hh+HPuzXxub6MDPH5h/wOuByoAO0emxzCfAD4LXAOuDbwJWjrn0RbfwrYE91fw/wsR7bPQZsGXW9i2hX334BbgK+ytw5CNcC94267iG0sQ18adS1LqONbwfeDHy3x+NF92HNNjbWh2M1As/MhzOz35mepX8FwC5gpro/A+weXSkDVadfdgGfzTn/AWyKiG1NF7oMpb/3+srMe4FnFtik9D6s08bGjFWA19TtKwBeMaJalmIiM48CVLdbe2yXwNcj4lD1dQYrXZ1+Kb3v6tb/1oj4dkR8NSJe30xpjSm9D+tqpA+LuyJPRPwL8PIuD30kM/fXeYku61bUsZQLtXERL3NdZj4ZEVuBeyLicDVyWKnq9MuK77s+6tT/LeA1mTkbETcB/wTsGHZhDSq9D+torA+LC/DMfOcyX2LFfwXAQm2MiGMRsS0zj1a/eh7v8RpPVrfHI+Ju5n59X8kBXqdfVnzf9dG3/sz8ybz7X4mIv42ILZk5Ll8CVXof9tVkH67GKZTSvwLgADBZ3Z8ELvqtIyI2RMRLnr8PvAvo+hfzFaROvxwAfqc6kuFa4MfPTycVom8bI+LlERHV/bcw9xl9uvFKh6f0PuyryT4sbgS+kIj4deB24GXAlyPigcz8lYj4ReauGHRTlvMVAL3sBe6MiNuAx4GbAea3EZgA7q7eQy8C7sjMr42o3lp69UtE/G71+N8BX2HuKIbvA6eA94+q3qWo2cbfAH4vIp4D/he4JatDG0oQEZ9n7iiMLRHxBPDnwFoYjz6EWm1srA89lV6SCrUap1AkaSwY4JJUKANckgplgEtSoQxwSSqUAS5JhTLAJalQ/wcgf8S26koeOwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "b_hats = pd.Series(b_hats)\n", "fig, ax = plt.subplots()\n", "b_hats.hist(ax=ax, bins=50, alpha=.4)\n", "ax.axvline(x=b, color=\"red\")" ] }, { "cell_type": "markdown", "id": "9e7ff372", "metadata": {}, "source": [ "### 4.) Manual OLS\n", "\n", "Run the regression from assignment 1, but now manually calculate the estimates for b, a and $\\sigma$" ] }, { "cell_type": "code", "execution_count": 25, "id": "3cf080dd", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABE9ElEQVR4nO3de3ycZZn4/88158mpObTpIWkphUKlnHoUAWsBFRQFpHUp6squrqDrV3H97q7o7pfdL/vd10/cVVfXVQvrrnjiYItQEeRcu/qilDalQKHYkpYmadq0OSdznrl/fzzPpJM0TTLJTJ6Z5Hq/XmGSe+Z55p6GXHPP/dzXdYsxBqWUUpPP5XQHlFJqutIArJRSDtEArJRSDtEArJRSDtEArJRSDvE43YFCcs0115jf/va3TndDKTX1yHCNOgLOcOLECae7oJSaRjQAK6WUQzQAK6WUQzQAK6WUQzQAK6WUQzQAK6WUQzQAK6WUQzQAK6WUQzQAK6WUQzQAK6WUQzQVuYhs3dfGxm2NNHWGmF9Vwm1rFrF2Sa3T3VJKjZOOgIvE1n1t3LllL229ESqDXtp6I9y5ZS9b97U53TWl1DhpAC4SG7c14nULJT4PItat1y1s3NbodNeUUuOkAbhINHWGCHrdg9qCXjfNnSGHeqSUmigNwEViflUJ4XhyUFs4nqS+qsShHimlJkoDcJG4bc0i4klDKJbAGOs2njTctmaR011TSo2TBuAisXZJLXddt5Ta8gDd4Ti15QHuum6proJQqojpMrQisnZJrQZcpaYQHQErpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDdEsiVdC27mtj47ZGmjpDzK8q4bY1i3RbJjVl6AhYFayt+9q4c8te2nojVAa9tPVGuHPLXrbua3O6a0rlhAZgVbA2bmvE6xZKfB5ErFuvW9i4rdHprimVExqAVcFq6gwR9LoHtQW9bpo7Qw71SKnc0gCsCtb8qhLC8eSgtnA8SX1ViUM9Uiq38hqAReS/RKRNRF7LaKsWkadFZL99W5Vx31dF5ICIvCkiV2e0rxCRV+37visiYrf7ReRBu/1FEVmYccwt9nPsF5Fb8vk6VX7ctmYR8aQhFEtgjHUbTxpuW7PI6a4plRP5HgH/GLhmSNsdwLPGmMXAs/bPiMh5wAZgqX3M90Uk/fnzB8CtwGL7K33OTwOdxpizgW8Dd9vnqgb+AXgnsBr4h8xAr4rD2iW13HXdUmrLA3SH49SWB7jruqW6CkJNGXldhmaM2ZY5KrVdD6y1v78P2Ap8xW5/wBgTBQ6KyAFgtYgcAiqMMS8AiMhPgBuAJ+xj/tE+1ybge/bo+GrgaWNMh33M01hB+/5cv0aVX2uX1GrAVVOWE3PAs40xrQD2bfqvqw5oynhcs91WZ38/tH3QMcaYBNAN1IxwrlOIyK0islNEdh4/fnwCL0sppbJTSBfhZJg2M0L7eI8Z3GjMPcaYlcaYlbNmzRpTR5VSKhecCMDHRGQugH2bXlXfDMzPeFw9cMRurx+mfdAxIuIBZgAdI5xLKaUKhhMBeAuQXpVwC/BoRvsGe2XDmVgX23bY0xS9InKJPb/7ySHHpM+1HnjOGGOAJ4H3i0iVffHt/XabUkoVjLxehBOR+7EuuM0UkWaslQlfBx4SkU8Dh4GPAhhj9orIQ8DrQAL4vDEmvQj0c1grKoJYF9+esNt/BPzUvmDXgbWKAmNMh4j8E/CS/bi70hfklFKqUIg1YFQAK1euNDt37nS6G0qpqWe461IFdRFOKaWmFQ3ASinlEA3ASinlEA3ASinlEA3ASinlEA3ASinlEA3ASinlEA3ASinlEA3ASinlEA3ASinlEA3ASinlEA3ASinlEA3ASinlEA3ASinlEA3ASinlEA3ASinlEA3ASinlEA3ASinlEA3ASinlkLxuyqmUUhO1dV8bG7c10tQZYn5VCbetWcTaJbVOdysndASslCpYW/e1ceeWvbT1RqgMemnrjXDnlr1s3dfmdNdyQgOwUqpgbdzWiNctlPg8iFi3XrewcVuj013LCQ3ASqmC1dQZIuh1D2oLet00d4Yc6lFuaQBWShWs+VUlhOPJQW3heJL6qhKHepRbGoCVUgXrtjWLiCcNoVgCY6zbeNJw25pFTnctJzQAK6UK1toltdx13VJqywN0h+PUlge467qlU2YVhC5DU0oVtLVLaqdMwB1KR8BKKeUQDcBKKeUQDcBKKeUQDcBKKeUQDcBKKeUQDcBKKeUQDcBKKeUQDcBKKeUQDcBKKeUQDcBKKeUQDcBKKeUQDcBKKeUQxwKwiPyViOwVkddE5H4RCYhItYg8LSL77duqjMd/VUQOiMibInJ1RvsKEXnVvu+7IiJ2u19EHrTbXxSRhQ68TKWUOi1HArCI1AFfBFYaY84H3MAG4A7gWWPMYuBZ+2dE5Dz7/qXANcD3RSRdJv8HwK3AYvvrGrv900CnMeZs4NvA3ZPw0pRSasycnILwAEER8QAlwBHgeuA++/77gBvs768HHjDGRI0xB4EDwGoRmQtUGGNeMMYY4CdDjkmfaxNwVXp0rJRShcCRAGyMaQH+FTgMtALdxpingNnGmFb7Ma1AughoHdCUcYpmu63O/n5o+6BjjDEJoBuoGdoXEblVRHaKyM7jx4/n5gUqpdQYODUFUYU1Qj0TmAeUisgnRjpkmDYzQvtIxwxuMOYeY8xKY8zKWbNmjdxxpZTKIaemIN4LHDTGHDfGxIGHgUuBY/a0AvZtm/34ZmB+xvH1WFMWzfb3Q9sHHWNPc8wAOvLyapRSahycCsCHgUtEpMSel70KeAPYAtxiP+YW4FH7+y3ABntlw5lYF9t22NMUvSJyiX2eTw45Jn2u9cBz9jyxUkoVBEf2hDPGvCgim4AGIAHsBu4ByoCHROTTWEH6o/bj94rIQ8Dr9uM/b4xJ71X9OeDHQBB4wv4C+BHwUxE5gDXy3TAJL00ppcZMdFB40sqVK83OnTud7oZSauoZdgWWZsIppZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDNAArpZRDHCnGoybH1n1tbNzWSFNniPlVJdy2ZhFrl9SOfqBSalLoCHiK2rqvjTu37KWtN0Jl0Etbb4Q7t+xl67620Q9WSk0KDcBT1MZtjXjdQonPg4h163ULG7c1Ot01pZRNA/AU1dQZIuh1D2oLet00d4Yc6pFSaqis5oBF5FqsreED6TZjzF257pSauPlVJbT1RijxnfwVh+NJ6qtKHOyVUirTmEfAIvJD4CbgC1jFhT8KnJGnfqkJum3NIuJJQyiWwBjrNp403LZmkdNdU0rZspmCuNQY80mg0xjzf4F3MXijTFVA1i6p5a7rllJbHqA7HKe2PMBd1y3VVRBKFZBspiDC9m1IROYB7VjbyqsCtXZJrQZcpQpYNgH4MRGpBP4FazNNA/xnPjqllFLTQTYB+BvGmCiwWUQew7oQF8lPt5RSaurLZg74hfQ3xpioMaY7s00ppVR2Rh0Bi8gcoA4IisgyTm6vXAHomiallBqnsUxBXA38GVAPfCujvRf4Wh76pJRS08KoAdgYcx9wn4isM8ZsnoQ+KaXUtDCWKYhPGGN+BiwUkS8Pvd8Y861hDlNKKTWKsUxBlNq3ZfnsiFJKTTdjmYLYaN/+3/x3Rymlpo9sakEsEpFfi8hxEWkTkUdFRAsLKKXUOGWzDvgXwEPAXGAe8Evg/nx0qpj0RxNOd0EpVaSyCcBijPmpMSZhf/0MKx15WjvWE6G1O0wskXK6K0qpIpNNKvLzInIH8ABW4L0J+I2IVAMYYzry0L+iEI4laYmHmRH0Uhn04nLJ6AepUemedmqqE2PGNogVkYMj3G2MMUU/H7xy5Uqzc+fOrI5pPN436GePy0VVqZfygDeXXZt20nvaed1C0OsmHE8STxotqamK1bCjsjGPgI0xWnpyiK372vi3Z/bT2hNmbkWQDavms3pRNcd7o3SH49SU+gn63KOfSJ0ic087gBKfh1AswcZtjRqA1ZSR7ZZE5wPnMXhLop/kulPFID1CM8ZQEfDQ3h/lO8/t53YWs3pRNbFEitbuMCU+D1WlXvweDcTZaOoMURkc/ClC97RTU002y9D+Afh3++sK4BvAdXnqV8FLj9CCXjeCdetxCQ+81DTocaFYgpbOMCf6oqRS0/6a5ZjNryohHE8OatM97dRUk80qiPXAVcBRY8yfAxcB/rz0qggMt+twwOviaE942Mf3hOM0d4YJxXTZ2ljonnZqOsgmAIeNMSkgISIVQBswbf8ahhuhReIp5lQET3tMIpXiaHeEYz0R4kldtjYS3dNOTQfZzAHvtLckuhfYBfQBO/LRqWJw25pF3LllLzGTIuB1EYmnSKQMG1aNvk9pfzRBKJakIuChssSHW5etDUv3tFNT3ZhHwMaYvzTGdBljfgi8D7jFnooYFxGpFJFNIrJPRN4QkXeJSLWIPC0i++3bqozHf1VEDojImyJydUb7ChF51b7vuyIidrtfRB60218UkYXj7etw0iO0mlI/vZEENaV+br/SugA3FsYYusNxmjpCdIVijHU5oFJq6hjzCFhEPgI8Z4zpNsYcsgPoDcaYR8b53N8BfmuMWS8iPqzdNb4GPGuM+bqd9HEH8BUROQ/YACzFSoN+RkTOMcYkgR8AtwLbgceBa4AngE8DncaYs0VkA3A3VvJIzqxdUsuCmoldFEoZQ0d/zAriZb6BZVdKqakvmzngf7D3gQPAGNMF/MN4ntSeQ14D/Mg+V8w+3/XAffbD7gNusL+/HnjA3ovuIHAAWC0ic4EKY8wLxhpC/mTIMelzbQKuSo+OC1E8ac0PH+2OEBkyt6xUPmzd18bN92zn8ruf4+Z7trN1X5vTXZp2sgnAwz12vMO1RcBx4L9FZLeI/KeIlAKzjTGtAPZtegKwDshc39Vst9XZ3w9tH3SMMSYBdAM1QzsiIreKyE4R2Xn8+PFxvpzcCcUSHOkKc6RLV0yo/EmvY2/rjVAZ9NLWG+HOLXs1CE+ybALwThH5loicZZem/DbWxbjx8ADLgR8YY5YB/VjTDacz3MjVjNA+0jGDG4y5xxiz0hizctasWSP3Osd2NHbw5Qf3cPO92/nyg3vY0XiynEYknuRod4TmzhB9WnFN5VhmpqGIdet1Cxu3NTrdtWklmwD8BSAGPIhVljIMfH6cz9sMNBtjXrR/3oQVkI/Z0wrYt20Zj89cXlAPHLHb64dpH3SMiHiAGUDBFAza0djBd57bT3t/dFAmXWYQBoglUrT1RGjpChOO6dSEyo3h1rFrpuHky2YVRL8x5o70aNEY8zVjTH/6fhH59yzOdRRoEpFz7aargNeBLcAtdtstwKP291uADfbKhjOBxcAOe5qiV0Qused3PznkmPS51mNdQCyYpQYPvNSExzV6Jl1aNJ6ktTvMsZ4ICV1DrCZIMw0LQy4vuV+W5eO/APzcXgHRCPw51hvCQyLyaeAw8FEAY8xeEXkIK0gngM/bKyAAPgf8GAhirX54wm7/EfBTETmANfLdMM7XlRetPWEqAoP/+UfKpEvrjyYIx5JUlfqoCFgfH5XKVnodeyiWGFRtTjMNJ5dja56MMS8DK4e566rTPP6fgX8epn0ncP4w7RHsAF6I5lYEae+PDvoYOFomXVrKGNr7ovSE48wo8VLu10CssrN2SS13Yc0FN3eGqNd6y47QRacO2bBqPt95bj/heDLrTLq0eDLFid4o3SENxCp7mmnovFwG4Gn1l5/eraHxRN+gWsCj2dHYwQMvNdHaE6bU5wFj6I0kmJPFOYZKB+Ku/vi4i8Hr7hNKTb4xBWARcQNfN8b8zQgP+05uulT4MndrGK4W8OmkVz54XNZx6VHvl646Z1yBd6hEKsXx3uhAVt1YaxCnX088maQ7FKe1O0zD4U4+v/YsvvjecybcL6XU8Ma0CsK+4LVipEwyY8yPc9WpQpdeQ9kVitMbSRDwuEZcwZCW7cqH8YrEkxzpitDRP7YaExu3NRJPJmnvi5M04HW7SBnDf2x9SxfmK5VH2UxB7AYeFZFfYiVOAGCMeTjnvSpwTZ0hAh4X7f0xAI67YlQGPURHSSEe78qH8TDG0BWK0R9NMKvcT8B7+tFwU2eI7lAcEXDZ77FugUTK6BZASuVRNokY1UA7cCXwYfvrQ/noVKGbX1VCNJFihr1lTjJlaO+Pc7wvxree/iOH2vuHPW5uRZBIfPAa3rGufBiveDLFka4wx3ujJE+zI8f8qhKiyRSZn2+MAb/HpQvzlcqjbDblHHfpyakmvYZyZpmPqhIP7X1x+qIJDPDYK6089korqxZWsX5FPSvPqBpYmZCLlQ/j1RuJE4olmFnmp9Q/+Nd+25pFNBzuJJkyuMUKvikM5QGvLsxXKo+y2RPuHBF5VkRes3++UET+Pn9dK1yZuzVE4inOmlXGP374PL5w5dnUVVqj2ZcOdfKVza/yqR/v5LFXjhCJJ1m9qJrbr1w87hrCE5VMGY71RGjriQwaDa9dUsvn156FS4REyuBxCzWlPnwety7MVyqPZKzZuSLyO+BvgI12AR1E5DVjzClJEMVq5cqVZufOnVkd03i8b9DPyZThxYPtbNrVwstNXQPtFQEPH7pwLtdfXMescue30nOJUFniZUbQOzBCTy9F04X5SuXcsAsYsgnALxljVonI7owA/LIx5uLc9dFZuQjAmd5q62NzQwvP7jtGPGn9O7tdwhXnzmLd8nrOnVM+of7mgtftoqrUR5lfc3KUyqNhA3A2f3UnROQs7JKOIrIeaM1Bx6ass2rL+NtrzuUv3n0mW/YcYcvLR+gKx3nmjTaeeaONC+oqWLeinsvOmunYvnDxpFVtrdvrpqbUN+JqCaVUbmUzAl4E3ANcCnQCB4GPG2Pezl/3JleuR8BDxRIpnn3jGJsbWmg8cXKlxJyKAB9ZXscHzp/j+Ei0IuilusSHSzcKVSqXJjYFMXCAtXOFyxjTm4teFZJsA/DWfW382zP7ae0JZ5WObIxhd1MXm3e1sL2xfaBKfInPzTXnz+HGZXXMq8zf0rTReFyucac0K6WGNeE54BqsPeAux5qG+D1wlzGmPVc9dFo2ATidvmuMGbSkLNtVDU0dIR7e3cKTrx0lkrDWCAtw6dk1rF9Rz4V1M/JSYCezJsXp3jz8eZyW0NoTapqZcAB+GtgG/Mxu+jiw1hjz3px0rwBkE4Bvvmc7bb0R3BnBMRxPUlPq51s3XZT1c/dG4vzmlVYeefkIbb3Rgfaza8tYv6KeK86dhdedTd7M6WXWpBjLm0d5wEt1qS9n89SZtTQya9Hedd3Sgg3C+oahJmjYP56sMuGMMf9kjDlof/0/oDInXStCw23pMpG04vKAlw2rF/CzT6/m/1z7Ds6ba62QONDWx9ef2MfN977IT7e/TVcoNuG+Z1uTojcSp7kzRGd/7LTZdNkotv3IdANLlS/ZBODnRWSDiLjsrz8BfpOvjhW64bZ0yUVascft4ooltXzvY8v53s3LuOLcWbgEOvpj/PcfDrHh3hf516fe5OCJ4dOdx6K1J0zAO/hXP9qbRzJl6AzFaOoI0dEfIzWBQFxs+5EV2xuGKh7ZXHK/Dfgy8FP7ZzfQLyJfBowxpiLXnStk6XTkmEnlLa34vHkVnDfvPG7tifDI7hYee7WV/miSx189yuOvHmXFGVWsX1HHqoXVA0V0xmKiu3F0hWL0RRJUlw2/fni0j+vzq0po641Q4jt5bCHvR9bUGaIyOPiCZCG/Yajikc2mnOXGGJcxxmt/uey2cmNMhYgszWdHC006HXky0opnVwS47T1n8dCt7+L2q86mvsoKlLve7uSrD7/Gp368ky17rHTnsdiwaj6JlCEcT2KwbrN980ikTu7WHIolBtrH8nH9tjWLiCcNoVgCY6zbQt6PrNA3sNy6r42b79nO5Xc/x833bNepkSKS9TK0055IpMEYszwnJ3NIvtcB50rKGF5s7GBzQzMNh7sG2isCHq69cC43jCHdOb0K4mhPeEK7caT5vW5mBL185r6dp4xuQ7EEteUB7r/1koG2Ykp7LuSLhoXcNzVIbtYBn/bsGSnKxapYAnCmt473sXnXqenO7zlnFutX1LFkzuTODH3s3u1UlfjwZKzYMMbQHY7zP1+5clL7kkuF+oaRXo0z2huectyEU5FHk5tIrrJy1qyT6c6/3nOELXuO0BmK89y+Np7b18b58ypYv6Key86enHTnOfb8conPg8cluFxSUB/Xx6tQN7DU+eniphVYpojqUh+3XLqQm1cv4Ll9bWxuaOat4/28dqSH1468zuwKPzcuq+MDF8wduHA2lmSMbKVrHodiCQJeF7FoipShYOd3i12xXdBUg+VyCmK7MaaoP/MU4xTE6RhjeLmpi01D0p2DXivd+ayZpfx8x+ExJ2NkY+j88s2r53PN+XOpCHryktU3nekccNGYcCacYGW/LTLG3CUiC4A5xpgdueujs6ZSAM7U3Bni4YYWfrv36KAtkQJeFzNLfVZChsiEMvnGQmtM5Eehzk+rQSYcgH8ApIArjTHvEJEq4CljzKrc9dFZUzUAp/VG4vzm1aM8srtlULqz3+OiMuilLOCmP5rkF58Z+weZ8Uxj+Dwuakr9BH1a+lJNGxNORX6nMebzQATAGNMJ+HLQMTVJygNeNqyaz8//4p0srC7BZ69UiCZSHOuNcvBECJAxpzuna0q090epCHho74/ynef2s6OxY8TjYokUrd1hjnZHxrx2WampKJsAHBcRNycLss/CGhGrIuN2CZ99z1nUlPmoLfdT5rdGoikDR3si3HTPdv71ydHTnbOtKTFUKJbgSFeY1u6wBmI1LWWzCuK7wK+A2SLyz8B6YFpuyumUXK5aWL2omttZPHCxbEFVKbPK/TQ0ddIbSfD4a0d5/LWjrFhQyboV9aw+89R059aeMBWBwf8LjacgUTiWJBwLE/C6qSrx6dSEmjay2Zb+5yKyC7jKbrrBGPNGfrpVPESEXK0kGUlmCcnMj/u3M/5VC6sXVZ9ybDie5Km9x9jc0ExzZ5hdh7vYdbiL+VVBblxez/uXzh6oITGRmhLDicSTtHaH8XvdVJV4By2tUmoqymoZmogs52RB9j8YYxry1TEnjOcinDGGaCJFNJ4ikkgSjadIpHI/M/PlB/ecEuzyuWohZQw7DnawadfgdOfygIdrL5jLR5bVcfB4f1Z1hbPldbuYUeKl3K/L11TRm/AqiDuBjwKb7ZPdAPzSrgs8JYwnAA8nnkwRiiWJxJPEEiniyYkH5Jvv3U5FwINk/B4Nht5IIqtVC+PReLyPhxtaePqNk+nOLoH3nDOL8+ZU8Ie32nNWU2I4bpcwI+ilIuDVvepUsZpwAH4DWGaMidg/B4EGY8w7ctZFh+UqAA+VShliSWuUHE0k7ZFidkF5skfAw+kMxfj1niM8+rKV7py2dF4F65bX8+7F+U13dolQEfQyI+h1bBdppcZpwrUgDgEB7GVogB94a2J9mh5cLiHgctt7q1lJCIlkikgiRTSeJJJIEUukRpxLTqf4huPJvNUfHk1ViY9PvmshG1Yt4Pk329i0y0p33nukh71HXqe23M9HltVx7QVzKQvkfv42XYu4OxynzO+hPODJy351oFsQqcmRzQj4EWAV8DTWHPD7sDbmbAMwxnwxP12cPPkaAY9F5lzy6UbJuS4hOVHGGPY0d7NpVzMvvHUy3TngdfGB8+dy47I66qryu7uz1+0aCMaeHO2Zp+m9Kg8mPAVxy0j3G2PuG0enCoqTAXg4iWQqY+rCGiXn4wJfLrR0hnl4dwtPvNY6kO4swLvOqmHd8jounl85oQtpoy3Bs7YKclMe8AykVo+XlnhUeZDfesBTQaEF4OFkBuX0qotUAf0O+yIJfvi7twZdsAM4a1Yp65bXU+H3sKmhJau1zNnu4ux2CeWB8c8VX373c1QGvYOC+FSoaawcNb45YBF5lRFq/RpjLpxAp1SWPG4XHreLkowk8GgiSTiWJBRLEh1lLjnfXj/Sw+6mLmaX+0mkDB39cWLJFG8d7+cbT76JS6DM76GmzDfmtcyZGXfAwLTAAy81DXtcMmXNFfeE41SWZL96IpclHnUuWY1kLFdKPmTfft6+TW/K+XFAqz4XAL/Hjd/jprLk1HXJkXgyJ1vJj9XQYFkR8NIVjhNLpOgKx0kZ6Ikk6I0krOkCn/u0gTQtM+OuP5agoz9GLJHiaE+EHY0dpz02ZQwd/dZFuwp7RDyWQJzecDUUSwyaA862pnHmXHLm/nh3gQZhBYyhFoQx5m1jzNvAZcaYvzXGvGp/3QFcPZEnFxG3iOwWkcfsn6tF5GkR2W/fVmU89qsickBE3hSRqzPaV4jIq/Z937XLZiIifhF50G5/UUQWTqSvxUJECHjdzCjxMrsiwBk1pdRVBakp81Pm9+R9+dZwW97PCHqs0pdlPqpKvLjE+kjVE0lwrCfK663dbG9sP+1UytyKIJF4iv5YgraeKImUwSUgwpiK/yRThs5QjKbOEF2hGKlR3pDSG67WlgfoDsepLQ+M6wKcbmevRpPNWqFSEbncGPN7ABG5FCid4PPfDrwBpDcuuwN41hjzdRG5w/75KyJyHrABWArMA54RkXOMMUngB8CtwHbgceAa4Ang00CnMeZsEdkA3A3cNMH+FqX0CBl765pYwhodP/d6G/e9cIgj3bnbESOdnpwOevFkCpdLqJ8RpCJoTTvUlProiSTs+w2xpOFrv3qN+qog65bX8f6lcwatd04vwWvviwIGjGCAWaV+3Hbxn7H0O5ka+4g4F1sQ6XZBajTZrNv5NPAfInJIRA4B3wc+Nd4nFpF64FrgPzOarwfSqynuw8q2S7c/YIyJGmMOAgeA1SIyF6gwxrxgrInPnww5Jn2uTcBV6dHxdOfzuGg41Mm3nvkjPZE4M0t9dIdj/PvzB9h5aOTR5Gg2rJpPXzTBsZ4IiWQKAZJJQ2c4zrL5M0ikrCmSGUEPs8v9VJV4WVxbBkBzZ5jvPHuAmzZu555tjbT1WEvOVy+q5vYrF2OwRs4et1BbHqDM7xlX8Z/0m8PhjhDtfdGcZCoOp9C3s1fOy6YYzy7gIhGpwFo90T3B5/434G+B8oy22caYVvv5WkUkPQSpwxrhpjXbbXH7+6Ht6WOa7HMlRKQbqAFOZHZCRG7FGkGzYMGCCb6k4pH58RigLOAiFEvwq91HuH5ZHX2RBP3RZNbL3lYvqqaqxEd/NEHKGLxuF9WlPlwi7G7q5vYrFw9ay/y/7FH3wRP9bG5o5unXj9EXTfDAS008tLPJ3t25ntWLqlk6d8Yp2YBdoRihWIqb792e9Sg+Za9s6A7HKfV7KPN7KPFNbAlbplzNJaupa8wBWES+PORngG5glzHm5WyeVEQ+BLQZY3aJyNqxHDJMmxmhfaRjBjcYcw9wD1jL0MbQlylhpI/Hfo8bf5mbmjKrVGRfNEE4NvZg3B9LcEZNySl1K472hIetwAZw5sxS/vr95/IXl5/Jr19p5dGXj9DRH+P5N4/z/JvHOW9uBRfXz+C5N63RbsDroisUo70/TnWJd8IV4vqjCfqjCVwilPo9zAh68Xkmltixdkktd4FuF6ROK5s54JX216/tn68FXgI+KyK/NMZ8I4tzXQZcJyIfxEpvrhCRnwHHRGSuPfqdi51lhzWyzcy5rQeO2O31w7RnHtMsIh5gBjCxz9dTyFiXWgV97oH6vJG4taoiXWjodCZSprKyxMefXnIGN62cz9Y329jU0MKBtj5eb+3h9dYea95WhO5wnFAsRXWJl+pSv9XXUZanjUXKGHojcXojccoCHioC3gmlOxfqdvaqMGTzFl8DLDfG/G9jzP/GCsazgDXAn2XzpMaYrxpj6o0xC7Eurj1njPkEsAVIZ9zdAjxqf78F2GCvbDgTWAzssKcrekXkEnt+95NDjkmfa739HNNmhDua29YsIp40hGIJjLFuR/t4HPC6qSzxMa8yyBk1VgH34Yqnb1g1n0TKEI4nMVi32dat8HlcvH/pHDZ+YjnfvukiLju7BgG6w3GOdEfoiVhTHKX+iReET9vR2MGXH9zDzfdu58sP7uG519s40hWmyZ4r1l07VK5lMwJeAGRuFhYHzjDGhEUkeppjsvV14CER+TRwGKv8JcaYvSLyEPA6kAA+b6+AAPgc8GMgiLX64Qm7/UfAT0XkANbId0OO+jglTPTjcTrbrDzgJZFM0Re11vbGk6lTdtuYSN0KEeGi+kouqq+kpSvMr3a38MSrRwcubr3dEcYlUFXipbrERyQxvoLwoxW87w6n6A7H8bhclPjd9gVA3blDTcyYUpHt0eX/wVphkB5hfhhrlPlN4B5jzMfz1MdJUwypyIUuHEvSG40TiibzliK97c3jfPOZP9oX+k62u11Wlt1fv+9cLls8M6tzjqfcp9ftotS+cOf3uLRovBrJ+MtRGmOMiFwPfAZrRwwBPmuMSUerog++KjfSc8amzJp66Ismch6MH3n5CJVBL3PK/bT3xwYy7JIp6A4n+NYzf+TA8T4+fNE8qkvHtnH3ePa3iydTdIVidIWsTwRBn5tSn1UMSAvHq7HIZgpiO+AyxnwnX51RU0c686vE58GUGSLxFOF4knA8SXSCc6npYCkIM8v8zCzzE44nON5nJXZ0huLc98Lb/GLHYa5aMpt1K+o4a1bZiOec6P52yZShL5KgL5KwsxFdlHg9BHwuKwlGqWFkE4CvAG4TkbeBfqxRsNFiPGo0IjJoNUXSvkDXH00QiiWzLh40XLAEYcnsCr76wSVs2XOEX+85Qk8kwW/3HuW3e4+ybEEl65bXccmimlN2d4bcFrw3xtg7PSehn4F54/KAR4OxGiSbesBnDNdu14mYEnQOePIlkil67eI8Y11nPJbylJF4kqdfP8bmhhYOd5xM/a2rDHLj8jquWTrnlBUck1Hw3ut2UWK/GU20brEqKloPeDQagJ0ViiXoCScIxRKjPnaswTJlDDsPdbJpVzM73+4caC/ze/jgBXP4yLI6ZlcEcvo6xsolMrCiwu9x6z53U5sG4NFoAC4M4xkVj8Wh9n4277J2d44lTp63MujlY6sXsH5l/QhH55/X7cLvcRHwuQl43BPOxFMFRQPwaDQAFxZjDH3RBD2RxIQv3GV6/o02vv3sfkKxwcvY5lcFueXShaxZPDNn+8tNhNftojzgodTvwVsA/VETogF4NBqAC1cknqQ3kqAvmpjwjh/pNb8Bj4veaILOUJxoxoi4ttzPDRfP49oL51Ie8I5wpsnj97oHigVpMC5KGoBHowF48ox3q55kyqpg1hOOj3tt8c33bh9Yxgb2qoV40l7GlhoYFQc8Lq5eOocbl9cxvzq/JSRH23Q0k8flIuC1tqbyumVg6kIv6BU0DcCj0QA8OXKx7XsqZeiNJuiNxAfN547FSFlvf331OVa682tHCcVOTntcsqia9cvrWbZgYrs7DyfbTUeHIyL4PC4CHpeusChMGoBHowE4t043ys31tu+ReJKeSJz+6NjWFI8l4PVHEzz+2lF+1dDCUbswPMCimaWsW17HVe+YnbOLZONJgx5NOhkk6LWWvOn6Y8dpAB6NBuDcGWmU+/ePvpaXbd+TKUNPOE5PJD7qRqRjXcaWTBn+8NYJNu9q5tWWnoH2yqCX6y6ax3UXjz3d+XSGTomAVT+5N5LgF5/J/g1pOOlU6RKfhxJNlXbC+GtBKJWtoTtulPg8hGIJNm5rzOm275ncLqGq1EdliZfeqJUWPFoJydGGH26XsGbxLNYsnsUfj/WyaVczz795nK5wnJ9sf5v7XzrMlUtqWb+8nrNqR053Pp2JpkGPxXCp0jo6dp5eTlV50dQZGpIqfHLHjfHUIs6GiFAR8DKvMkh9VQnlgcGj7fQURHt/dFDpydF2Vz5ndjlf++A7uP8z7+Tj71xARcBDPGl4cu8xPvPTXXz5oT384cCJrC8O5qJ+cjbSqdId/TFaOsMcOtHP0e4I3aE40YTWPJ5MOgWRQacgcme0ed70/PBkbdVjVS6L0xdN8FcPvJyTOddIPMkzbxxj864W3h6S7vyRZXV84PxT051PZzLSoMfK7ZKB0XGJz6MZermhc8Cj0QCcO7lY6ZAPyZTh3Xc/R6nfM+gvYiJzrsYYdr7dyeZdzew4dDLdudTv5toL5nLDsjrmOJTunAvpQFzqcxdEgkqR0jlgNXkKdUNKt0s4o6aUYz1hAl6PlepsJjbnKiKsWljNqoXVvN3ez8MNLTz1+jH6o0ke2tnMpl3NXL54JuuX17N0XkXRLQ9LV3Zrx9qWqtTnocSvCSG5oCPgDDoCnh4yR+cBj4v+WJJYIsUXs1h3O5ruUJzHXj3CI7uP0N5/cievJXPKWbe8nvecUxjpzhPh87go9Vmp0lq3YlQ6BTEaDcDFYbxZdMOdI3N0bu39Frfq+E5AZlbb7PIA75hTzu7mLv54rG/gMbPK/NywbB7XXjCXiuCp6c7ZZMYVgsztmXSvvGFpAB6NBuDCNxlzy9FEkm77gl22Tpfk8cUrzqbE72FTQ7O9UsJ6fMDe/fnG5XUssNOdc5EZ56R0zeNSv0dTpE/SADwaDcCFL9dZdCOJJpJ0heL0ZxGIx5LV1tpt7e78+KuD053feWY165bX8fPth+kIxXKaGZcL4xmVp9cc++3ymn6Pa7rOHetFOFX8mjpDVA75yJ5eX5wpF9MUfo+b2RVuInFrzexoSR0wts09584I8pdrz+aWdy3kt3uP8nBDC63dEV482MGLBzvwuITqUi9+j2tg+6TRNgjNt8xReeba6dsZeVQ+aHsmm8uuW5H+8ntc+NzTc6Q8Ld+KVPGaX1VCeEggHJpFl56maOuNUBn00tYb4c4te9m6r21czxnwuplXGWTujCD+UeY351YEicQHFwc63QqLUr+Hdcvr+cmnVnPXdUu5sH4GAImUoa03xsETIU70RUkkUznPjMvWAy814bHXBwvWrcclPPBSU9bnShlj1e8IxznRG7WSQdpDHOkK02m/0U2XT+YagFVRGUsW3cZtjcQSSY52R3jzWC9HuyPEEkk2bmuc0HMHfW7qKoPMmRE4bYLFeLLa3C7h8sUz+bebLuaHn1jOsvmVACSNoSMUp7HdCsRrFs+cUP8norUnTMA7OFzkclRu7KDcGYpxpCvM2+2hgey8bKvdFRMNwKqorF1Sy13XLaW2PEB3OE5teeCUC3B/PNZDe3+MRNLgFiGRNLT3x9h/rGeEM49dic/D3BlWmnOpf/B0w+pF1dx+5WJqSv30RhLUlPqzunh2zuxyvvknF/G1a5Ywu9xP+lN5OJ7iu88f4K8efJk/HDgxarGhXMtmZJ8LKfvNtb0/SnNniMPtIdp6I/RFE5P+2vNJL8Jl0ItwU8OF//gk4XgSj+vk+CKRShH0unnlH6/O+fOFY0k6QrGcbpuUFo0nefqNNjY3NPN2+8l57nmVAW5cVsc1588ZdEEyXwptZYbP46LE5yHodRPwFsX8sa6CGI0G4Klh5f97mu5QHJdLEAFjrALuM0q87Pz79+XtecOxJN3h+Jh2dc7WQLpzQws7Dp4sGlTqc/PBC+bykWV1zJmR33TnQqpXkUnErl3hdRPwuQq1upsG4NFoAJ4abr5nO4fa++gJJ4glU/jcLiqCHhbWlOV8qdpwYokU3eF4TvavG87b7f08vLuFp/YeG9jLziVw+dkzWb+iONOdcym9ZZPf48bvLZjtmjQAj0YD8NRQKIWAkilDbyROTzhh1ZzIse5wnN+80sojL7dwou9kuvO5c8pZv7yO95wzq+jTnXNBRPC6ZVBAdmCUrAF4NBqAp47JLnc5EmMMPZEEXaFYXi4gJZIpfvfHE2za1cybx3oH2meW+bjh4jquvXAuM4ZJd57O0muRA173wF56eX6z0gA8Gg3AKp9S9o7O3RPY0Xkkxhj2Hulh065mfp+R7uz3uHj/ebNZt7yeBTX53d25mLld1ig54E0niLhzWQtZA/BoNACr4eQiqy5TMmXoDMXojeRnjhjgaHfETndupT8jC231mdWsX17HijOqCmFetOB53faUhddtT12Mez5ZA/BoNACrofI5nxxLpGjvj064+tpIQrEEv33tKA/vbuFI18ndnRfWlLBueT3vfUftqNl96qR0Pelx0AA8Gg3AI8v1SLAYTEbxn/5ogva+WF4u1KUlU4btje1s2tXMnubugfaKgIcPXzSPGy6eR02ZP2/PP1W4RFg4UwNwXkz3ADxSgC2UlQWT7fK7n6MyOHhTT2Osudz/+cqVOXueVMrQZc8P5/tvcv+xXjY3tPDcvjYS9kSxxyVcsaSW9cvrWDy7PK/PX8xyHYB1jYoCRi9gk7nNvIh163XLhOsrFLqxFP/JBZdLqC71UVcZzHtB88Wzy7njA0u4/zPv5E8vWcCMoJdEyvD068e47WcNfOnBl/n9/slPd56OtBzlFJCLqYHMAAtWvYNQLMHGbY2sXVI75jKQU81taxZx55a9hGKJQSP/zOI/ueTzuJhXGaQ3EqezP57XaYmaMj9/ftmZfGz1Ap7d18amXc0cag/xSnM3rzR3M3dGgBuX13HN0jmn1LxQuaEj4CKXq9KLTZ2hQQXAYXCAnayRYKEZS/GffCgPeJlfHaS61DdQEzhf/F4rnflHt6zkG+su4J1nWunFrd0R/uP5t9hwz3a+v/UArd3O1SOeqvRtrciNNnIdq/lVJadcbMoMsJM9Eiwka5fUOjLPLSJUlvgoD3jpCsXoyeOytfTzrVxYzcqF1RxuD7F5dzNP7T1GfyzJpl0tPNzQwmVnW7s7n183vdOdc8WREbCIzBeR50XkDRHZKyK32+3VIvK0iOy3b6syjvmqiBwQkTdF5OqM9hUi8qp933fF/r9CRPwi8qDd/qKILJz0FzoJRhu5jtVodXadGgkqa+lTTZmf+qrgpE0FLKgp4a/eew4P3noJn3n3mcws85Ey8D/7T3D7gy/zlz/fzTNvHCOenLq1eieDI6sgRGQuMNcY0yAi5cAu4Abgz4AOY8zXReQOoMoY8xUROQ+4H1gNzAOeAc4xxiRFZAdwO7AdeBz4rjHmCRH5S+BCY8xnRWQD8BFjzE0j9asYV0HkcplUIaXvqtMLx5Kc6ItOavBLpztvbmhm39GT6c41ZT5uuHgeH7pw3rRId56Sy9BE5FHge/bXWmNMqx2ktxpjzhWRrwIYY/4/+/FPAv8IHAKeN8Yssdtvto+/Lf0YY8wLIuIBjgKzzAgvuBgD8HRdHjadbd3Xxg9/9xZvd4SYXR7IuizkRLa8H0h3bmjm9/sHpzu/77zZrFteN95EhaKQ6wDs+BywPTWwDHgRmG2MaQWwg3A6gtRhjXDTmu22uP390Pb0MU32uRIi0g3UACeGPP+twK0ACxYsyNnrmixrl9RyF+jIdZrIfMOtKfXRHY6NaXPMtPFurpkmIpxfN4Pz62ZwtCfCrxpaePy1VvqjSR57pZXHXmll9cIq1q2oZ6WmO4/K0QAsImXAZuBLxpieEX5Zw91hRmgf6ZjBDcbcA9wD1gh4tD4XIqcuEqnJN/Sia1nAi0TjPLSraUwBNHNzTWDgU9MDL43t+ExzKgJ8bu1Z3HLpGTy59xgPN7TQ0hVmx6FOdhzq5IzqEtatqON975it6c6n4dgyNBHxYgXfnxtjHrabj9lTD+l54vRaqmYgc1fDeuCI3V4/TPugY+wpiBlAB0oVseEuupb4PBzvjVJbERi0DdNw8rG5ZonPw0eW1fHjP1/FP12/lIvnW7s7v90R4ltP7+eme7bzo98f5ERfdNzPMVU5MgK2Vyr8CHjDGPOtjLu2ALcAX7dvH81o/4WIfAvrItxiYId9Ea5XRC7BmsL4JPDvQ871ArAeeG6k+V+lisFIywV3Huzgh797i8MdIWorAmxYeerc7tyKIO39UVLG0NEfI55M4RLJyXput0u47OyZXHb2TA609bG5oZnn9rXRE0nw8xcP8+BLTaw9dxbrV9RzjqY7A86tgrgc+B/gVSB9KfdrWEH0IWABcBj4qDGmwz7m74BPAQmsKYsn7PaVwI+BIPAE8AVjjBGRAPBTrPnlDmCDMWbEvNlivAinppfTXXRdv7yOTQ0tA+2hWIJowvDFK89m1Zkng/COxg7ufnIfPeE46VK3KQMVQS9fuXpJzvd46+iPseXlI2zZc4SucHyg/YK6GaxbUcdlZ83MZc3dvJuSqyAKhQZgVQyGWy64cVvjsMsRZ5X5+f4nVtCTUQT+L378Es3dYVIpg9ftoqrEZ601LvXzrZsuykufY4kUz75xjE0NLRw80T/QPqciwEeW1/HB84sj3VkDcB5pAFbFarSqbZnV1jbc8wIVAQ+SERMMht5Igl98Jr+blhpj2H24i00NzWxvPHlJpsTn5gPnz+HG5XXMnRHMax8mYsotQ1NKTdxoqeTpamtlfg91lUGO90YHXcyLxFPMqch/4BMRlp9RxfIzqjjcEeJXu1t48rWjhGJJNje08KvdLVx61kzWr6jjgroZU34ZmxbjUWoKGC2VPM3ncfHFKxdjgEgiicEQjidJpAwbVs0f/uR5sqC6hNuvWsyDt13Cre8+k1llflIGfn/gBF96cA+f/VkDT70+tdOddQoig05BqGKWTSp5OpvucEeI2nFk0+VDIpnif/afYFNDM2+0ZqQ7l/q4/uJ5fPjCecwocTbdWeeA80gDcP5Mx+2MikUknqS9P0Y0nr+96bK190g3m3e1sG3/8YF0Z5/HxfveMZsbl9dx5viC4IRpAM4jDcD5ofUqikN3OE5nf2xgtUQhONYT4ZHdLTz2qpXunLbyjCrWrahj1cLqvNdLzqQBOI80AOfHZGxsWWiKdcT/7OvH+P7Wt2juCmVdqCefwrEkT+49ymY73TltQXUJ65bX8b7zZud9KyfQAJxXGoDzY7I2tiwUxTriz+y33+2iz76Qd/uVYyvUMxlSxtrdeXNDC7sPdw20VwQ8fOjCuVx/cR2zyvO3u7NuyqmKznTbzqhYNzDN7Lfb7aIi4MXvcfHAzianuzbAJcKlZ83kmx+9iHv/dAXXLJ2D1y30RBL8YkcTH/vPF/nn37zBmxk1iwuZrgNWeTdVtjMa67RCsW5gOrTfIkKZ30NTex9/88tXCm5a4qzaMv72mnP5i3efyZY9R/j1niN0huI8u6+NZ/e1cUFdBeuW13PZ2YWb7qwBWOXdVKhZnPnxPHPz07vglNcxWlKE0073RjJcv0/0RemLpegKx6gp8dERyq5+8GSoLvXxZ5cuHNjdeXNDM43H+3m1pYdXW1630p2XzeMDF8ylrMDSnXUOOIPOAavTyeZCYiHPAY/UN+CU+5o7w1SXeplZFgCsufu+aIKqEh/f/JP81I2YKGMMu5u62LRrcLpz0OvmAxfM4cZldcyrHF/Wn6YiqzEp1qvwhSqbaYXxjPgn6/c10i7a9996ySn97g7HqSk9eVErPS1xvDfCzHI/Xf1xEqnCylQTEZYvqGL5giqaOkI8bKc7h+NJHm5o4VcNLVx6dg3rV9RzocPpzjoCzjBVRsCFPAIrVvlcSjeZv69sV6SM9rozi/wUcizpjcT5zatHeWR3C229JwvDn11bxvrldVyxpBave/Q1CboKQo2qWK/CF7Kx1loYj8n8fWW7ImW0150u8lNfFSzocpLlAS8bVs3nZ59ezf+59h2cN9cqCH+grY+v//ZNbr73RX66/W26Q/FRzpRbGoCnoOG2rSmGq/CFbO2SWu66bim15QG6w3FqywM5G6FO5u8r2zeSsb5ur9vF7IoAc2aMvi2SkzxuF1csqeV7H1vO925exhXnzsIlVuH4//7DIW66dzv/+tSbg2oW55NOQWSYKlMQ0zHzrJhN9u8rm6I945FKGU70R+mLJHJ2znxq64nwyMtHeOyVVvqiJ/u84owq1i2vY/WZJ9OdNRMuj6ZKANY54OIyVX9foViC9r5Y0ZSTDMeTPGWnOzd3nkx3nl8VZN2Ket533mxKfR4NwPkyVQIw5H+Uo3Jrqv6+jDH0hBN0hgqryM9IUsaw42AHm3Y105CR7lwe8PDhC+fxhavOHs+uHRqARzOVArBShSSZsnZh7o1M7kWuiWo83sfmhhaeeeMY8aQVKz0u4f5bL2HVwqwSUTQAj0YDsJru8r0eOZpIcrw3SixRHNMSaZ2hGL/ec4RHXz6C1+3i91+5Ep8nq4uNGoBHowFYTXUjBdjJmos2xtAZitMViuXsnJMlkTSkMCyZU5HtoboOWKnpLB1g23ojg+pZbN3XBkzeemQRa+3wvMrgmJIfConP4xpP8D2t4nr1SqlxGy3ATvb68YDXTX1VkMoSX17OXwwKN3VFKZUzW/e10XC4k2Qqhd/jZla5n/KAd1CAdaKKW3o0HPS6Od4bLbi6EvmmI2Clprj01IOIlUiQSBmOdEXojcQHBdh8pluPJuhzU1cVLLhykfmmAVipKS499TC7PIAB0v852h0ZFGDzmW49Fm6XUGunMxfb3PB4Ta+3G6WmoXQpTfFZF+JP9EWJJQ0GTgmwa5fUOp4AUuLzEKxy0xkq/CprE6UBWKk8KKR6zJlzuxVBLxVB70CtCaeD7emk54bL/B7a+6OEY8nRD8rCjsYOHnipidae8CnbLI10X67pOuAMug5Y5UKh1XYotP6MR28kTkd/jGRq4vFqR2MH33luPx6XEPC6iMRTJFLW7s/AiPc9sLOJE33R8bypaiLGaDQAq1woxGp0U6HWRK7Smb/84B7a+6ODltyF48mBnT+Gu8/rdtm3QkXAO543Md2SSKnJUIi7IhfC3O5EuV1iL5/zWPPY40xnbu0JUxEYHPoCXhdHe8IYGPa+Q+0h5lQECHjdA2uo01s5TeTfdXpcalRqEmW764TKjpXAUUJNmX+gTm825lYEicQHB+9IPMWciuBp77Oed3C4zMWbqgZgpXLMyfW008mMoJf6cawd3rBqPomUIRxPYrBuEynDhlXzT3vf/MpTA3Mu3lQ1ACuVY06vp51OPG4XtRUB5s4Ye12J1Yuquf3KxdSU+umNJKgp9XP7lYtZvaj6tPfduuYsEilDJJ7M6ZuqXoTLoBfhlCpe6Spr+Vo7vKOxgwd2NtHeFx3PhUy9CKeUmroy1w6f6IsSied27fDqRdVcclbNeLckGpZOQSilphSfx8W8yiAzy/24XdlfpJtMGoCVUlNSRcBLfVVJQRf4mfIBWESuEZE3ReSAiNzhdH+UUpMns8CPx1V44a7wepRDIuIG/gP4AHAecLOInOdsr5RSk63E56G+KkjFkAQZp03pAAysBg4YYxqNMTHgAeB6h/uklHKAyyXMLPMX1FZIhdGL/KkDmjJ+brbbBojIrSKyU0R2Hj9+fFI7p5SafOmtkKpKfMg4MulyaaoH4OH+dQctEDTG3GOMWWmMWTlr1qxJ6pZSykkiQlWpj7rKIIEh++BNpqkegJuB+Rk/1wNHHOqLUqrApJesjbeuxERN9QD8ErBYRM4UER+wAdjicJ+UUgUmXVcis4ToZCjcBXI5YIxJiMj/Ap4E3MB/GWP2OtwtpVQB8rhdzJkRoC+aoL0vmpPi76M+Z96fwWHGmMeBx53uh1KqOJT5PQS9btr7o/RFEnl9rqk+BaGUUllzu4Ta8kDel6xpAFZKqdNIL1mrKc3PRbopPwWhlFITISLMKPFSFvDQHZ7YfnRD6QhYKaXGwO2yyl3mkgZgpZRyiAZgpZRyiAZgpZRyiAZgpZRyiAZgpZRyiAZgpZRyiAZgpZRyiAZgpZRyiAZgpZRyiAZgpZRyiAZgpZRyiAZgpZRyiAZgpZRyiBiT/203ioWIHAfenuBpZgInctAdJxRr37Xfk6tY+w3O9f2EMeaaoY0agHNMRHYaY1Y63Y/xKNa+a78nV7H2Gwqv7zoFoZRSDtEArJRSDtEAnHv3ON2BCSjWvmu/J1ex9hsKrO86B6yUUg7REbBSSjlEA7BSSjlEA/AYiMh/iUibiLw2zH1/LSJGRGZmtH1VRA6IyJsicnVG+woRedW+77siIk70W0S+YPdtr4h8oxj6LSIXi8h2EXlZRHaKyOoC7Pd8EXleRN6w/21vt9urReRpEdlv31YVUt9H6Pe/iMg+EXlFRH4lIpWF1O+R+p5xf8H+fQJgjNGvUb6ANcBy4LUh7fOBJ7GSN2babecBewA/cCbwFuC279sBvAsQ4AngA5Pdb+AK4BnAb/9cWyT9fir9vMAHga0F2O+5wHL7+3Lgj3b/vgHcYbffAdxdSH0fod/vBzx2+92F1u+R+m7/XNB/n8YYHQGPhTFmG9AxzF3fBv4WyLySeT3wgDEmaow5CBwAVovIXKDCGPOCsX7bPwFucKDfnwO+boyJ2o9pK5J+G6DC/n4GcKQA+91qjGmwv+8F3gDq7D7eZz/svox+FETfT9dvY8xTxpiE/bDtQH0h9Xukvtt3F/TfJ+gUxLiJyHVAizFmz5C76oCmjJ+b7bY6+/uh7ZPtHODdIvKiiPxORFbZ7YXe7y8B/yIiTcC/Al+12wuy3yKyEFgGvAjMNsa0ghUwgFr7YQXX9yH9zvQprFEhFGC/YXDfi+Xv05PvJ5iKRKQE+Dusj2in3D1MmxmhfbJ5gCrgEmAV8JCILKLw+/054K+MMZtF5E+AHwHvpQD7LSJlwGbgS8aYnhGmEguq70P7ndH+d0AC+Hm66TT9K4h/c6y+FsXfp46Ax+csrPmjPSJyCOujWYOIzMF655yf8dh6rI/LzZz8CJfZPtmagYeNZQeQwipQUuj9vgV42P7+l0D6IlxB9VtEvFiB4OfGmHR/j9kfcbFv09M+BdP30/QbEbkF+BDwcfujeUH12+7j0L4Xz99nvieZp8oXsJAhF+Ey7jvEyUn+pQye5G/k5CT/S1gjz/Qk/wcnu9/AZ4G77O/Pwfo4JkXQ7zeAtfb3VwG7Cu3f236enwD/NqT9Xxh8Ee4bhdT3Efp9DfA6MGtIe0H0e6S+D3lM4f595vsJpsIXcD/QCsSx3ik/fbpfsP3z32FdXX2TjCupwErgNfu+72FnIk5mvwEf8DO7Hw3AlUXS78uBXfYfz4vAigLs9+VYH1tfAV62vz4I1ADPAvvt2+pC6vsI/T6A9QadbvthIfV7pL4PeUxB/n0aYzQVWSmlnKJzwEop5RANwEop5RANwEop5RANwEop5RANwEop5RANwEpNkIgsHFpxTqmx0ACslFIO0VoQSp2GiDyClbYaAL4DuIEzjTF/a9//Z8AK4JuAW0TuBS4FWoDrjTFhB7qtiogmYih1GiJSbYzpEJEgVprqVcAfjDFn2/c/AfwzVrbeAWClMeZlEXkI2GKM+ZlTfVfFQacglDq9L4rIHqxauPOxaweIyCUiUgOcC/zBfuxBY8zL9ve7sGpZKDUinYJQahgishar3OW7jDEhEdmKNRXxIPAnwD7gV8YYY5ebjGYcngSCk9lfVZx0BKzU8GYAnXbwXYJVJQuskpg3ADdjBWOlxk0DsFLD+y3gEZFXgH/CmobAGNOJVaLxDGPVU1Zq3PQinFJKOURHwEop5RANwEop5RANwEop5RANwEop5RANwEop5RANwEop5RANwEop5ZD/H4SwxXzTm/42AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sns.lmplot(x=\"avh\", y=\"rgdpe_per_capita\", data=pwt2019)" ] }, { "cell_type": "code", "execution_count": 26, "id": "348430ed", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "a-hat 114221.65253311764\n", "b-hat -42.361747336876576\n", "sigma-hat, sigma_hat\n" ] } ], "source": [ "n = pwt2019[\"avh\"].count()\n", "\n", "pwt2019 = pwt2019.loc[(pwt2019.avh.notna()) & (pwt2019.rgdpe_per_capita.notna()), ]\n", "\n", "b_hat = np.cov(pwt2019[\"avh\"],pwt2019[\"rgdpe_per_capita\"], ddof=1)[0,1]/np.var(pwt2019[\"avh\"], ddof=1)\n", "\n", "\n", "a_hat = np.mean(pwt2019[\"rgdpe_per_capita\"]) - b_hat*np.mean(pwt2019[\"avh\"])\n", "\n", "\n", "\n", "y_hat = a_hat + b_hat*np.mean(pwt2019[\"avh\"])\n", "\n", "sigma_hat = np.sqrt(np.sum((pwt2019[\"rgdpe_per_capita\"]-y_hat)**2)/n-2)\n", " \n", "print(\"a-hat\", a_hat)\n", "print(\"b-hat\", b_hat) \n", "print(\"sigma-hat, sigma_hat\")\n", " " ] }, { "cell_type": "code", "execution_count": null, "id": "8fa70624", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.12" } }, "nbformat": 4, "nbformat_minor": 5 }