1-Introduction-To-Causality.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why Bother?\n",
    "\n",
    "First and foremost, you might be wondering: what's in it for me? Here is what:\n",
    "\n",
    "## Data Science is Not What it Used to (or it Finally Is).\n",
    "\n",
    "Data Scientist has been labeled [The Sexiest Job of the 21st Century](https://hbr.org/2012/10/data-scientist-the-sexiest-job-of-the-21st-century) by Harvard Business Review. This was no empty statement. For a decade now, Data Scientist has been at the spotlight. AI experts had [salaries that rivaled those of sports superstars](https://www.economist.com/business/2016/04/02/million-dollar-babies). In the search for fame and fortune, hundreds of young professionals entered into what seemed a frenetic golden rush to get the Data Science title as quickly as possible. Whole new industries sprang around the hype. Miraculous teaching methods could make you a Data Scientist without requiring you to look at a single math formula. Consulting specialists promised millions if your company could unlock the potential of data. AI or, Machine Learning, has been called the new electricity and data, the new oil. \n",
    "\n",
    "Meanwhile, we kind of forgot about those that have been doing \"old fashioned\" science with data all along. During all this time, economists were trying to answer what is the true impact of education on one's earnings, biostatisticians were trying to understand if saturated fat led to higher chance of heart attack and psychologists were trying to understand if words of affirmation led indeed to a happier marriage. If we were to be completely honest, data scientist is not a recent field. We are simply aware of it just now due to the massive amount of free marketing the media has provided.\n",
    "\n",
    "To use a Jim Collins analogy, think about pouring yourself an ice cold cup of your favorite beer. If you do this the right way, most of the cup will be beer but there will be a 1 finger thick layer of foam at the top. This cup is just like Data Science. \n",
    "\n",
    "1. It's the bear. The statistical foundations, the scientific curiosity, the passion for difficult problems. All of this was proven very valuable throughout hundreds of years.\n",
    "2. It's the foam. The fluffy stuff built on unrealistic expectations that will eventually go away. \n",
    "\n",
    "This foam might come down crashing faster than you think. As The Economist puts it\n",
    "\n",
    "> The same consultants who predict that AI will have a world-altering impact also report that real managers in real companies are finding AI hard to implement, and that enthusiasm for it is cooling. Svetlana Sicular of Gartner, a research firm, says that 2020 could be the year AI falls onto the downslope of her firm’s well-publicised “hype cycle”. Investors are beginning to wake up to bandwagon-jumping: a survey of European AI startups by MMC, a venture-capital fund, found that 40% did not seem to be using any AI at all.\n",
    "\n",
    "In the midst of all this craze, what should we, as Data Scientists - or better yet, as \"just\" Scientists - do? As a starter, if you are smart, you will learn to ignore the foam. We are in it for the beer. Math and statistics has been useful since forever and it is unlikely it will stop now. Second, learn what makes your work valuable and useful, not the latest shining tool that no one figured out how to use. \n",
    "\n",
    "Last but not least, remember that there are no shortcuts. Knowledge in Math and Statistics are valuable precisely because they are hard to acquire. If everyone could do it, excess supply would drive its price down. So **toughen up**! Learn them as well as you can. And heck, why not? have fun along the way as we embark on this quest only **for the Brave and True**.\n",
    "\n",
    "![img](./data/img/intro/tougher-up-cupcake1.jpg)\n",
    "\n",
    "## Answering a Different Kind of Question\n",
    "\n",
    "The type of question Machine Learning is currently very good at answering is of the prediction kind. As Ajay Agrawal, Joshua Gans and Avi Goldfarb puts it in the book Prediction Machines, \"the new wave of artificial intelligence does not actually bring us intelligence but instead a critical component of intelligence - prediction\". We can do all sorts of wonderful things with machine learning. The only requirement is that we frame our problems as prediction ones. Want to translate from english to portuguese? Then build a ML model that predicts portuguese sentences when given english sentences. Want to recognize faces? Then build a ML model that predicts the presence of a face in a subsection of a picture. Want to build a self driving car? Then build one ML model to predict the direction of the wheel and the pressure on the brakes and accelerator when presented with images and sensors from the surroundings of a car.\n",
    "\n",
    "However, ML is not a panacea. It can perform wonders under very strict boundaries and still fail miserably if the data it's using deviates a little from what the model is accustomed to. To give another example from Prediction Machines, \"in many industries, low prices are associated with low sales. For example, in the hotel industry, prices are low outside the tourist season, and prices are high when demand is highest and hotels are full. Given that data, a naive prediction might suggest that increasing the price would lead to more rooms sold.” \n",
    "\n",
    "ML is notoriously bad at this inverse causality type of problems. They require us to answer \"what if\" questions, what Economists call counterfactuals. What would happen if instead of this price I'm currently asking for my merchandise, I use another price? What would happen if instead of this low fat diet I'm in, I do a low sugar one? If you work in a bank, giving credit, you will have to figure out how changing the customer line changes you'r revenue. Or if you work at the local government, you might be asked to figure out how to make the schooling system better. Should you give tablets to every kid because the era of digital knowledge tells you to? Or should you build an old fashioned library? \n",
    "\n",
    "At the heart of these questions there is a causal inquiry we wish to know the answer. Causal questions permeate everyday problems, like figuring out how to make sales go up, but they also play an important role on dilemmas that are very personal and dear to us: do I have to go to an expensive school to be successful in life (does education cause earnings)? Does imigration lowers my chances of getting a job (does imigration causes unemployment to go up)? Does money transfer to the poor lower the crime rate? It doesn't matter the field you are in, it is very likely you had or will have to answer some type of causal question. Unfortunately for ML, we can't rely on correlation type predictions to tackle them.\n",
    "\n",
    "Answering this kind of question is tougher than most people appreciate. Your parents have probably repeated to you that \"association is not causation\", \"association is not causation\". But actually explaining why that is the case is a bit more involved. This is what this introduction to causal inference is all about. The rest of this series is dedicated to **figuring how to make association be causation**.\n",
    "\n",
    "## When Association IS Causation\n",
    "\n",
    "Intuitively, we kind of know why association is not causation. If someone tells you that schools that give tablets to its students perform better than those who don't, you can quickly point out that it is probably the case that those schools with the tablets are richer. As such, they would do better than average even without the tablets. Because of this, we can't conclude that giving tablets to kids during classes will cause an increase in their academic performance. We can only say that tablets in school are associated with high academic performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "jupyter": {
     "outputs_hidden": true,
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "In /Users/matheus.facure/miniconda3/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n",
      "The text.latex.unicode rcparam was deprecated in Matplotlib 3.0 and will be removed in 3.2.\n",
      "In /Users/matheus.facure/miniconda3/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n",
      "The savefig.frameon rcparam was deprecated in Matplotlib 3.1 and will be removed in 3.3.\n",
      "In /Users/matheus.facure/miniconda3/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n",
      "The pgf.debug rcparam was deprecated in Matplotlib 3.0 and will be removed in 3.2.\n",
      "In /Users/matheus.facure/miniconda3/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n",
      "The verbose.level rcparam was deprecated in Matplotlib 3.1 and will be removed in 3.3.\n",
      "In /Users/matheus.facure/miniconda3/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n",
      "The verbose.fileo rcparam was deprecated in Matplotlib 3.1 and will be removed in 3.3.\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy.special import expit\n",
    "import seaborn as sns\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import style\n",
    "\n",
    "style.use(\"fivethirtyeight\")\n",
    "\n",
    "np.random.seed(123)\n",
    "n = 100\n",
    "tuition = np.random.normal(1000, 300, n).round()\n",
    "tablet = np.random.binomial(1, expit((tuition - tuition.mean()) / tuition.std())).astype(bool)\n",
    "enem_score = np.random.normal(200 - 50 * tablet + 0.7 * tuition, 200)\n",
    "enem_score = (enem_score - enem_score.min()) / enem_score.max()\n",
    "enem_score *= 1000\n",
    "\n",
    "data = pd.DataFrame(dict(enem_score=enem_score, Tuition=tuition, Tablet=tablet))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAHwCAYAAAAVRjAeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3debgkdX3v8fcHRhTC4hI1BBAwjAu5PjEsBuMCETRXjTImbqiIBjVRL66JQa+J5gZRNBHjEzUuRBENEbcBEiOyiKhRUAxxQxxUhBEUFBhBR1nme/+oGmyaPksz53Sf35z363n6OV1Vv6r+dp0+/elf1e9Up6qQJKkFW0y7AEmS5svQkiQ1w9CSJDXD0JIkNcPQkiQ1w9CSJDXD0JIWWZKzk7xn2nUMSvKsJDfN0eaAJJVk50WqoZI8YzG2PeKxdusf76GTeDwtHkNrmUvyvv6Pefh2/Yg2x45Y/1ZvPEkumWF7/z6izREjtveWftkZi/F8Nxez7OdbblOq6xljPPaOwEcW6HEfk+STSX6SZH2Si5L8c5L7LMT2tXQYWgL4LN0byODt3kNt1gMvnOebwDEjtjf8ifpS4LmDM5LcCTgU+P6Y9S+6JFtNu4Yh+/KrfbtXP+9PuPU+X9Kq6odV9YtN3U6SvwFOBb4L/DFwP+BPgRuAozZ1+1paDC0B3NC/gQzerhxq81/A+cCb5rG960ds79qhNv8G3DvJ7w3MeyJwDfCZuR4gyXOSXJjkF/2n63MGD2Ml2bv/5P3TJNcnOW/wsZIcluSbSX6ZZG2So5KsGFh+dpLjkvxdkiuAH/Tz90jy0STXJrkmyaeSPGAe+2SLJG9I8uO+pvck2brf5rP77W0z9Bxfk+R7STK8saq6auO+Ba7qZ189MO+Ifv/8PMllfa9jhxH78aAk3+j343lJ9hpuM9R+xuef5ADghP7+xh7f+2bZ1nAvvZK8IMkJSa7r637FHPXsDfwt8H+r6gVV9Zmq+n5Vfb6qXgT82Szrvm62fZRk+yTvTfLD/nVyWZI3Dyx/aJLP97Vel+R/kvzhbPVq0xlaGsdLgccl+YMF2NZ1dME12Nt6HvAeYNbDS/0b1T8DrwfuCxwAvH9g+W8D59AF4COA3wWOpX+9J3ks8C90b7APAF4OvBB4zdBDPRm4O3Ag8Igk9wQ+B1wJPAzYD7gIODvJ3ed4vk8E7tav93Tg8XQ9Uvr9UMCTBp7DFsCzgffU7bvW2nq6/bkn8Cy6ffTWoTZbAG8EXgA8qH9e/zEcngM1zfX8/wv4P33zjb29F49Z92vofncPpPuAdMwcr7dDgZ8Dbx61sKqumWXdufbRUXS92IOBlcBTgAsBkmwJnAKc27fZC3htX4sWU1V5W8Y34H3ATcD1Q7dTh9qc0d8/EfhvYIt+uoBnDLS9BPjliO29YqjNq+neKK8HtqM7pHMDcM/Bx5uh5icA64DtZ1h+AvA/G2scsfyzwElD815M9ya2VT99NvDtwW3QvSl9cWi9AN8BXjJLvWf3z3nLgXnP6/fTr/XTbwU+N7D8D4EbgR3n8Tvcuf89HDDHPvvlwO/tWf06Bw60uUv/+3hOP31A32bn+T5/usPANc/X3vBrp4C3DrX5FvD6WbbxCeCr83is3frtP3SMfXQy8L4Z2t5lrn3ubXFutxwO0bJ2LnDY0LyZPjEeSfdG8iy63soobwPePjTvJ8ONquq8JGuAQ+hC69Sq+tGIo2HDTqc7f/G9JKcDZwEfq6of98v3Bj5ZVRtmWP+3gQ8NzfsMcCfgt+g/TQPnD21jX2DvDAxS6W1N90l8NudV1c0D058Htuof76vAO4GvJ9mzqr5J1wP9j6q6Yo7tjpTkj4GXAHsA29P1qrYCfgO4fKDpFzbeqaprklxI1/MYZVOe/3xdMDT9A7oPMjMJc/TMZ1xx7n30duCjSfYBzgQ+CZxWVRv6ffUe4LQkZ9G9fj5eVRfdnlo0fx4eFMD6qrp46Hb5qIZV9X26Q21HJdl2hu1dPWJ7Mx2meTfwfOCZwLvmU2xVXQ/sQ/fJ+NvAnwMX94cNb2k212aGpjNi/s+G2mxB9+b1wKHbfel6IeO4VTJX1TfoDr09J8k96A4fzmt/3GbD3bm7D9MdZnsC3aGrP+8XzzWgZLZPDAv5/Gdyw9B0Mfv71EXAHuMOlJnPPqqq04B7Aa+j+0DzAeCs/tAgVfVcug9IpwP7033omPEcmhaGoaXb4/V0r52/WoBtfYDuU/r1dH/881JVN1fVOVX1N3RvHFcAT+sXnw8c1J8XGuUbdG8ygx5Od3jwu7M87Jfpemk/GBHKV82yHsC+G9/seg+me4P+zsC8d9KF9/OAH9J9sr89Hgr8uKpeXVXnVtW36Q4hjrLfxjtJ7kzX471whrbzef439NvacoZtLLQPANsALxu1MMldZlhvXvuoqq6uqhOr6s+Ax9K9bvYcWP71qnpzVT0aOI7ud6dFZGgJYKskvzHiNvJTd1VdB/w13QCGUbYdsa27zbCtnwI7AQ+Y5XDerSQ5OMlL040QvBewCtgF+Gbf5I10QfjBJPsk+a0kT0ry4H7564E/SXJkkvskeTJdT+Efqmr4k/6gfwK2BFYneVi6f1h9aD8K7ffnKPtuwNuS3L8fCPJ3wLurarA3t/F/lv4aOG6++2OEi4C7Jzk8yb2TPJNusMWwAt6Y5OH9CMD30/Uu/3WG7c7n+X+v//n4JHefpTe+IKrqy8D/A45O8rb+ueya5MHp/q/wnTOsOuc+6p/XHye5b5KVdANorgcu7UdRHtM//13719bD+NVrUItl2ifVvE33RjfooWa4/fpAmzOG1tuCbrDDqIEYo7b19aE2r56jptkGYjyc7jzWVcAvgDV059oy0OZBwBl0b8LX0Z23e9DA8sPoehQ30J03eR2wYmD52XQj94Yfe1fgg/1j/5Luf8o+AOw+S71n053/exPdub3r+ultRrQ9FriZfvDDPH+HtxmIQReKP+qf/yfozhsWsFu//Fl0A3Ae1e+HXwJfAvYZ2MYBDAzEmO/zB97SP/YGZhjI0LcbNRDjGUNtzphtGwPtHgd8im7E6C/oQuntwB798t0YGogxj33018DX6YJqHd15q4f2y3YEPgas7ffD5XSHuneY9t/05n5L/wuQtAQkOQnYuqoeN+1apKXI0YPSEtCfe3kY3aCAR065HGnJMrSkpeG/6c57vbGqzp5yLdKS5eFBSVIzmu1prVu3zrSVpM3YDjvscJsRzA55lyQ1w9CSJDXD0JIkNcPQkiQ1w9CSJDXD0JIkNcPQkiQ1w9CSJDXD0JIkNcPQkiQ1w9CSJDXD0JIkNcPQkiQ1w9CSJDXD0JIkNcPQkiQ1w9CSJDXD0JIkNWPFtAuQJIBVq1ZNu4RmrF69etolTI09LUlSMwwtSVIzDC1JUjNSVdOu4XZZt25dm4VLWvKGz68t53NI07TDDjtkeJ49LUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMyYWWklemuQbSb6e5MQkd0qye5Jzk6xJ8qEkW/Vt79hPX9wv321SdUqSlq6JhFaSnYAXAftU1f8CtgSeChwDHFtVK4FrgMP7VQ4HrqmqPYBj+3aSpGVukocHVwBbJ1kBbANcATwC+Ei//HhgVX//4H6afvmBSTLBWiVJS9CKSTxIVf0gyd8DlwLrgU8B5wPXVtVNfbO1wE79/Z2Ay/p1b0qyDrgb8ONR21+zZs0iVi9pufM9ZnJWrlw56/KJhFaSu9D1nnYHrgU+DDx6RNPauMosy25jricpSZvC95ilY1KHBw8CvldVV1XVjcDHgN8H7twfLgTYGbi8v78W2AWgX74DcPWEapUkLVGTCq1Lgf2SbNOfmzoQ+CbwaeCJfZvDgJP7+6f00/TLz6qqGXtakqTlYSKhVVXn0g2o+Arwtf5x3wX8FfCyJBfTnbM6rl/lOOBu/fyXAUdOok5J0tI2kXNaAFX1GuA1Q7O/CzxoRNtfAE+aRF2SpHZ4RQxJUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzVky7AM1s1apV0y6hGatXr552CZImwJ6WJKkZhpYkqRmGliSpGZ7TWsKW2nma4XNsS60+SZs/e1qSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmTCy0ktw5yUeSfCvJhUkenOSuSU5Psqb/eZe+bZK8NcnFSb6aZK9J1SlJWrom2dP6R+CTVXU/4HeAC4EjgTOraiVwZj8N8GhgZX97HvCOCdYpSVqiJhJaSbYHHg4cB1BVN1TVtcDBwPF9s+OBVf39g4H3V+eLwJ2T7DiJWiVJS9eKCT3OvYGrgPcm+R3gfODFwD2r6gqAqroiyT369jsBlw2sv7afd8Woja9Zs2ax6tYs3O9aLnytT87KlStnXT6p0FoB7AUcUVXnJvlHfnUocJSMmFczNZ7rSWpxuN+1XPhaXzomdU5rLbC2qs7tpz9CF2I/2njYr/955UD7XQbW3xm4fEK1SpKWqImEVlX9ELgsyX37WQcC3wROAQ7r5x0GnNzfPwV4Zj+KcD9g3cbDiJKk5WtShwcBjgA+mGQr4LvAs+lC86QkhwOXAk/q234CeAxwMfDzvq2kTbBq1aq5G2kk9914Vq9evWjbnlhoVdUFwD4jFh04om0BL1z0oiRJTfGKGJKkZhhakqRmTPKclqQl5ITdt552CdpMHPq99RN7LHtakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZowVWkkemeS4JKf20/skecTilCZJ0q3NO7SSHAG8A1gDPLyfvR44ahHqkiTpNsbpab0EOKiq3gBs6Od9C7jvglclSdII44TWdsBl/f3qf94BuGFBK5IkaQbjhNY5wJFD814EfHrhypEkaWYrxmh7BHBqkucC2yW5CPgp8LhFqUySpCHzCq0kWwD3Bx4GPADYle5Q4XlVtWG2dSVJWijzCq2q2pDk5KraDjivv0mSNFFjndNKst+iVSJJ0hzGOaf1feA/k5xMd2hw4whCqupvFrowSZKGjRNaWwOr+/s7D8yvEW0lSVpw8w6tqnr2YhYiSdJcxulpkWQlcAiwE/AD4MSqWrMYhUmSNGycaw8+DjgfuB9wNd3lm76c5PGLVJskSbcyTk/raODgqrrlChhJDgD+CThlgeuSJOk2xhnyvjPw2aF5n+PWgzIkSVo044TWBcDLh+a9rJ8vSdKiG+fw4PPprj34Yrr/09oF+BngOS1J0kSMM+T9W0nuD+wH/CZwOXBuVd24WMVJkjRo3qGV5IHAT6rqcwPzdkly16r6n0WpTpKkAeOc0/oA3Zc+DtoKOGHhypEkaWbjhNa9quq7gzOq6jvAbgtakSRJMxgntNYm2WtwRj99+cKWJEnSaOOMHjwWODnJG4HvAL8F/AXwusUoTJKkYeOMHnx3kmuBw+mGu18GvLyqPrJYxUmSNGisC+ZW1YeBDy9SLZIm6NDvrZ92CdLYxrlg7iH9/2mR5D5JPpPkrCT3W7zyJEn6lXEGYhxFd3V3gH8AvgScA7x9oYuSJGmUcQ4P3r2qfpTkTsBDgScCNwI/XpTKJEkaMk5oXZVkD+ABwJeq6pdJtgGyOKUtvlWrVk27hKa5/+Zv9erV0y7hNk7Yfetpl6DNxCTPj44TWn9H9yWQNwNP6ecdCHgJJ0nSRMz7nFZVvQ/YEdi5qk7vZ58LPHVjmyQPWdDqJEkaMO6Q958PTV851OQ/ge03tShJkkYZK7TmodnzWwBrHvWqaZegzcTKTx097RKkzdI4Q97noxZ4e5Ik3WKhQ0uSpEVjaEmSmrHQodX0OS1J0tI29kCMJNsD2w7Oq6rL+5/bLVBdkiTdxrxDK8lBwLuAXbl1j6qALRe4LkmSbmOcw4PHAUcDOwB3GLhttQh1SZJ0G+McHrwT8N6qunmxipEkaTbj9LSOBV6RxMEWkqSpGKen9VHgNOCVSW71dSRVde8FrUqSpBHGCa2PAJ8FPgz4Pd2SpIkbJ7R2B363qjYsVjGSJM1mnHNaJwOPWKxCJEmayzg9rTsCpyT5LPCjwQVV9cwFrUqSpBHGCa1v9DdJkqZi3qFVVX+7qQ+WZEvgy8APquqPkuwO/BtwV+ArwKFVdUOSOwLvB/YGfgI8paou2dTHlyS1bawL5iZ5ZJLjkpzaT++TZJzzXC8GLhyYPgY4tqpWAtcAh/fzDweuqao96P4/7Jhx6pQkbZ7mHVpJjgDeAawBHt7PXg8cNc/1dwYeC7ynnw7dwI6P9E2OB1b19w/up+mXH+g/NUuSxulpvQQ4qKreAGwc9v4t4L7zXP8twCsG1r0bcG1V3dRPrwV26u/vBFwG0C9f17eXJC1j4wzE2I4+SOiu7A7dBXNvmGvFJH8EXFlV5yc5YOPsEU1rHstuY82aNXOVIE2Vr1EtJ5vyel+5cuWsy8cJrXOAI4HXDcx7EfDpeaz7EODxSR5Dd+Hd7el6XndOsqLvTe0MXN63XwvsAqxNsoLuyvJXz7TxuZ6kNG2+RrWcLObrfZzDg0cAT0hyCbBdkouAJwEvm2vFqnplVe1cVbsBTwXOqqqn0wXeE/tmh9H9AzPAKf00/fKzqmrGnpYkaXkYZ8j7FUn2BR4E3IvuUOF5m3hZp78C/i3JUcB/031nF/3PE5JcTNfDeuomPIYkaTMxzuFB+t7OuUm+tHFeki3GCa6qOhs4u7//XboQHG7zC7penCRJtxhnyPteSb6Q5GfAjf3tpv6nJEmLbpye1vHAqcCfAj9fnHIkSZrZOKG1K/B/HRAhSZqWcUYPfhx41GIVIknSXMbpad0J+HiSzwE/HFzgV5NIkiZhnND6Zn+TJGkqxvpqkiSPBA4B7tF/tcg+dFe3kCRp0d2eq7x/G3hYP3veV3mXJGlTTfIq75IkbZJxQut2X+VdkqSFME5obbzK+6D5XuVdkqRNNs7owSOAU5M8l19d5f2nwOMWpTJJkobcnqu870t3dYyFuMq7JEnzdnuu8n5ef5MkaaLGOaclSdJUGVqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZhhakqRmGFqSpGYYWpKkZox1RYzN3cpPHT3tEiRJs7CnJUlqhqElSWqGoSVJaobntAasedSrpl2CNhOeH5UWhz0tSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzvIyTtEwd+r310y5BGps9LUlSMwwtSVIzDC1JUjM8pyUtE6tXr552Cc1YtWrVrabdd0uHPS1JUjMMLUlSMwwtSVIzDC1JUjMMLUlSMwwtSVIzHPI+YOWnjp52CZKkWdjTkiQ1w9CSJDXD0JIkNWNZn9Py0izj8dI2kqbNnpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGRMJrSS7JPl0kguTfCPJi/v5d01yepI1/c+79POT5K1JLk7y1SR7TaJOSdLSNqme1k3Ay6vq/sB+wAuT7AkcCZxZVSuBM/tpgEcDK/vb84B3TKhOSdISNpHQqqorquor/f3rgAuBnYCDgeP7ZscDq/r7BwPvr84XgTsn2XEStUqSlq6Jn9NKshvwu8C5wD2r6grogg24R99sJ+CygdXW9vMkScvYikk+WJJtgY8CL6mqnyaZsemIeTVT4zVr1ixAdRqX+13Lha/1yVm5cuWsyycWWknuQBdYH6yqj/Wzf5Rkx6q6oj/8d2U/fy2wy8DqOwOXz7TtuZ6kFof7XcuFr/WlY1KjBwMcB1xYVW8eWHQKcFh//zDg5IH5z+xHEe4HrNt4GFGStHxNqqf1EOBQ4GtJLujnvQp4A3BSksOBS4En9cs+ATwGuBj4OfDsCdUpSVrCJhJaVfU5Rp+nAjhwRPsCXrioRUmSmuMVMSRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc0wtCRJzTC0JEnNMLQkSc1YMe0CNLNVq1ZNu4RZLaX6Vq9ePe0SJE2APS1JUjMMLUlSMwwtSVIzPKe1hHmeRpJuzZ6WJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpbGcvXVV/OqV72Ka665ZtqlSFqGDC2N5aSTTuLCCy/kQx/60LRLkbQMLdnQSvK/k1yU5OIkR067HnW9rLPOOouq4qyzzrK3JWnilmRoJdkSeBvwaGBP4JAke063Kp100kls2LABgA0bNtjbkjRxqapp13AbSR4MvLaq/rCffiVAVb1+Y5t169YtvcI3c4cccgjr16+/ZXrrrbfmxBNPnGJF2pysWrVq2iU0Y7l8bdEOO+yQ4XlLsqcF7ARcNjC9tp+nKdp///1ZsaL7CrYVK1aw//77T7kiScvNUv0SyNukKzBjz2rNmjWLWIo22nfffTnjjDMASMK+++7rvpemYHP+u1u5cuWsy5dqaK0FdhmY3hm4fKbGcz1JLZyDDjqI0047jQMPPJC999572uVIy9Jyfs9bqqH1JWBlkt2BHwBPBZ423ZIE8OQnP5lLL72UpzzlKdMuRZuZ5XKeRptmSQ7EAEjyGOAtwJbAv1TV6waXOxBDkjZvowZiLNnQmouhJUmbt5ZGD0qSdBuGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGliSpGYaWJKkZhpYkqRnNfgmkJGn5saclSWqGoSVJasaKaReg6UtyM/C1gVmrquqSGdpeX1XbTqQwaZEkuRtwZj/5G8DNwFX99IOq6oapFKY5eU5LYwWRoaXNTZLXAtdX1d8PzFtRVTdNryrNxMODuo0k2yY5M8lXknwtycEj2uyY5JwkFyT5epKH9fMfleQL/bofTmLAqQlJ3pfkzUk+DRyT5LVJ/mJg+deT7Nbff0aS8/rX/zuTbDmlspcdQ0sAW/d/fBck+TjwC+AJVbUX8AfAPyTJ0DpPA06rqgcCvwNckOTXgVcDB/Xrfhl42eSehrTJ7kP3+n35TA2S3B94CvCQ/vV/M/D0CdW37HlOSwDr+z8+AJLcATg6ycOBDcBOwD2BHw6s8yXgX/q2q6vqgiT7A3sCn+8zbivgCxN6DtJC+HBV3TxHmwOBvYEv9a/zrYErF7swdQwtjfJ04O7A3lV1Y5JLgDsNNqiqc/pQeyxwQpI3AdcAp1fVIZMuWFogPxu4fxO3Phq18W8gwPFV9cqJVaVbeHhQo+wAXNkH1h8Auw43SLJr3+bdwHHAXsAXgYck2aNvs02S+0ywbmkhXUL3uibJXsDu/fwzgScmuUe/7K7934MmwJ6WRvkgcGqSLwMXAN8a0eYA4C+T3AhcDzyzqq5K8izgxCR37Nu9Gvj24pcsLbiPAs9McgHd4fBvA1TVN5O8GvhUki2AG4EXAt+fWqXLiEPeJUnN8PCgJKkZhpYkqRmGliSpGYaWJKkZhpYkqRmGlrQEJamN/+82Ytmzknxu0jVJS4GhJS2iJNcP3DYkWT8wvejXq+svAnvUYj+ONCn+c7G0iAa/xqW/HNZzquqM6VUktc2eljQFSR7Uf4XLtUmuSPJPSbYaavaYJN9N8uMkb+qvvjBqW/dLcnqSq5NclOTJ/fzn0V1H8hV9z+7URX5a0qIztKTpuBl4KfDrwIPprhz+gqE2TwD2obv+3cHAnw5vJMmvAacD/wrcAzgEeHuS366qd9FdkuuNVbVtVT1ukZ6LNDGGljQFVXV+VX2xqm6qqkuAdwL7DzU7pqqurqpLgbfQBdKwPwIuqar39tv6Ct018564mPVL0+I5LWkK+qvfv5muJ7UN3d/i+UPNLhu4/33gN0dsalfg95JcOzBvBXDCwlUrLR32tKTpeAfd1fNXVhiuq2EAAADJSURBVNX2wKvovqdp0C4D9+8FXD5iO5cBn6mqOw/ctq2q5/fLvSK2NiuGljQd2wE/Ba5Pcj/g+SPa/GWSuyTZBXgx8KERbf4duE+SQ5Pcob/t238lPMCPgHsvxhOQpsHQkqbjL4CnAdcB72Z0IJ1Md8jwAuA/6L5s81aq6jrgUcBT6XpiPwSOATZ+n9lxwJ79KMXVC/wcpInz+7QkSc2wpyVJaoahJUlqhqElSWqGoSVJaoahJUlqhqElSWqGoSVJaoahJUlqxv8HBtJMGRDg98UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(6,8))\n",
    "sns.boxplot(y=\"enem_score\", x=\"Tablet\", data=data).set_title('ENEM score by Tablet in Class')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get beyond simple intuition, let's first establish some notation. This will be our common language to speak about causality. Think of it as the common tongue we will use to identify other brave and true causal warriors and that will compose our cry in the many battles to come.\n",
    "\n",
    "Let's call \\\\(T_i\\\\) the treatment intake for unit i. \n",
    "\n",
    "$\n",
    "T_i=\\begin{cases}\n",
    "1 \\ \\text{if unit i received the treatment}\\\\\n",
    "0 \\ \\text{otherwise}\\\\\n",
    "\\end{cases}\n",
    "$\n",
    "\n",
    "The treatment here doesn't need to be a medicine or anything from the medical field. Istead, it is just a term we will use to denote some intervention for which we want to know the effect. In our case, the treatment is giving tablets to students. As a side note, you might sometimes see \\\\(D\\\\) instead of \\\\(T\\\\) to denote the treatment.\n",
    "\n",
    "Now, let's call \\\\(Y_i\\\\) the observed outcome variable for unit i.\n",
    "\n",
    "$\n",
    "Y_i\n",
    "$\n",
    "\n",
    "The outcome is our variable of interest. We want to know if the treatment has any influence in it. In our tablet example, it would be the academic performance.\n",
    "\n",
    "Here is where things get interesting. The **fundamental problem of causal inference** is that we can never observe the same unit with and without treatment. It is as if we have two diverging roads and we can only know what lies ahead of the one we take. As in Robert Frost poem:\n",
    "\n",
    "```\n",
    "Two roads diverged in a yellow wood,\n",
    "And sorry I could not travel both\n",
    "And be one traveler, long I stood\n",
    "And looked down one as far as I could\n",
    "To where it bent in the undergrowth;\n",
    "```\n",
    "\n",
    "To wrap our heads around this, we will talk a lot in term of **potential outcomes**. They are potential because they didn't actually happen. Instead they denote **what would have happened** in the case some treatment was taken. We sometimes call the potential outcome that happened, factual, and the one that didn't happen, counterfactual.\n",
    "\n",
    "As for the notation, we use an additional subscript:\n",
    "\n",
    "\\\\(Y_{0i}\\\\) is the potential outcome for unit i without the treatment. \n",
    "\n",
    "\\\\(Y_{1i}\\\\) is the potential outcome for **the same** unit i with the treatment.\n",
    "\n",
    "Sometimes you might see potential outcomes represented as functions \\\\(Y_i(t)\\\\), so beware. \\\\(Y_{0i}\\\\) could be \\\\(Y_i(0)\\\\) and \\\\(Y_{1i}\\\\) could be \\\\(Y_i(1)\\\\). Here, we will use the subscript notation most of the time.\n",
    "\n",
    "![img](./data/img/intro/potential_outcomes.png)\n",
    "\n",
    "Back to our example, \\\\(Y_{1i}\\\\) is the academic performance for student i if he or she has tablets at the classroom. Whether this is or not the case, it doesn't matter for \\\\(Y_{1i}\\\\). It is the same regardless. If student i gets the tablet, we can observe \\\\(Y_{1i}\\\\). If not, we can observe \\\\(Y_{0i}\\\\). Notice how in this last case, \\\\(Y_{1i}\\\\) is still defined, we just can't see it. In this case, it is a counterfactual potential outcome.\n",
    "\n",
    "With potential outcomes, we can define the individual treatment effect:\n",
    "\n",
    " \\\\(Y_{1i} - Y_{0i}\\\\)\n",
    " \n",
    "Of course, due to the fundamental problem of causal inference, we can never know the individual treatment effect because we only observe one of the potential outcomes. For the time being, let's focus on something easier than estimating the individual treatment effect. Instead, lets focus on the **average treatment effect**, which is defined as follows.\n",
    "\n",
    "\\\\(ATE = E[Y_1 - Y_0]\\\\)\n",
    "\n",
    "Another easier quantity to estimate is the **average treatment effect on the treated**:\n",
    "\n",
    "\\\\(ATET = E[Y_1 - Y_0 | T=1]\\\\)\n",
    "\n",
    "Now, I know we can't see both potential outcomes, but just for the sake of the argument, let's suppose we could. Pretend that the causal inference deity is pleased with the many statistical battles that we fought and has rewarded us with godlike powers to see the alternative potential outcomes. With that power, say we collect data on 4 schools. We know if they gave tablets to its students and their score on some annual academic test. Here, tablets are the treatment, so \\\\(T=1\\\\) if the school gives tablets to its kids. \\\\(Y\\\\) will be the test score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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>i</th>\n",
       "      <th>y0</th>\n",
       "      <th>y1</th>\n",
       "      <th>t</th>\n",
       "      <th>y</th>\n",
       "      <th>te</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>500</td>\n",
       "      <td>450</td>\n",
       "      <td>0</td>\n",
       "      <td>500</td>\n",
       "      <td>-50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>600</td>\n",
       "      <td>600</td>\n",
       "      <td>0</td>\n",
       "      <td>600</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>800</td>\n",
       "      <td>600</td>\n",
       "      <td>1</td>\n",
       "      <td>600</td>\n",
       "      <td>-200</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>700</td>\n",
       "      <td>750</td>\n",
       "      <td>1</td>\n",
       "      <td>750</td>\n",
       "      <td>50</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   i   y0   y1  t    y   te\n",
       "0  1  500  450  0  500  -50\n",
       "1  2  600  600  0  600    0\n",
       "2  3  800  600  1  600 -200\n",
       "3  4  700  750  1  750   50"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(dict(\n",
    "    i= [1,2,3,4],\n",
    "    y0=[500,600,800,700],\n",
    "    y1=[450,600,600,750],\n",
    "    t= [0,0,1,1],\n",
    "    y= [500,600,600,750],\n",
    "    te=[-50,0,-200,50],\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \\\\(ATE\\\\) here would be the mean of the last column, that is, of the treatment effect:\n",
    "\n",
    "\\\\(ATE=(-50 + 0 - 200 + 50)/4 = -50\\\\)\n",
    "\n",
    "This would mean that tablets reduced the academic performance of students, on average, by 50 points. The \\\\(ATET\\\\) here would be the mean of the last column when \\\\(T=1\\\\):\n",
    "\n",
    "\\\\(ATET=(- 200 + 50)/2 = -75\\\\)\n",
    "\n",
    "This is saying that, for the schools that were treated, the tablets reduced the academic performance of students, on average, by 75 points. Of course we can never know this. In reality, the table above would look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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>i</th>\n",
       "      <th>y0</th>\n",
       "      <th>y1</th>\n",
       "      <th>t</th>\n",
       "      <th>y</th>\n",
       "      <th>te</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>500.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0</td>\n",
       "      <td>500</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>600.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0</td>\n",
       "      <td>600</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>NaN</td>\n",
       "      <td>600.0</td>\n",
       "      <td>1</td>\n",
       "      <td>600</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>NaN</td>\n",
       "      <td>750.0</td>\n",
       "      <td>1</td>\n",
       "      <td>750</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   i     y0     y1  t    y  te\n",
       "0  1  500.0    NaN  0  500 NaN\n",
       "1  2  600.0    NaN  0  600 NaN\n",
       "2  3    NaN  600.0  1  600 NaN\n",
       "3  4    NaN  750.0  1  750 NaN"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(dict(\n",
    "    i= [1,2,3,4],\n",
    "    y0=[500,600,np.nan,np.nan],\n",
    "    y1=[np.nan,np.nan,600,750],\n",
    "    t= [0,0,1,1],\n",
    "    y= [500,600,600,750],\n",
    "    te=[np.nan,np.nan,np.nan,np.nan],\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is surely not ideal, you might say, but can't I still take the mean of the treated and compare it to the mean of the untreated? In other words, can't I just do \\\\(ATE=(600+750)/2  - (500 + 600)/2  = -75\\\\)? Well, no! Notice how different the results are. That's because you've just committed the gravest sin of mistaking association for causation. To understand why, let's take a look into the main enemy of causal inference.\n",
    "\n",
    "## Bias\n",
    "\n",
    "Bias is what makes association different from causation. Fortunately, it too can be easily understood with our intuition. Let's recap our tablets in the classroom example. When confronted with the claim that schools that give tablets to their kids do better on tests score, we can rebut it by saying those schools will probably do better anyway, even without the tables. That is because they probably have more money than the other schools, and hence can pay better teachers, afford better classrooms, and so on. In other words, it is the case that treated schools (with tablets) are not comparable with untreated schools. \n",
    "\n",
    "To say this in potential outcome notation is to say that \\\\(Y_0\\\\) of the treated is different from the \\\\(Y_0\\\\) of the untreated. Remember that the \\\\(Y_0\\\\) of the treated **is counterfactual**. We can't observe it, but we can reason about it. In this particular case, we can even leverage our understanding of how the world works to go even further. We can say that, probably, \\\\(Y_0\\\\) of the treated is bigger than \\\\(Y_0\\\\) of the untreated schools. That is because schools that can afford to give tablets to their kids can also afford other factors that contribute to better test scores. Let this sink in for a moment. It takes some time to get used to talking about potential outcomes. Read this paragraph again and make sure you understand it.\n",
    "\n",
    "With this in mind, we can show with very simple math why is it the case that association is not causation. Association is measured by \\\\(E[Y|T=1] - E[Y|T=0]\\\\). In our example, this is the average test score for the schools with tablets minus the average test score for those without it. On the other hand, causation is measured by \\\\(E[Y_1 - Y_0]\\\\).\n",
    "\n",
    "To see how they relate, let's take the association measurement and replace the observed outcomes with the potential outcomes. For the treated, the observed outcome is \\\\(Y_1\\\\). For the untreated, the observed outcome is \\\\(Y_0\\\\).\n",
    "\n",
    "$\n",
    "E[Y|T=1] - E[Y|T=0] = E[Y_1|T=1] - E[Y_0|T=0]\n",
    "$\n",
    "\n",
    "Now, let's add and subtract \\\\(E[Y_0|T=1]\\\\). This is a counterfactual outcome. It tells what would have been the outcome of the treated, had they not received the treatment.\n",
    "\n",
    "$\n",
    "E[Y|T=1] - E[Y|T=0] = E[Y_1|T=1] - E[Y_0|T=0] + E[Y_0|T=1] - E[Y_0|T=1]\n",
    "$\n",
    "\n",
    "Finally, we reorder the terms, merge some expectations, and lo and behold:\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",
    "This simple piece of math encompasses all the problems we will encounter in causal questions. I cannot stress how important it is that you understand every aspect of it. If you're ever forced to tattoo something on your arm, this equation should be a good candidate for it. It's something to hold onto very dearly and really understand what is telling us, like some sacred text that can be interpreted 100 different ways. In fact, let's take a deeper look. Let's break it down into some of its implications. First, this equation tells why association is not causation. As we can see, association is equal to the treatment effect on the treated plus a bias term. **The bias is given by how the treated and control group differ before the treatment, that is, in case neither of them has received the treatment**. We can now say precisely why we are suspicious when someone tells us that tablets in the classroom boost academic performance. We think that, is this example, \\\\(E[Y_0|T=0] < E[Y_0|T=1]\\\\), that is, schools that can afford to give tablets to their kids are better than those that can't, **regardless of the tablets treatment**.\n",
    "\n",
    "Why does this happen? We will talk more about that once we enter confounding, but for now you can think of bias arising because many things we can't control are changing together with the treatment. As a result, the treated and untreated schools don't differ only on the tablets. They also differ on the tuition cost, location, teachers...\n",
    "For us to say that tablets in the classroom increase academic performance, we would need for schools with and without them to be, on average, similar to each other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqgAAAGECAYAAAD6PtKMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd5hU5dnH8e99ztRtdFBAsQC2aOwVEMTYW+wmMRpNYnyNqRpjNDHFFlvM66tGY4kaEyNqRGNvqKRYY4kYREQpgnTYNvU87x8zwJbZXRZmZ2d3f5/rmoud8zxzzj1ndtl7n2rOOUREREREyoXX3QGIiIiIiDSlBFVEREREyooSVBEREREpK0pQRURERKSsKEEVERERkbKiBFVEREREyooSVBGRbmJm08zstu6OoyUz+4OZPbse9T42s4tLEZOI9C1KUEWkoHyS4go86grU+U2B1zsz+0qT5x+3cb6/FahzboHzXZ8v6zBx6svauc9rH+txmu8CJzQ5521mNq1AvT2AVp99VzGzr5jZS2a2yszqzew/ZnaVmY0o4jUyZnZ6sc4nIhtGCaqItOdlYNMWj61a1GkEzjGzsetxvl8XON9XWtSZC3yj6QEziwGnAp90Mv4uZ2aR7o6hhT1Yd293zR87jub3vF3OuVXOuRXrUW+Jc65+I2Jdb2Z2O3A78BJwKLA98B1gE+CHpYhBREpHCaqItCflnFvU4rG4RZ1/AG8AV6/H+eoKnG9lizr3AVuZ2V5Njh0PrABe7OgCZvZ1M3vfzBJmtizf4jaySfluZvakma02szoze7XptczsNDObYWZJM5tvZpeaWahJ+TQzu93MfmVmC4EF+eOjzexBM1tpZivM7Gkz23E97olnZlea2dJ8TLeZWTx/zq/lz1fR4j1eYmZzzMxaniyfNC5yzi0CluQPL29y7F8tu+VbtpA27eI3s58DZwL7N2mFPT1f1qyL38yqzewWM1uSv/+vm9lBTcq3yL/+RDN71MwazOwjMzu1vRtkZscBZwCnOecuds79wzn3iXPueefcV4FfNal7mJm9kf/8FpvZTWZW2aR8BzN7Kn9f6/PfK6eueT+AD9zZidZmEekCSlBFpBi+DxxpZpOKcK5acklq01bUbwK3Ae0mDGa2G/A74ApgG2AicHeT8h3ItcCtAA4AdiHXRe3lyw8H7gDuAXYk1zJ3DnBJi0udCAwBJgMHmNkwYDqwGBgP7A3MBKaZ2ZAO3u/xwKD8674MHEWupZn8fXA07273gK8Bt7nS7FV9DfAn4J+sa4H9Sxt17wAOJtcqvgvwd+BvZrZti3pXkrvHOwH3k0sIx7QTw6nAh865+woVrmntNbOdgEfIfcY7A6cBR5D7nljjz8AyYF9yn/EPyH0/QK71OQt8j/VsbRaRLuKc00MPPfRo9QD+AGSAuhaPR1vUeTb/9Z+BfwNe/rkDvtKk7sdAssD5ftSizsXAnvmyamBbIAUMa3q9NmL+IrAKqGmj/B7g7TUxFih/Gbi/xbHvkhvGEMk/nwZ80PQcwM+Bf7V4nQGzge+1E++0/Hv2mxz7Zv4+Veaf/y8wvUn5wUAa2HQ9PsOR+c9hYst73KLebcC0Qp9rofJC5wJG5691WIs6bwJ35L/eIl/nB03KQ/nP+qx23scM4JH1eL/3AK+2OHY0EACj8s9XAae3c45Me+V66KFHaR5ru61ERAp4hVwrVFMNbdT9MfBf4HRyLWmF3Ajc1OLYspaVnHOvmtks4BRyCeqjzrnPCvRot/QM8BEwx8yeAZ4HHnLOLc2X7wY86ZwL2nj9DrRuHXwRiAFbA+/nj73R4hx7ALtZkwlkeXGgvZZByCVU2SbP/w5E8td7B7gF+I+Zbe+cm0GuZfkx59zCDs5batvn/32pxfGXgH1aHHtrzRfOuYyZfUbuD5C2GB20nuftQO4zb+rF/Ou3JzeG+RrgtvwwhWnkEt831+PcIlJCSlBFpD2NzrkP16eic+4Ty83mv9TM7m+j2vL1PR/we+BsYDNyXd/rE0Odme0O7AccCHwLuMrMJjvn3lhTraPTtHhuBY63nBjkAc8B3y5wvlUdBl74ermLOveemU0Hvm5mV5IbAnBMJ8/ZVNDyGkB4I87XkULJZarFc0f7Q85mkks+10dbn2+uWd+5X5nZvcAh5IZ5/MTMrnLOabkskTKiMagiUkxXkPt/5YIinOuP5Fof68i1jK4X51zWOfeSc+5n5FpMFwJfyhe/ARyYH8dZyHvA/i2OTSDXxf9RO5d9nVwCtcA592GLx5J2Xgewh5n5TZ7vQy6Bm93k2C3AV8l1/y8CnuzgnO1ZDAxvcWyXDl6TIjd5qD3v5f+d0OL4+CZlG+qPwGgzO7lQoZkNaBJDy89vf3LJ6Yw1B5xzHznnbnLOHQ/8jNwfQmusz3sVkS6mBFVE2hMxs00KPAr2tTvnaoGf0vayP1UFzjWojXOtBkYAO7bTJd+MmR1tZt/Pz9TfnFxL42asS06uIpf03mtmu5vZ1mZ2gpmt6YK+AjjOzH5sZmPN7ERy40uvdc61bPVr6v/IJTUPm9n4/Gz1cWZ2mZnt20HYg4AbzWy7/CStXwG/d82Xb3og/+9PgdvX93604VngJDM7yMy2ybd6j+rgNXOAbfMz4AebWbRlBefcbGAKcJOZHWxm25rZb4HPsX4rPLTJOfcAuclud1lu9YR9zGxzM9vfzO4kd1/IX2dXM7suf/1DgBuAe51zc82sysxuNLMDzGxLM9uFXEvqjCaXmwNMMrPhZjZ4Y+IWkQ2nBFVE2jOeXAtky0fBpDLvdmBWG2UXFDhXm0tHudx6nLWdiHcFcCS5FsYPyCWkl5IfE+uce5fczP4h+eu+BZxHbuY2zrnHyS9nBPyH3Az/m4BftHdR59xn5Fo+lwIPkeuSvpdc4tfRWNEHyK1cMJ3crP3HgR+1OH+C3ASgELn7uzF+DTxGbqzty+SGIEzp4DW3A6+RW1JsCbmxwYV8HXiKXIvn2+SGWhzhnPvvRsaMc+40ci3Ik/LXeJ/cmObPyCfAzrl3yA2B2D9//XvIvddv5U+TAQbk38/7+fN8xroWdsj9cbUbuUS1o9ZvEeki5pyWeRMRKXf5cb1x59yR3R2LiEhX0yQpEZEylh9fOZ7cElpf6OZwRERKQgmqiEh5+ze5IRVXOeemdXMsIiIloS5+ERERESkrmiQlIiIiImWlx3bxr1q1Sk2/IiIiIj1cv379Wi1dqBZUERERESkrSlBFREREpKwoQS2RWbPaWrdc1tA96pjuUcd0jzqme9Qx3aOO6R51TPdowylBFREREZGyogRVRERERMpKj53F3xbnHHV1dQRB0N2hNBOLxVi1alWXX8fzPKqqqjBrNSFOREREpEfodQlqXV0d0WiUSCTS3aE0E41GicViXX6dVCpFXV0d1dXVXX4tERERka7Q67r4gyAou+S0lCKRSNm1HouIiIh0Rq9LUEVERESkZ1OCWkTLly9n3LhxjBs3jrFjx7LddtutfZ5KpZrVPfvss5k6dWqrc7z88sucdNJJ7V7nnXfe4emnny5q7CIiIiLloteNQe1OAwcOZPr06QBcccUVVFVVce655wKQSCSKdp13332Xt956i4MOOqho5xQRkZ7B933CyQbMOTLRCtIa1iW9kFpQu9hdd93FpEmTOOCAAzj11FNpaGhYWzZt2jQOPfRQdtttN5588slWr62vr+ecc85h0qRJjB8/nscee4xUKsUVV1zBQw89xLhx43jooYdK+XZERKQbVWZTxF97gdBV5+Nf9l2iT/6FqkwSz9Ovc+ld1ILaxY488khOO+00EokE11xzDffccw9nnXUWAHPnzuWxxx5jzpw5HHnkkUycOLHZa6+99lomTJjAjTfeyMqVK5k8eTITJ07kwgsv5K233uLqq6/uhnckIiLdoTKTwn57Mbzz6rqDc2bC0w9RccWd1IW7fqUYkVJRgtrFZsyYwWWXXcaKFStoaGhg8uTJa8uOOeYYPM9j6623ZtSoUXzwwQfNXvv888/zxBNPcMMNNwCQTCaZP39+SeMXEZHu53keNvfD5snpGksWwhP3Ez76NNLOlT44kS6gBLWL/c///A/33nsvY8aM4cEHH1w7RhVotZh+y+fOOe6++27GjBnT7Pjrr7/edQGLiEjZCeGw56bSZvr58hOEDzmBdCReyrBEuowGrXSxuro6NtlkE9LpNFOmTGlWNnXqVIIgYM6cOXzyySetEtHJkydz66234vJ/Eb/99tsAVFVVUVtbW5o3ICIiZcBw7Y0z9fzShSJSAkpQu9hFF13E5MmTOfHEE1sloKNHj+bwww/n+OOP57rrrmu109T5559POp1mv/32Y5999uHyyy8HYMKECcycOVOTpERE+ogMwEHHtV1h0pGkYxWlCkeky6mLv4tceOGFa78+88wzSSQSzRLQm2++ueDrxo8fz/jx4wGIx+Ncf/31reoMGDCAF154ocgRi4hIuQqCgGD45nj7TIZ/Pte8cOSWcOAxpAONP5XeQwmqiIhID9DgR6j4xgV4B34RHr8P0ikYfwjsNo76UBQ0QUp6ESWoIiIiPUSDH8EbuxPhrbbDcKTDUbLZrJJT6XWUoIqIiPQgQRCQXDMpKpvt3mBEuogmSYmIiIhIWVGCKiIiIiJlRQmqiIiIiJQVjUHtAgMHDmT77bdf+/zee+9l2LBhBeuOGDGCBQsWlCo0ERERkbKnBLULxOPxZluaAiQSiW6KRkRERKRnURd/CdTV1XH88cczYcIE9t13Xx577LFWdRYtWsShhx7KuHHj2GefffjHP/4BwPPPP88XvvAFJkyYwGmnnUZdXV2pwxcREREpqT7fgvpJbZpL36xlYUOWTSt8Lt61mlHV4Y06Z2NjI+PGjQNg1KhR3HXXXdx5550MGTKEZcuWceCBB3LYYYdhZmtfM2XKFCZPnsx5551HNpuloaGBZcuWcfXVV/Pwww9TWVnJ9ddfz4033sgFF1ywUfGJiIiIlLM+naB+UpvmmKeWMad23Tpyry9J8fDBgzYqSW3ZxZ9Op7n88st55ZVX8DyPhQsXsnjx4mbjUnfddVe+/e1vk06nOfzww9lpp5148sknmTlzJgcffPDa8+yxxx4bHJeIiIiUD9/38bMZ8DzSDpw2XFirTyeol75Z2yw5BZhTm+XSN2v5/f4Di3ad+++/n2XLlvHiiy8SDofZcccdW41J3W+//Xj88cd5+umn+da3vsW5555L//79mTRpErfffnvRYhEREZHuZWZUpJPYe2/DS09ALE7kkBPIDh1Oox/p7vDKQp9OUBc2FN6BY1EbxzfU6tWrGTx4MOFwmJdeeol58+a1qjN37lyGDx/OaaedRn19PW+//TbnnXce559/Ph999BFbbbUVDQ0NfPrpp4wePbqo8YmIiEjpVKYTcOm5MOeDdQenPYZ/8HHETzpLSSp9PEHdtMIveHyTNo5vqBNPPJETTzyRiRMnsuOOOzJ27NhWdaZPn84NN9xAKBSiqqqKm2++mcGDB3PjjTdy5plnkkwmAbj44ouVoIqIiPRQYc/guanNk9M1nnoQ/8AvYsM26/Pd/X06Qb1412peX5Jq1s2/ZXVuotTGaLmu6aBBg3jssceIxWJt1v3Sl77El770pVbl+++/Py+88MJGxSMiIiLlIZxoyCWobXnqAUKn/5B0tri9uT1Nn05QR1WHefjgQVz6Zi2LGrJsUqRZ/CIiIiJtSqXaLLJUAqNvt55CH09QIZekFnNClIiIiEhbMpEYkb0nwhNTCpa7SUeSUX6qhfpFRERESiWFYV/8GvQb0KrMbbszbvMxBEHQDZGVlz7fgioiIiJSSg2xSiqu+RPukXvglWkQjcFBx2H7fYH6UBT6+AQpUIIqIiIiUlJZ56iNxImc8E1CR52Kw0jHK8hkAyWneUpQRURERLpBykEqEs89yapbvyklqEW0fPlyjjrqKAAWL16M7/sMGjQIgMcff7zgMlMiIiIi0pwS1CIaOHAg06dPB+CKK66gqqqKc889F4BEIkEmkyEU0i0XERERaU/JsiUz+z7wdcAB7wJfAzYF7gMGAm8CpzrnUmYWBe4GdgOWASc55z4uVazFdPbZZzNgwADeeustdtllF6qqqpolrvvssw/33Xcfo0aN4i9/+Qu33HILqVSK3XffnWuvvRbfL+6uViIiIiLlriTLTJnZCOA7wO7Ouc8BPnAy8GvgN865McAK4Mz8S84EVjjnRgO/ydfrmtgWf0r0d5cSu+J7RH93Kbb406JfY/bs2UyZMoXLLruszTozZ87koYce4qmnnmL69On4vs/9999f9FhERKRn8X2feDpJPJ0g5Gl1SOkbStnfHALiZpYGKoCFwAHAmv097wJ+DtwMHJ3/GuAB4P/MzFyRN6a1xZ8Sv/o8vCZJqT97Bo3nX4MbOrxo1zn66KM7bAl98cUXefvtt5k0aRKQGxIwePDgosUgIiI9i5lRkUlir/8DnnoQ0in8cQfDxMNpiMS1Vqb0aiVJUJ1zC8zsGmAu0Ag8DbwBrHTOZfLV5gMj8l+PAOblX5sxs1XAIGBpMeOKPHRHs+QUwFv8KZGH7iD5rYuLdp3Kysq1X4dCoWb/qSQSCQCcc5xyyilccsklRbuuiIj0XBXpJFx3Ibz3xrqDc2bC0w9Scfkd1IWi3RecSBcrSYJqZgPItYpuCawEpgCHFqi6poXU2ilrZdasWWu/jsViRKPr90MbWba44HG3fMnaxHFDZTIZ0uk02WyWVH7P3UQiwaabbsozzzxDIpHgnXfe4ZNPPiGZTLL33ntz+umnc8YZZzBkyBBWrFhBXV0dm222WaevvXr1ahYvLvzeyl3Tz1IK0z3qmO5Rx3SPOtad9ygejzOmdjF+0+R0jc8WEDz1AMv3/ALLVq4qfXBN6PuoY7pHhY0ZM6bd8lJ18R8IzHHOLQEws4eAfYH+ZhbKt6KOBNY0Z84HNgPmm1kI6Acsb+vkTd/kqlWr1ns5Jxs0tPDxgUM2ekmoUChEOBzG930ikQiQS56PO+44HnzwQQ488EB23XVXRo8eTTQaZezYsfz0pz/llFNOIQgCwuEw11xzTYcfYCE1NTUblNh2t1mzZm3Q++1LdI86pnvUMd2jjnX3PYqYw7/pjjbLvelPMfILxzJwSOHfY6XQ3feoJ9A92nClSlDnAnubWQW5Lv7JwOvAC8Dx5GbynwZMzdd/JP/8n/ny54s9/hQgdewZ+LNnNOvmD4YOJ3XsGRt97gsvvLDZ8zUtsvF4nL/+9a8FX3Psscdy7LHHbvS1RUSkpzNob0KU5+OsUGejSO9QqjGor5jZA+SWksoA/wZuBR4D7jOzS/PHbs+/5HbgHjP7kFzL6cldEtfQ4TSefw2Rh+7AVi7D9R9E6tgzijpBSkREpLMy5hE5+Hj4x7OFKxxwJOloBWiilPRSJZvF75y7BGg5A+gjYM8CdRPACSWJa+jwok6IEhER2VhBEBCM3Apvj/3htRebF26+Nex/BBklp9KLaVsjERGRMtQQilB59kXYwcfBE/dDOgX7HwY770t9OArFH/kmUjaUoIqIiJSpej+Ct83OhLfeAcORDkfJZrNKTvs4Lz8+uTevhdvrElTP80ilUmtnzvc1qVRq7TeuiIj0fEEQkPTym71ks90bjHSriMsSyaThwxngAmz0DqRCYZJer0vnel+CWlVVRV1dHY2Njd0dSjOrV6+mpqamy6/jeR5VVVVdfh0REREpnViQIfTaNLj115BJA+D8EOEvn4NNOoKEF+7eAIus1yWoZkZ1dXV3h9HK4sWLe+TapCIiItK9zIzQiqVw06XNC7IZuPu3hLbZCW/zMb2qy199wSIiIiJlLOwCmHp32xUeuoNIJlW6gEpACaqIiIhIGfOyGVj6WdsVli7GMpnSBVQCSlBFREREylg2HIHtdm67wjY7ku1lk8OVoIqIiEivZWZEgyxVqUaqEnVUpBP4vt/dYXVKOnDYQcdCrKJ1YTiCHXMaqV6W0vWudyMiIiKSZ2ZUphoI33E1fOtIOOsIvJ+dRfyDt4kGPatLvDFWiV15J4zZYd3BLcZil91OY2XXrxJUar1uFr+IiIgIQEU6AT8/B+Z/tO7gp5/AL79N5Fe3kt5i2x4z8z3joH7QpkR/fB1+Jg3OEYQjJCLxHvMeOkMtqCIiItLreJ6HzZ/TPDltwv3hN0RT5bVmekeCIKAxFKUuVkVdvJqGULRXJqegBFVERER6Ic/z4J1X264w+318bRlbtpSgioiISK/jnIOhm7Zdobo/Sk/LlxJUERER6XWy2SzsOg5CbWwBesQppKLx0gYl600JqoiIiPRKyUgUu/gGiMaaF+wxAXfQsaQDtaGWK83iFxERkV4pjYfbchtiNz4Ms/4Dq1fAdruQrayh0WujZVXKghJUERER6bUyGHXhGN7n9sTMCIIgNz5VypoSVBEREen1eutyTL2VxqCKiIiISFlRgioiIiIiZUUJqoiIiIiUFSWoIiIiIlJWlKCKiIiISFlRgioiIiIiZUUJqoiIiIiUFa2DKiIiItLDmBnRbJpQOgXpJETjJCMx0r1kDwIlqCIiItJj+b6P7xkEjoxzfWJBfjOjMtUIN18Gb04H56CiiugJZxLe/wga/J6/jasSVBEREelxkhamMp3AFs0n/OxDeKkEkYlH4DYfTUM42qu3M61IJ+DS78DHH6w72FAHd/0WLxojNP5wMj08UVeCKiIiIj1KwouwcnktVY/eSvTZBwFwAC8+DtvvSsV5V1LvR7o1xq5iZtiyxc2T06b+8nuie04kE46XNrAi0yQpERERWStsUJFqJJ5OEPLKL03wPI959Y6a5QvWJqdN2Yw3sX8+V5axF4PneTBnZtsVVi3HMtnSBdRF1IIqIiIi+GZUJBtwj/4JXp0G0Rj+Icfj9jqAhlCkbLrMk/i8vbCWnV+Y0nalpx4gvNckMuFY6QIrEeccbLpZ2xVicfD90gXURZSgioiICBWJetz5X4ZVK9YdvOUKmP40FT+4gvpQeXSZpwIjkcliyca2KyUTWJkk1MUWBAEM3wIGDIYVS1tXOPh4UpGe3b0P6uIXERHp8yI43IN3NE9O8+y9N7D5H+W6lstAhR+Q8aMs2+OgtivtOZFMtOcnaW1pjMaxX94CQzZtXrDvF3BHf5VU94RVVGpBFRER6ePCqQS88kKb5fbCo/hb70A5zAu3bIajt65iSf/PM2jUGPhkVvMK/QZgR32FVO9sQAUgGwTU9x9C/Io7sRVLYdVyGLEFmWicRC9YYgqUoIqIiPR5DrBQ24mNi0TBShdPRwZ4aRL9BrLwvP+l5sWHqXzpEUilcHsdgHfc12iIVuTWBu3FgiCgPhzDhm0GwzYrmzHCxaIEVUSkDzM/xOqsjwPiXkDEZXrdLzrpWDpWQXTyUfCXWwtXOOg4MtlyaD/NC7KMiDjqwv1IHn064YOPxcPhojGSjl6fnDbVW39ey2NAiYiIlJSZ0eDF+N37DYx7cB6f/9PHnPvyMj7LRnBez58BLJ2TDhwcdBxsvnXrwoOPIztgSNklQkEQUOFSRIIkqUiURCSfnEqvoBZUEZE+qJ4IJz+5kJcWNKw9dt8HtTzyUR1vfmlLNg0FZZeQSNeqD8eovOQmeO9NmP4URONw2Ilkh42gsYcveu95HtFUI34qCY0NUNOPdChC0vTHWLlSgioi0sdEo1E+XJ1plpyu0ZBxXPj3JdwycTAR1xvmAsv6cs5RF4ri7zoO//N7gXmksR7/h4rneVQ2rMZdcwHMei93MBQifOhJhL54OvW9ZFJRb6MufhGRPqayspL7Z9W2Wf7Ex3WknH499FXZbJYUHinXdeMbPc+jzqLMT0f4oDFErcXIeF2TKMaTjbhLzl6XnAJkMvDovdizfyVsZTT7S9ZSC6qISB8TBAH9om0noFURD0fPbjWTMub5fNjocezf5vHx6jQA8ZDx870Gc9q21cSCZNEuZWbY4gWwaF7hCo/8kegBR5HuhTtO9XT6E1lEpI+pra3l1G37tVl+5g79qPJ6/l7eUp5qXZj9H5i7NjkFaMw4Lvj7El5emMQv4jadnufBvI/arlC3OteaKmVHCaqISB+TzWbpHwq4cr8hrcp2Ghzlu58fiGX1S1uKz/d9Hvu4jtpU4SWrfvHKUlZni5egBkEAI7dsu0JlNYTUmVyO9KmIiPRBMZfm9G0qOXKram77z0qWJrKcPKaanYdEqXQJdfBLl/A8jzcWt92FP2tlisA8ivUN6JzDDRuJDRsBny1oXeGIL5GMxPvUuqk9hVpQRUT6qJhLMzKc4tI9arh53ADGDzEqg0SPn7Ut5SsIAvYc1vZ4z20GRPBccTcEaIzGsV/8DrbcZt1BPwSHnog7+HjS+n4vS2pBFRHpw5xzpNPpjiuKFEE2m+WQUZX0i3qsSrZORH+59xBq/CzZIg6BDoKA+sp+xC7+X/xkI66+DvoPyq2Dqk0pypZaUEVERKRkqi3Ny8ePYuyAdYv/V4U9frv/UPYeGiZbzOw0LwgCGkJRaiv7Uzd0JHWRuJLTMqcWVBERESmdIMsWMcfzXxxJfQYSWcegqEeFZfACbQ4hOUpQRUREpKSCIKCKJFU+4JObFKWhoNKEElQREZFezsyIZNP4mTTO80hG4rklmETKlBJUERGRXizsskRXLIM/3wyz34dBQ6k47gyC0TvQoH3opUxpkpSIiEgv5fs+sfkfwQ9OgVdegKWLYOY7cPn38J74C1FtyCBlSgmqiIhILxVN1ONuuhSCAjPjH7yDSKZ4+96LFJMSVBERkV7Ky6Th008KFwYB7sMZuf3qRcqMxqCKiIj0VtZB8hmJliYO2ShhHNF0ElYtxyJRgngViUiMbC/eBUsJqoiISC+VDYXxt/k8zHy7dWEkBluM1Wz+MhcP0vjTn4Z7/w8SjTjAhg6n4oJraBwynEwvzVFL1q5vZv3N7AEz+6+ZvW9m+5jZQDN7xsxm5f8dkK9rZva/Zvahmb1jZruWKk4REZHeIhmJYd++BKpqmhd4Hnz3lyRDakEtZ77v4384A26/GhKN6woWf4q76OvEko1tv7iHK+XAk98CTzrntgU+D7wP/Bh4zjk3Bngu/xzgUNt4KH8AACAASURBVGBM/vFN4OYSxikiItIrBEFAQ/9BcP39cMYPYa+JcOSXsRseIrXDbqTNujtEaUck0QD3/a5wYaIBXnmBUKh3doaX5F2ZWQ0wATgdwDmXAlJmdjQwMV/tLmAacAFwNHC3c84B/8q3vm7qnFtYinhFRER6i2zgqAvHCB1wDN7+R+D8EOku2O9eis9zASxoY5IbYLPfxyYdVcKISqdULahbAUuAO83s32Z2m5lVAsPWJJ35f4fm648A5jV5/fz8MREREdkAmWyWFKbktAcJzIMRW7RZ7kZv32t3iDVXghlgZrY78C9gP+fcK2b2W2A1cK5zrn+TeiuccwPM7DHgCufc9Pzx54AfOefeWFN31apVawOfNWtWl78HERGRYqquriYSiZBKpaitre3ucJrxfZ/Kykqcc9TX12siVTepqKhgdO1i/Cu+37owXknmuvt4d8Gi0gdWBGPGjFn7db9+/VqNNSnVwIX5wHzn3Cv55w+QG2/62ZquezPbFFjcpP5mTV4/Evi0rZM3fZPlatasWT0izu6ke9Qx3aOO6R51TPeoY115j0JmxFIN8Pp0+Oh9GL0D7LofjZEY2TJoDotl04QbanGvvQSeD3tMIBOrINFiW1R9H3WsGPfIBUPhmxfA3Tfkxp0CbDIS+9HVZCprGDOmugiRlp+SJKjOuUVmNs/MtnHOzQQmAzPyj9OAK/P/Ts2/5BHg22Z2H7AXsErjT0VEpKfzPY/4kk9x1/4Yli+Bxnp4+iGorKLisjuoH7RJt7ZWVmTTePffintyyrqDd15L6MRvED/kRBpbJKnS9RJemPB+hxDd8wBYvQKLRAhiFTRGK8j24pbtUk79Ohe418wiwEfA18iNgb3fzM4E5gIn5Os+DhwGfAg05OuKiEiZ8X0f55y6gNdTMmssC/dn9lnXM7QixOCVnzLk7l/D3Nm4q84n9stbaAjHuiU2z/PwPpgBTZPTNe7/Pf5u47BNt6AUQwOluTRGOhyFQZusO9jLf+ZKlqA6594Cdi9QNLlAXQec0+VBiYjIBolnUvj1tTDjzdwam9vuRDIUJd3RzkWd5HkeDc4n7Yy45wi7TI9NkBq9GD+YvoT7Zq5eO7FlmwED+NsPbmDL334P5nyAn0xANyWokXQSpt7TdoVH7iX8zQtJoaWppOv1zsWzRESky1Rm09htv4Z/PrfuYDhC9EdXwegdSXt+Ua6TtAjvr8xy6WtLmFebZvehMS7aczBDw1lCLlOUa5RK4IW56s0V/Hnm6mbHZ65IMfk5x7/O+ClDf3oqLp3spgjBshmoXdl2hVXL8bIZUDe/lEApF+oXEZEeLuR52EtPNE9OAdIpuPI8YqlEUa6TthB3f1DPxAfn8uzcemauSHHvzNXsdO9HvLfa4Xnl+evL933C4TC+3zxJr3c+t7y7ouBr5tamme31gxFbQlW/UoRZUCYag532bLvCrvuSCWvnKSmN8vwJFxGRshRNNsDf7i1cmM3g/vFMUXa2aSTET/6xuNXxTABff3YRta68WvHM86izGA/PTfOTV1fzwCdpai2WmwUPpAJHQzubps+uzcLJZ5HsxgQw7Qw74stQUdW6sN9AmHAYmSKOewyFQoRDoVbJvAioi19ERDrBAFYsa7vCkoXYRm6f6YXCvLM4RbqNXGjWyhQNWagskyYWz/OYnwozfsrHLEusWwS/X9TjxeNGsVXcEXHGwJjP8kThRfI/N3IA2ao9SHfz+M6GWCUVV92Du/NaePPvYF5ue9TTvktjtKIoE3NCuFxL+7RnsI9n4bbZCXYfT2M0TjbomeOLpfiUoIqIyHrLmoe/zY4w49+FK+y8D9mN2Kmo0YvxwKw6Bsbab1XzymieTq0Lc/xjC5olpwCrkgFf/Nt8Xj5+M6otwwW7DeSCvy9p9fqdBkcZ3i9GY1Cc4REbIwvU1Qwk+u2fE8qmASPth0j5YVwRklPfIPbpXPjZWZBK5CaLvfAoVNVQcdnt1A8cphUhBFAXv4iIdEIyEoPTvgeFWkk33Ry23m6DE4yUhbn0teWcO+0ztqiJEA8VzkJ3HhIjVka/vRqyMGN54clNc1anWZlyEGQ4bdtqfrbnICqavK+DNq/kb0eNpNJ13+SolpxzJLwQdeE4deEYSS9UtJUTYqkEXPkDaDlWuW417roLiaUai3Id6fnUgioiIustCAJSQ0cS/cUtuN9fCfM+Aj8Eex8Ap3+fhkh8g7uBG/H5XX4i0c3vrOC6CcP4n+cXNdtrvCrscecXNqHa0pRLb3BjO2NLAerSAYQhFiT53o7VfP1z/alPO+IhI2pZKkj22KWzOstbtRy3so0hIp98iJdKQkgTsUQJqoiIdFLK88lutR3RS27OLTvkeaTDEVJeaIO7gc2MebWZtVt9/vG/qxha4fP88aO4e8Yq5tWl2WfTOF/foR/9LE0QBIQ9w8+kCfwQaazbkrwBUY+qsJdLRFuI+saIyhDkl8XygjQ1pKlZ89vXQd9ITXOf8dqtOtuqk+1Zy4dJ1ymjThIREekpstksDeEodbFK6iJxkuZvVILonGNAi3Gn1725nOP+Np+wD3ttEufsHfvTzyUIBVmqkw1Ep95N6LoLidxzPVUrlxANNnzs68aosAy/3GdwwbKL9hhEHCVdkPuMbchwCLWxAkN1fywaL21QUraUoIqISFnoF4ZtB0SaHVueyHLruyt5bl49BFk8zyO+bCHuO8fBA7fDu6/BMw/jvnsC4XdeIexKP8HGDzJ8eUwlUw4dzvYDo3gGYwdE+ONBm3LWDjU9blOBrmJmuMZ67OhTC5effFZRJmJJ76AEVUREykKNpZl65EhGVDUffbZVvzB/Ong4laSJpRpx/3sJJFpMpnEObvwl0UyqhBHn+GYMWjCLL/7pQp7begHzDwnx4phFnDzlYgZ//C6hjVx2q7dwzuWGM2wyEjvnZ7DlWIhXwtgdsfN/DY313R2ilBGNQRUR6UJmRtpCNASGD9T42Y1ahqk3C4KATcNpXjlxFDNWpJixLMnnh8QY2z9MNSmCwOGlU/DxB4VPkErAvI+wrXco6XjUWCoB1/0YFi9kyNv/albm5swkes29ZMKxksVTzlxlNTx8N1T1ww49CRswCLd4Ie6B26Gymsyko7o7RCkTSlBFRLpI4IVYmva59LWlvLigkaFxn/N2HcT+w6PEgvJZVqicBEFAFQn2Huixz6AYQRDgggRrO347SjyDbK4ruYQJqpdowC1eWLhw5TKsdiUM3KRk8ZSzxmgFlRddj/vZt3A3/WrdBLGRW2Hf+VVuJ60+sqKBtK9TCaqZfQE4GRjqnDvSzHYHapxzz3dJdCIiPZTnecxpNPb5yxwS+anpn6xOc9ITC/jGDv351V4DiLnSd0f3FG2tpRqEI3jDN4dP57YuDIVh89ElX+jddTQ5K6MxqGsEQUB99UDiV96FLZwLn34CW4zBDd6EhkhcY1BlrfUeg2pm5wI3A7OACfnDjcClXRCXiEiPVu9CfHvaZ2uT06Z+/95KVmU0LnFDJKMV2DmXQKhA+8pXv0MqFGl9vKtVVEF1/8Jl8UoYUHiGf18VBAH14Sj1o8bSuN/B1G+6BfWhqHaQkmY6M0nqe8CBzrkrYW1vy3+BbYoelYhID5dyHn//tO1dcZ6eW4/vt7+dp7SWDQISI7bErp8Ck46AkVvBrvvB5XeQGX8IKa/09zQRiWHf+knhwjPPz3VbSyvOObLZbJ/ZpEA6pzNd/NXAvPzXa76bwoD6qEREWjDL7TteoAEVgJivFtQNlTaPTM1Awqf/AD+VwoVCpEKRLm+BMzNSFiIR5Np2avwsLpsh4yC1/W5Erv4j7r7f5XbXGj4KTvommU02I40+a5HO6kyC+hLwY+CyJse+A7xQ1IhERHqBKFmO2LKKqR/VtSrzDL6weSXZbKLAK2V9OOdI4UEkPzu+i5NTZz5LsiEu+ddSHv+4jqqwxzc/15+v79CPSpck6fmkhm1G9JxL8LIZAj9EKqxua5EN1ZkE9VzgUTP7BlBtZjOB1cCRXRKZiEgPFiPDbyYM45VFCRY1NJ8k89sJw6iwbN/Z47KAUCiEmREEQdkvu2VmLMmG2OVPH6/dznRlMuBn/1rKo3PqePjwEVSSwDlHwg+Dn98pScmpyAZbrwTVzDxgO2A8sCMwilx3/6vOdcO2HSIiZc45x2A/xeunbMHU2bU8Obee4ZUhztlpAMOijohLd3eI3SIExJMNuBeeweZ+iNtuF9hlXxojcbJlOhYxZbmW0zXJaVOvfZbgP8tT7DPIU2upSBGtV4LqnAvMbKpzrhp4Nf8QEZF25Nb0bOSro6OcvHWckAFBus9OCgmZEZ8/G/fzsyGVzDUgP/8oVPen4vI7qO8/uCyTvETg8fjHrYdqrPHnD1az37gBJY/dzIhm04TSKSyVxEVjpCIxUhrzKp0U9gwvnSIIR0gH5fH/U6fGoJrZ3s65f3VcVURE1shms/hk+3KPPkBum9IrfwipFpsU1K7EXX8R0Z9cT2OoPGe8V4U9ViYLJ6D9o17JU0IzozKdgN//Gl6dlls/NBIjcvjJhI/8MvVrhhmItCPsHNHGWnj8Pvh4FmwxhthhJ5OIV5G2ziz0VHydSVA/AZ4ws6nkuvfX/l/rnPtZsQMTEZGuEQ6HabAISedRnw6oiXhUWJZQ0HXDDswMli2G1SsKV5j9Pn46BWWYoNb4Wb75uf787F9LC5afuX1/MpnSDtmIp5O4ay/EZry57mAqAX/9A+Z7RI48VS2p0q6QGbFPZuJ+ec66zSTefQ33xP3EfnojbsttyXRjb09n0uM48DC5xHQksFn+MbIL4hIR6TPMjLQXpo4IaS+cS+a68Fo1m43mS09/xpZ3zuZzf5zD6Ls+4vI3V9Pod/F+8YmG9mPLlueOSy6b4euf688ew1rfn0v2GsyQKCUftuE31DZPTpv625+JpLRChLQvmmrE/eai1judZTK5Ho1U2+s4l8J6t6A6577WlYGIiPRFKQszvxEufW0ZM1ek2GFghIv2HMzwWEA4KH7CVkeEYx5fwDtL13Wzp7KO6/69nKqw8f0dq7AuuK5zLrc2qB+CQolov4EE0XjRr1sslUGCqUeM4N1lKf78wWoGRH3O2L4fQ6KUfMtaz/NwhbZ6XSPRmBtGESnf+yndz+prYUXhXgFWLM2V9+/iP1rb0ZkufsxsDHAKMAJYAPzZOTerKwITEentAi/EC5+mOfmJBWvHTL23LMmUWbX89YgRTBzmQ0f7vHeCmbGwMWiWnDZ1/Vsr+Mbn+lNN17RkpkMRwseeDlNua114xnkkoxVluzSTc44Kl2CfQR77jRuAAZlM10548zwPM8M512wClnMOBm/S9gtDIYh0w5av0rN09L3bzT+L693Fb2ZHAm8A2wLLyW1x+rqZHdVFsYmI9Gr1LsRZzy9sNXnKAd94bhF1rrgTXcyMmcvbbu1bnQpItrX1VREkPZ/gsJPgh1fCqNEQi8M2O2G/upXMTnuSKdPktKkgCMik06TTXZecep5HZSZJxZz3iT/3MBWz3qUqncT3cr+ynXO4fgNhxBaFTzD+ENLaXlU6UlUD1f0Kl1X3a7usRDrTgno5cLRzbu3OUWY2Efg/4JEixyUi0ustbsy2OTN8SWOWlamAqiLmqM45tu7f9gkrQkbUty7dQKDBC+Pvsh/R7XbGcATmk4zEynJ5qe7geR6V9atwPzsLFn+6rmDAECouvZX6mkEEQUBjNE7lT2/AXfbd3Naqa+w2Dk79DknzSx+89CiJSIzYty6Cay5o3ppqhp31ExrD3de9D51LUEcCL7c4Nh1NkhIR2SBeB5OhOirvLOccIyt9RveP8OHK1i2pZ+04oCQ7XGWzWRqa/vJTcrpWNNWIu/bHzZNTgBVLcFf8gNgvfkdDKLeFan1lP2I/vxmvdlVuLOGwEWSjcRp9de9LxzIO0tvvSuSae3FTboP5c2DkltgJXyc1cGgXDfRZf51JUN8Cfgj8usmxH+SPi4hIJw2OeQyJ+yxpbD3OdHhliJowRU8WayzN40dsyvFPLlo7FtWAL29bw492G4Cf1ezv7hRKp3Cz3itcOH8OfkMd1OS674MgoCEUhQFDsYHD+uwGELLhkl6I1JARRM66CC+TIghFSPmhsvhe6kyCejbwqJl9l9w6qJsB9YDGoIqIbIBqS3PPwcM5fOo8mg79DHvwx4OHU22ZDucxdFYQBCQWfswTR41hdRqWJbKMrAoRt4CoktNu55LtL+3j6muxfoNbJRDlkFBIz+ScI+n561Z9KJPvpc4sM/VfM9sO2BsYDnwKvOJcH91QWkRkYwVZdh8YYsapW/Hbt1bw3vIkOw2K8p2dBzAwlMV1wXJPAMlkkkqXpDIEw6sN5zJd3q0v66myBiKx3KL7LfkhGKSWUukb1jtBNbOdgWXOuelNjm1mZgOdc293SXQiIr1cyGUY5me4bI8aEs6IeQ7LJqFEwzKV7JSXVDhC5Jivwv23ti485HjSYY0vlb6hMztJ/RFoOf0zAtxTvHBERPomCzLEXbpsd1OS0kjhERxyAnz9RzBgcO5gzQA49Tu4487U7HzpMzozBnVz59xHTQ8452ab2RZFjUhERKQPa/DDhPY/guheB2Aui/N8UtE46UCt3dJ3dCZBnW9muzrn1m7+a2a7khuLKiIiIkWSCQIykaZLcSk5lb6lMwnqb4CpZnYVMBvYGjgPuKwrAhMREekpPM8jFGQwB5lQmGy2eFvUivRFnZnF/3szWwmcSW6JqXnAD51zD3RVcCIiIuWuIpvCmz0LnpoCmQzhSUfittuFhnBUk9BENlBnWlBxzk0BpnRRLCIiIj1KRTaFd/Nl8NqL6w6+MR3G7EDlhb+hbj12dfJ9H8/zcM6RyWiSnAh0Yha/mZ2SXwcVMxtrZi+a2fNmtm3XhSciIlKePM/D+/iD5slpns16D/75HCGv7V+znudR58V4bEGGb09fwXXv1rPMRUlbp9qORHqlzvwUXArsm//6WuA1oA64CTigyHGJiIh0m6qqKuqJEGBEvICYy7Tqrg8FGXiynVFuzz5MeK8DyISjrYrMjOVBhP0fnMsnq9ftd/OLV5Yy5bARTBgWIuRat6Z6nkc4mwEcmVBEY12l1+pMgjrEOfeZmcWAccDxQBpY2iWRiYiIdIOEF+VtN5SrHlvI/Lo0uw+Nc9m+g9ks7nJJaZ45oL11a7OZNreNTBDiey991iw5Bcg6OOXJT/nwtK2opsm1zKjIJLH//gee+Ss4R/jAY3DbfJ7PwtWkso6wBzVeFqe1dKUX6EyCusTMRgM7Aq8555JmVgFY14QmIn2N53l4nkcQBARBibZS6iPW3FtN2mlf0sJc/e+VXPvm8rXHHv+4jic/qeOZL27OngP9ta2WmVCY8KQjcmNOCxl/MJlYvOASUUl8HvmornAMWcc/FzZyyHB/7c9BRSaJ/eYiePe1dRXfmA7b78qyMy/l84+uYFiFz8/2GsxxW1USC5IbeAdEykNndpL6FfAGcDtwdf7YZEDbnIrIRvE8j3ovxhsrjRtmNPDEp1nqLAaeds3ZWGvu7T+XwTXv1vPI/AyDt9oOp3tbUH3gc12T5HSNwMHZLyxiVbCuXSebzeK22xU3evvWJxo2Ag44us3F9bNB+0ubrkoFmOXaf3zfx2a82Tw5zbMZbzJ81muMH1HBZw1ZznnhM+79oA7naRyr9GydWWbqD2Z2f/7rhvzhV4CT19Qxs/2cc38vbogi0puZGStchIP+Oo+ZK1Jrj9dEPJ4/dnNGVwKBxtltiDXjHCc/NJfZq9Z1JVfn7+0Y3dtmPM/j9YWNtJU3frAiRTKAqib9hg3hKJU/uR7++Rw8+1fIZGD8ITD5aBoicWijJyDsOXYYFOW9ZYVbOieMqCCbzZWFUwl46sE24x7w0sN85YideHlB7vmvXl3KCWO2pIbWXf1hM8LJBgxIR6KkzVerupSlzi4z1dDi+eIWVZ4AajY2KBHpOxoJ8z8vfNYsOQVYnQo4dOo8/n3KFlShJGpDNBLmm88tapacAtSmAg6bOo83T9mCat3bZipC7Xcs+tZ8VJtzjjo/QmjC4YT3mgQOMrH8tqTtDFPp52W4ceIwDnhobquW1FO37UdNyLE2U3YO2psMlcngNwlrZTKgMeOoadJIbmZUppPwxP3w/FRIp4nsuT/RE79JQ7yKrJJUKTOd6eJfHxqPKiKdknQej88pPBZvSWOW2avTa7s6pXMSzuOZufUFy5Y0ZvlY97aZIAj4/OAo8VDhezJhRAVRK5x0ZoKAxnCMxkiszW79prLZLDv0M149eQsO3KySqrDH6P4Rbp40jKvHDSbu1v3BlonGYeLhbZ5r9b6H8dCidc89g5jf/D1UphPwi7Phgdtg+RKoXQnPTcWd92UqEoW/R0S6U7ETVP0JJiKdkgxcu/9xLG3MKonaQKkO7u2yhO5tS3FL86dDhuO1uC1D4j6/n7wJlVa8GfJhl2FsRZY/fmEo75+6BdOOHcmXt44Syyaa1ctks7D7eBg1pvVJRm7Fks/vzxOfrOvgPHLLKmK2rsXV932Y8SbMnd369bUrcVPvpuPtBERKS6OoRaRbVfgwtMJncUPhLsydBkcJNCN5g8R9a/fe7qh720ooyDJ+WIgZXx7FH/5bx+xVaSaNjHPkVtXUkCz66hJBEBAjRWzN8zZ68utDUSovuRH+/jQ8/wg4cBOPYN4uB3LQc3Vr/xD5/JAo/zdxGDGSa4/52QxMf6rtIF57idAxp5GKxIv0rkQ2nhJUEelW1V6WK/YdwpnPLmpVdtLYaqp9p76ZDVTjZbhy36Gc8ezCVmWnbdePSk/3tpBQkKF+7mwu2mUbMs4ImyOTaXvyVCmsHes66SjC+xyYG+saryQShLltchVzVqfZeUiMERUeVSSbT3zyPIhXtn3yeAVOLelSZjQGVUS6lctmOHzzOH85dDhb9wsDMCDq8fO9BnP9hKHEXKqDM0hbXDbD4aNiPHDYCMYOyHXiDon7XLbPYK7cb7DubTuCICCbyWDZNJlM+Sx8n8k2GeuazVLlEuw1EE7ZMsI2FRkqXbLVrPy0Aw49qe2THnwC6WhF1wYu0kmdbkE1sxqgqukx59yn+X+rixSXiPQhMZfisBEh9jtuM7IYHlBtGWgxFk86LxqkOHi4z97HjCQww3A0LF5APNt6+03pmToaduCcIxi+Od45l8D9t8KSdS3q7nO7w94HkNHGGFJm1jtBNbMDgVuBUTRvKXWAVnwWkY2SzWabLyelrueiWXtv8/d00coVDBsyuHuDkpJIWZjVWY8pM2tZEdubYy6ZzFaJxQx+fgrsMxmGj6IhFGlzS1aR7tKZFtTbye0mdR/Q2DXhiIiISDEkLcyUjxr59rTP1h676o3lTBxRwX1f+gEVLj/pS8mplKHOjEGNAXc65+qcc9mmj64KTkRERDbMyozXLDldY9qCBu6YsQpnxZ6GIlI8nfnu/A3wI9OieSIiImUtFApx1/ur2iy/6Z0VrM5qdJ6Ur8508T8IPAVcaGZLmxY457YqalQiIiKywcyMxY1td3CuSGZzs0m6qHc/5PtEEvV4QNaMZCRe9DVkpXfrTIL6APAyMAWNQRURkU7I+mEanQ8OQgTESbdaDkmKJ5vN8sWtq7jtPysLlh+4WSURuiZhrAjSeK9PgwfvgMUL8bccS+VXziW9+WgSnpZfl/XTme+ULYFdnHMb/B1tZj7wOrDAOXeEmW1JbtLVQOBN4FTnXMrMosDdwG7AMuAk59zHG3pdERHpHmbGaoty+avLuWvGShozjgM3r+Q3E4YyPBLgufJZY7Q3CYKAnQfH2HFwlHeXNt8tLOobl+87hBjpojegRoMM3sN3w9R71h2c9R7ukm8R+uEVhHbel0ygP0ykY50ZgzoVOGAjr/dd4P0mz38N/MY5NwZYAZyZP34msMI5N5rc2Ndfb+R1RUSkG9RZlAMfmset765gx8Exxo+o4NXPGtn9zx+zOOOjaQ1dp8olefLokfxgl4HURDw8g4M2r+S1k7dg02i2S1qwI5kUPHpv4cI7riWaVAesrJ/OtKBGgf9n777j26rOBo7/nnuvtmeGsweEJBA2JIS9N5TRsqGslrdldZfZAaWD8pbRTd9SKKMtpZRdKHuUAmGFFQI4eydO4hEPjXvvef+Q7Nix5BHbsmw/38/Hn0jnXMlHJ8fSozMfE5H/AG2WBRpjzu3swSIyHjgO+Anwrcxiq0OBszKX3A1cB/weODFzG9JTC34jImJ0PEgppQYM27Z5bVWc48bYfOOgEsIfv4MVbyB52F68VBfkhjnruXX/YQRNqr+LOigZY4iZON/bo5hv7V4OQAA/3XPaB72YIoJZthByzTWtXo8k4hAI9/rvVoNPdwLUeZmfrXUbcAXQfNrUcKDGmJbxnRXAuMztccByAGOMKyK1mevbLM5SSg1tIoKI6OKLQiUWxDfxPfMBJd/+MXibh/NP2v8YtjvxcuK+RbAfizgUWL5LEZvrvk97eqJFHec7OgdVdU2XW4ox5vqt/SUicjywzhjzjogc3Jyc7dd0Ia+dysrKrS1aXg2UcvYnraPOaR1BJFZEdOQ4Kus8VjW47D4yTMzE2bRmOZCuo3A4jGVZJBIJPE+3a95SPtrRsGHDOaJEKPrh9e02gw+++hRTd9iT+mGHU7m4MNu0/q11bss62nHsKILFZbApy+Ks6btSG0+waPnq9nmDmLaj7KZOndphfre+yojIEcAZQIUx5nMiMhMoMca80MlD9wNOEJFjSW/4X0K6R7VMRJxML+p4YFXm+hXABGCFiDhAKbAx15N39iILQWVl5YAoZ3/SOuqc1hFg2Sxqsjj4n8tZ32obnSMnxvjzETOoXb2c4tETeG99gnWNLrMnRCgPQtQk+7HQhSVf7SjgOAQfvCPnSUXFT/2For0PoLgA27T+rXUuWx35AnLNrZgfXgzJa1sUmAAAIABJREFU+OaMsuHI139EoLicqWXD81zS/qPtaOt1OUAVkctJL3K6Azglk9wE/ArYt6PHGmOuBq7OPM/BwHeMMWeLyD8yz3U/cB7phVgAj2Xuv57Jf0HnnyqlAOpNgMMfWsLGeNte0WeWNfDTtzfylZ3GssN9i9mU3Dzsf9w2Me44bDQRL77l0/WIsR3qPQsLKLG9PumpDVgWgXgDAG4oQrJA3gmdzFCt6+ZehS/GYG1Yk3v4q3Zj7vmKqk+ICI5tgQHPmF6fHuMaaBq7DeHfPgxzXoTlC2HGnrDTnjTqXqiqG7rTg/oN4DBjzBIRuTKT9gkwvQe//0rgfhH5MTAX+FMm/U/AvSKygHTP6Rk9+B1KqUHCsiw+WJ9oF5w2u3NeDZ/frqRNcArwr8UN/Pr9Gq7YtRi8ni/IEcui1gT51dxqHltUT1HA4tJdyjl2crTXgmARIZaKw7MPw4uPg+cS3OdwgiecQ2M41m8f9E0SZGNKePiTOiwRPr9dMWWOIZylh9oDzB77wyv/zv5k2++G5wT6tsCqRcRNYldXwbMPIYk45uDjMROm0BgI9eqKfheoD0ZwDjmhZY6453n6ZUR1S3cC1GIyC5fYPB80AHRr3MwY8xLwUub2ImCvLNfEgVO787xKqcFPRFhalzvAbHQNVo5di373QTUX71zaskpza9m2zXovyF73L6aq1RSDLz23mqMnxbjzsFFETaKDZ+iamJuA6y+BZQs3Jz7+F3j1GaL/ey/1/bASuskKcf2bG7n9w83zC695rYpv7T6MK/YoI+y3fd2e58HOs2DkGKjaYt6h7SDnXEbCDuScAqB6T8RLYj94Bzz1AJD5EH/pXzBjD6LfuZEGu/eXqnXUu65UZ7qzD+orwFVbpH0NeLH3iqOUUrn5vs/MUbkDs/FFDjWJ7L00NQkfk3X9ZdeICA1WmLeq4fo5VW2C02b/XtrAgjoXy+rOW2t7tm3D+3PaBqfNqqvg3w8SyPP+oZZlMXd9qk1w2uyWuRv5tDb7624MRpCf3AH7HQF2pk9k6o7Iz+4iXjpcT5PKAxHBrlrTEpy2yfv4XeT153F62GaV6m3daZGXAyeLyBKgWEQ+Jd3L+a2+KJhSSm3JGMO4mMWsHEHqj/YZyR0fVWfN22FYEKsHG+zUS5jDH15OTcLnsUX1Oa+755O6dIDZA46bhFefzn3BG8/jJJvyusl9o7G5+d2ca1W55d2NJGn/un3fpz5cROrLV8EfnoD/exLvqltpGDOJlGhQlA8BAZ76e+4Lnn6QQKIxb+VRqiu6/O5gjFkNzAJOJ725/nnAbGPMmj4qm1JKtRMzSR45fjxnTy8hkHkHGx11+NPhozlqYozFOaYA3LR/BaXW1g05iu1w7/xa5m9M4hlafm82YbvnQaMRC4Kh7Jk77Ib1latxaquJrVhIUaqJoOn7bbQ8Iznn/gKsj3ukcsT/xhgSlk19IEx9KEqTE9TFMvnk+9DUQQAab0K0J1sVmG5tM5VZST9HRN5qThMRyxij7zRKqbwwxlBEnNv2H8bP9htJ0jeELKHUcjF+nCeOG801c6r5R2Udrg+TSwL8Yv8KZo1w8LZygdQmz+LuT2oBeHJxPadOLeE372fvqb1wRmmP5965lk3w2NPh9efbZuy5P9ZxZ+Df9j3YsC6d5gQInvIlnCO/QKPddwuOorbPUZOizK3KvgjsuMlFRC1DHmJl1U2e7cABR6VX1Wez18G4wUh+C6VUJ7rcgyoie4jI6yLSAKQyP27mX6WUyhtjDAE/RYmJM0ISFJs4vudijKF66Wf8av9hLDp/CgvP35ZXT5nAUWNtwj04TtMAzSdD3v9ZLWdvX8q08vaLSi7euYwxEenxvErf9/HHTYaDj9ucaFlYp12Ef9N3NwenAG4K7r8da97bPZ5a0BHxXC7epZxh4fa/oyJqc+4OpRhPF8UUIs/zYIc9YFKW/ThLypETz+neamel8qA7Pah3A48DFwI6WUUp1accyyKYaEQAz7JJBkJdGhb2PI+AnyJAKn0mnYGebk9aYvucPb2EH7yxnibXcNZTK/nzkWN5tyrOM0sbKA5YXL5bOduVWET83vmob7SDRM/9OtbRp8HTDyJjJ2Iq50G8KfsDHvgjoR12p7EPV/eXkWTO6ZO44tUqHlu0CRH4wnbF/HTfCkpIoENphashECL2g9/Cs/+E5x+DZAL2PhT5wgU0hmK6k4IqON0JUCcB1+qG+UqpvhbzU8hLT8ETf4PaapwZuxM47+vEyyv6ZWGN8Vy+tGMZ93xSx4KaJIvrUhz04FIOGh/lzGklnDyliLAXx/i924PYaAeRsZNxvvRdHNvBvusXuS9es6IHexR0jTE+FXaSPxw8gl8fVAECQTxCJoHv60dDITPGUO8ECRx3FoHDTkKAVDBMEtHgVBWk7gSoDwNHAh0sLVVKqZ6Juknk9p/AWy9vTnzvdcyHbxH66Z14Yyf3ywKbYhK8duokVja4xD1DyjcYA1NLnXRw2kcf8sYYUp7BxyMyfRd45qHsF47fBr/PQ9R0eYImSRBadsTW8GbgSPmGlM43VQNAdwLUMPCwiLwKtFm5b4w5t1dLpZQakkQEq666bXDazHPhjpsIXXULTU7vbyremUYC/HNBHT9+cwOrG1yGhW2+sXs500pL87KXZyCVwJqxO35xKWyqbX/B2ZeSDEd7Pp9BKaUKQHcC1I8zP0qpQcq2bWp9h4QPlggh8SkSN289lpZlwXuv576g8iNs3wXyG6D6lsOd8zdxzWtVLWkb4x4/eH09Kza53LBXGaEeLMLqjG3bOB+/jf/YfVhX3oz/xxth6YJ0ZqwYOfdrpKbskF4Mo5RSg0CXA1RjzPUicgRwJlBhjDleRGYCJX1WOqVU3rji8EGN4X+eX8En1emFPnuPiXDn4WMYF0xh8jWsXtTBW0ogCP0wB7XBOPzkrfVZ8+6YV8M1s4YT6sMNTULxBrjnV7B6Gf7//Qz5/PnI6Ingpv+fTOkwEnZQ5xIqpQaNLgeoInI58HXgDuALmeQm4FfAvr1fNKVUvogIa5MWhz60GLdVHPrG6ib2fWAJH5y9LSXkWD3eizzPQ3bbG2PZ4GfpDTzgKJKthvdFhJCbJJBKYGo3IqXDmDxiGLZt42Mhkl7g1NMh+PqUoSHHLvS+gcV1KfYo7fn2UrlYvg+rl6XvLFuIue376XmfklngcsG3sA47WXtQc7AsiwbjkDIWCITwCNPzdqGU6jvd6Yr4BnC4MeZGaNlN5BNgeq+XSimVV0lxuG7O+jbBabOahM/9n9X26R6brSWcEHz9hnTw1drYSXDWJaSsdDlEhFgqgfPbH2G+cjxccS7mii/ipQyvVvl89T8b+eZrNXza6NAkPZsSEHY6Xnw0LGz3bbBjWRCJtU9v/p0jRmuwlYOxbJbEHU5/ei0T71rI5LsW8j8vrafKC6anlCilClJ35qAWA8szt5vfCQOg+/sqNdDFfeGttdlPCAJ4eWUTF06P4dD3PXRJy4ZdZhP8/aPw8pNQtRZmHQjbzaAxGGmZahByk3D7T+Cd/6QfaNmsv/YPnPQ2vLF2Rcvz/WleDZfsUsb3Z5YTMVv3dhUSn/3HRnl1VfstoKeUBhgekj5dyp4Ihggd9QV45J72meEoTN9Fjw7NYaMXYO+/L6LRTf8H+QYeXVTP62uaePfMbSjKw8iAUqr7uvP18RXgqi3SvgbkODtNKTVQOALji3Ifk7lNSYBA3+9g1CJp2dRHSogffw6JC79D44w9qXfabtQfcJPw9ist9/1ZB/HXqhBvrE20e77ffVDDsgaDbNkr20UxUtx71BimlLato1FRm8dPmECJ9O2BeikjcMI5sPsWs6liRXD974n34eb8A5lnBfj52xtagtPW1jV6PPBZHVaeRgaUUt3TnR7Uy4HHReQioFhEPgXqgM/1ScmUUnkTE5erZw3juEfb9xAKcOku5RivfeDX1zo6095sqmuzKGjDfsfzhwW5e0h//2E1v9y3HLzuB5PGGIZbSV7+wkQW1KaYWxVnenmQXYaHKJFkXnov6+0g0cuvw6qvg8qPoGwETNqOeDBMlvhrwBARPMsh6UPEMr16XGqjL7y4IvfBh/9e1sA506IE8zAyoJTqnu6s4l8tIrOAWaRPlVoOvGmM0XElpQY43/fZc0SIH+w1gh+/tb7l3PmwLdx15BhGBPyC241dikswsvkUHN8J0phtEm3GpqSPD1u9lb3v+xQRZ/dSYc/yMMYYfD9OPg9QarSDUDoCe/Zhmd9feP8v3ZGQAKvjwq1zN7J0U4oDx0X50o6llEoq+yK5bnIEKqI2lTXZ80dHHDqZXqyU6ifd6UElc8zpm5kfpVQPeFaABmPjG3DEUGq5/boKO+wnuGynIr60YylvrY0TdoRdR4SI4OL04R6fWyvlBHFmHwJvvABA2bzXOWnSqfz64+w9cOdsX4JtvB6fF2+M6dX/p+aFOt3phR0Mq/VT4vDIkjhfeWHzuS8vrWjktrkbee20SUwMmR73TBdZHt/dYxj/XbUya/43di/H9t2BHOMrNWjpEkY1ZFiWVRCrdkWEBivCNXNqmHb3IibftZDDH17Bf6oMCck9DzQfgiZFKXGOHGNx4Eih2MRxTO+eL99b4nYALroK9jkMRAi98Ajf3T7AiEj7OYW7jwwze1S4JeAphHaQlAC1EubFdYb/rjdskjCe1a0+gwGtCYfLXlrTLr0u6XPR82uoNz2vC8/z2HtUiEt2KWuTLsD/7l/BuEjfbQ2mlOqZofNuqIaspARoNDavr2xCBPYZEyGCR7CfegUbCHLSEyvarJr/pDrJsY8u598nTWC/EXa/95ANlBXh9U6Q8EVXEbjg25j6OsYUObx12kT+d24NDy/cRNgWLtqpjPN2KKFYUtRJiPUJw5LaJNuWBhkWEopI5j1IiVtB7vi4nuvmVLVs7RVxhLuOGMOhYwIEC/RLQW+xLIs318RJ5Whmr69uImEsor3wu8J+gh/MLOfruw3j6aUNBG04cmIRMcsjtJW7Oiil+p4GqGpQi0uQP39azzWvVbXMFbQFfnFABWdNjRHy8/sBJSIsrfdybun07f+s4+kTx1E0RBZt2LZNrWdjRAhgiG7FsapxO5DuTQ2lw5nqRZ9ww6zpXL1nOQAltgd+gg1+iM89voIP129e7DVrVJiHjhtHiSTyFqRalsWHG12+93pVm/Qm13DmU6v4+IvbMsYZ3AEqgNfJ5F1j2PoJw1sImyRhGy6YGkr/bi8+oOfuKjUU9P84l1J9xLIsFm7yueq/VW0WsngGvvnKOpb2YNuhrWXbNi+tzL2qeN6GBF4XP5Vt28ZxnF5/DbZt9/kQuG3bxO0w9y9OctgjK9nxviV88bl1VDbauNKz782e52H7KYpJUEwC47nUE+DUJ1e2CU4B3lob5/xn19BI/qZWNBiHH7+1IWueAX7zfjVi9+9Uj77m+z6zR0ewcjTd3UeGCVm9H0F6ntfvoxNKqa7RAFUNWglsbnw7eyAAcNM7G0j1MBjqLmMM42O5f2dRwOr0j9ITmzoJ89CyFL/7uJFPGx0aJdTjQLVJgqx0g9yzIMETKz3qJNzjYHFLCQlQbcLM22Rx5X/Xc9Hza/i0Okld0ufZZQ3Mun8JH9aaXj+1qt7NfRDB88sbiJuuvRU6jkPQsnCcra+XlA/LN+WeXrKwLkVyCPTuRcXj+r1HtEsP2cIfDhtNcR/vLauUKmw6xK8GraQvrGzIPVS6ot4l6ef3j8DzPA4eHyVsC3GvfRTyPzuVEbNcco3wx0qH8dZGw4mPLyLR8vgq9h8b4YFjxhE1W3cqTtwO87WX1/GPyk0taUFbuP/osRwwyiHQjTmRnpWe85v0DUFLKLFcjOcSlyD3fFbPj+as58mTJnH3/Np2j/UNXPriWp49aRyxLkxzEBFS4pDwhYhlsPz2QY2IsL6p4+falPQp7+A0VAdDJNGIefpJZGklZtpOsO8RxENR3G5ODQhbhpkVYRbVZg/ADhgbISSmx7sNFLqgSXHRjGL2HxvlZ29vYGW9y96jw1w1czjDne5P9VBKDS7ag6oGrYjts/+YSM78A8dGiPbDITIxSfHkiROIbLEB44Hjonxnj2FYHWxUHh4+ihMeX9EqOE17dVUTN8/diNmKVeBiO/zts01tglOApGc49cmVXV5NLSI0WhGufbOGqXcvYsqfFzHz/iX8dWGCuBViWaPhiler2H5YiDfX5A6kP96YINmF2CQhARY0OVz8ygZOfHI117xZQ5UfIlY2rM11xhjGFjk5J04ELCgL5X4rtAUiKxdjLj0Z/vpbzH+fgbtuga+dQmT96pbpEFamZ7WznuwQLj+cPQIny68sCVqcs30Jfi9uVl/Iwn6SPcsM9x42kqdPGMtNe5dTYSexemEPVKXUwKYBqhq0LM/lG3sMI5plJ+6igMUlu2zdqUI9LpfvsVs5fHbuttx/zFhuO7CCd86czP1HjSZmsg9DQ3re5vPLm0hm6XkF+L+PamjYiq15Nnk2t71XnTXPM3D/p3VdGtJuIMB5z67m9g9rWnqH1zZ6XPziGv65qLHlHPu4a4h1cG6qJWDnmpyY4VoOTyyLs+fflvCPyk28vTbO7R/WsON9i1kfHNFuikBUfE6fVpz1ub60YxkRyR0QhZNxzE3fhdQWC+qaGjA3X0UklaQoFSc67y3CzzxIbOlnxNxEznm8xhhGhXye//xEprXqtp05Ksx/T51E2RAb2vZ9n5BJUUQS20/ptk9KKUCH+NUgV26l+O9pk/jKC2tbeu32GRPhD4eOptRK0V/jqJbvUYzHsWNtRBx838UY0+HCYhFhVWPuQKr5pKTuMiKs6WAqxJK6VJfmt1Yn4dllDVnzrnujioc+NwFILwTbdWSYkC3teoIBjp1cREi8DldZNxmHy15a1i496Rm+/OI6Xjh5fJudEEImya0HVlAStPnz/FqSniHiCF/duZwr9iwn4OX+YiC1G6Emx1zmolLsmir8H3wF6us2P2bsJGLX3059KJo14HJ8lz3KbF48eTxNHtgihCyfoq3YxUAppQYjDVDVoCa+x5SI4eFjRpPKDBgEpXACge6eHnTIuDA/zJG/+8gw9laEqA4++4yJ8FKOM8uPnhzrdOWzZVl8sMUK+dY2xD1GR2yeOXkiQVuIOMKdR4zli0+vbLPDwvgih18dNIqwSeSMT0WET6uTNOU4gP6z6iRxD4q26MCMeHF+OruM7+01nCY3HaDGxMPuKDgVgWTu12WdfhH+DV9rE5wCsGop5nc3ELr8+vQWWFl4nkcMj1hzOQ15PTZVKaUKmQaoatDzfZ8YrYZnB2ggYIxhXNgwa1S43Yp0AW4+oIISy6O7u+jExOWm/SvY++9L2tXLlNIAe40K4/u5g7jmso0tyv52Ygvcf+x4HlywiZveXk91wifiCL8/dDQfnzuFuz+uYUW9y1ETohwyIUYxCfzO/oO2csMCx09RQooSm3QPbSe/xhgDI8dAMNQ+UC0tB8+D6qrsD37vdRw3BYN8yyillOoLOgdVqQGkbtUSHjl+PJfuUt4yt3aXESGe+/xEdiqztmqPR9/32TZmePELE9llRHojc8eCM6YV88IXJlJM54cZGGOYUhJgdLR9kPrlncqZuy7O1f9dR3Ui3cPb5Bq+8591hBtquKFoKXduepoznv8lo174O9FkU4f7sBpjmF4WbLfIrNn08iDhXlz8lnSCcNal7TMiMcyW81JbM6b9vFWllFJdoj2oSg0grusS85v40axSrp45DM+kh+hLbA+vBwu+mudEPnXCWFLGwhIIi0fIdKE3M6NYkjxz8gQOf3gZ61rNlb1wxzIO/eeSdtf/+5Bixv36Svj0/bYdok8/SPSGP1LvhHL+riguvztkNBc8u7pNesgW7ji0glLL7XZPci5JsbAPPAZ7wjbw9/+DNStg/DZwzmVQ3n4fzxbFZRAK904heomIkBSH4tET8KwAjnF1UZJSqiBpgKrUAGT7KYrYHJD2RjDWPCcSaBn+7k7oYnyfSWGXd86YzGc1SZbVpditIkzYsWhItX2mmaPCTFjxEfLp++2faNUyeOlfhI85jXgq++It27gcOyHE3LMmc9PbG1lUl2JWRZhv7DGMZNUKPK+DwHErNNkBrGm7ErryZizfw7NtkoEwwVQC56Bj4eUn2z/orItJhCIFM58kKQGWN8FP39rAotoUu410uXLmcEY4Hk439rlVSql80ABVKdVrjO9TTJyZZcJew4L4vksNNkLbYPeEsTbDX3k89xO98iSB7WYQGDGaZLSIRJb9XUMmxbZh4TcHDCPuQ8Q2WF6cytpqqOjdABXSUyGanGDrBOJ2gNj530TGToTH/5peLFUxBs68BG/XvUkVSHDqWQ5PLUtwXqse57lVce79pJZnTprIzHKrIBYNKqVUM52D2gW2bZOQAHEJ9PoRjEoNRsYYPM/DGENYfI7bpqhNvm8ErA7+lsSClUswl51M4N1XCfrZe/iMMdh+ihipDg846EsNVoDEMWdibnsA/u9JzE/vIj7zIJoKaHFUk3G4+MU17dJdHy54bjV1pnDKqpRSoAFqh9In44R5aqXHWc9WccYz6/jHkhSNVrjH554rNVSETYrfHzKKPSo2z8d8eGWK6oNOyvkY2ecwzNzX0guNfnsDIbewFxuljKEhEKY+FKUhEMYtoN5IEWFRXYrGHNtyLa1L0agj/EqpAqND/B1okBBffGY1zy/fvD/kSysa2XVkiKdOGE+0g1N/lFJpxhhKJMHjx49lbZPPRxsSbFsaoDhSBrvvA3Nfb/uAydOQnWdi/v6H9H3fw7z/Jvbeh23VLgUqfTpXRwrp+7Zt24TjDYjrgm2TCgRJiI5cKTXUaICag23bvLE63iY4bfZ+VYIHF9Zz/nYhfP3AVAXEsqz0iVQFtjLbGEOMBFMiwnYTHIzxiRuH2GXXI5UfwdMPgucisw5EJk/Dv/lqaN0L2ZT9hCrVOWMMk4sDFAUs6lPte3a3LQ0Qad4Xtp+FfRfnk3fhrlvSuyU4DoH9jyZwzmU0BCMF166VUn1Hh/hzaPIt/jivNmf+XfNqqfO0+lT/ExEaJcjKVJCnVnl8uMmiXsId7iXaX4wx+L7fEmg0OEGadpqFuew6ZN/DMXNewv/+/8CGdW0eJ7vto72nPRAhxZ8OH93ufIOgLdx95FhKZOu3KOsttm3jLJwPP/tWOjgFcF146Qn4ydeIpnKf6KWUGny0B7UDHX1XT+cV0LiYGrLqM1NRXmjV2z+uyOHZkycyPpgq+NXZnufRFAwTmzID7rip/QUHHI0bK2qfrrrMNh6HjAnywdnb8It3N1JZk2TPijBf262c4XZhHPsbijfCn2/Nnrn4M2TtCmT8FO1FVWqIKLwulgIRsXwu2KE0Z/6525dQYvf/m7oa2lwrwPffWN8mOAVYWe9yxMPLqDPBHI8sLL7vEx85Frnlfph9CJSUw/ht4ZLvYy78Nk32wHgdzUSEoPEpSsUpSjYRcZP93qMdNCkmhlLcul85d+4T5kczSxhpJbBMz3umRQTHcXAcZ6sXkFrGhxWLcl/w8bu6i4pSQ4j2oObgeR4HjA2zz5gIr69uapO3w7Agp08rwfeacjxaqfxoNDb3fVKXNW9lvcviuhS7lMiA6HVKiYVbNpLQV6/FdlMYEZLhGO4AG9q3LItovAH++lv4z9PgprCn70rsoiuIjxhDSvovUDXGYHkpNq1byejSaK88Z5MEWZ+Ef1Sm2+GpU0sYEYSI6ebOCyIQLYLG+uz5I0YPiHaslOodGqB2IGYS/PPYsfx7aSN/nFeD6xvO2b6UU7YrosjEC2FNgRrikp4h6eVuicvrXXYttQfMB7sxhrjlQDDz1jQAglPLsgim4ojn4QXDBJNN8IOvwOplmy/69H3MVecRvvlvuOUVA+b/ozNNVojr39zI7R/WtKT9aM4GvrpzGT/caxgRv+vzRhPBCKGjT4GH/tw+MxiGnWbqPGSlhhANUDtgjCFq4pw6OcDRE0ZhkMyZ5xqcqsIQsoXhYZsN8ewf3DsND+H7hb2H6EAW9VJYlfPh4buhej32Iccj4yZjWgenzVwXc++vCV3yfeIFtIn/1rIsiw82pNoEp81u/7CGk6cUs8/wrp9QlTKG0PFnweLPYO5rmzPCUeT7v6YpEOqtoiulBgANULvA8zyimTPK9Qu8KiTFlsv39hrON19Z1y7vgLFRRoSkILYPGozCXgrrkbvhsfta0mTtSli5NPeD5r+L7bkwCALUJmNzy9yqnPm3zN3IPYeNJETX5+rX20Gil12HVV8Ln7wPZcNhyg40BUK4uihVqSFFA1SlBjLP5Yyp6RXuN7y5gY1xj4AFp08r4ef7VRDTqSh9xkk0tQlOAUxjPTJyTO4HFZdjBkmg5RmoSeQOPqsTPh3MPsmp0Qki5RVY+x3Vsi2ZUmro0QBVqQEu4ie4YGqEU7abTNwzhGwhIh4Bv0mD0z5i2zZ88Gb7jHdfQ37wG8xj96WPad3S8WeSjMQGxVBMzDYcOznGnDXZF4sePzlGzDZszSYBxhidb6rUEKfbTCk1GPguxSbOSCtBiYkT8Pt/4/VBz8kyTB9vxMx5EbnwO7DltlKzD8Hse8SA25UgF+O5XDijjIpo+62fKqI2F8wow3huP5RMKTUYaA+qUkp1k+d5yM4zMZbV9khWwDx6L/LlK5DbH8fMfQ0a6mGP/fBLymgcBHNPWyuVBG+cNpnvvV7Fgws2AXDKdsX8eJ+RlEoCX7vwlVJbSQNUpVQLsWzqTIAG12ALRCwotV0dbs0iGQgRuPDbcMf/ts0oLoM99qM+XIS139GICJ7nDZqtpVrzfZ/hVoJf7j+Mm/YbCUBIPIImga/RqVKqBzRAVUoBkBSHN9Z5fOWFpaxuSA/NTi8P8pejxzIlCuJrkNpawnKQfY/E2XFPePwvsLEKdtkbDjyaxmBkyMz1tJ4FAAAbiklEQVSjNMYQMCkCZKaVGN04QinVcxqgKqUQEZY3CSc+vqJNcPFpdZKDHlzKx+dsSymDP9jqrrgdQIaPIXD+txHPxXOCuL6ffYGUUltJRBAR3dFADSm6SEopRUICXDdnfdaer4aU4ZllDYSNT8hz9Tz0LRhjSJp0j6qrAYTqRZZlEUslKKpaSfTT9yhqqCXi6sEbamjQHlSlFAkfPlzf/ljK8pDFvw8tZtrad3Eeehh8j8ChJ2B2mkVjIDQo51UqVQgsyyJWtwFzw+WYNSta0u1ZB1J08feot4P9WDql+p4GqEopAgLblARYUtd2e6p/H1rMbvdej/3RW5sT358DU3ckdvWt+iGpVB+JJBox110CVavbZrz1Cgz/I8EzLiYpOgiqBi9t3UopouLyg9kj2qTNHBVm8qqP2wanGVI5D956BUeH+5XqdSKCVK1uH5w2e+Exgqn2Ix5KDSYaoCql8H2fHUptfn3QKCJO+ijO08bbDH/lkdwPeu5hAvHGPJVQqaFDRGDtqtwXJBPg6mEcanDTIX6lFABhk+TMKWFO2m4Ki2qSbBdKwnMdrNzXBUFK9QljDEyakvuC4lII6PQaNbhpD6pSqoVjXIr9JnYr9SmNheCQz+W++KDjSIUj+SucUkOEMQa/pBym7pj9gs9fQCKof3tqcMtLgCoiE0TkRRGZLyLzROTrmfRhIvKsiFRm/i3PpIuI/EpEFojIByKyRz7KqdRgI7ZNvYRY64XYaMI0SAhryzPiszDG4LouZpe9YMqM9heM3xb2PxLX015UpfpCUyAMV94Msw+B5r/ZSAzOugT/4ONI6Q4aapDL1xC/C3zbGPOuiBQD74jIs8D5wPPGmBtF5CrgKuBK4BhgauZnNvD7zL9KqS5KSICXV6W4/KWVrG1MD9XvPjLMfUeNYXzIxXRhiL7RCRG79jZ462V47pH0sP7Bx8F+R9LghPplQ3rLsjDG6BZXalAzxtAQjBD6yjU4X/oupJIQCpMIhkgZ6e/iKdXn8hKgGmNWA6sztzeJyHxgHHAicHDmsruBl0gHqCcC95j0J9AbIlImImMyz6OUasWyLOLYeD7EbIPvuViWxWe1Pqc9ubLNtXOr4uz/j6V8cPY2FBPv9LmNMdTbQZz9jiYw80AwhlQ4mu45zWOAKCI0EGBjEj5Yn2BckcOUkgDFkuxSoK3UQGSMIW4HwA5AKJpJ7N8yKZUveV8kJSKTgd2BOcCo5qDTGLNaRCoyl40Dlrd62IpMmgaoSrXSZIV4f0OKX763nrqkx0nbFnHm9FJC4nPN62uyPqY64fPE4nrOmRLq8lnxrufhOqH0nX4Y1m+QEF98ZjXPL9+8a8DoqMMzJ09gUrhrvcFKKaUGDsnnMJmIFAEvAz8xxjwkIjXGmLJW+dXGmHIR+RfwM2PMq5n054ErjDHvNF9bW1vbUvDKysq8vQalCkXZ2En8+IMEd39S1yZ9fJHDf06dxEEPLmXZJjfrY8+aXsIN27ts3LA+H0XtkZIRo/hFpcUf59W2y6uI2rz++bFsXPJpP5RMKaXU1po6dWrL7dLS0nbzVvLWgyoiAeCfwF+MMQ9lktc2D92LyBhgXSZ9BTCh1cPHAzk3hWv9IgtVZWXlgChnf9I66lxzHYkIC5sc7v6kqt01K+pdfv72Br66cznXvNY+H2Cn4SFGVZQzfFh5Xxe5xzZJmHs+WZQ1b12jx5JGYea0aS1zUrUddU7rqHNaR53TOuqc1tHWy9cqfgH+BMw3xtzSKusx4LzM7fOAR1uln5tZzb83UKvzT5XazHEc7t2i57S1v31Wx4lTirPmBW3hzOkluG723tVCk/QNCS/3SM/yTaku7UyglFJq4MjXu/p+wBeBQ0XkvczPscCNwBEiUgkckbkP8CSwCFgA/BG4JE/lVGrA8DuYnmMMFAcsrps9gkCrv/JhYZtnTppAiQycU2iCljAqmvtI1d1GhvF1DqpSSg0q+VrF/yqQa1+Mw7Jcb4BL+7RQSg1gruvyxe1L+eV71VnzT59aQkxcLt2xiAt3LOWT6iRFAYuJRQ4xSSF+1xZHFYISy+W62SO5+MX2i772HROhIiy65ZRSSg0yOi6m1ABkjGFc1OLMaSXt8kZHHb6313BsP0XApCgxcWaXG3Ys8igy8QEVnAIYz+XkbaP89uBRjIyke1IDFpyzfQn/OHYcMZL9XEKllFK9Le/bTCmlekfEJLjlgBGcMb2E2+ZuZFPK54RtirhgRhmlkqD1qPdA72EM+wnO2S7MCdtOIuGl59FGxCPgN+m2kEopNQhpgKrUABbxExxSYTH7yAo8IxTZPsZrwh+EUZvxXIpxKbZIb1Y+CF+jUkqpNA1QlRrgfN8nTLq71Ays0XullFIqKw1QleoG33Ko9y0EKLF9jDcwtmpSSimlBhINUJXqAsuyqDFBbn23mocWbiJsC1/esZSzp5cQM4kBP8dTKaWUKiQaoCrVBTUmyKz7l7C2cfMY+ndfreKByk08cvw4Yibej6UbWESEpDgkfSFiGSx/4OzJqpRSKj90mymlOuHbDrfOrW4TnDZ7a22ct9cl9CSjLkpIgMomh4tf2cAJ/1rNNW/WUOWH8CX3RvxKKaWGHu1BVaoTDZ7NQws35cy/Z34tB44ajqCnGXXEFYfHlsa56PnNG+6/sy7OnR/X8uLnJ7JLqY3n6SovpZRS2oOqVBcYQnaug9Ag6lg5j0lTmzXhcPlLa9ulJz3DBc+tptbT78tKKaXSNEBVqhOlts9FO5bmzL9klzLE19X8HRER5lcniXvZF5N9Vp2kSTuglVJKZWiAqnqNiBCwLYICtj145hT6nsvZ25cya1S4Xd6XdixlUpGlq/i7oLNeZu2FVkop1UzH1FSvCHspnLpq+PcDULMRZh0Eu+1NYzCC7w/8rrGYH+fR48fx1roE98yvJepYXLJLGZOKLKJGz4LvjDGG7ctDRByhyW0fzG9fHiRso6dDKaWUAjRAVb0g5Lk4L/8L7r5tc+IbL0DFGKI/+RMNodiA72E0xhA1cQ6psDILokB8F2N0aL+rIuLyu0NGc8Gzq9ukh2zhz0eOodRy0TVSSimlQIf4VS8IJBvbBqfN1q2Gv/6OgD94og7f9xEvBV5qwAfd+eb4LsdNCDH3rG04c1oJe4+JcPmu5cz74rZMKzK6gl8ppVQL7UFVPeI4Dsx5OfcF/3ma4NmXkQy0n7+php6gSbFtWPjNAcNI+BCxDOLH0R26lFJKtaYBquq5ZFPuvCF6Vr2xHOqNQ23SJ2gLURtKJDUo5uP2lDEG26SIggamSimlstIAVfWI53kw8yD4y++yX7DrbFxr8Kzo74qEFeCpZQm++cpyqhPpCGxaeZC/HzOWbSICeZryICLYto2I4LquTklQSik1YOgc1EHCsSwiqQSRVAInj1s8GWPwSsth38PbZwbDyIXfIeEE81ae/mZZFh/X+Jz/7OqW4BTS+3we+OAy6kzX6sK2beolRL2EqZcQlt2975KNEuTTRofvvFHLN16rYV69TaOEuvUcSimlVH/RHtQBTkSIuQl49Xl46QlAsA87EfY6iAYnlJdesyY7SOyiK5E994fH/wr1tbDzLOTUi2iKlQypnrsG4/C919dkzduU9Pnngjoumh7BdXNPfUhIkP+uTnLt62uYvzHJuCKHq/YczqnbFSHS+W6hTVaIq/67gXs+qW1J++NHNZw8pYjfHTyKiB/v8PFBfIKpBHgexnaIh2O6gEkppVReaYA6wMXcBFx/KSyt3JxY+RE8+xCxa39JvZ2f3ssGK4A9+zCCu+6NGB83ECZl2UMqOAVwjTB/Y+59UeesjXPh9GjuJ7AdHl3cxP+0Oq9+Zb3L5S+v5dOaJJdPGdPh77csiw82pNoEp80eXljPmdNLOWZs9jPvRYRYKg5/+z28/CSkksjoCUTO/ybu9F2IW/p2oZRSKj90iH8AcywL5rzYNjhttnA+vPdGXk908jyPpkCYxmCUpAzN05UsMWxXlvtLwa4jQtgd7Ea/yXe46r9VWfN+90E1diz3kasACWx++V51zvxfv19NvZ+9TURTccyN34LnHoFUJshesxxu/BbOgnmD6nQwpZRShU0D1AHMSTbBi0/kvuCFRwkkOx7OVb2rxPL44ezhWfMijnDW9JIOh/cbXMPGePbhdN/Ap9XJDof5PQP1qdxL4+uTPn6W+FhEkPVrkcp52R/451sJxRtyPq9SSinVmzRAHeikg//CIbZ6vhB4nsfMEQFuObCCiLM5kBwbc3j+8xMplo633QrZHc8xLQt1/H8atXxOmVqcM//kKUUU2e0DWMuy4OO5uZ945RIs3SJLKaVUnuiksgEsFYriHH4SfPp+9guOOJlUMIyeH5lfIT/JudtFOG3qFJZvShFxLIaHLYol1ekWUyHx2XdMhNdWt99bdnTUYXSYDqdO+J7HF6YUceu7G1lcl2qTNzbmcOGOZRiv/XMbY2DEqNwFC0fB0u+zSiml8kM/cQYwz/Mwe+6HmbZzuzwzYw/Mjnvq6ut+4hiXIr+JHWIuk0NJik28S/ufxkhxz5FjGF/U9rtjcdDiiRPGE1+/qtPnKCbJK6dM4pu7D6MiajM8bHPpLuW8fvokSkhkfYzv+7DDbhCOZH/SI07W08CUUkrljfagDnCNTojoVTfDR2/Dc4+CAEd8HmbsTqMTgiG4UGkgM8ZQ4aR44/RJvLM2zmtrmpgxLMSh46MUS4rKFfWdPofv+xTTxPf3KOZbu5UB6aF/249nnX/aLB4IEfnBbzE/ugzijZvLtNNM5OTzyL03gVJKKdW7NEAd4IwxNNhB7D0OILDzXgCkAqF0z6kGpwOS7/sUEeeQURaHjSlOH4bgxelg8X9Wlu9STGbOaxemj7oITeO2Ifybh2DBPKheD9vviikZpl92lFJK5ZUGqIOE53l4zftU6rD+oOD7fnroPY9chPpAGGvHWYgIvu+n56dqcKqUUiqPNEBVSrWT78BYKaWUak0XSSmllFJKqYKiAapSSimllCooGqAqpZRSSqmCogGqUkoppZQqKBqgKqWUUkqpgqIBqlJKKaWUKigaoCqllFJKqYKiAapSSimllCooGqAqpZRSSqmCogGqUkoppZQqKBqgKqWUUkqpgqIBqlJKKaWUKigaoCqllFJKqYLi9HcBlOoPnhWg0dhsTPjEAhYRy6dYXHzf7++iKaWUUkOeBqhqyIlbIf4wr44b395Ak2sA2G1kmAeOGcvoQEqDVKWUUqqf6RC/GlIs2+bRxY388I31LcEpwHtVcQ7+5zLqCPZj6ZRSSikFGqCqIabOd/jRm+uz5q1qcHmvKoFl6Z+FUkop1Z/0k1gNKT7Cyno3Z/57VXENUJVSSql+pp/EakgRDGNiuade7zwipHNQlVJKqX6mAaoaUkotj6tnDs+aVxG12XNkWANUpZRSqp9pgKqGFN9zOW27GFfsOYxAq9Y/vTzIy1+YRIkk+69wSimllAJ0myk1BIX9BN/dtYTLdi1nVYNLSdCmJAAloltMKaWUUoVAA1Q1JAVMigApSqIALhjwTWePUkoppVQ+6BC/UkoppZQqKBqgKqWUUkqpgqIBqlIDhOM4lI8Yie3ozByllFKDm37SKVXgxLKoM0Ee/qye/6yymV5ezwUzyii1XByT+9ABpZRSaqDSAFWpAmZZFiuSAfZ7YDHVic07DNz49gYe/9wEZg93sDRIHdKCxifoJmHjOgiGMUUlxIMRPKOr/pRSA1fBBqgicjTwS8AG7jDG3NjPRVIq7+qNw7nPrGoTnAK4Ppz+1Eo+PmcbitAAdaiKeCnsl5+Av90OyQQAMmIU0StvpmnUeFyNUZVSA1RBzkEVERv4LXAMMAM4U0Rm9G+plMq/hLF4e208a15d0qeyNoWI5LlUqhDYto392Qdw9y9bglMA1q/FfP8iwsmm/iucUkr1UEEGqMBewAJjzCJjTBK4Hzixn8ukVN65nWzO2pjSgwWGqmC8Af7+h+yZ8SZ440UcXVCnlBqgCvXdaxywvNX9FcDsXBdXVlb2eYF6w0ApZ3/SOmqrfOJUJhQ7LN/UfhjfsWCHMpvPPpvfDyUrbEOhHe0weiT26uU582XRfNbtuDcbNmzImj8U6qintI46p3XUOa2j7KZOndphfqEGqNnGLHN2JXX2IgtBZWXlgChnf9I6as+ybX53yGhOeGxFuz+Aq2eOoMgRhmmdtTFU2lEglYDx28CCj7Pmm2m7MHLkSIYNG9Yub6jUUU9oHXVO66hzWkdbr1CH+FcAE1rdHw+s6qeyKNVvfM9j9gibOWdM5oiJMUZGbGaOCvPQceO4dOcSgibV30VU/SQZjsJZl2TPjBXBzAPwPC+/hVJKqV5SqD2obwFTRWQbYCVwBnBW/xZJqf4RNC7bxyzuO7yC2sY4sUiYEtvD8xKdP1gNWp7n4U6ejnPxtXDPL6GhPp0xbjLy3Z/TGAx3MO6klFKFrSADVGOMKyKXAU+T3mbqTmPMvH4ullL9xvd9QiSpXbGQiqlT0Y4xBRC3AwT2OYLQngdAQz0SCOAFw8RDEXxfF9AppQauggxQAYwxTwJP9nc5lFKqkKUMpAJhKAtvTtTgVCk1wBXqHFSllFJKKTVEaYCqlFJKKaUKigaoSimllFKqoGiAqpRSSimlCooGqEoppZRSqqBogKqUUkoppQqKBqhKKaWUUqqgaICqlFJKKaUKihgzMM/Cq62tHZgFV0oppZRSLUpLS2XLNO1BVUoppZRSBUUDVKWUUkopVVAG7BC/UkoppZQanLQHVSmllFJKFRQNUHuBiEwQkRdFZL6IzBORr2fSh4nIsyJSmfm3PJMuIvIrEVkgIh+IyB79+wryR0RsEZkrIk9k7m8jInMydfR3EQlm0kOZ+wsy+ZP7s9z5IiJlIvKgiHySaU/7aDtqS0S+mfk7+0hE/iYiYW1HICJ3isg6EfmoVVq3246InJe5vlJEzuuP19JXctTR/2b+3j4QkYdFpKxV3tWZOvpURI5qlX50Jm2BiFyV79fRl7LVUau874iIEZERmfvajtqmX55pF/NE5KZW6UOuHfUKY4z+9PAHGAPskbldDHwGzABuAq7KpF8F/Dxz+1jgKUCAvYE5/f0a8lhX3wL+CjyRuf8AcEbm9u3AxZnblwC3Z26fAfy9v8uep/q5G/hy5nYQKNN21KZ+xgGLgUir9nO+tiMDcCCwB/BRq7RutR1gGLAo82955nZ5f7+2Pq6jIwEnc/vnrepoBvA+EAK2ARYCduZnIbBt5m/0fWBGf7+2vqyjTPoE4GlgKTBC21G7dnQI8BwQytyvGMrtqDd+tAe1FxhjVhtj3s3c3gTMJ/1BeiLpgIPMvydlbp8I3GPS3gDKRGRMnouddyIyHjgOuCNzX4BDgQczl2xZR8119yBwWOb6QUtESki/8f0JwBiTNMbUoO1oSw4QEREHiAKr0XaEMeYVYOMWyd1tO0cBzxpjNhpjqoFngaP7vvT5ka2OjDHPGGPczN03gPGZ2ycC9xtjEsaYxcACYK/MzwJjzCJjTBK4P3PtoJCjHQHcClwBtF64ou1os4uBG40xicw16zLpQ7Id9QYNUHtZZghxd2AOMMoYsxrSQSxQkblsHLC81cNWZNIGu9tIv8H5mfvDgZpWHw6t66GljjL5tZnrB7NtgSrgLklPg7hDRGJoO2phjFkJ/AJYRjowrQXeQdtRLt1tO0OuTW3hQtI9gqB11EJETgBWGmPe3yJL62izacABmalEL4vIrEy61tFW0gC1F4lIEfBP4BvGmLqOLs2SNqi3UxCR44F1xph3WidnudR0IW+wckgPG/3eGLM70EB6WDaXIVdHmTmUJ5IeKhsLxIBjslw6lNtRV+SqlyFbXyJyLeACf2lOynLZkKsjEYkC1wI/yJadJW3I1VGGQ3o6w97Ad4EHMqM1WkdbSQPUXiIiAdLB6V+MMQ9lktc2D7lm/m3u8l9Bej5Ps/HAqnyVtZ/sB5wgIktID2UcSrpHtSwzVAtt66GljjL5pWQfdhpMVgArjDFzMvcfJB2wajva7HBgsTGmyhiTAh4C9kXbUS7dbTtDsU2RWcRzPHC2MaY5SNA6SptC+gvh+5n37/HAuyIyGq2j1lYAD2WmO7xJeqRwBFpHW00D1F6Q+Zb0J2C+MeaWVlmPAc2rF88DHm2Vfm5mBeTeQG3zMNxgZYy52hgz3hgzmfRilReMMWcDLwKnZC7bso6a6+6UzPWD+tulMWYNsFxEpmeSDgM+RttRa8uAvUUkmvm7a64jbUfZdbftPA0cKSLlmd7qIzNpg5aIHA1cCZxgjGlslfUYcIakd4LYBpgKvAm8BUyV9M4RQdLvZ4/lu9z5Yoz50BhTYYyZnHn/XkF6UfAatB219gjpjhdEZBrphU/r0Xa09fp7ldZg+AH2J901/wHwXubnWNJz3Z4HKjP/DstcL8BvSa/g+xCY2d+vIc/1dTCbV/FvS/qPdQHwDzavgAxn7i/I5G/b3+XOU93sBrydaUuPkB4y0nbUto6uBz4BPgLuJb06dsi3I+BvpOflpkgHEV/amrZDeh7mgszPBf39uvJQRwtIzwVsfu++vdX112bq6FPgmFbpx5LerWUhcG1/v66+rqMt8peweRW/tqPN7SgI3Jd5X3oXOHQot6Pe+NGTpJRSSimlVEHRIX6llFJKKVVQNEBVSimllFIFRQNUpZRSSilVUDRAVUoppZRSBUUDVKWUUkopVVA0QFVKqX4mIk9lNovPlX+7iHw/n2VSSqn+pNtMKaVULxGR+lZ3o0AC8DL3v2KM+Uv7R7V7jvOBLxtj9u/9Eiql1MDgdH6JUkqprjDGFDXfzhwL+WVjzHP9VyKllBqYdIhfKaX6mIhcJyL3tbo/WUSMiDiZ+y+JyJdFZAfgdmAfEakXkZpM/p9F5MetHn+RiCwQkY0i8piIjG2VZ0TkqyJSKSLVIvLbzLGwSik1YGiAqpRSBcIYMx/4KvC6MabIGFO25TUicijwM+A0YAywFLh/i8uOB2YBu2auO6ovy62UUr1NA1SllBpYzgbuNMa8a4xJAFeT7nGd3OqaG40xNcaYZcCLwG75L6ZSSm09DVCVUmpgGUu61xQAY0w9sAEY1+qaNa1uNwJFKKXUAKIBqlJK9b0G0qv6m43u4NrOtlZZBUxqviMiMWA4sHKrS6eUUgVGA1SllOp77wEHishEESklPSyfy1pgvIgEc+T/FbhARHYTkRDwU2COMWZJr5ZYKaX6kQaoSinVx4wxzwJ/Bz4A3gGe6ODyF4B5wBoRWZ/luZ4Hvg/8E1gNTAHO6O0yK6VUf9KN+pVSSimlVEHRHlSllFJKKVVQNEBVSimllFIFRQNUpZRSSilVUDRAVUr9f7t1LAAAAAAwyN96GjuKIgBYEVQAAFYEFQCAFUEFAGBFUAEAWBFUAABWAuNO2gprZMExAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,6))\n",
    "sns.scatterplot(x=\"Tuition\", y=\"enem_score\", hue=\"Tablet\", data=data, s=70).set_title('ENEM score by Tuition Cost')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we understand the problem, let's look at the solution. We can also say what would be necessary to make association equal to causation. **If \\\\(E[Y_0|T=0] = E[Y_0|T=1]\\\\), then, association IS CAUSATION!** Understanding this is not just remembering the equation. There is a strong intuitive argument here. To say that \\\\(E[Y_0|T=0] = E[Y_0|T=1]\\\\) is to say that treatment and control group are comparable before the treatment. Or, in the case that the treated had not been treated, if we could observe its \\\\(Y_0\\\\), then its outcome would be the same as the untreated. Mathematically, the bias term would vanish:\n",
    "\n",
    "$\n",
    "E[Y|T=1] - E[Y|T=0] = E[Y_1 - Y_0|T=1]\n",
    "$\n",
    "\n",
    "Also, if the treated and the untreated only differ on the treatment itself, that is,  \\\\(E[Y_0|T=0] = E[Y_0|T=1]\\\\)\n",
    " \n",
    "We have that the causal impact on the treated is the same as in the untreated (because they are very similar).\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "E[Y_1 - Y_0|T=1] &= E[Y_1|T=1] - E[Y_0|T=1] \\\\\n",
    "&= E[Y_1|T=1] - E[Y_0|T=0] \\\\\n",
    "&= E[Y|T=1] - E[Y|T=0]\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "Hence, in this case, the **difference in means BECOMES the causal effect**:\n",
    "\n",
    "$\n",
    "E[Y|T=1] - E[Y|T=0] = ATE = ATET\n",
    "$\n",
    "\n",
    "Once again, this is so important that I think it is worth going over it again, now with pretty pictures. If we do a simple average comparison between the treatment and the untreated group, this is what we get (blue dots didn't receive the treatment, that is, the tablet):\n",
    "\n",
    "![img](./data/img/intro/anatomy1.png)\n",
    "\n",
    "Notice how the difference in outcomes between the two groups can have two causes:\n",
    "\n",
    "1. The treatment effect. The increase in test score that is caused by giving kids tablets.\n",
    "2. Other differences between the treatment and untreated that are NOT the treatment itself. In this case, treated and untreated differ in the sense that the treated have a much higher tuition price. Some of the difference in test scores can be due to the effect of tuition price on better education. \n",
    "\n",
    "The true treatment effect can only be obtained if we had godlike powers to observe the potential outcome, like in the left figure below. The individual treatment effect is the difference between the unit's outcome and another theoretical outcome that the same unit would have in case it got the alternative treatment. These are the counterfactual outcomes and are denoted in light color.\n",
    "\n",
    "![img](./data/img/intro/anatomy2.png)\n",
    "\n",
    "In the left plot, we depicted what is the bias that we've talked about before. We get the bias if we set everyone to not receive the treatment. In this case, we are only left with the \\\\(T_0\\\\) potential outcome. Then, we see how the treated and untreated groups differ. If they do, it means that something other than the treatment is causing the treated and untreated to be different. This something is the bias and is what shadows the true treatment effect.\n",
    "\n",
    "Now, contrast this with a hypothetical situation where there is no bias. Suppose that tablets are randomly assigned to schools. In this situation, rich and poor schools have the same chance of receiving the treatment. Treatment would be well distributed across all the tuition spectrum.\n",
    "\n",
    "![img](./data/img/intro/anatomy3.png)\n",
    "\n",
    "In this case, the difference in the outcome between treated and untreated IS the average causal effect. This happens because there is no other source of difference between treatment and untreated other than the treatment itself. All the differences we see must be attributed to it. Another way to say this is that there is no bias.\n",
    "\n",
    "![img](./data/img/intro/anatomy4.png)\n",
    "\n",
    "If we set everyone to not receive the treatment in such a way that we only observe the \\\\(Y_0\\\\)s, we would find no difference between the treated and untreated groups. \n",
    "\n",
    "This is the herculean task causal inference is all about. It is about finding clever ways of removing bias, of making the treated and the untreated comparable so that all the difference we see between them is only the average treatment effect. Ultimately, causal inference is about figuring out how the world really works, stripped of all delusions and misinterpretations. And now that we understand this, we can move forward to mastering some of the most powerful methods to remove bias, the weapons of the Brave and True to identify the causal effect. \n",
    "\n",
    "\n",
    "# Key Ideas\n",
    "\n",
    "So far, we've seen that association is not causation. Most importantly, we've seen precisely why it isn't and how can we make association be causation. We've also introduced the potential outcome notation as a way to wrap our head around causal reasoning. With it, we saw statistics as two potential realities: one in which the treatment is given and another in which it is not. But, unfortunately, we can only measure one of them, and that is where the fundamental problem of causal inference lies.\n",
    "\n",
    "Moving forward, we will see some of the basic techniques to estimate causal effect, starting with the golden standard of a randomised trial. I'll also review some statistical concepts as we go. I'll end with a quote often used in causal inference classes, taken from a kung-fu series:\n",
    "\n",
    "```\n",
    "'What happens in a man's life is already written. A man must move through life as his destiny wills.' -Caine\n",
    "'Yes, yet each man is free to live as he chooses. Though they seem opposite, both are true.' -Old Man\n",
    "```\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's 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",
    "The beer analogy was taken from the awesome [Stock Series](https://jlcollinsnh.com/2012/04/15/stocks-part-1-theres-a-major-market-crash-coming-and-dr-lo-cant-save-you/), by JL Colins. This is an absolute must read for all of those wanting to learn how to productively invest their money.\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": 4
}