5-The-Unreasonable-Effectiveness-of-Linear-Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## All You Need is Regression\n",
    "\n",
    "When dealing with causal inference, we saw how there are two potential outcomes for each individual: \\\\(Y_0\\\\) is the outcome the individual would have if he or she didn't take the treatment and \\\\(Y_1\\\\) is the outcome if he or she took the treatment. The act of setting the treatment \\\\(T\\\\) to 0 or 1 materializes one of the potential outcomes and makes it impossible for us to ever know the other one. This leads to the fact that the individual treatment effect \\\\(\\tau_i = Y_{1i} - Y_{0i}\\\\) is unknowable. \n",
    "\n",
    "$\n",
    "Y_i = Y_{0i} + T_i(Y_{1i} - Y_{0i}) = Y_{0i}(1-T_i) + T_i Y_{1i}\n",
    "$\n",
    "\n",
    "So, for now, let's focus, on the simpler task of estimating the average causal effect. With this in mind, we are accepting the fact that some people respond better than others to the treatment, but we are also accepting that we can't know who they are. Instead, we will just try to see if the treatment works, **on average**. \n",
    "\n",
    "$\n",
    "ATE = E[Y_1 - Y_0]\n",
    "$\n",
    "\n",
    "This will give us a simplified model, with a constant treatment effect \\\\(Y_{1i} = Y_{0i} + \\kappa\\\\). If \\\\(\\kappa\\\\) is positive, we will say that the treatment has, on average, a positive effect. Even if some people will respond badly to it, on average, the impact will be positive.\n",
    "\n",
    "Let's also recall that we can't simply estimate \\\\(E[Y_1 - Y_0]\\\\) with the difference in mean \\\\(E[Y|T=1] - E[Y|T=0]\\\\) due to bias. Bias often arises when the treated and untreated are different for reasons other than the treatment itself. One way to see this is on how they differ in the potential outcome \\\\(Y_0\\\\)\n",
    "\n",
    "$\n",
    "E[Y|T=1] - E[Y|T=0] = \\underbrace{E[Y_1 - Y_0|T=1]}_{ATET} + \\underbrace{\\{ E[Y_0|T=1] - E[Y_0|T=0]\\}}_{BIAS}\n",
    "$\n",
    "\n",
    "Previously, we saw how we can eliminate bias with Random Experiments, or **Randomised Controlled Trial** (RCT) as they are sometimes called. RCT forces the treated and the untreated to be equal and that's why the bias vanishes. We also saw how to place uncertainty levels around our estimates for the treatment effect. Namely, we looked at the case of online versus face-to-face classrooms, where \\\\(T=0\\\\) represent face-to-face lectures and \\\\(T=1\\\\) represent online ones. Students were randomly assigned to one of those 2 types of lectures and then their performance on an exam was evaluated. We've built an A/B testing function that could compare both groups, provide the average treatment effect and even place a confidence interval around it.\n",
    "\n",
    "Now, it's time to see that we can do all of that with the workhorse of causal inference: **Linear Regression**! Think of it this way. If comparing treated and untreated means where an apple for dessert, linear regression would be cold and creamy tiramisu. Or if comparing treated and untreated is a sad and old loaf of white wonder bread, linear regression would be a crusty, soft crumb country loaf sourdough baked by Chad Robertson himself.\n",
    "\n",
    "![img](./data/img/linear-regression/you_vs.png)\n",
    "\n",
    "Lets see how this beauty works. In the code below, we want to run the exact same analysis of comparing online vs face-to-face classes. But instead of doing all that math of confidence intervals, we just run a regression. More specifically, we estimate the following model:\n",
    "\n",
    "$\n",
    "exam_i = \\beta_0 + \\kappa \\ Online_i + u_i\n",
    "$\n",
    "\n",
    "Notice that \\\\(Online\\\\) is our treatment indication and hence, a dummy variable. It is zero when the treatment is face to face and 1 if it's online. With that in mind, we can see that linear regression will recover \\\\(E[Y|T=0] = \\beta_0\\\\) and \\\\(E[Y|T=1] = \\beta_0 + \\kappa \\\\). \\\\(\\kappa\\\\) will be our ATE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   78.5475</td> <td>    1.113</td> <td>   70.563</td> <td> 0.000</td> <td>   76.353</td> <td>   80.742</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>format_ol</th> <td>   -4.9122</td> <td>    1.680</td> <td>   -2.925</td> <td> 0.004</td> <td>   -8.223</td> <td>   -1.601</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import statsmodels.formula.api as smf\n",
    "import graphviz as gr\n",
    "%matplotlib inline\n",
    "\n",
    "data = pd.read_csv(\"data/online_classroom.csv\").query(\"format_blended==0\")\n",
    "result = smf.ols('falsexam ~ format_ol', data=data).fit()\n",
    "result.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's quite amazing. We are not only able to estimate the ATE, but we also get, for free, confidence intervals and P-Values out of it! More than that, we can see that regression is doing exactly what it supposed to do: comparing \\\\(E[Y|T=0]\\\\) and \\\\(E[Y|T=1]\\\\). The intercept is exactly the sample mean when \\\\(T=0\\\\), \\\\(E[Y|T=0]\\\\), and the coefficient of the online format is exactly the sample difference in means \\\\(E[Y|T=1] - E[Y|T=0]\\\\). Don't trust me? No problem. You can see for yourself:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "format_ol\n",
       "0    78.547485\n",
       "1    73.635263\n",
       "Name: falsexam, dtype: float64"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(data\n",
    " .groupby(\"format_ol\")\n",
    " [\"falsexam\"]\n",
    " .mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected. If you add to the intercept the ATE, that is, the parameter estimate of online format, you get the sample mean for the treated: \\\\(78.5475 + (-4.9122) = 73.635263\\\\).\n",
    "\n",
    "## Regression Theory\n",
    "\n",
    "I don't intend to dive too deep into how linear regression is constructed and estimated. However, a little bit of theory will go a long way in explaining its power in causal inference. First of all, regression solves a theoretical best linear prediction problem. Let \\\\(\\beta^*\\\\) be a vector of parameters:\n",
    "\n",
    "$\n",
    "\\beta^* =\\underset{\\beta}{argmin} \\ E[(Y_i - X_i'\\beta)^2]\n",
    "$\n",
    "\n",
    "Linear regression finds the parameters that minimise the mean squared error (MSE). \n",
    "\n",
    "If you differentiate it and set it to zero, you will find that the linear solution to this problem is given by\n",
    "\n",
    "$\n",
    "\\beta^* = E[X_i'X_i]^{-1}E[X_i' Y_i]\n",
    "$\n",
    "\n",
    "We can estimate this beta using the sample equivalent:\n",
    "\n",
    "$\n",
    "\\hat{\\beta} = (X'X)^{-1}X' Y\n",
    "$\n",
    "\n",
    "But don't take my word for it. If you are one of those that understand code better than formulas, try for yourself:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-4.9122215 , 78.54748458])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = data[[\"format_ol\"]].assign(intercep=1)\n",
    "y = data[\"falsexam\"]\n",
    "\n",
    "def regress(y, X): \n",
    "    return np.linalg.inv(X.T.dot(X)).dot(X.T.dot(y))\n",
    "\n",
    "beta = regress(y, X)\n",
    "beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The formulas above are pretty general. However, it pays off to study the case where we only have one regressor. In causal inference, we often want to estimate the causal impact of a variable \\\\(T\\\\) on an outcome \\\\(y\\\\). So, we use regression with this single variable to estimate this effect. Even if we include other variables in the model, they are usually just auxiliary. Adding other variables can help us estimate the causal effect of the treatment, but we are not very interested in estimating their parameters. \n",
    "\n",
    "With a single regressor variable \\\\(T\\\\), the parameter associated to it will be given by\n",
    "\n",
    "$\n",
    "\\beta_1 = \\dfrac{Cov(Y_i, T_i)}{Var(T_i)} \n",
    "$\n",
    "\n",
    "If \\\\(T\\\\) is randomly assigned, \\\\(\\beta_1\\\\) is the ATE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-4.91222149822695"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kapa = data[\"falsexam\"].cov(data[\"format_ol\"]) / data[\"format_ol\"].var()\n",
    "kapa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we have more than one regressor, we can extend the following formula to accommodate that. Let's say those other variables are just auxiliary and that we are truly interested only in estimating the parameter \\\\(\\kappa\\\\) associated to \\\\(T\\\\).\n",
    "\n",
    "$\n",
    "y_i = \\beta_0 + \\kappa T_i + \\beta_1 X_{1i} + ... +\\beta_k X_{ki} + u_i\n",
    "$\n",
    "\n",
    "\\\\(\\kappa\\\\) can be obtained with the following formula\n",
    "\n",
    "$\n",
    "\\kappa = \\dfrac{Cov(Y_i, \\tilde{T_i})}{Var(\\tilde{T_i})} \n",
    "$\n",
    "\n",
    "where \\\\(\\tilde{T_i}\\\\) is the residual from a regression of all other covariates \\\\(X_{1i} + ... + X_{ki}\\\\) on \\\\(T_i\\\\). Now, let's appreciate how cool this is. It means that the coefficient of a multivariate regression is the bivariate coefficient of the same regressor **after accounting for the effect of other variables in the model**. In causal inference terms, \\\\(\\kappa\\\\) is the bivariate coeficiente of \\\\(T\\\\) after having used all other variables to predict it.\n",
    "\n",
    "This has a nice intuition behind it. If we can predict \\\\(T\\\\) using other variables, it means it's not random. However, we can make it so that \\\\(T\\\\) is as good as random once we control for other available variables. To do so, we use linear regression to predict it from the other variables and then we take the residuals of that regression \\\\(\\tilde{T}\\\\). By definition, \\\\(\\tilde{T}\\\\) cannot be predicted by the other variables \\\\(X\\\\) that we've already used to predict \\\\(T\\\\). Quite elegantly, \\\\(\\tilde{T}\\\\) is a version of the treatment that is not associated with any other variable in \\\\(X\\\\).\n",
    "\n",
    "By the way, this is also a property of linear regression. The residual always orthogonal or uncorrelated with any of the variables in the model that created it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Orthogonality imply that the dot product is zero: [-5.54223334e-13 -1.84741111e-12]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>format_ol</th>\n",
       "      <th>e</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>format_ol</th>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>1.963734e-16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>e</th>\n",
       "      <td>1.963734e-16</td>\n",
       "      <td>1.000000e+00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              format_ol             e\n",
       "format_ol  1.000000e+00  1.963734e-16\n",
       "e          1.963734e-16  1.000000e+00"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e = y - X.dot(beta)\n",
    "print(\"Orthogonality imply that the dot product is zero:\", np.dot(e, X))\n",
    "X[[\"format_ol\"]].assign(e=e).corr()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And what is even cooler is that these properties don't depend on nothing! They are mathematical truths, regardless of how your data look like.\n",
    "\n",
    "## Regression For Non-Random Data\n",
    "\n",
    "So far, we worked with random experiment data but, as we know, those are hard to come by. Experiments are very expensive to conduct or simply infeasible. It's very hard to convince McKinsey & Co. to randomly provide their services free of charge so that we can, once and for all, distinguish the value their consulting services brings from the fact that those firms that can afford to pay them are already very well off.\n",
    "\n",
    "\n",
    "For this reason, we shall now delve into non random or observational data. In the following example, we will try to estimate the impact of an additional year of education in hourly wage. As you might have guessed, it is extremely hard to conduct an experiment with education. You can't simply randomize people to 4, 8 or 12 years of education. In this case observational data is all we have.\n",
    "\n",
    "First, let's estimate a very simple model. We will regress log hourly wages on years of education. We use logs here so that our parameter estimates have a percentage interpretation. With it, we will be able to say that 1 extra year of education yields a wage increase of x%. \n",
    "\n",
    "$\n",
    "log(hwage)_i = \\beta_0 + \\beta_1 educ_i + u_i\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>    5.9995</td> <td>    0.094</td> <td>   63.640</td> <td> 0.000</td> <td>    5.814</td> <td>    6.185</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>educ</th>      <td>    0.0596</td> <td>    0.007</td> <td>    8.757</td> <td> 0.000</td> <td>    0.046</td> <td>    0.073</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wage = pd.read_csv(\"./data/wage.csv\").dropna()\n",
    "model_1 = smf.ols('np.log(wage) ~ educ', data=wage).fit()\n",
    "model_1.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimate of \\\\(\\beta_1\\\\) is 0.0529, with a 95% confidence interval of (0.040, 0.066). This means that this model predicts that wages will increase about 5.3% for every additional year of education. This percentage increase is inline with the belief that education impacts wages in an exponential fashion: we expect that going from 8 to 9 years of education to be more rewarding than going from 2 to 3 years."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAE0CAYAAABeuCYQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3gU5fr/8femkr4hJIFAElpooSMkdIIQhEAoUkQU5SBS9GtBqUcQOCigFEGlC6KggBAIHZSSSEnkSO+hg9S0Ta+7vz9ysj+W3YQNJNmU+3VdXJqZ2Z07k9397DPzPM8o4uPjNQghhBBljJmpCxBCCCGKggScEEKIMkkCTgghRJkkASeEEKJMkoATQghRJknACSGEKJMk4MQLSUhIYOLEiTRp0oRKlSqhVCr5888/i7WGRo0aoVQqi3WfRUGpVBIUFGTqMkQhWrduHUqlknXr1pm6lHKpTAecUqksEx98BXHr1q1i/aD8/PPPWbp0KZ6ennz00UdMmDABLy+vfB8TFBSk/dvk9W/ixInFUn9xUiqVNGrUyNRllBq5r+VnHbPcLzi3bt0qpspKhkGDBqFUKtmzZ4/B9Z07d0apVDJo0CCD63/++WeUSiWjRo0qyjJNysLUBYjSbe/evdjb27NlyxYsLS0L9NjBgwfnGYYtW7YsjPJKlb/++gsbGxtTlyFKiU6dOrF3717CwsJ45ZVXdNapVCpOnz6NQqHgyJEjZGVlYWGh+3Gfe6alY8eOxVZzcZOAEy/k/v37VKtWrcDhBvD666/Tvn37IqiqdKpTp46pSxClSG4whYWF6a07fPgw2dnZ9O3bly1btvD333/j5+ens014eLjO85RFZfoUpSFPnsJ79OgR7733Hj4+Pnh4eBAYGMiRI0cASEpKYvLkyTRs2BA3Nzf8/PzYunWr3vPlnmOfNWsWf/31F71798bT0xNPT0/69+/PqVOn9B5z//59Zs+eTWBgIHXq1MHV1ZV69eoxfPhwLl68mGftp06d4p133sHX1xc3Nzd8fHzo3r07P/zwg7aWJk2aAHDkyBGdU36zZs0y6vhcv36dMWPG0KBBA1xdXfHx8eHtt9/m7NmzOtvlnmbUaDTcuXNHu5+iOjWq0WhYvnw5/v7+uLu7U79+fT799FNUKpXB7Z917aNRo0Z5nvoKDQ2lX79+1KxZEzc3N3x9fXn99dc5dOiQdpuMjAyWL19O//79ta8Rb29vgoOD2bt3r87z/fnnn9pT5U8eK6VSyejRo7Xb5XX8EhIS+M9//kPLli1xd3fHy8uLnj17sn37dr1tn3x9x8TE8OGHH1K3bl3c3Nzw9/fnp59+Mvg75+f06dO89dZb+Pj44Orqiq+vL++99x43b97U23bWrFna4x4eHk5QUBDVqlXD09OTAQMG5Pv6Lgrh4eEMGDCAGjVq4ObmRpMmTZgwYQKPHz/W2zb3NW1I7t/w6fdR7mNu3rzJt99+q319vv766wafJysri3r16uHp6UlSUpLBbWbOnIlSqWT16tX5/m4NGjTAzc2Nixcv6v0+4eHhmJmZaU/1Px2Cly9f5sGDB9SpUwcPDw8AduzYwYgRI2jevDkeHh5UrVqVDh06sHjxYrKzsw3WEBUVxRtvvIG3t7f2M3Tv3r06n4tPU6lUfPHFF7Ru3ZoqVapQrVo1XnnlFYOfry+q3LbgVCoV3bp1w9nZmQEDBnDv3j1CQ0N59dVX2bdvHx999BEpKSn06NGDxMRENm/ezLBhw6hatarB02d///03CxYsICAggBEjRnDt2jW2b9/OkSNH2Lp1q863p6NHj7Jw4ULat29PcHAwtra2XLt2jdDQUHbv3s2ePXto3LixzvP//PPPfPzxxwAEBgZSt25d4uLiOHfuHAsXLmT48OE0atSIUaNGaa+JPfkma9eu3TOPycmTJ+nduzcJCQl069YNX19fbty4wfbt29m9ezdr166la9euQE7rq127dsyZMwdHR0ftB/Wzrr89r4kTJ7Js2TLc3d0ZOnQo1tbW7Nq1i7///pvMzMxC28/777/P2rVrcXJyokePHlSpUoV79+4RGRnJhg0b6NSpEwBxcXFMnDgRPz8/AgICqFSpEg8ePGDXrl0MGjSIb775hrfffhvIOSYTJkzQO1bAM68vxcfH88orr3Dp0iUaN27MqFGjUKlUbN26lTfffJPx48czefJkvcflvr6trKwIDg4mPT2d0NBQPvjgA8zMzHjjjTeMOh579uxh6NChqNVqevXqRY0aNTh//jzr1q1jx44dbNu2Tful6kl79+5l9+7ddOnShWHDhnH58mX27dvHiRMniIyMpFKlSkbt/0WsXr2asWPHYmNjQ+/evalcuTKRkZEsW7aMnTt3snv3bjw9PQtlX+PHjycyMpJu3boRGBiIvb29we0sLCx46623mDNnDr/99hvDhg3TWZ+VlcXatWtxcHBgwIABz9xvhw4d2LRpE+Hh4bz66qva5eHh4TRu3Ji6detSp04dwsLCGD9+vHZ9buA92XqbPn06ZmZmvPTSS3h4eKBSqQgLC2Py5MmcOHGClStX6uz78uXLBAYGolKpCAwMpGHDhty6dYs33nhD+znxtHv37tGrVy+uXbtG69atefvtt0lJSWHfvn28/fbbTJgwgUmTJj3z9zZWuQ24c+fOMXLkSGbPno1CoQBg/vz5zJgxg549exIQEMDKlSu1p946d+7MiBEj+Oabbwy2Cv744w++/vprRowYoV0WGhrKW2+9xfvvv89ff/2l3U+HDh24cuUKDg4OOs9x6tQpevTowfTp09m8ebN2+aVLl/j444+pUKECO3bsoGnTpjqPu3v3LgCNGzfGycmJpUuX4uXlVaAXikajYdSoUSQkJLB48WKdcDx06BB9+/Zl1KhRnD17FltbW4YMGQLAnDlzcHJyeq4X5S+//MLhw4cNrnv99dfx9vYG0H4oeXl5cfDgQVxcXACYMmUKvXv35sGDBwXetyFr1qxh7dq1+Pr6sm3bNu1+IOf43Lt3T/uzUqnk7NmzVK1aVec54uPj6datG9OmTWPQoEHY2Njg7e3NpEmTnutYTZs2jUuXLjFkyBC+++477Wto3LhxdO7cma+//ppu3brRokULncedO3eOt99+m3nz5mFubg7AmDFjaNu2LYsWLTIq4JKSkhgzZgyZmZmEhobSoUMH7bqffvqJDz74gFGjRnH06FFtXbl27tzJ1q1bdU5BT58+nQULFrB27Vo++ugjo4+BSqXK9wyEoVb87du3mTBhAra2tvzxxx/Ur19fu27mzJnMnTuXTz75hI0bNxpdR37Onj1LeHi49jWbn9y/y6pVq/QCbufOnTx48IDhw4fnGZJPyg24sLAwbcA9evSIixcv8sEHHwDQvn17fv75Z1JSUrC1tQUMB9zGjRupUaOGzvOr1WpGjRrFxo0bGTlypM6X+08++QSVSsWcOXMYOXKkdvnBgwfp27evwXpHjx7N9evXWblyJf3799cuT0hIoGfPnnz11VcEBQXpfcF/XuXuFGUuOzs7pk6dqvPGHDhwIJBzsGfOnKlzXalfv35YWlrqnarLVbNmTYYPH66zrHfv3vj5+REVFUVkZKR2uaurq164ATRt2pT27dtz+PBhnVbJDz/8QFZWFp988oleuAFUq1bNyN86b5GRkVy+fJnmzZvrnV7p1KkTPXv2JCYmhp07d77wvnL9+uuvzJkzx+C/27dva7fL/ULxySef6ISOtbU1U6ZMKbR6li1bBsA333yjsx8AhUKhE2bW1tZ64QY5wffGG28QHx/PiRMnXqiezMxMNm7ciK2tLdOnT9d5rVatWpWxY8ei0WgMnna0tbVl5syZ2nADqFevHv7+/ly5coXExMRn7n/Xrl3ExsbSu3dvnXADGDp0KE2bNuXixYv89ddfeo/t37+/3vXV3BZtQY9LQkJCnq+TOXPmkJCQoPeYjRs3kpGRwfDhw3XCDXK+HFSpUoV9+/bpfGl5Ef/3f/9nVLgBVKlShZ49e3L27Fn++9//6qzLPS35dPDlxdB1uKevrbVv356MjAwiIiKAnNA6cuQI5ubmOmd2ng43ADMzM8aMGQPAgQMHtMvv3LnD4cOH8fb25p133tF5TEBAAAEBAXrPdf78ecLCwggKCtIJNwBHR0cmTpyIRqPht99+M+p3N0a5bcHVqlULOzs7nWWVK1cGcj6knj51YW5ujqura55viNatW2Nmpv99oU2bNkRGRnLmzBn8/f21y/fu3cuqVas4deoUMTExZGVl6TwuJiZGW0/umyAwMLCAv6XxTp8+DaD3QZarU6dObN++ndOnTxt16sQY27dvN6qTSW5tbdu21Vvn7++PhYWF3vErqKSkJC5cuICzs7PRPTgvXrzIokWLOHr0KA8ePCA9PV1n/f3791+opitXrpCSksJLL71k8JRe7unS3OPzpFq1ahlsAeSGskqlMvgl60nPek107NiRU6dOcfr0ab0ODIa+iOXuOz4+Pt/9Ps3T0zPPL5aQc5r3zp07RtdubW2Nv78/W7Zs4cyZM9prUC/ipZdeKtD277zzDlu3bmXVqlXax964cYOwsDD8/Pxo2LChUc/j7e1N9erVuXnzJjdv3qR69eqEhYVhaWmp/bxp164dCoWCsLAwOnfuzOnTp4mPj6dFixY61xxjY2NZtGgR+/bt49atWyQnJ+vs68nXc+7fo2XLljpfonL5+flx8OBBnWW5X/ITExMNtshjYmKAnNd9YSm3AWfozZ3bjTavN765uXmeH6Rubm4Gl7u6ugLofMtcunQpEydORKlUEhAQgKenJxUqVEChULBz507OnTun82GZewrGUIuhsOTWl9fv4e7urrNdccrdZ+6xfJK5uTkVK1bk0aNHhbIPYz/sjh8/TnBwMFlZWXTs2JHu3bvj4OCAmZkZZ8+eZdeuXXqB97w1Pc/fxNHR0eBjcj+M8uo0UJT7z31/GbPvF1Xcr+e89pOXdu3a0aBBA7Zs2cKXX36p7VSi0WiMbr3l6tSpEz/++CPh4eFUr16d8PBwXnrpJe0X+EqVKlG/fn1tK89Q78n4+HgCAgK4desWLVq04LXXXsPZ2Rlzc3NUKhVLly7VeT3nngEw9J4Ew8cjNjYWyGltGur5mevpYH0R5TbgClteH7C5vZty3/BZWVnMmjULd3d3wsLCtK20XMePH9d7DicnJyDnAm1RDVzPrS+v3+Phw4c62xWn3H0+fvxYeyxyZWdna984T8ptTef1YapSqXSeK/f/jW11zZ07l9TUVIOt0Pnz57Nr1y6jnic/pv6bmHr/L+J5as99zRgaM5ZXb91cT1+DNMbw4cP55JNP+PXXXxk+fDi//PILzs7O9OnTp0DP07FjR3788UfCwsLo0KEDt27d0hvc3b59e1asWEFcXJzB628///wzt27dMtjJ46+//mLp0qU6y3IbAYZ6o4Lh4557rGfOnMn7779foN/xeZXba3CFLSIiArVarbf86NGjANqLpjExMahUKlq1aqUXbklJSQZPN+WeMtu3b98z68j9hm6olvzk9oTLa5qt3DeFoVNPRe3JoQ9Pi4iIMNiqzv0ikNsB50nXrl3T++ZuZ2dHgwYNiI2N1bsuYsj169dxdnY2eIrVUJ2Q8wFakL9LnTp1sLW15cKFC9rTN08q6r/Js14TuS0BU7wmniW/2tPT07Wny57sAZrfa+bkyZOFXuOgQYNwdHTkxx9/ZNu2bURHRzNkyBAqVKhQoOfp0KEDCoWC8PBw7Wvi6VOz7du3R61Wc+DAASIiIqhQoYLOaeXr168DEBwcrPf8hl7PuZ9nx48fN/gl8sk+B7latWoFwLFjx4z91V6YBFwhuXbtmnY8Wq7Q0FAiIyPx8fHRvphcXV2xtbXl5MmTOuNgMjMzmThxosEPsuHDh2Npacm8efMMXov4559/tP/v7OyMQqHQuybxLH5+ftStW5e///6bDRs26KwLCwtj+/btuLi40KNHjwI9b2HI7fQyb948ndZaeno6//nPfww+pnnz5piZmbFx40ad45ycnMy4ceMMPiZ3yqKPP/6YuLg4vfVPXn/18vLSDtN40k8//cT+/fsNPr+LiwvR0dGkpqYaXP80S0tLBg0aREpKCtOnT0ej0WjX3b9/nwULFqBQKIzu8l9QQUFBVKxYkdDQUL0PuXXr1nHy5Enq169fImedGThwIFZWVvzwww9613Tmz5/PvXv3CAwMpEqVKtrlub/HqlWrdLY/c+aMXgumMNjb2zNo0CAuX77MZ599hkKhKPDpSch5Xfn6+vL48WOWLVuGjY2NNkxytWvXDjMzM+bNm0dKSgp+fn46QZo7vOfpLwSnT59mwYIFevusVq0a7dq149atW3rDBw4ePKh3/Q1yvgi1bduWXbt2sWbNGp3Xc66rV68W+LMrP3KKspB06dKFzz77jD/++ANfX1/tODgbGxu+/fZb7SkMMzMzRo4cyYIFC2jTpg09evQgMzOTP//8k7i4ONq3b6/3Iqtbty7z58/no48+IiAggG7dulG3bl1UKhXnz5/n3r17nDlzBshpifj7+3Ps2DEGDRpE06ZNsbCwoE2bNgY7aeRSKBQsWbKEPn36MGrUKLZs2aIdB7dt2zasrKxYunSptptxYchvmECdOnW03Z79/f159913Wb58Oa1btyY4OFg7Ds7JyYnKlSvrDRVwd3dnyJAh/Pzzz7Rv357AwEDS0tLYv38/Xl5eOh9suYYOHUpkZCTr1q2jWbNmBAUFUblyZR4+fEhERAQtW7ZkyZIlQE535/3799O9e3f69OmDo6MjJ0+eJCIigt69exMaGqr3/AEBAWzcuJFXX32VNm3aYG1tTcOGDenevXuex+jzzz/n2LFj/PTTT5w5c4ZOnTppx8HFxcUxfvz4AndwMJadnR2LFy9m6NCh9OnTh+DgYKpXr865c+fYt28fTk5OLFmy5LlOzxU1Ly8v5syZw9ixYwkICKBPnz64u7sTGRnJkSNHqFq1KvPmzdN5TO5QjEWLFnH+/Hl8fX25efMmu3fvJjg4WGfoTmF55513WLFiBQ8ePKBTp07UqlXruZ6nY8eOnDt3jgsXLtCpUyesrKx01iuVSho2bKj9nMjtoJTrtddeY9GiRUyePJnDhw9Tq1Ytrl27xt69e+nVqxchISF6+5w7dy7dunVjwoQJ7N+/n0aNGnHr1i1CQ0Pp0aMHu3bt0ut4t3LlSnr37s2HH37IsmXLaNmyJc7Ozty7d49Lly5x5swZ1q5dW2jjEyXgCkmLFi0YP348M2fOZPny5UDOB9qUKVP0TuH8+9//xsXFhZ9//pkff/wRR0dHOnXqxGeffZbneJ8333yTBg0a8O2333L06FH27duHs7MzPj4+jB07VmfbpUuX8u9//5ujR4/y+++/o1armTBhQr4BBzmtnkOHDvH1119z6NAh9u/fj5OTE0FBQXzyySeFNjYl16+//prnuh49eugMXJ0zZw61a9dm5cqVrFmzhooVK9KzZ0+mTJmS5yD2+fPn4+bmxoYNG1i1ahXu7u4MGDCA8ePH633DzfX999/TuXNnfvzxR3bs2EFqaipubm40bdqU1157Tbtdly5dWL9+PXPnzmXLli2YmZnRokULtm/fzs2bNw0G3OzZszEzM+PgwYNERkaSnZ3N4MGD8w04pVLJ3r17WbhwIdu2bWPx4sVYW1vTuHFjRo4cafCUUmF65ZVX2LdvH/PnzycsLIzQ0FBcXV0ZPHgw48ePp3r16kW6/xcxbNgwatasybfffsvOnTtJTk6mSpUqvPvuu3z66ad6HSFcXFzYuXMnU6dO5ciRIxw9ehRfX19Wr16Nk5NTkQRc3bp1admyJcePH3+u1luujh078v333wN593pt3769NuCenp6rSpUq7N69m2nTphEREcGBAwfw8fFh3rx5dOzY0WDA1atXj99//50ZM2YQHh7O4cOH8fX1Ze3atVy5coVdu3bpXZ+tUqUKBw8eZMWKFYSGhrJ582YyMzNxc3Ojdu3azJ4926hJKYyliI+P128nCqOtW7eO9957r9BH4Ashyr7k5GQaNGiAjY0N586d0+vcUlqNGDGC3377jZCQEDp37myyOuQanBBCmMiqVatQqVQMGzas1IWbRqMxOItQWFgYISEhuLi4PPOsUVErXUdUCCFKOZVKxfLly3nw4AFr167Fzc2tVN6TLTs7G19fXzp06ICPjw8WFhZcunSJgwcPaju0WFtbm7RGCTghhChG8fHxfPHFF1hbW9OkSRPtHKWljbm5OSNGjODw4cPaXuFKpZLu3bvzwQcf6MzcZCpyDU4IIUSZJNfghBBClEkScEIIIcokCTghhBBlUpkMuKioKFOXUCrIcTKOHCfjybEyjhwn47zocSqTASeEEEJIwAkhhCiTJOCEEEKUSRJwQgghyiQJOCGEEGWSBJwQQogyyWQBt2LFCtq0aYOnpyeenp507dqVvXv3atePHj0apVKp869Lly6mKlcIIUQBLL2QxFVVpklrMNlkyx4eHkyfPp1atWqhVqv59ddfGTJkCIcOHaJhw4ZAzl1nly1bpn3M03epFUIIUfJ8dy6Rz44n4GZjxtZulWjgbGmSOkwWcEFBQTo/T5kyhR9++IHjx49rA87a2hp3d3dTlCeEEOI5rLqUzGfHEwB4lKomaPdjtgRWomml4m+glIhrcNnZ2WzevJnk5GRatWqlXX7s2DFq165NixYt+OCDD3j8+LEJqxRCCJGfX6+mMPZYvM6yjGxIzzbNTWtMej+48+fPExgYSFpaGnZ2dqxduxZfX18AunTpQq9evfD29ub27dvMnDmT4OBgDh06ZPKb6AkhhNC19UYq7x2O01lmbQ6/vOyCn7tpPrNNej+4jIwM7t69i0qlYtu2baxZs4YdO3bQoEEDvW3v379Po0aNWLVqFcHBwXk+p8zxJoQQxevPWDPGXbQmW6PQLrNQaJhbP522FdUv9Nw+Pj7P/ViTtuCsrKyoWbMmAM2aNePEiRMsXryY7777Tm/bKlWq4OHhwfXr1/N9Th8fH6Kiol7ooJQXcpyMI8fJeHKsjFOWjtOhe2lMOhbDk2chzRTwQycXele3eaHnLlOTLavVajIyMgyui4mJ4f79+9LpRAghSohjD9N5fX8s6dn/f5kCWNLe+YXDrTCYrAU3bdo0AgMDqVq1KklJSWzatInDhw+zceNGkpKSmD17NsHBwbi7u3P79m1mzJiBq6srPXv2NFXJQggh/udkdAYDf48hJUv3KteCNkoG1bI1UVW6TBZwDx8+5N133+XRo0c4Ojri6+vLpk2bePnll0lNTeXChQusX78elUqFu7s77du3Z/Xq1Tg4OJiqZCGEEMD52Ez67YsmMVM33L5s5cTbde1MVJU+kwXckiVL8lxnY2NDSEhIMVYjhBDCGFGqTPrsjSYuXTfcPmvuyBhfexNVZViJugYnhBCi5LqZmEXvPdE8TtPtGTm2sT2fNil5Z9ck4IQQQjzTP8nZ9N4Tzb0U3XAbWd+OKc0dTVRV/iTghBBC5OtRajZ99kZzKylbZ/nQOrbM9nNCoVDk8UjTkoATQgiRp7h0NX32RhOlytJZPqCmDQtaK0tsuIEEnBBCiDwkZKh5dV80F+J0w62nVwWWtHfG3KzkhhtIwAkhhDBAlaGm375oTkTr3tOtS1VrfuhUEYsSHm5g4qm6hBBClDzx6Tktt7+fCrd2la34ubML1uYlP9xAAk4IIcQT4tPV9N0Xzcmnwq2VqxW/dnHBxqJ0hBtIwAkhhPif3A4lp2N0w83PzYrfurrgYFm6rmpJwAkhhCAuXU3vPdGcidUNt9buVmwsheEGEnBCCFHuxaZl03tvDGcNhNtvXV2wL4XhBhJwQghRrsWkZRO8J5rzTw0FaFvZig1dSm+4gQwTEEKIcis6LZteBsKtXWUrNpbycANpwQkhRLn0ODVnbskL8brh1qGKNeu7VMTWonSHG0gLTgghyp1HqTktt6fDrZNH2Qk3kBacEEKUKw9Tcq65XX5qbskAD2t+ebl0jXN7lrIR00IIIZ7pQUpOy+3pcHu5atkLN5AWnBBClAv3U7LptTuaqwm64dalqjVrO7tQoYyFG0jACSFEmXcvOZteex5zLUH3fm6B1az5KaBshhtIwAkhRJl2JymL3nuiuZ6oG27dqlnzUymaOPl5SMAJIUQZdVWVSZ+9MdxN1g23VzwrsCagYpkON5CAE0KIMulsbCb99kbzOE2ts7yHVwV+7FQRqzIebiABJ4QQZU7kw3QG/hGDKkOjs7x39Qqs6FA+wg1kmIAQQpQph+6l0Xeffri94WPLqo7lJ9xAAk4IIcqMHbdSGfh7DClZuuE2uoEdi9oqMTcrP+EGcopSCCHKhPVXU3jvcBzZutnGxKYOTGjqgEJRvsINJOCEEKLUW3ExiXERKr3lX7ZyYoyvvQkqKhkk4IQQopTSaDTMP5PEf04k6Cw3U8DCNkrerGNnospKBgk4IYQohTQaDdP+m8DCc0k6yy3NYEWHivSpYWOiykoOCTghhChlstUaPo2IZ/XlFJ3lNuYKfu5ckS7VKpiospJFAk4IIUqRTLWG0X/Gsel6qs5yR0sF67u40KaytYkqK3kk4IQQopRIzdLw9qFY9t5J01nuYm3G5kAXmlayMlFlJZPJxsGtWLGCNm3a4OnpiaenJ127dmXv3r3a9RqNhlmzZlGvXj0qV65MUFAQFy9eNFW5QghhUomZagb+Hq0XblVszdjVo5KEmwEmCzgPDw+mT59OWFgYBw8epEOHDgwZMoRz584BsHDhQr7//nvmzJnDgQMHcHV1pW/fviQmJpqqZCGEMInotGx674nmzwcZOstrOJizu4crdZWWJqqsZDNZwAUFBdG1a1dq1qxJ7dq1mTJlCvb29hw/fhyNRsOSJUv46KOP6N27Nw0aNGDJkiUkJSWxadMmU5UshBDF7mZiFt12PuZEdKbO8vpKC3b3cKW6g1xpykuJmKorOzubzZs3k5ycTKtWrbh16xYPHz6kc+fO2m1sbGxo06YNkZGRJqxUCCGKz6noDLru0L9RafNKluzsXonKtuYmqqx0MGn0nz9/nsDAQNLS0rCzs2Pt2toMxqEAACAASURBVLX4+vpqQ8zV1VVne1dXV+7fv5/vc0ZFRen8V+RPjpNx5DgZT46VcZ51nI7FmTHhojWpat0ptlo6ZfN17RRi7qiIKcoCywCTBpyPjw9//vknKpWKbdu2MXr0aHbs2KFd//TcaRqN5pnzqfn4+BAVFYWPj0+R1FyWyHEyjhwn48mxMs6zjtOvV1MYeyGOp+ZMZmBNG75r51xu7gjwol+WTHqK0srKipo1a9KsWTM+//xzGjVqxOLFi3F3dwfg0aNHOttHR0frteqEEKKsyJl6K5HRf+qH24cN7VnaofyEW2EoEdfgcqnVajIyMvD29sbd3Z2DBw9q16WlpXHs2DH8/PxMWKEQQhSNbLWGcREqZvytO6+kApjt58T0lk6YlcM7ArwIk52inDZtGoGBgVStWlXbO/Lw4cNs3LgRhULB6NGjmTdvHj4+PtSuXZu5c+diZ2dH//79TVWyEEIUidQsDSPCYtlxW3eMm5UZLJd5JZ+byQLu4cOHvPvuuzx69AhHR0d8fX3ZtGkTL7/8MgAffvghqampjBs3jvj4eFq0aEFISAgODg6mKlkIIQpdXLqawX/EEPFId4ybo5WCX152oZ1MvfXcTBZwS5YsyXe9QqFg0qRJTJo0qZgqEkKI4nUnKYv++2K4rMrSWV7V1pzfAl1o4CwDuF+EjBAUQggTOBebyYDfo7mfotZZXl9pwW9dXahmLx/PL0qOoBBCFLPj8WZMjHxMQqZuV8k27lb88rILSusS1f+v1CpwwN25c4cjR47w+PFj+vbtS7Vq1cjKyiIuLg5nZ2csLCQzhRAiL5uvp/DBeWuyNLrh1qe6DUvbO1PBQnpKFpYCpdHkyZNZvnw52dnZKBQKGjduTLVq1UhJSaF58+ZMnDiR9957r6hqFUKIUkuj0fDduSSm/DeBnM7//9/I+nbM8pNhAIXN6HbwokWLWLJkCe+99x5bt25F88S3D0dHR4KCgnRmIRFCCJEjU63h46Px/ws3XTNecmS2hFuRMLoFt2bNGgYOHMj06dOJjY3VW+/r68uBAwcKtTghhCjt4tPVvHUwlrD76TrLLc3g+3bODKxla6LKyj6jW3B3796lTZs2ea53cHBApVIVSlFCCFEWXE/IouvOx3rhZmeu4beuLhJuRczoFlzFihV58OBBnuvPnz9PlSpVCqUoIYQo7Y48SOeNAzHEpet2JvGyN+crnyQ6eVQwUWXlh9EtuMDAQNasWUNMjP4NGk6fPs3atWsJCgoq1OKEEKI0+iUqmT57o/XCzc/NigO9XKllp8njkaIwGR1wkydPxszMjDZt2jBt2jQUCgXr1q3jX//6F127dsXDw4Nx48YVZa1CCFGiqTUaZvytYszheDJ1x28zsKYNod0qUamC3KS0uBgdcO7u7hw6dIhXXnmF7du3o9Fo+O233/jjjz8YNGgQ+/btQ6lUFmWtQghRYqVk5XQmmX8mSW/d5GYOLOsgY9yKW4HGwVWqVImFCxeycOFCoqOjUavVVKpUCTMzGXUvhCi/7qdkM/iPGE7FZOostzaHJe2c6VdTOpOYwnNPO1KpUqXCrEMIIUqlMzEZDP4jln9SsnWWu1Yw45eXXWjpZmWiyoTRATdnzpx81ysUCipUqICHhwdt27aVHpVCiDJv1+1URoTFkfzU7bcbOFuwvosLXjJhskkZffRnz56N4n8j7TVPzaH29HJzc3OGDRvGV199pV0nhBBlRe60W1P/m8DT/SEDq1mzsmNFHK3k0o2pGf0XuHTpEo0aNWLgwIEcPHiQ27dvc/v2bQ4cOMCAAQNo1KgRJ06c4NChQ/Tr148ffviBhQsXFmXtQghR7DLVGj7837RbT4fbqAZ2/PKyi4RbCWH0X2HChAlUr16dpUuX0rRpUxwcHHBwcKBZs2YsW7YMb29vZsyYQZMmTVi+fDkdO3Zk3bp1RVm7EEIUq+i0bPrsjeanKyk6y80VMK+1E7P9lFiYyVmrksLogDtw4AAdO3bMc33Hjh3Zv3+/9udu3bpx+/btF6tOCCFKiFPRGXTa9pgjDzJ0ljtaKvitqwvD69mbqDKRF6OvwSkUCs6fP5/n+vPnz+tcb9NoNNjY2LxYdUIIUQKsv5rCR0fjSNPtKIm3vTkburpQT2lpmsJEvoxuwQUFBfHjjz/y7bffkpycrF2enJzMokWLWLNmjc5UXZGRkdSqVatwqxVCiGKUpdYwKTKeUX/qh1u7ylbs7+Uq4VaCGd2C+/LLL7lx4wZTp05l+vTpuLu7A/Dw4UOys7Np1aoVX375JQBpaWk4OjoyevTooqlaCCGKWExaNsMOxRH+1J0AIKczyX9aOmEp19tKNKMDTqlUsmvXLnbs2MH+/fu5c+cOAF27dqVLly4EBQVpT1FWqFCBRYsWFU3FQghRxE7HZPDGgVjuJOk226zN4Zs2zgyuLTOTlAYFGoWoUCjo1asXvXr1Kqp6hBDCpDZdT+H/DseTmq07CKCqrTlrX65Is0oyM0lpIcPshRCCnOtt0/6bwHfn9SdLbuNuxZqAirjayJ0ASpMCBdylS5dYunQpp06dQqVSoVbr3g9CoVBw6tSpQi1QCCGKWmxaNv8Ki+PQPf3rbSPq2/FlK7neVhoZ3YsyMjKSgIAAdu7cibu7Ozdv3qR69epUqVKFO3fuYGdnR5s2bYqyViGEKHRnYzMJ2P5YL9yszOC7dkq+9ldKuJVSRrfgZs6ciYeHB/v37yc7O5vatWszduxYOnbsSGRkJAMHDmTmzJlFWasQQhSqkOspvGfgepuHrRk/d3ahhatcbyvNjG7BnTx5kqFDh6JUKrX3f8s9Renn58dbb73FF198UTRVCiFEIcpWa/j8uIp/hcXphVtrdysO9nKTcCsDjA44hUKBk5MTALa2OV1kY2Njtetr167NxYsXC7k8IYQoXI9Ss+m3L4aF5/Q7k7xTz47QbpVwt5XOJGWB0QHn5eXF9evXAbC2tsbb25uDBw9q1x89epSKFSsWfoVCCFFIjjxIp0PoI8Lu619vW9RWydzWSqzM5XpbWWF0wAUEBBAaGqq959tbb73FunXrCA4OplevXmzYsIEBAwYUWaFCCPG81BoN35xJJHhPNA9SdXt/V7E1Y2d3V4bWsTNRdaKoGB1wn376KWvWrCErKwuAjz76iClTphAXF0diYiITJ05k8uTJRu94/vz5BAQE4OnpSa1atRg0aBAXLlzQ2Wb06NEolUqdf126dDF6H0IIEZeuZvD+WKb9ncBTl9toWznneltLN7neVhbl24ty0qRJtGnTBn9/f1xdXWnatKl2nUKhYOzYsYwdO/a5dnz48GGGDx9O8+bN0Wg0fPnll/Tp04fIyEicnZ2123Xq1Illy5Zpf7aykheiEMI4Jx5n8NYh/Sm3AD5pbM+kZo5y/7YyLN+AW7p0qTZcatasSevWrfH396dNmzbUqFHjhXYcEhKi8/OyZcvw8vIiIiKC7t27a5dbW1trJ3YWQghjaDQaVlxM5t/HVWTqnpHE2VrBsvYVCfSsYJriRLHJN+DOnTvHsWPHOHbsGJGRkfzyyy+sXbsWhUKBm5sb/v7+2tBr3Lixzv3gCiopKQm1Wo1SqdRZfuzYMWrXro2TkxNt27ZlypQpuLq6Pvd+hBBlW2Kmmg8Ox7PlZqreupdcLVndqSKe9jJLYXmQ71+5atWq9O/fn/79+wOQkJDAX3/9xbFjx4iIiGDfvn2EhoaiUCiwt7enVatWbNq06bkKmThxIo0aNaJVq1baZV26dKFXr154e3tz+/ZtZs6cSXBwMIcOHcLa2vq59iOEKLvOx2by1sFYriZk6a0bWT/nFjfSS7L8UMTHx2uevZlhWVlZ/P777yxcuJDIyEgUCoXO2DhjTZ48mZCQEPbs2UP16tXz3O7+/fs0atSIVatWERwcbHCbqKioAu9fCFH6bX9ozpxrVqSrdQPMzlzDFJ8MXq6kfx1OlHw+Pj7P/dgCt9OvX7+ubcFFRERw7do1FAoFDRs2xN/fv8AFTJo0iZCQELZv355vuAFUqVIFDw8P7Xg8Q3x8fIiKinqhg1JeyHEyjhwn45niWKVkqRkXoWJdVIreOl9nC34KcKGWU8k6JSmvKeO8aIMl37+6Wq3m9OnT2kCLjIzk8ePH2NnZ0aJFC/r164e/vz8vvfQSDg4OBd75hAkTCAkJYceOHdSpU+eZ28fExHD//n3pdCKEAOCqKpOhB2O5EKd/SvJNH1u+8ldiYyGnJMurfAPO29ub5ORkqlatip+fH59++imtWrWiUaNG2vkon9enn37Khg0bWLt2LUqlkocPHwJgZ2eHvb09SUlJzJ49m+DgYNzd3bl9+zYzZszA1dWVnj17vtC+hRCl32/XUhh7LJ7ETN2rLDbmCua1duJ1Hxm4Xd7lG3BJSUlYWFjg7e1N9erVqVGjBtWrV3/hcANYuXIlAL1799ZZPmHCBCZNmoS5uTkXLlxg/fr1qFQq3N3dad++PatXr36u1qIQomxIzFTz6bF4NlzT7yXp42TBj50q4lvR0gSViZLGqGECERER7NmzhwULFgBQt25d/P398fPzw9/fH29v7wLvOD4+Pt/1NjY2emPlhBDl29+PM3gnLJYbifodRl6tYcM3bZU4WL74F3BRNhR4mEBkZCSRkZFERESwfv160tLScHd314bdqFGjiqVwIUT5odZoWHg2iS9OJJD1VL9va3P4oqUTw+vZvdBYXFH2FKhrkaOjI127dqVr166A7jCB0NBQtm3bJgEnhChU91OyGRkeR/hTdwAAqK+0YGVHOSUpDCuUYQK5dxjw8vIq9AKFEOXXrtupvH84nth0td66d+rlDNyWXpIiL881TECj0WBmZkaDBg145513aN26Na1bt6Zy5crFVbcQogxLzdIw5biKlZeS9dY5Wyv4rq0zQd42JqhMlCb5BpyXlxcpKSloNBpsbGxo1qwZb775Jq1bt6ZVq1bSm1EIUeguxGXyzqFYLsTrj21rX9mKZR0q4mEnd9wWz5ZvwLVr1057u5xmzZphaSnnuYUQRUOj0fDDpWQ+O64i7alOkhYK+HdzRz5oaI+53N5GGCnfgFu/fn1x1SGEKMdi0rJ5/3A8u++k6a2r7mDODx0r0sJV7gUpCqZkTdAmhCh3wu6lM+rPWO6n6HckGVjLhrn+ShytZGybKDgJOCGESaRmafjPCRVLzifz9C1NHCwVzG2tZFAtW5PUJsoGCTghRLE7GZ3BqPA4Lqv0O5K0qGTJyo4VqeEoH0/ixcgrSAhRbDLVGuaeTmTu6USyn2q2KYCPG9szqZkjltKRRBQCCTghRLG4FJ/JqPA4TsVk6q3ztDdncTtn2lexNkFloqwy+srt1KlTuXDhQlHWIoQog7LVGr49l0jHbY8MhtsbPrYc6e0m4SYKndEBt2zZMtq1a0e7du347rvvePDgQVHWJYQoA24mZtFrTzRTjieQ/tTYNjcbM359uSLftXOWXpKiSBj9qrp8+TLz5s3DwcGBqVOn0rBhQ/r168fGjRtJSdG/VbwQovzSaDSsuZxMu62POPowQ2997+oVONbHje5eMt2WKDpGB5xSqWTYsGHs3r2bU6dOMWHCBP755x9GjhxJ3bp1GT16NGFhYdqJl4UQ5dODlGwG/RHDh0fjSXrq3jZOVgpWdHDmx04Vcakg022JovVc5wW8vLwYN24ckZGRHDhwgJdffpkNGzbQt29ffH19mT59Onfv3i3sWoUQJdy+x+a03vqQfXf1b23zclVrjvVxZ0AtW7lvmygWz92LMjk5me3bt7NhwwbCw8OxtLSkW7duWFlZsXTpUpYuXcrixYvp27dvYdYrhCiBYtOy+TRCRcgNa3hq2LadhYKZLZ14u64EmyheBQo4tVrNgQMH2LBhA7t27SIlJYXmzZszZ84c+vfvj1KpBCAmJoZ//etfTJ06VQJOiDJu281UxkXE8zBVf6qt1u5WLG7nLIO2hUkY/aqbNGkSISEhPH78mCpVqvDuu+8yePBg6tSpo7eti4sLr7/+utzdW4gy7GFKNuMi4tl2S3+CZCszmNLckTG+Mvu/MB2jA27NmjX07NmTwYMH06lTp2eeavD39+f7779/4QKFECWLRqPh16spTP5LRXyGfqeyxhUtWdrBmQbOcnstYVpGB9zly5cLdINTb29vvL29n6soIUTJdCsxi4+PxnPgnn4nEgsFvF0tky8DPLAyl1abMD2jA07u3i1E+aXWaFhxMZkZfyeQnKXfamtWyZJv2zpjHXNTwk2UGHkG3HvvvVfgJ1MoFHz33XcvVJAQomS5Ep/JB0fiiXikP2C7gnnOnbZHN7DHwkxBVIwJChQiD3kGXHh4eIG79EoXYCHKjky1hkVnk5hzKoEM/Q6StK1sxbdtnakpPSRFCZXnK/Ps2bPFWYcQogQ5FZ3B+0fiORerPzmyg6WC/7R0YmgdW8zkS60owYz66pWenk5ISAh16tShRYsWRV2TEMJEUrM0zDmVwLfnkvTu1wbQzbMC81srqWon02yJks+oqbqsra358MMPpVUnRBl25EE67UMf8c1Z/XBzsTZjZUdn1r9cUcJNlBpGnzz38fHh4cOHRVmLEMIEYtKymfrfBNZFGb4ryICaNszyc6KSTI4sShmjJ1seP348K1as4Pz580VZjxCimKg1Gn66ksxLIQ8NhpuHbc792lZ0rCjhJkolo1tw4eHhuLq60qFDB1q1akWNGjWwsdG9l5NCoWDu3LmFXqQQonCdj83kk2OGu/4DDKtry7SXnHCSG5GKUszogFu1apX2/yMiIoiIiNDbRgJOiJItOVPNnFOJLD6fhIHx2tRTWjCvtZK2la2LvzghCpnRX8/i4uKe+S82NtboHc+fP5+AgAA8PT2pVasWgwYN4sKFCzrbaDQaZs2aRb169ahcuTJBQUFcvHjR+N9OCKG163YqflseseicfrjZmCv4vIUj4cFuEm6izDDZ+YfDhw8zfPhw9u7dy7Zt27CwsKBPnz7ExcVpt1m4cCHff/89c+bM4cCBA7i6utK3b18SExNNVbYQpc7tpCwG/xHD6/tjuZucrbe+WzVrIvq68XFjB5lmS5QpJpuCICQkROfnZcuW4eXlRUREBN27d0ej0bBkyRI++ugjevfuDcCSJUvw8fFh06ZNDBs2zBRlC1FqZKo1LDmfxOxTiaQYOB9Z1dac2f5O9PSqILMQiTLJ6IBzdnY26k1QkNOUT0pKSkKtVmtvmnrr1i0ePnxI586dtdvY2NjQpk0bIiMjJeCEyEfEw3TGHo3nQnyW3jpzBYxuYM/EZg7YW0onElF2GR1w48eP1wu47Oxsbt26xe7du6lduzbdunV77kImTpxIo0aNaNWqFYB2zJ2rq6vOdq6urty/fz/P54mKitL5r8ifHCfjlJbjFJ8J3920IvSh4bd2I4dsJtXOwMcuhfs3HxVJDaXlWJmaHKeiV6A7eufl3r17dOnSxeDdvY0xefJkIiIi2LNnD+bmuuNtng5VjUaTb0vSx8eHqKgofHx8nquW8kSOk3FKw3HKUmtYfTmZL04mGLwJqdJKwfSXnHiziOePLA3HqiSQ42ScF/0SUCjnJzw8PBg2bBhfffVVgR87adIkNm/ezLZt26hevbp2ubu7OwCPHul+y4yOjtZr1QlRnoXfT6fDtkeMizB8h+3BtW053s+dt+rayeTIolwptBPwSqWSGzduFOgxEyZMYNOmTWzbtk2v9eft7Y27uzsHDx7ULktLS+PYsWP4+fkVSs1ClGa3k7J462AMwXuiuRCnf62trpMFO7pXYkl7Z1xtZCYSUf4USi/K6Oho1qxZg5eXl9GP+fTTT9mwYQNr165FqVRqr7nZ2dlhb2+PQqFg9OjRzJs3Dx8fH2rXrs3cuXOxs7Ojf//+hVG2EKVSSpaahWeTWHg2kTT9Xv/YWyj4tIkDY3ztpdu/KNeMDrhevXoZXK5Sqbhy5QqZmZk6s508y8qVKwG0QwByTZgwQXu978MPPyQ1NZVx48YRHx9PixYtCAkJwcHBwej9CFFWaDQaQm+m8dlxlcHxbACDatkw7SUnqthKi00IowNOrVbrde5QKBR4e3sTEBDA0KFDqVWrltE7jo+Pf+Y2CoWCSZMm5dvBRYjy4HxsJhMi4zn8wPDckc0qWTLHz4lWbjILiRC5jA64nTt3FmUdQggD4tLVfHkigR8uJ6M2MHekawUzprZwZIiP3F1biKeZbCYTIUTestUafrySzMwTCcSl6yebhQJGNrBnfFMHmfFfiDwUKOASEhL49ttv2bdvH7dv3wbAy8uLbt268f777+Po6FgkRQpRnoTfT2dSZDznDfSMBHi5qjWzWjlRR2lZzJUJUboYHXAPHjzglVde4datW/j4+NC2bVs0Gg1RUVF8/fXX/Pbbb+zevZvKlSsXZb1ClFmX4jP5/LiKvXfTDa6v4WDOl62ceMVT5o4UwhhGB9y0adN4+PAh69ato0ePHjrrdu/ezb/+9S9mzJjB4sWLC71IIcqyhynZzD6VwJorKQavs9k90e3fWrr9C2E0owNu//79vPvuu3rhBtC9e3dGjBjBL7/8UqjFCVGWJWeq+f58EgvPJpFs6O6jwMBaNkxr4YSHnXT7F6KgjA64xMREqlWrluf6atWqkZSUVChFCVGWZas1/HothS9OJHA/RW1wm9buVvynpRMvuVoVc3VClB1Gd7+qVasW27ZtQ63Wf0Oq1Wq2b99eoHFwQpRHB/9Jo8O2R7x/ON5guNV0MOfnzhXZ1b2ShJsQL8joFtzIkSP58MMP6du3L2PGjNHOhH3lyhWWLl3KkSNHWLhwYZEVKkRpdj42k8//q+KPfwx3IKlobcaEpg4Mq2sn02sJUUiMDrihQ4cSExPDnDlz+PPPP7XLNRoN1tbWTJ06lTfffLNIihSitLqfks2XJxJYd9VwBxJr85ybj37UyAGltYxnE6IwFWgc3Mcff8xbb73FoUOHdMbBderUiYoVKxZJgUKURomZar47l8S355JIyasDSU0bPmvhiJe9zLcgRFEo8DurYsWK9OvXryhqEaLUS8vS8MPlZOafTiQm3XAHkraVrZjZ0olmleQamxBFSb46ClEIstQafrmawpyTifyTYnimfx8nC6a/5Eh3GagtRLHIN+CcnZ0L9EZUKBTExMS8cFFClBZqjYatN1L54mQC1xIMB1ulCmZMaubA0Dp2WJpJsAlRXPINuLffflsn4FJTU1m/fj3dunXDw8OjyIsToqTSaDT8fjed/5xI4GxspsFt7CwUjPG15/8a2uMoEyILUezyDbgFCxbo/BwTE8P69esZPXo0HTt2LNLChCipjj7ICbZjDw3fm83aHIbXs+PjRg642sgMJEKYSoGuwcl1A1GenY7JYObfCfyex1g2cwUM8bFlfBMHqknPSCFMTt6FQjzDzRQFXx6MZcvN1Dy36VvdhsnNHfBxklvYCFFSSMAJkYebiVnMO53IL1EVyMZwuHWtas1nLRxp4iJd/oUoaSTghHhKbrD9ejWFnDHa+qfmW7tbMaW5I20qWxd7fUII4+QbcH///bfOzwkJCQBERUVhb29v8DEtWrQopNKEKF76waavUUVLprZwpEtVa7kmLUQJl2/AdenSxeCbePz48XrLNBoNCoWC2NjYwqtOiGJwMzGLuf8Ltuw8gq22owX/bu5A7+o2mEmwCVEq5Btw33//fXHVIUSxu5GQxdwziazPJ9hqOpjzZuUU/q+NBxYySFuIUiXfgHv99deLqw4his31hJwW24ZreQdbLUdzPm3iyICaNty4dlXCTYhSSDqZiHLD2GAb18SR/jVtJNSEKOUk4ESZdz0hi69PJ7Ixn2Cr7WjBuKYOvFpDgk2IskICTpRZZ2Iy+OZsEltvphq82SjkBNv4/wWbuQSbEGWKBJwoc44+SGfBmcQ8p9SCnFvXjG/iQD8JNiHKLAk4USZoNBr23c0JtohHhidBBqjjlNNi61tdgk2Isk4CTpRqWWoNW2+msuBMIufjsvLcroHSgrFNJNiEKE8k4ESplJal4derKSw8l8jNRMM3GgXwc7Pi48b2BFarIAO0hShnJOBEqZKYqWb1pWS+P5/Ew1R1ntu9XNWasY0daONuJVNqCVFOmfQ2w0eOHOG1116jfv36KJVK1q1bp7N+9OjRKJVKnX9dunQxUbXClKLTspl5IoGGGx8w9b8JBsNNAfSpbsOhXq5sDqxE28oyX6QQ5ZlJW3DJyck0aNCAwYMHM2rUKIPbdOrUiWXLlml/trKS25KUJ1GqTBafT+LXqymk5XEm0tIMXqtly4eN7Kkt92MTQvyPSQMuMDCQwMBAAMaMGWNwG2tra9zd3YuzLGFiGo2Gww8y+O58EnvvpOW5nZ2Fgrfq2vKerwNV7cyLsUIhRGlQ4q/BHTt2jNq1a+Pk5ETbtm2ZMmUKrq6upi5LFIFMtYYtN1L5/nwSp2My89xOaaVgZAN7Rta3o2IFCTYhhGGK+Pj4POZ4KF5Vq1blq6++YsiQIdplmzdvxsbGBm9vb27fvs3MmTNRq9UcOnQIa2vDN5qMiooqrpJFIUnMgq0PLFh/z4JHGXlfFna3VjPYI4u+lbOwlVwTolzw8fF57seW6Bbcq6++qv1/X19fmjZtSqNGjdi7dy/BwcEGH+Pj40NUVNQLHZTywtTH6VZiFksvJPHzlRSS8rrDKNCskiXv+9oTXN0GSxOMYTP1cSpN5FgZR46TcV60wVKiA+5pVapUwcPDg+vXr5u6FPEC/vs4g+/OJbHtVt5zRCqA7l4VeN/XntbS1V8I8RxKVcDFxMRw//596XRSCmWpNey8ncaS80n5TqVlY67gdR9bRjewkx6RQogXYtKAS0pK0rbG1Go1d+/e5cyZMzg7O+Ps7Mzs2bMJDg7G3d2d27dvM2PGDFxdXenZs6cpyxYFEJ2WzZrLKay6lMw/KXnPF8yaOwAAGCtJREFUOOJmY8a79e35V11b6TgihCgUJg24kydP0qtXL+3Ps2bNYtasWQwePJj58+dz4cIF1q9fj0qlwt3dnfbt27N69WocHBxMWLUwxsnoDJZfTCbkRgrpeecaDZQWjGloz4Catliby2lIIUThMWnAtW/fnvj4+DzXh4SEFGM14kVlZGsIvZnK8otJHH+cdzd/gM4e1rzf0J4AD5ltRAhRNErVNThRMj1IyWb15WR+vJyc7/yQthYKBta04d0G9jRwlutrQoiiJQEnnotGo+H445zTkFtvpJJPL3+qO5jzTj073vCxQ2lt0ulPhRDliAScKJDULA0hN1JYfjE539lGIOc05LsN7OhatYLcg00IUewk4IRRrsRnsvpyMr9eTSE+I+/mmr2FgsE+trxb3w4f6eYvhDAhCTiRp/RsDdtvpbL6cjJHHuQ9dg3Ax8mCEfXseK22LY5WchpSCGF6EnBCzzVVFj9eSeaXqBRi0vPuNKIAAj0rMLK+HZ08rOWO2UKIEkUCTgA5Xfx33U5j9eVkwu6n57uts7WCIbXtGF7PjhqO8hISQpRM8ulUzt1MzOKnK8msjUrhUT5d/AFau1sxrK4dwd42VLCQ1poQomSTgCuHMtUa9txJ4/tz1kQefkh+90tyslLwWi1bhtWzo55SOo0IIUoPCbhy5FJ8JuuiUlh/NYXHaWog7zkfW7la8XZdW/rUsMHWQjqNCCFKHwm4Mi4hQ82WG6msjUp+5vRZDpYKBtWy5e26djSsKK01IUTpJgFXBmk0Go48zGDtlWRCb6aRmp3/TdubVbJkWF07Xq1hg52ltNaEEGWDBFwZ8k9yNr9eTWFdVDI3EvOZwh9wtFTQxSWDD1p60LSSVTFVKIQQxUcCrpRLz9aw+3Yaa6OSOXAvPc87ZOfqUMWaIT629PKuwD83ruEj4SaEKKMk4EohjUbDqZhMfr2awm/XU4hLzz/VqtmZ87qPLa/XtqW6g/zJhRDlg3zalSK3k7LYeC2VjddSuKLKyndba3Po6WXDGz62dKhiLZMdCyHKHQm4Ei4+Xc22W6msv5rC0Yf5zwcJ0MTFkjd8bOlf0xZnuTWNEKIck4ArgTLVGv64m8aGa6nsvpNKev79RXC2VjCwpi1DfGxp7CLX1IQQAiTgSgyNRsPf0ZlsuJZCyPXUfCc5BrA0g67VKjColi2veFbA2lxOQQohxJMk4EzsZmIWG6+lsPFaKlcT8r+uBtDS1ZJBtWzpV8OGihXynolECCHKOwk4E7ifks2WG6mE3Ejhv8+YXQSguoM5A2vZMqimLbWc5E8mhBDGkE/LYhKTls22m2lsupHC0QcZ+U5wDKC0UtCvhi0Da9ng52aFQu61JoQQBSIBV4RUGWp23kol5EYqB++l84wZs7A0g27VKjCoti2B1eS6mhBCvAgJuEKWnKlm7500Nt9I5fe7aWTk31cEAH83KwbWsqVvDRvp2i+EEIVEAq4QpGdr2P9PGiE3Utl9O43krGedgMyZ4LhfDRv6Vrehmr38GYQQorDJJ+tzSs3S8Mc/aWy/mcqeO2kkZD471BooLehXM6cHZE1HOfRCCFGU5FO2AJIy1fx+N43Qm2n8fte4llpNB3NtqDVwlnusCSFEcZGAewZVRs41tdCbqez/J420Z8wqAlDV1px+NW14tYYNTVwspQekEP+vvXuPirrMHzj+RvCCF5wFYTCU8TKjMIlZgp0FNRQDwQt5ycuyraus68HMQtQwWDHNRrBUOhXlbdmE1g7IUczM1gVSE93Tbhp7WovFBWMzCAQEQkmY3x/8ZnIYFLywM4yf1zlziuf7nXk+85zH+Xy/z/f5fh8hLEASXBuqrjdz5FIDh4tbZj92ZKKI0rEbM1WOzBnmyDi3HnSTpCaEEBYlCe7/3WjWk1b4I9nFDZy4fJ0OjD7i0dueGUN6MVPVcq+aPLFfCCGshyS4/2dvB8kFte2uhO3Z157wIY7MVDky1rW7nKkJIYSVsuhNV5999hkLFizA29sbhUJBenq6yXa9Xo9Op8PLywt3d3emTZvGv/71r06Jxc7Ojpkqxza3DXeyZ9XovuTNcOX8XCWb/PrjJ8OQQghh1Sya4Orr69FqtWzZsgVHR/PkkpyczFtvvUViYiI5OTm4uroya9YsamtrOyWe8CE/x+CtcGDtmH58Fu7G57OVrB/bnzED5JFZQgjRVVh0iDI4OJjg4GAAli9fbrJNr9eTkpLCCy+8QHh4OAApKSloNBoyMzNZvHjxfY/n0QHdecXPiScH9WKkQqb0CyFEV2a1z4UqKSmhrKyMyZMnG8scHR3x9/fn7NmznVKnnZ0dK0b1k+QmhBA2wGoTXFlZGQCurq4m5a6urpSXl1siJCGEEF2I1c+ibH3NS6/X3/Y6WGFhocl/xe1JO3WMtFPHSVt1jLRT57PaBKdUKgEoLy9n0KBBxvKKigqzs7qbaTQaCgsL0Wg0nR5jVyft1DHSTh0nbdUx0k4dc68HAVY7RKlSqVAqleTm5hrLrl27Rn5+Po8//rgFIxNCCNEVWPQMrq6ujosXLwLQ3NxMaWkpX375Jb/4xS8YPHgwUVFRvP7662g0GtRqNa+99hp9+vRh7ty5lgxbCCFEF2DRBPfFF18wY8YM4986nQ6dTsfChQtJSUnh+eefp6GhgTVr1lBdXc3YsWPJysqiX79+FoxaCCFEV2BXXV3dgacuCiGEEF2L1V6DE0IIIe6FJDghhBA2SRKcEEIImyQJTgghhE2SBCeEEMIm2VSC2717N6NHj0apVPLEE09w+vRpS4dkdXQ6HQqFwuQ1YsQIS4dlcda0NqE1a6+doqKizPrXlClTLBSt5Wzbto1JkyYxePBghg8fzvz58/nqq69M9pE+1bF2upc+ZTMJLisri9jYWGJiYjhx4gTjxo3j6aef5ttvv7V0aFZHo9Hw9ddfG19yIGB9axNaq/baCSAwMNCkf2VkZPyPo7S8U6dOERkZybFjx8jOzsbBwYGnnnqKqqoq4z7SpzrWTnD3fcpm7oMLCgri4Ycf5o033jCWPfbYY4SHh5OQkGDByKyLTqcjOzub/Px8S4ditTw8PEhKSiIiIgJoOdL28vJi6dKlrF69GoCGhgY0Gg2bNm3qlLUJu4LW7QQtR9tXrlzhgw8+sGBk1qeurg5PT0/S09MJDQ2VPnULrdsJ7q1P2cQZXGNjI+fOnTNZOw5g8uTJnbZ2XFdWXFyMt7c3o0ePZsmSJRQXF1s6JKtmibUJu7L8/HzUajVjx45l5cqV/PDDD5YOyeLq6upobm5GoVAA0qdupXU7Gdxtn7La1QTuRGVlJU1NTbJ2XAf4+vry9ttvo9FoqKioYOvWrQQHB3PmzBmcnZ0tHZ5Vut3ahJcvX7ZESFZrypQpzJgxA5VKxaVLl3jllVeYOXMmeXl59OzZ09LhWUxsbCw+Pj6MGzcOkD51K63bCe6tT9lEgjO407XjHkRPPvmkyd++vr6MGTOG999/nxUrVlgoqq5B+lf75syZY/z/hx9+mDFjxuDj48OxY8eYOXOmBSOznJdeeokzZ87w8ccfY29vb7JN+tTPbtVO99KnbGKI0sXFBXt7e7OztfbWjhPQt29fvLy8jKs6CHM3r014M+lf7Rs4cCAPPfTQA9u/1q1bx4EDB8jOzmbIkCHGculTpm7VTm25kz5lEwmuR48ejBkzxmTtOIDc3FxZO64d165do7Cw0PgPTpiTtQnvXmVlJZcvX34g+9eLL75IZmYm2dnZZrfiSJ/62e3aqS130qfsY2NjN9yHGC2uX79+6HQ63N3d6dWrF1u3buX06dO8+eab9O/f39LhWY34+Hh69OhBc3Mz//73v1mzZg0XL15k+/btD3Q71dXVceHCBcrKyti3bx9arRYnJycaGxvp378/TU1NbN++HbVaTVNTE3FxcZSVlbFjx44H6trS7drJ3t6ejRs30rdvX27cuEFBQQHPPfccTU1NbN269YFqp9WrV7N//35SU1MZNGgQ9fX11NfXAy0H5HZ2dtKnaL+d6urq7qlP2cxtAtByo3dycjJlZWV4e3vz6quvEhAQYOmwrMqSJUs4ffo0lZWVDBgwAF9fX+Li4vDy8rJ0aBZ18uRJk7UJDQxrE+r1erZs2UJqaqpxbcLXXnsNrVZrgWgt53bttG3bNiIiIvjyyy+pqalBqVQyYcIE4uLiGDRokAWitZzWswANXnzxRdatWwcgfYr226mhoeGe+pRNJTghhBDCwCauwQkhhBCtSYITQghhkyTBCSGEsEmS4IQQQtgkSXBCCCFskiQ4IYQQNkkSnBAWkpubyxNPPIG7uzsKhYKSkpJOqysqKgofH59O+/z7IT09vdPbQTxYJMEJi9Pr9YSFhTF06NA2l8Gor6/Hx8eHX/7yl/z0008WiPD+q6mp4be//S16vZ6kpCTeffddBgwY0Oa+JSUlZisa3/yKj4//H0d/9+rq6tDpdJw8edLSoYgHgE2tJiC6Jjs7O5KTkxk/fjxxcXHs3LnTZPuWLVsoLS3l448/pnv37haK8v4yPJlh7dq1TJ8+vUPvmT17NiEhIWblI0eOvN/hdZr6+noSExMBmDBhgsm2BQsWMGfOnAfmMVWi80mCE1ZBo9GwatUqdDodCxYsMC4EWVBQQEpKCpGRkZ3+EFq9Xs/169fp1atXp9YDLU+NB+7o+Z8+Pj7Mnz+/s0KyOHt7e7PlZIS4FzJEKaxGdHQ0Xl5erFq1ioaGBpqbm4mOjsbNzY3169dTVFTEkiVLGD58OG5ubvj7+5OWlmbyGY2NjWzevJnAwEBUKhXu7u4EBQXx0UcfmdWnUCiIjo7m4MGD+Pv74+bmxoEDBwD49NNPCQ0NRaVS4eHhga+vLzExMR36Hu+99x7+/v4olUrUajXLli0zWcRy2rRpLF68GIAZM2agUCiYNm3a3TabmbS0NMaOHYtSqSQgIICjR4+a7XPy5EkUCkWbQ4UKhQKdTmdS9v333/PCCy+g1Wpxc3PDx8eHlStXUltbC0BVVRXx8fH4+/szaNAgPDw8mD59OmfOnDF+RklJifFsMzEx0TjEGhUVBdz6GtxHH31EUFAQAwcORKVSERERwTfffGOyj06nQ6FQUFRURHR0NEOHDsXDw4NFixZx5cqVu2hFYQvkDE5YjR49epCcnMzUqVNJSkrCw8ODzz//nLS0NC5fvkxISAguLi48++yz9O/fn08++YQVK1Zw9epVli9fDkBtbS1//OMfmTVrFr/+9a9paGggIyODiIgIMjMzCQoKMqkzPz+fQ4cOsXTpUpRKJSNGjODChQvMmzcPrVZLbGwsvXv3pri4mGPHjrX7HbZv387LL7+Mv78/GzdupLS0lF27dpGfn8+JEydQKBSsXr0arVbLrl27iImJYcSIEbi5ubX72T/++COVlZVm5U5OTsahW8PCtY899hi/+93v+OGHH1i2bNk9Pey4rKyMoKAgKioq+M1vfoNWq+X777/nww8/5MqVK/Tr14/i4mIOHTpEeHg4w4YNo6amhvfee4/w8HByc3PRarUMGDCArVu3smbNGqZPn258aPPQoUNvWXdmZiZLly5l1KhRxMXFcfXqVXbu3ElwcDB5eXlma4dFRkaiVCqJi4ujqKiInTt30r17d3bv3n3X3190XfKwZWF1YmJi+NOf/oSjoyOBgYHs27ePWbNm8d1335Gbm0vv3r2N+y5evJjjx49z4cIF+vTpQ1NTEzdu3DC5jtPY2MiECRMYOHAgBw8eNJYrFArs7Oz49NNPGT16tLE8JSWFdevWUVRUhIuLS4fjrqysRKvV4ufnx8GDB3FwaDl+PHLkCBEREaxevdo4IeTAgQNERkZy+PBhs2tRrZWUlPDII4/ccvv+/fuZOnUqN27cQKvV4uzsTF5ennGoNScnh9mzZzN48GAKCgqAn1cFaKt+hUJh8tT7qKgo9u/fz5EjR/D39zfZ17AC9fXr13FwcDAZYqyqqsLPz4+wsDDeeOMNoCVZjhw50uTzDdLT03n22Wc5f/48KpWKn376iVGjRuHk5EReXh59+vQB4Pz580yaNIm5c+car9fqdDoSExOZN2+eyTXc2NhYdu3axcWLFx/o5aAeVDJEKaxOQkICLi4uxhmG1dXV5OXl8dRTT9HQ0EBlZaXxNWXKFGpra/niiy+Alus4huTW2NhIVVUVtbW1BAQEcO7cObO6Hn/8cZPkBi1rC0JLYmpubu5w3Hl5eVy/fp3ly5cbkxu0DElqNJoOnQHezjPPPMPBgwfNXn5+fgD84x//oLy8nMWLF5tcR5w8efJdL4fU3NzMkSNHmDJlillyg5YJQgA9e/Y0Jrdr165x5coVmpubGTt2bJvt3hHnzp2jrKyMyMhIY3IDeOSRRwgMDOSTTz5Brzc9Po+MjDT5OyAggKamJkpLS+8qBtG1yRClsDpOTk6o1WrKy8sZOHAgf//739Hr9SQmJhpn4LVmmLQBLdfA3n77bb7++muTH0DDj/HNWg9xAcyZM4d9+/axcuVKNmzYwMSJEwkLC2PWrFm3ncV56dIlgDZXJR4xYgSnTp265Xs7YtiwYQQGBt5y+7fffgu0TNhpTa1Wc/78+Tuus6KigqtXr7a7RllzczPJycmkpqaaXUNTqVR3XC/cvj1HjhxJTk4OV69eNTkzGzx4sMl+hvXGqqqq7ioG0bVJghNWz3AWtXz5coKDg9vcx/ADnJmZycqVKwkNDeX555/H1dUVBwcH0tPTycjIMHufo6Njm2VHjx7l1KlTHD9+nL/+9a/8/ve/58033+TYsWNtvqc9rc80OoOhjrYSeev629oHoKmpqcOfebMdO3awceNGFi5cSHx8PM7Oztjb27Nt2zb+85//dPg7dNSt2vNWszD/F+0vrI8kOGH1DGdZDg4Otz2DAcjKymLIkCG8//77Jj/K6enpd1Rnt27dmDhxIhMnTmTjxo3s2bOHmJgYDh8+zLx589p8j6enJwDffPMNarXaZFthYaFxe2e5uf5JkyaZbCsqKjL523BmU1NTY1JuOGsycHV1xcnJia+++uq2dWdlZTF+/HhSUlJMylvPxmwvUd7s5u9juG3EoLCwEIVCgZOTU4c/Tzx45BqcsHqurq5MnDiR1NTUNq+l3Dw8aTiCv/mIvbi4mA8//LDD9bU1rdwwyaO6uvqW7wsMDKRnz5688847JmdCR48epbCwsM2btO+nRx99FFdXV1JTU7l27ZqxPCcnhwsXLpjs6+npib29vdltAq1vsu/WrRvTpk3jL3/5C2fPnjWr09DO9vb2ZmdJZ8+e5W9/+5tJmWGC0O3a0WDMmDEolUr27t1LQ0ODsbygoIDc3FyCg4PvKGGKB4+cwYkuYdu2bYSEhBAQEMCiRYsYPnw4lZWVnD9/npycHOP1p9DQUA4fPszChQsJDQ3lu+++Y8+ePQwfPpx//vOfHaorKSmJU6dOERISgqenJ9XV1ezdu5c+ffowderUW77PxcWF2NhYXn75ZcLDw5kxYwb//e9/2blzJ56enqxYseKe2qCgoIAPPvjArPyhhx5iwoQJdO/enfXr1/Pcc88RFhbG008/TUVFBbt27cLb25u6ujrje5ycnJgzZw67d+/Gzs4OjUbDyZMnKS4uNvv8hIQE8vLyCA8PZ9GiRXh7e1NeXs7hw4dJS0tDpVIRGhrKli1bWLZsGf7+/hQVFZGamoqXl5dJvX379kWj0ZCVlYVarcbZ2RmVSoWvr69Zvd27d2fz5s0sXbqUkJAQ5s+fb7xNwMnJiZdeeume2lPYPklwoktQq9Xk5eWRlJRERkYGFRUVuLi4MHLkSDZt2mTc71e/+hUVFRXs2bOHvLw8hg0bxquvvsrFixc7nODCwsIoLS3lz3/+MxUVFTg7O+Pn58fatWvbHWaMjo7GxcWFd955hz/84Q/07duX8PBwEhISjMOCdysrK4usrCyz8pCQEONU/2eeeQa9Xs+OHTtISEhArVbz7rvvkp2dbTbJJTExkRs3bpCWlka3bt0IDg4mMzPTbHjV3d2d48ePs3nzZrKysqipqcHd3Z3Jkycbb6Mw3JyfkZHBoUOH8Pb2Zu/evRw4cMCs3rfeeot169YRHx/P9evXWbhwYZsJDmDu3Lk4Ojry+uuvs2nTJnr06MH48ePZsGFDmxOEhLiZ3AcnhBDCJsk1OCGEEDZJEpwQQgibJAlOCCGETZIEJ4QQwiZJghNCCGGTJMEJIYSwSZLghBBC2CRJcEIIIWySJDghhBA2SRKcEEIIm/R/QluLg/4PXIsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import style\n",
    "style.use(\"fivethirtyeight\")\n",
    "\n",
    "x = np.array(range(1, 25))\n",
    "plt.plot(x, np.exp(2.2954 + 0.0529 * x))\n",
    "plt.xlabel(\"Years of Education\")\n",
    "plt.ylabel(\"Hourly Wage\")\n",
    "plt.title(\"Impact of Education on Hourly Wage\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, it is not because we can estimate this simple model that it's correct. Notice how I was carefully with my words saying it **predicts** wage from education. I never said that this prediction was causal. In fact, by now, you probably have very serious reasons to believe this model is biased. Since our data didn't come from a random experiment, we don't know if those that got more education are comparable to those who got less. Going even further, from our understanding of how the world works, we are very certain that they are not comparable. Namely, we can argue that those with more years of education probably have richer parents, and that the increase we are seeing in wages as we increase education is just a reflection of how the family wealth is associated with more years of education. Putting it in math terms, we think that \\\\(E[Y_0|T=0] < E[Y_0|T=1]\\\\), that is, those with more education would have higher income anyway, even without so many years of education. If you are really grim about education, you can argue that it can even *reduce* wages by keeping people out of the workforce and lowering their experience.\n",
    "\n",
    "Fortunately, in our data, we have access to lots of other variables. We can see the parents' education `meduc`, `feduc`, the `IQ` score for that person, the number of years of experience `exper` and the tenure of the person in his or her current company `tenure`. We even have some dummy variables for marriage and black ethnicity. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>wage</th>\n",
       "      <th>hours</th>\n",
       "      <th>lhwage</th>\n",
       "      <th>IQ</th>\n",
       "      <th>educ</th>\n",
       "      <th>exper</th>\n",
       "      <th>tenure</th>\n",
       "      <th>age</th>\n",
       "      <th>married</th>\n",
       "      <th>black</th>\n",
       "      <th>south</th>\n",
       "      <th>urban</th>\n",
       "      <th>sibs</th>\n",
       "      <th>brthord</th>\n",
       "      <th>meduc</th>\n",
       "      <th>feduc</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>769</td>\n",
       "      <td>40</td>\n",
       "      <td>2.956212</td>\n",
       "      <td>93</td>\n",
       "      <td>12</td>\n",
       "      <td>11</td>\n",
       "      <td>2</td>\n",
       "      <td>31</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>8.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>825</td>\n",
       "      <td>40</td>\n",
       "      <td>3.026504</td>\n",
       "      <td>108</td>\n",
       "      <td>14</td>\n",
       "      <td>11</td>\n",
       "      <td>9</td>\n",
       "      <td>33</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>14.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>650</td>\n",
       "      <td>40</td>\n",
       "      <td>2.788093</td>\n",
       "      <td>96</td>\n",
       "      <td>12</td>\n",
       "      <td>13</td>\n",
       "      <td>7</td>\n",
       "      <td>32</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>3.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>12.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>562</td>\n",
       "      <td>40</td>\n",
       "      <td>2.642622</td>\n",
       "      <td>74</td>\n",
       "      <td>11</td>\n",
       "      <td>14</td>\n",
       "      <td>5</td>\n",
       "      <td>34</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>10</td>\n",
       "      <td>6.0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>11.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>600</td>\n",
       "      <td>40</td>\n",
       "      <td>2.708050</td>\n",
       "      <td>91</td>\n",
       "      <td>10</td>\n",
       "      <td>13</td>\n",
       "      <td>0</td>\n",
       "      <td>30</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>8.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   wage  hours    lhwage   IQ  educ  exper  tenure  age  married  black  \\\n",
       "0   769     40  2.956212   93    12     11       2   31        1      0   \n",
       "2   825     40  3.026504  108    14     11       9   33        1      0   \n",
       "3   650     40  2.788093   96    12     13       7   32        1      0   \n",
       "4   562     40  2.642622   74    11     14       5   34        1      0   \n",
       "6   600     40  2.708050   91    10     13       0   30        0      0   \n",
       "\n",
       "   south  urban  sibs  brthord  meduc  feduc  \n",
       "0      0      1     1      2.0    8.0    8.0  \n",
       "2      0      1     1      2.0   14.0   14.0  \n",
       "3      0      1     4      3.0   12.0   12.0  \n",
       "4      0      1    10      6.0    6.0   11.0  \n",
       "6      0      1     1      2.0    8.0    8.0  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wage.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can include all those extra variables in a model and estimate it:\n",
    "\n",
    "$\n",
    "log(hwage)_i = \\beta_0 + \\kappa \\ educ_i + \\pmb{\\beta}X_i + u_i\n",
    "$\n",
    "\n",
    "To understand how this helps with the bias problem, let's recap the bivariate breakdown of multivariate linear regression.\n",
    "\n",
    "$\n",
    "\\kappa = \\dfrac{Cov(Y_i, \\tilde{T_i})}{Var(\\tilde{T_i})} \n",
    "$\n",
    "\n",
    "This formula says that we can predict `educ` from mother and father education, from IQ, from experience and so on. After we do that, we'll be left with a version of `educ`, \\\\(\\tilde{educ}\\\\), which is uncorrelated with all the variables included previously. This will break down arguments such as \"people that have more years of education have it because they have higher IQ. It is not the case that education leads to higher wages. It is just the case that it is correlated with IQ, which is what drives waged\". Well, if we include IQ in our model, then \\\\(\\kappa\\\\) becomes the return of an additional year of education while keeping IQ fixed. Pause a little bit to understand what this implies. Even if we can't use randomised controlled trials to keep other factors equal between treated and untreated, regression can do this by including those same factors in the model, even if the data is not random!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.04114719101005262"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "controls = ['IQ', 'exper', 'tenure', 'age', 'married', 'black',\n",
    "            'south', 'urban', 'sibs', 'brthord', 'meduc', 'feduc']\n",
    "\n",
    "X = wage[controls].assign(intercep=1)\n",
    "t = wage[\"educ\"]\n",
    "y = wage[\"lhwage\"]\n",
    "\n",
    "beta_aux = regress(t, X)\n",
    "t_tilde = t - X.dot(beta_aux)\n",
    "\n",
    "kappa = t_tilde.cov(y) / t_tilde.var()\n",
    "kappa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This coefficient we've just estimated tells us that, for people with the same IQ, experience, tenure, age and so on, we should expect an additional year of education to be associated with a 4.11% increase in hourly wage. This confirms our suspicion that the first simple model with only `educ` was biases. It also confirms that this bias was overestimating the impact of education. Once we controlled for other factors, the estimated impact of education fell.\n",
    "\n",
    "If we are wiser and use software that other people wrote instead of coding everything yourself, we can even place a confidence interval around this estimate. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>    1.1156</td> <td>    0.232</td> <td>    4.802</td> <td> 0.000</td> <td>    0.659</td> <td>    1.572</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>educ</th>      <td>    0.0411</td> <td>    0.010</td> <td>    4.075</td> <td> 0.000</td> <td>    0.021</td> <td>    0.061</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>IQ</th>        <td>    0.0038</td> <td>    0.001</td> <td>    2.794</td> <td> 0.005</td> <td>    0.001</td> <td>    0.006</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>exper</th>     <td>    0.0153</td> <td>    0.005</td> <td>    3.032</td> <td> 0.003</td> <td>    0.005</td> <td>    0.025</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>tenure</th>    <td>    0.0094</td> <td>    0.003</td> <td>    2.836</td> <td> 0.005</td> <td>    0.003</td> <td>    0.016</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age</th>       <td>    0.0086</td> <td>    0.006</td> <td>    1.364</td> <td> 0.173</td> <td>   -0.004</td> <td>    0.021</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>married</th>   <td>    0.1795</td> <td>    0.053</td> <td>    3.415</td> <td> 0.001</td> <td>    0.076</td> <td>    0.283</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>black</th>     <td>   -0.0801</td> <td>    0.063</td> <td>   -1.263</td> <td> 0.207</td> <td>   -0.205</td> <td>    0.044</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>south</th>     <td>   -0.0397</td> <td>    0.035</td> <td>   -1.129</td> <td> 0.259</td> <td>   -0.109</td> <td>    0.029</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>urban</th>     <td>    0.1926</td> <td>    0.036</td> <td>    5.418</td> <td> 0.000</td> <td>    0.123</td> <td>    0.262</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>sibs</th>      <td>    0.0065</td> <td>    0.009</td> <td>    0.722</td> <td> 0.470</td> <td>   -0.011</td> <td>    0.024</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>brthord</th>   <td>   -0.0080</td> <td>    0.013</td> <td>   -0.604</td> <td> 0.546</td> <td>   -0.034</td> <td>    0.018</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>meduc</th>     <td>    0.0089</td> <td>    0.007</td> <td>    1.265</td> <td> 0.206</td> <td>   -0.005</td> <td>    0.023</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>feduc</th>     <td>    0.0069</td> <td>    0.006</td> <td>    1.113</td> <td> 0.266</td> <td>   -0.005</td> <td>    0.019</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_2 = smf.ols('lhwage ~ educ +'+\"+\".join(controls), data=wage).fit()\n",
    "model_2.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Omitted Variable or Confounding Bias\n",
    "\n",
    "The question that remains is: this parameter we've estimated is causal? Unfortunately, we can't say for sure. We can argue that the first simple model that regress wage on education probably isn't. It omits important variables that are correlated both with education and with wages. Without controlling for them, the estimated impact of education is also capturing the impact of those other variables that were not included in the model.\n",
    "\n",
    "To better understand how this bias work, let's suppose the true model for how education affect wage looks a bit like this\n",
    "\n",
    "$\n",
    "Wage_i = \\alpha + \\kappa \\ Educ_i + A_i'\\beta + u_i\n",
    "$\n",
    "\n",
    "wage is affected by education, which is measured by the size of \\\\(\\kappa\\\\) and by additional ability factors, denoted as the vector \\\\(A\\\\). If we omit ability from our model, our estimate for \\\\(\\kappa\\\\) will look like this:\n",
    "\n",
    "$\n",
    "\\dfrac{Cov(Wage_i, Educ_i)}{Var(Educ_i)} = \\kappa + \\beta'\\delta_{Ability}\n",
    "$\n",
    "\n",
    "where \\\\(\\delta_{A}\\\\) is the vector of coefficients from the regression of \\\\(A\\\\) on \\\\(Educ\\\\)\n",
    "\n",
    "The key point here is that it won't be exactly the \\\\(\\kappa\\\\) that we want. Instead, it comes with this extra annoying term \\\\(\\beta'\\delta_{A}\\\\). This term is the impact of the omitted \\\\(A\\\\) on \\\\(Wage\\\\), \\\\(\\beta\\\\) times the impact of the omitted on the included \\\\(Educ\\\\). This is important for economists that Joshua Angrist made a mantra out of it so that the students can recite it in meditation:\n",
    "\n",
    "```\n",
    "\"Short equals long \n",
    "plus the effect of omitted \n",
    "times the regression of omitted on included\"\n",
    "```\n",
    "\n",
    "Here, the short regression is the one that omits variables, while the long is the one that includes them. This formula or mantra gives us further insight into the nature of bias. First, the bias term will be zero if the omitted variables have no impact on the dependent variable \\\\(Y\\\\). This makes total sense. I don't need to control for stuff that is irrelevant for wages when trying to understand the impact of education on it (like how tall the lilies of the field). Second, the bias term will also be zero if the omitted variables have no impact on the treatment variable. This also makes intuitive sense. If everything that impacts education has been included in the model, there is no way the estimated impact of education is mixed with a correlation from education on something else that also impacts wages. \n",
    "\n",
    "![img](data/img/linear-regression/confused_cat.png)\n",
    "\n",
    "To put it more succinctly, we say that **there is no OVB if all the confounding variables are accounted for in the model**. We can also leverage our knowledge about causal graphs here. A confounding variable is one that **causes both the treatment and the outcome**. In the wage example, IQ is a confounder. People with high IQ tend to complete more years of education because it's easier for them, so we can say that IQ causes education. People with high IQ also tend to be naturally more productive and consequently have higher wages, so IQ also causes wage. Since confounders are variables that affect both the treatment and the outcome, we mark them with an arrow going to T and Y. Here, I've denoted them with \\\\(W\\\\). I've also marked positive causation with red and negative causation with blue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"342pt\" height=\"188pt\"\n",
       " viewBox=\"0.00 0.00 342.50 188.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 184)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-184 338.4994,-184 338.4994,4 -4,4\"/>\n",
       "<!-- W -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>W</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"69\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"69\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">W</text>\n",
       "</g>\n",
       "<!-- T -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>T</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">T</text>\n",
       "</g>\n",
       "<!-- W&#45;&gt;T -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>W&#45;&gt;T</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M59.0466,-144.937C53.9729,-136.2393 47.6966,-125.4799 42.047,-115.7948\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"44.944,-113.8148 36.882,-106.9405 38.8975,-117.3419 44.944,-113.8148\"/>\n",
       "</g>\n",
       "<!-- Y -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>Y</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"63\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"63\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Y</text>\n",
       "</g>\n",
       "<!-- W&#45;&gt;Y -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>W&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M68.2401,-143.7623C67.2167,-119.201 65.3853,-75.2474 64.1814,-46.3541\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"67.6671,-45.9352 63.7537,-36.0896 60.6732,-46.2267 67.6671,-45.9352\"/>\n",
       "</g>\n",
       "<!-- T&#45;&gt;Y -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>T&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M35.7146,-72.5708C39.9597,-64.0807 45.1536,-53.6929 49.8663,-44.2674\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"53.024,-45.7782 54.3657,-35.2687 46.763,-42.6477 53.024,-45.7782\"/>\n",
       "</g>\n",
       "<!-- IQ -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>IQ</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"160\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"160\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">IQ</text>\n",
       "</g>\n",
       "<!-- Educ -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>Educ</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"147\" cy=\"-90\" rx=\"28.0565\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"147\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Educ</text>\n",
       "</g>\n",
       "<!-- IQ&#45;&gt;Educ -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>IQ&#45;&gt;Educ</title>\n",
       "<path fill=\"none\" stroke=\"#ff0000\" d=\"M156.7196,-143.8314C155.3078,-136.0125 153.625,-126.6923 152.0595,-118.0221\"/>\n",
       "<polygon fill=\"#ff0000\" stroke=\"#ff0000\" points=\"155.4556,-117.1324 150.2344,-107.9134 148.567,-118.3762 155.4556,-117.1324\"/>\n",
       "</g>\n",
       "<!-- Wage -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>Wage</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"160\" cy=\"-18\" rx=\"30.0438\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"160\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Wage</text>\n",
       "</g>\n",
       "<!-- IQ&#45;&gt;Wage -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>IQ&#45;&gt;Wage</title>\n",
       "<path fill=\"none\" stroke=\"#ff0000\" d=\"M169.7067,-145.1142C175.0856,-134.7678 181.2314,-121.0465 184,-108 187.3215,-92.3486 187.3215,-87.6514 184,-72 182.0194,-62.6668 178.3104,-52.9882 174.3903,-44.4468\"/>\n",
       "<polygon fill=\"#ff0000\" stroke=\"#ff0000\" points=\"177.4888,-42.8155 169.9555,-35.3657 171.1987,-45.8873 177.4888,-42.8155\"/>\n",
       "</g>\n",
       "<!-- Educ&#45;&gt;Wage -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>Educ&#45;&gt;Wage</title>\n",
       "<path fill=\"none\" stroke=\"#ff0000\" d=\"M150.2804,-71.8314C151.6922,-64.0125 153.375,-54.6923 154.9405,-46.0221\"/>\n",
       "<polygon fill=\"#ff0000\" stroke=\"#ff0000\" points=\"158.433,-46.3762 156.7656,-35.9134 151.5444,-45.1324 158.433,-46.3762\"/>\n",
       "</g>\n",
       "<!-- Crime -->\n",
       "<g id=\"node7\" class=\"node\">\n",
       "<title>Crime</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"293\" cy=\"-162\" rx=\"32.4237\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"293\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Crime</text>\n",
       "</g>\n",
       "<!-- Police -->\n",
       "<g id=\"node8\" class=\"node\">\n",
       "<title>Police</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"263\" cy=\"-90\" rx=\"32.4237\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"263\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Police</text>\n",
       "</g>\n",
       "<!-- Crime&#45;&gt;Police -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>Crime&#45;&gt;Police</title>\n",
       "<path fill=\"none\" stroke=\"#ff0000\" d=\"M285.5843,-144.2022C282.1693,-136.0064 278.0427,-126.1024 274.2561,-117.0145\"/>\n",
       "<polygon fill=\"#ff0000\" stroke=\"#ff0000\" points=\"277.4074,-115.4776 270.3304,-107.593 270.9459,-118.17 277.4074,-115.4776\"/>\n",
       "</g>\n",
       "<!-- Violence -->\n",
       "<g id=\"node9\" class=\"node\">\n",
       "<title>Violence</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"293\" cy=\"-18\" rx=\"41.4988\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"293\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Violence</text>\n",
       "</g>\n",
       "<!-- Crime&#45;&gt;Violence -->\n",
       "<g id=\"edge8\" class=\"edge\">\n",
       "<title>Crime&#45;&gt;Violence</title>\n",
       "<path fill=\"none\" stroke=\"#ff0000\" d=\"M297.5787,-144.1385C300.0216,-133.6847 302.7716,-120.1847 304,-108 305.6049,-92.0807 305.6049,-87.9193 304,-72 303.1498,-63.5665 301.5706,-54.5029 299.8689,-46.2554\"/>\n",
       "<polygon fill=\"#ff0000\" stroke=\"#ff0000\" points=\"303.2578,-45.3697 297.6934,-36.3537 296.4209,-46.8719 303.2578,-45.3697\"/>\n",
       "</g>\n",
       "<!-- Police&#45;&gt;Violence -->\n",
       "<g id=\"edge9\" class=\"edge\">\n",
       "<title>Police&#45;&gt;Violence</title>\n",
       "<path fill=\"none\" stroke=\"#0000ff\" d=\"M270.4157,-72.2022C273.8052,-64.0675 277.8958,-54.2501 281.6591,-45.2181\"/>\n",
       "<polygon fill=\"#0000ff\" stroke=\"#0000ff\" points=\"284.9502,-46.4195 285.5656,-35.8425 278.4886,-43.7271 284.9502,-46.4195\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a19407828>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = gr.Digraph()\n",
    "\n",
    "g.edge(\"W\", \"T\"), g.edge(\"W\", \"Y\"), g.edge(\"T\", \"Y\")\n",
    "\n",
    "g.edge(\"IQ\", \"Educ\", color=\"red\"), g.edge(\"IQ\", \"Wage\", color=\"red\"), g.edge(\"Educ\", \"Wage\", color=\"red\")\n",
    "\n",
    "g.edge(\"Crime\", \"Police\", color=\"red\"), g.edge(\"Crime\", \"Violence\", color=\"red\"), \n",
    "g.edge(\"Police\", \"Violence\", color=\"blue\")\n",
    "\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Causal graphs are excellent to depict our understanding of the world and understand how confounding bias works. In our first example, we have a graph where education causes wage: more education leads to higher wages. However, IQ also causes wage and it also causes education: high IQ causes both more education and wage. If we don't account for IQ in our model, some of its effect on wage will flow through the correlation with education. That will make the impact of education look higher than it actually is. This is an example of positive bias.\n",
    "\n",
    "Just to give another example, but with negative bias, consider the causal graph about the effect of police on city violence. What we usually see in the world is that cities with higher police force also have more violence. Does this mean that the police are causing the violence? Well, it could be, I don't think it's worth getting into that discussion here. But, there is also a strong possibility that there is a confounding variable causing us to see a biased version of the impact of police on violence. It could be that increasing police force decreases violence. But, a third variable crime, causes both more violence and more police force. If we don't account for it, the impact of crime on violence will flow through police force, making it look like it increases violence. This is an example of negative bias.\n",
    "\n",
    "Causal graphs can also show us how both regression and randomized control trials are correct for confounding bias. RCT does so by severing the connection of the confounder to the treatment variable. By making \\\\(T\\\\) random, by definition, nothing can cause it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"280pt\" height=\"116pt\"\n",
       " viewBox=\"0.00 0.00 280.28 116.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 112)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-112 276.2774,-112 276.2774,4 -4,4\"/>\n",
       "<!-- W -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>W</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">W</text>\n",
       "</g>\n",
       "<!-- Y -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>Y</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"61\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"61\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Y</text>\n",
       "</g>\n",
       "<!-- W&#45;&gt;Y -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>W&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M35.2305,-72.5708C39.1979,-64.1691 44.0429,-53.9091 48.4567,-44.5623\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"51.7401,-45.8057 52.8453,-35.2687 45.4104,-42.8166 51.7401,-45.8057\"/>\n",
       "</g>\n",
       "<!-- T -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>T</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"99\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">T</text>\n",
       "</g>\n",
       "<!-- T&#45;&gt;Y -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>T&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M89.9945,-72.937C85.5172,-64.4537 80.0048,-54.0092 74.993,-44.513\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"77.9678,-42.6509 70.2048,-35.4407 71.7771,-45.9182 77.9678,-42.6509\"/>\n",
       "</g>\n",
       "<!-- IQ -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>IQ</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"171\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"171\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">IQ</text>\n",
       "</g>\n",
       "<!-- Wage -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>Wage</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"207\" cy=\"-18\" rx=\"30.0438\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"207\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Wage</text>\n",
       "</g>\n",
       "<!-- IQ&#45;&gt;Wage -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>IQ&#45;&gt;Wage</title>\n",
       "<path fill=\"none\" stroke=\"#ff0000\" d=\"M179.7146,-72.5708C183.9597,-64.0807 189.1536,-53.6929 193.8663,-44.2674\"/>\n",
       "<polygon fill=\"#ff0000\" stroke=\"#ff0000\" points=\"197.024,-45.7782 198.3657,-35.2687 190.763,-42.6477 197.024,-45.7782\"/>\n",
       "</g>\n",
       "<!-- Educ -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>Educ</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"244\" cy=\"-90\" rx=\"28.0565\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"244\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Educ</text>\n",
       "</g>\n",
       "<!-- Educ&#45;&gt;Wage -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>Educ&#45;&gt;Wage</title>\n",
       "<path fill=\"none\" stroke=\"#ff0000\" d=\"M235.0433,-72.5708C230.6803,-64.0807 225.3422,-53.6929 220.4985,-44.2674\"/>\n",
       "<polygon fill=\"#ff0000\" stroke=\"#ff0000\" points=\"223.5579,-42.5633 215.8742,-35.2687 217.3319,-45.7628 223.5579,-42.5633\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1047f2e48>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = gr.Digraph()\n",
    "\n",
    "g.edge(\"W\", \"Y\"), g.edge(\"T\", \"Y\")\n",
    "\n",
    "g.edge(\"IQ\", \"Wage\", color=\"red\"), g.edge(\"Educ\", \"Wage\", color=\"red\")\n",
    "\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regression, on the other hand, does so by comparing the effect of \\\\(T\\\\) while maintaining the confounder \\\\(W\\\\) set to a fixed level. With regression, it is not the case that W cease to cause T and Y. It is just that it is held fixed, so it can't influence changes on T and Y."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"292pt\" height=\"116pt\"\n",
       " viewBox=\"0.00 0.00 291.85 116.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 112)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-112 287.8495,-112 287.8495,4 -4,4\"/>\n",
       "<!-- W=w -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>W=w</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"29.8276\" cy=\"-90\" rx=\"29.6561\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"29.8276\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">W=w</text>\n",
       "</g>\n",
       "<!-- T -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>T</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"104.8276\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"104.8276\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">T</text>\n",
       "</g>\n",
       "<!-- Y -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>Y</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"104.8276\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"104.8276\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Y</text>\n",
       "</g>\n",
       "<!-- T&#45;&gt;Y -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>T&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M104.8276,-71.8314C104.8276,-64.131 104.8276,-54.9743 104.8276,-46.4166\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"108.3277,-46.4132 104.8276,-36.4133 101.3277,-46.4133 108.3277,-46.4132\"/>\n",
       "</g>\n",
       "<!-- IQ=x -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>IQ=x</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"178.8276\" cy=\"-90\" rx=\"28.6955\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"178.8276\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">IQ=x</text>\n",
       "</g>\n",
       "<!-- Educ -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>Educ</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"253.8276\" cy=\"-90\" rx=\"28.0565\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"253.8276\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Educ</text>\n",
       "</g>\n",
       "<!-- Wage -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>Wage</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"253.8276\" cy=\"-18\" rx=\"30.0438\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"253.8276\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Wage</text>\n",
       "</g>\n",
       "<!-- Educ&#45;&gt;Wage -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>Educ&#45;&gt;Wage</title>\n",
       "<path fill=\"none\" stroke=\"#ff0000\" d=\"M253.8276,-71.8314C253.8276,-64.131 253.8276,-54.9743 253.8276,-46.4166\"/>\n",
       "<polygon fill=\"#ff0000\" stroke=\"#ff0000\" points=\"257.3277,-46.4132 253.8276,-36.4133 250.3277,-46.4133 257.3277,-46.4132\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a194070f0>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = gr.Digraph()\n",
    "\n",
    "g.node(\"W=w\"), g.edge(\"T\", \"Y\")\n",
    "g.node(\"IQ=x\"), g.edge(\"Educ\", \"Wage\", color=\"red\")\n",
    "\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, back to our question, is the parameter we've estimated for the impact of `educ` on wage casual? I'm sorry to bring it to you, but that will depend on our ability to argue in favor or against that fact that all confounders have been included in the model. Personally, I think they haven't. For instance, we haven't included family wealth. Even if we included family education, that can only be seen as a proxy for wealth. We've also not accounted for factors like personal ambition. It could be that ambition is what causes both more years of education and higher wage, so it is a confounder. This is to show that **causal inference with non-random or observational data should always be taken with a grain of salt**. We can never be sure that all confounders were accounted for.\n",
    "\n",
    "## Key Ideas\n",
    "\n",
    "We've covered a lot of ground with regression. We saw how regression can be used to perform A/B testing and how it conveniently gives us confidence intervals. Then, we moved to study how regression solves a prediction problem and it is the best linear approximation to the CEF. We've also discussed how, in the bivariate case, the regression treatment coefficient is the covariance between the treatment and the outcome divided by the variance of the treatment. Expanding to the multivariate case, we figured out how regression gives us a partialling out interpretation of the treatment coefficient: it can be interpreted as how the outcome would change with the treatment while keeping all other included variables constant. This is what economists love to refer as *ceteris paribus*.\n",
    "\n",
    "Finally, we took a turn to understanding bias. We saw how `Short equals long plus the effect of omitted times the regression of omitted on included`. This shed some light to how bias comes to be. We discovered that the source of omitted variable bias is confounding: a variable that affects both the treatment and the outcome. Lastly, we used causal graphs to see how RCT and regression fixes confounding.\n",
    "\n",
    "## References\n",
    "\n",
    "I like to think of this entire series as a tribute to Joshua Angrist, Alberto Abadie and Christopher Walters for their amazing Econometrics class. Most of the ideas here are taken from their classes at the American Economic Association. Watching them is what is keeping me sane during this tough year of 2020.\n",
    "* [Cross-Section Econometrics](https://www.aeaweb.org/conference/cont-ed/2017-webcasts)\n",
    "* [Mastering Mostly Harmless Econometrics](https://www.aeaweb.org/conference/cont-ed/2020-webcasts)\n",
    "\n",
    "I'll also like to reference the amazing books from Angrist. They have shown me that Econometrics, or 'Metrics as they call it, is not only extremely useful but also profoundly fun.\n",
    "\n",
    "* [Mostly Harmless Econometrics](https://www.mostlyharmlesseconometrics.com/)\n",
    "* [Mastering 'Metrics](https://www.masteringmetrics.com/)\n",
    "\n",
    "My final reference is Miguel Hernan and Jamie Robins' book. It has been my trustworthy companion in the most thorny causal questions I had to answer.\n",
    "\n",
    "* [Causal Inference Book](https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/)\n",
    "\n",
    "\n",
    "![img](./data/img/poetry.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}