14-Difference-in-Difference.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Three Billboards in the South of Brazil\n",
    "\n",
    "I remember when I worked with marketing and a great way to do it was internet marketing. Not because it is very efficient (although it is), but because it is very easy to know if its effective or not. With online marketing, you have a way of knowing which customers saw the ad and you can track them with cookies to see if they ended up on your landing page. You can also use machine learning to find prospects that are very similar to your customers and present the ad only to them. In this sense, online marketing is very precise: you target only those you want to and you can see if they respond as you would like them to. \n",
    "\n",
    "But not everyone is susceptible to online marketing. Sometimes you have to resort to less precise techniques, like a TV campaign or placing a billboard down the street. Usually, diversity of marketing channels is something marketing departments look for. But if online marketing is a professional fishing rod to catch that specific type of tuna, billboard and TV are giant nets you throw at a fish shoal and hope you catch at least some good fish. But another problem with billboard and TV ads is that it is much harder to know how effective they are. Sure, you could measure the purchase volume, or whatever you want to drive, before and after placing a billboard somewhere. If there is an increase, there is some evidence that the marketing is effective. But how would you know if this increase is not just some natural trend in the awareness of your product? In other words, how would you know the counterfactual \\\\(Y_0\\\\) of what would have happened if you didn't set up the billboards? \n",
    "\n",
    "![img](./data/img/diff-in-diff/secrets.png)\n",
    "\n",
    "One technique to answer these types of questions is simple Difference-in-Difference, or diff-in-diff for close friends. Diff-in-diff is commonly used to access the effect of macro interventions, like the effect of immigrants on unemployment, the effect of gun law changes in crime rates or simply the difference in user engagement due to a marketing campaign. In all these cases, you have a period before and after the intervention and whsh to assess what untangles the impact of the intervention from a general trend. As a motivating example, let's look at a question similar to the one I had to answer.\n",
    "\n",
    "In order to figure out how good billboards were as a marketing channel, we've placed 3 billboards in the city of Porto Alegre, the capital of the state of Rio Grande do Sul. As a note for those not very familiar with Brazilian geography, the south of the country is one of the most developed regions, with lower poverty rates when compared to the rest of the county. Having this in mind, we decided to also look at data from Florianopolis, the capital city of the state of Santa Catarina, another state from the south region. The idea is that we could use Florianopolis as a control sample to estimate the counterfactual \\\\(Y_0\\\\). What we were trying to boost with this particular campaign was deposits into our savings account (By the way, this was not the true experiment, which is confidential, but the idea is very similar). We've placed the billboard in Porto Alegre for the entire month of June. The data we have looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import style\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "style.use(\"fivethirtyeight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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>deposits</th>\n",
       "      <th>poa</th>\n",
       "      <th>jul</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>42</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>52</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>119</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>21</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   deposits  poa  jul\n",
       "0        42    1    0\n",
       "1         0    1    0\n",
       "2        52    1    0\n",
       "3       119    1    0\n",
       "4        21    1    0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_csv(\"data/billboard_impact.csv\")\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "deposits is our outcome variable. POA is a dummy indicator for the city of Porto Alegre. When it is zero, it means the samples are from Florianopolis. Jul is a dummy for the month of July, or for the post intervention period. When it is zero it refers to samples from May, the pre-intervention period.\n",
    "\n",
    "## DID Estimator\n",
    "\n",
    "To avoid confusion between Time and Treatment, I'll use D to denote treatment and T to denote time from now on. Let \\\\(Y_D(T)\\\\) be the potential outcome for treatment D on period T. In an ideal world where we have the ability to observe the counterfactual, we would estimate the treatment effect of an intervention the following way:\n",
    "\n",
    "$\n",
    "\\hat{ATET} = E[Y_1(1) - Y_0(1)|D=1]\n",
    "$\n",
    "\n",
    "In words, the causal effect is the outcome in the period post intervention in case of a treatment minus the outcome in also in the period after the intervention, but in the case of no treatment. Of course, we can't measure this because \\\\(Y_0(1)\\\\) is counterfactual. \n",
    "\n",
    "One way to solve this is a before and after comparison.\n",
    "\n",
    "$\n",
    "\\hat{ATET} = E[Y(1)|D=1] - E[Y(0)|D=1]\n",
    "$\n",
    "\n",
    "In our example, we would compare the average deposits from POA before and after the billboard was placed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "41.04775"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "poa_before = data.query(\"poa==1 & jul==0\")[\"deposits\"].mean()\n",
    "\n",
    "poa_after = data.query(\"poa==1 & jul==1\")[\"deposits\"].mean()\n",
    "\n",
    "poa_after - poa_before"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This estimator is telling us that we should expect deposits to increase R$ 41,04 after the intervention. But can we trust this?\n",
    "\n",
    "Notice that \\\\(E[Y(0)|D=1]=E[Y_0(0)|D=1]\\\\), so this estimation above is assuming \\\\(E[Y_0(1)|D=1] = E[Y_0(0)|D=1]\\\\). It is saying that in the case of no intervention, the outcome in the latter period would be the same as the outcome from the starting period. This would obviously be false if your outcome variable follows any kind of trend. For example, if deposits are going up in POA, \\\\(E[Y_0(1)|D=1] > E[Y_0(0)|D=1]\\\\), i.e. the outcome of the latter period would be greater than that of the starting period even in the absence of the intervention. With a similar argument, if the trend in Y is going down, \\\\(E[Y_0(1)|D=1] < E[Y_0(0)|D=1]\\\\). \n",
    "\n",
    "So this didn't work. Another idea is to compare the treated group with an untreated group that didn't get the intervention:\n",
    "\n",
    "$\n",
    "\\hat{ATET} = E[Y(1)|D=1] - E[Y(1)|D=0]\n",
    "$\n",
    "\n",
    "In our example, it would be to compare the deposits from POA to that of Florianopolis in the post intervention period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-119.10175000000001"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fl_after = data.query(\"poa==0 & jul==1\")[\"deposits\"].mean()\n",
    "poa_after - fl_after"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This estimator is telling us that the campaign is detrimental and that customers will decrease deposits by R$ 119.10. \n",
    "\n",
    "Notice that \\\\(E[Y(1)|D=0]=E[Y_0(1)|D=0]\\\\), so we are assuming we can replace the missing counterfactual like \\\\(E[Y_0(1)|D=0] = E[Y_0(1)|D=1]\\\\). But notice that this would only be true if both groups have a very similar baseline level. For instance, if Florianopolis has way more deposits than Porto Alegre, this would not be true because \\\\(E[Y_0(1)|D=0] > E[Y_0(1)|D=1]\\\\). On the other hand, if the level of deposits are lower in Florianopolis, we would have \\\\(E[Y_0(1)|D=0] < E[Y_0(1)|D=1]\\\\). \n",
    "\n",
    "So this didn't work as well. To solve this, we can use both space and time comparison. This is the idea of the difference in difference approach. It works by replacing the missing counterfactual counterfactual the following way:\n",
    "\n",
    "$\n",
    "E[Y_0(1)|D=1] = E[Y_1(0)|D=1] + (E[Y_0(1)|D=0] - E[Y_0(0)|D=0])\n",
    "$\n",
    "\n",
    "What this does is take the treated unit before the treatment and it adds a trend component estimated using the control \\\\(E[Y_0(1)|T=0] - E[Y_0(0)|T=0]\\\\). In words, it is saying that the treated, had it not been treated, would look like the treated before the treatment plus a growth factor that is the same as the growth of the control. It is important to notice that this assumes that the trends in the treatment and control are the same:\n",
    "\n",
    "$\n",
    "E[Y_0(1) − Y_0(0)|D=1] = E[Y_0(1) − Y_0(0)|D=0]\n",
    "$\n",
    "\n",
    "where the left hand side is the counterfactual trend. Now, we can replace the estimated counterfactual in the treatment effect definition \\\\(E[Y_1(1)|D=1] - E[Y_0(1)|D=1]\\\\)\n",
    "\n",
    "$\n",
    "\\hat{ATET} = E[Y(1)|D=1] - (E[Y(0)|D=1] + (E[Y(1)|D=0] - E[Y(0)|D=0])\n",
    "$\n",
    "\n",
    "If we rearrange the terms, we get the classical Diff-in-Diff estimator.\n",
    "\n",
    "$\n",
    "\\hat{ATET} = (E[Y(1)|D=1] - E[Y(1)|D=0]) - (E[Y(0)|D=1] - E[Y(0)|D=0])\n",
    "$\n",
    "\n",
    "It gets that name because it gets the difference between the difference between treatment and control after and before the treatment. \n",
    "\n",
    "Here is what that looks in code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.524557692307688"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fl_before = data.query(\"poa==0 & jul==0\")[\"deposits\"].mean()\n",
    "\n",
    "diff_in_diff = (poa_after-poa_before)-(fl_after-fl_before)\n",
    "diff_in_diff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Diff-in-Diff is telling us that we should expect deposits to increase by R$ 6.52 per customer. Notice that the assumption that diff-in-diff makes is much more plausible than the other 2 estimators. It just assumes that the growth pattern between the 2 cities are the same. But it doesn't require them to have the same base level nor does it require the trend to be zero. \n",
    "\n",
    "To visualize what diff-in-diff is doing, we can project the growth trend from the untreated into the treated to see the counterfactual, that is, the number of deposits we should expect if there were no intervention."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApcAAAErCAYAAACLjKgCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de3hkVZ3v//d31zXp3NNJ03Y3NEq4i42jyE1AQAVvjOPlACqIeBnhKKPDDDI6oyMiyPiTUTnjFQSkf4gKas8chnNwFIFxgFF0AGmYINduoJPu3Dup1GWv80ddUlWpJJVkJ6kkn9fz1LOrdu3atap5gE+vtb5rmXMOEREREZEgeEvdABERERFZORQuRURERCQwCpciIiIiEhiFSxEREREJjMKliIiIiAQmvBRfOjg4qBJ1ERERkWWuubnZys+p51JEREREAqNwKSIiIiKBWTXhsru7e6mbICIiIrIgainnrJpwKSIiIiILT+FSRERERAKzJNXiU3HOMTIygu/7gd87Ho8zODgY+H1XO8/zaGhowGxSsZiIiIisQjUVLkdGRojFYkSj0cDvHYvFiMfjgd93tUsmk4yMjNDY2LjUTREREZEaUFPD4r7vL0iwlIUTjUYXpKdZRERElqeaCpciIiIisrzV1LB4LWhra+PQQw8tvN66dSvPPvss11xzDbfccssStkxERESk9ilclqmrq+Pee+8tOffss88uUWtEREREpjaS8knX2Ow0hUsRERGRGpLMOHoTPj1jGXaNZegZ89k1mjuOlR73ph3nHVTPBR1L3eoJNR0uW763M9D7DZy3YcZrxsbGOP744wHYb7/92Lp1a6BtEBERkdXHd46+cZ9do/nQOPnYmzv2jVffFRkPge8WsOFzUNPhcilUGhYXERERKeecYzjl6J3Uozg5PPaM+WSqDIGeQWfco7MuxLo6j8763LFu4tiZOzZFDDOju3v3wv7YWajpcFlNT2O1EolEYPcSERGRlWs84wqBsHwYunx4eqzaxAi0xozOeDYYrqvPHetCJaFxXZ1HW8wj5C3fzUlqOlyKiIiIBCHjO/aM+xM9iqMZehOVQ+NAsvrAWBcy1tXnQ+LUx466ELHQ8g2Ms6FwWaVf/epXJUsUXX/99Rx11FFL2CIREZHVzTnHUMpNDD+PVp7L2DOWDZLVzk0MG4Vh5/Jh6JLQWO/REDZtgVxG4bLMzp2Ti4he+9rX8uKLLy5Ba0RERFafsfRUw9LZ58WhMZGp/r5tMW/S3MVKcxpbYx6eAuOczRguzWwTcCOwD+AD33bOfdXM2oBbgM3A08C7nXP9lo3vXwXeBIwC73fOPbgwzRcREZHlIOM7dicqzGGscByaxbB0Q9hK5jBO6l3MneuIe0RXybD0Uqum5zIN/KVz7kEzawR+a2Z3Au8H/s05d6WZfQr4FHAJcDrQlXu8BvhG7igiIiIriHOOwaRjV/kw9GimJCz2Jnx2z2JYOuKRLXypr9DLmDuuqwvRUefRENFO1rVmxnDpnHsBeCH3fNjMtgMbgDOAk3KX3QDcRTZcngHc6JxzwH1m1mJm63P3ERERkRo3mvYnDT9PNacxWeWSjAasjXvTFL2ECoUxLVHNY1zOZjXn0sw2A0cC9wPr8oHROfeCmXXmLtsAPFf0sR25cxXDZXd3d+F5PB4nFovNpkmzouWIFsbQ0BA9PT1L3QwREZlG2kF/0tiTgj1Jyz5SZcfc872Z6oPdmpCjPepoj0xxjDraI9AacYSn6mT0gb3Zx26yD5m94ky1kLq6uqZ9v+pwaWYNwK3AXzjnhqb5G0WlN6bsCC9u4ODgIPF4vNomzUoikViwe692TU1NbNq0aambISKy6jjn6C9eXic/DF1hLuOehD/1/4zLRD1mLHrJV1DXT5kYZTF1d3fPGPoWS1Xh0swiZIPlVufcbbnTu/LD3Wa2Hsh3Xe0AipPGRuD5oBosIiKy0u1N+YVgONXSOj2jPj2JDKlZDEuXL69TfOwomsvYrGFpmYdqqsUNuBbY7pz7StFb24BzgStzx58Vnf+fZvYDsoU8g8tpvmVbWxuHHnoomUyGAw88kG984xvU19ezc+dOLr74Yh5//HF83+eNb3wjl112GdFotPDZs846i927d3PnnXcu4S8QEZFalMw4ehP5gFh5t5f8Ujsj6eqrpZuiNuMC3uvqQrTHPcLLeNcXWT6q6bk8Dngf8LCZ/T537m/Ihsofmtn5wLPAu3Lv3U52GaInyC5FdF6gLV5gxXuLf+hDH+K6667jwgsv5H3vex8f+MAHuPnmm8lkMlx00UVcdtllXHbZZQAMDAzw0EMPsWbNGp5++mk2b968hL9CREQWg19hWHqqhbz7xqvsYgTiofJh6bKwmF92Jx4iHlZglNpSTbX4vVSeRwlwSoXrHXDhPNtVE4455hj+8Ic/cPfddxOLxXjve98LQCgU4otf/CKveMUruPTSS6mvr2fbtm2cdtppdHR0cNttt/HJT35yiVsvIiJz4ZxjJO3oGZ28gPekXV/GfKrtZPQMOuOTh6ErzWlsimhYWpavmt6hp+Hck4K7FzByw11VX59Op7nzzjs59dRT2b59O1u2bCl5v6mpiY0bN/Lkk09y+OGHc+utt3LJJZfQ2dnJOeeco3ApIlJjxjOlu75UKnrJH0dnMSzdUjwsXV+2tE5ReGyLeYQ0LC2rQE2Hy6UwNjbG8ccfD2R7Lt/3vvdx7bXXVvwbpHMOM6Onp4cnn3ySY445BjMjHA7z6KOPluxFLiIiwcv4jr7yYekp5jQOzGLXl7qQFdZcnGoOY0cuOMa064tIiZoOl7PpaZxJIpGgmoWIiudc5h1yyCFs27at5NzQ0BA7d+5k//3358Ybb2RwcJAjjjgCgOHhYW677TaFSxGROXDOMZRyU85h7C0elk74ZKrMjGHLVktPGo4un9NY79EQ1rC0yFzVdLisFSeeeCKf+9znuPnmmznrrLPIZDJ85jOf4eyzz6a+vp5bb72VH//4xxx11FEAPP3007z97W/nM5/5zBK3XESkdiTSjp5E5Srp8jmNiUz1922LeRWX1ilfl7E15uEpMIosOIXLKpgZN910E3/5l3/JP/zDP+D7Pq9//ev5u7/7O5555hl27NjBq1/96sL1mzdvprGxkd/85je86lWvWsKWi4gsrIzv2J2Yeu5i8XFoFsPSDWErrYquC9EZnzynsSPuEdWwtEhNsWxx9+IaHBys+KWDg4M0NzcvyHdqh56Fs5D/3ERk8TnnGEy6itXR5cfdCR9/FsPS+WHnqRbyzs9lbIho1xeR2ViqHXqam5sn/e1OPZciIqvEaNovLNQ907qMyeqXZGRt3JtmAe/8c48WDUuLrAoKlyIiy1jaz+76UjyHsfx1PlAOpWax60vECtXQU1ZM14dYG/eIaHkdESmicCkiUmOccwwkXXZf6dGypXXKjnsSPtVGxqhXvuvL5KKXfE9jfVjD0iIyNwqXIiKLZG/KLwTDSnMYe4uCY6rKYWmDQsHLuqLCl/LQuK4uRHNUy+uIyMJTuBQRmYeU7+gtCojlS+sUL7szMotdX5qKd32Z5tge9whrWFpEaojCpYhIGd85+st3falQ9LJrzKdvvPrKl1goWy1dceHuoiHpzroQdWEFRhFZnhQuy+zatYtLL72UBx98kFgsxr777ssVV1zBAQccEMj977nnHqLRKK95zWtm/dnzzz+f7du38573vIcLL7yw6s8NDAzw4x//mA9+8IOz/s68j370o5x22mmcccYZc76HyFJyzjGSdvSMlha5VFrAu2fMp9pORs+gIz710jrFw9NNEQ1Li8jKp3BZxDnHe9/7Xs466yyuu+46AB566CF6enoCC5f33nsvDQ0NswqX6XSaPXv2cP/99/PII4/M+jsHBwe59tpr5xUuRWpVMuPKQuLUx9FZDEu3FA9Lly3cXRwe22IeIQ1Li4gUKFwWufvuuwmHw3zgAx8onDviiCNwzvG3f/u3/PznP8fMuPjii/mzP/sz7rnnHq655hpuueUWAP7qr/6KLVu28J73vIeXv/zlnHXWWdxxxx2k02muv/56YrEY3/ve9wiFQtxyyy1cddVVHHjggXziE59gx44dAFxxxRUcffTRXHHFFbz44os8++yztLe38+ijj7J7926OP/54rrrqKrq7u7n++utJJpO89KUv5Vvf+hb19fX09PTwiU98gqeffhqAr3zlK3zrW9/iqaee4vjjj+d1r3sdb3jDG6Zs95e+9CXuuOMOEokERx11FP/4j/+onhZZdL5z7En4k4ahJ6/RmKF/vPrAWBcy1tVPHoauNDwd064vIiJzUtPhcu8vTpvV9V7jAdS9+ppJn19z8h1VfX779u1s2bJl0vlt27bx8MMPc++997Jnzx5OPvlkjj322Bnv197ezt133813v/tdvv71r/P1r3+d8847j4aGBj72sY8B8MEPfpALLriAY445hueee453vOMdPPDAAwD8/ve/54477qCuro5nnnmGM888k3vvvReAgw8+mHPPPReAL3zhC3z/+9/nIx/5CJdccgnHHXccW7duJZPJMDIywmc/+1m2b99e+Ow999wzZZs//OEPc8kllxSe33HHHZx++ulV/fmJTMc5x3Bqml1fiuY09iZ8MlVmxpBNrpauOKex3qMhrGFpEZGFVtPhslbcd999vOMd7yAUCtHZ2cmxxx7Lgw8+SGNj47Sfe+tb3wrAli1b+Od//ueK19x111089thjhdfDw8MMDw8DcPrpp1NXV1fxc48++iiXX345g4ODjIyMcMoppwDZ3tdvfvObAIRCIZqbmxkYGKj6t95999187WtfY2xsjP7+fg455BCFS5lWIu3oSZRWRVeqmO4Z8xmrNjECbTFvyjmMnfk5jvXZYWnt+iIiUjtqOlxW2+MY1OcPOeQQfvazn006P9X+6+FwGN+fqBRNJBIl78diMSAb8tLpdMV7+L7PnXfeWTFErlmzZsq2XnDBBWzdupWXv/zlbN26tdArWY2p2p1IJLj44ov55S9/ycaNG7niiism/SZZHTK+Y3fCpyeR61msGBqzx8Fk9YFxTdim3OmleE5jR9wjqmFpEZFlqabD5WI74YQT+PznP88NN9xQGHJ+8MEHaWlp4Sc/+Qlnn302/f39/PrXv+ayyy4jlUrx2GOPMT4+TiKR4Fe/+hVHH330tN/R0NBQ6JkEOPnkk/nOd77Dxz/+cSBbQHTEEUfM2NaRkRH22WcfUqkUP/rRj1i/fj0AJ554Itdeey0XXHABmUyGvXv30tjYWPKdmzZtqtjufJBsb29nZGSEbdu28ba3vW12f4hSs5xzDCanGZYuOu5O+PhVZsawURh2nqpiel1diI46j4aIdn0REVnpFC6LmBk33XQTl156KVdffTXxeLywFNHevXs5/vjjMTM+//nPs27dOgDe/va3c9xxx/Gyl72sqlB4+umnc84553D77bdz1VVX8aUvfYmLL76YY489lkwmw7HHHsvVV189430+/elPc8opp7Bp0yYOPfRQRkZGALjyyiu56KKLuOmmm/A8j6985SscddRRHH300RxzzDGceuqpXHbZZRXb3dLSwrnnnsuxxx7Lvvvuy5FHHjmPP01ZLGNpV3EpnUqhcTxT/X3Xxr0pi16Kg2OLhqVFRKSITTXku5AGBwcrfung4CDNzc0L8p2JRIJ4PL4g917tFvKf22qVzg1LTxqGLhue7hnLMJSq/t/hxogVFb9MsetLfYi1cY+IltcREVk2uru76erqWvTvbW5unvQ/C/VciiwS5xwDyWwv467RiaV0Ji/m7bM74VNtZIx6TLtwd2c8Gxg74h5rNCwtIiILbMZwaWbXAW8Bepxzh+fObQG+CcSBNHCBc+4By67x8VXgTcAo8H7n3IML1XiRWrA35dObyPYqTjWHMR8cU1XuFGjkd30p61msnzyXsTmq5XVERKR2VNNzeT1wDXBj0bmrgL93zv2rmb0p9/ok4HSgK/d4DfCN3FFkWUn5jt6igFhpaZ1do9n1GIdnMSzdVLzryzTH9rhHWMPSIiKyDM0YLp1zd5vZ5vLTQFPueTPwfO75GcCNLjuR8z4zazGz9c65FwJqr8ic+c7RP57d9aW3YmicWMh7z3iVXYxALERJkctUcxo760LUhRUYRURkZZvrnMu/AP6PmX0Z8ID8djUbgOeKrtuROzdluOzu7p5oTDiMmRGNRufYrOlpzcbgJZNJ+vr66OnpWbI2jGZgT9Kyj1TZMWnsSVE4l3HVhTsPR2sE2qMu+4iUHXPP10Yda0Iw5ah0MvtIDmb/ZRAREVkoxZlqIc1UODTXcPlR4BPOuVvN7N3AtcCpZKeKlZt2zLC4gc45RkZGGB8fn2OzpjY0NERTU9PMF8qseJ7H5s2bA5/zl8y4SUUuu8Yy9Fbobdybrn5YuqV4WLps4e7iLQTbYx4hDUuLiMgysVTV4pXMNVyeC1yUe/4j4Lu55zuATUXXbWRiyHxGZjbjlopz1dPTw6ZNm2a+UBaM7xx9435RpfRUxwz949UHxrqQTVpKp9Jcxs66EDHt+iIiIrKg5hounwdOBO4CTgby/bDbgP9pZj8gW8gzqPmWK5tzjuHUNLu+5OYw9ub2nq52a+mQUdKTWHosrp72aAirWlpERKRWVLMU0c1kK8HXmtkO4LPAh4CvmlkYSAAfzl1+O9lliJ4guxTReQvQZlkE41MMS1dayHus2sQItMW8CoUvZcd6jzbt+iIiIrIsVVMtftYUb/1JhWsdcOF8GyULI+M79uSqpXvGMpNCYvFxMFl9YKwPW2HNxamHp7OLeEc1LC0iIrKiaYeeZc45x2DSTVqwu9IwdW/Cx68yM4aNwrBzpd7F4hDZoF1fREREJEfhskaNpacalp4cGscz1d+3PT8sXaHopThAtmhYWkREROZA4XIRpX3H7sTUcxh7EhPnh2YxLN0YsSkX7i4eol4b94hoeR0RERFZQAqX85Qflt5VsUq6NETuTvjTL/pZJOpR6EXsqFgxnQ2MHXGPNRqWFhERkRqhcDmF0bRf6FWcaj3G/NzGZJU7BRrQEfcm9yzWl4XGuhDNUS2vIyIiIsvPqgqXKd/RWxQMK1VJ598bTlU/LN1UvOtL0bGj7PXauEdYw9IiIiKygq2acHn27+J031v1ZkHEQlReg7HCri91YQVGEREREVhF4TJsDs/yw9KV5zDmh6c74hqWFhEREZmLVRMu/+nwcY446ABCGpYWERERWTCrpsy4IYyCpYiIiMgCWzXhUkREREQWnsKliIiIiARG4VJEREREAqNwKSIiIiKBUbgUERERkcAoXIqIiIhIYBQuRURERCQwCpciIiIiEhiFSxEREREJjMKliIiIiARG4VJEREREAqNwKSIiIiKBmTFcmtl1ZtZjZo+Unf+YmT1uZn8ws6uKzl9qZk/k3nvjQjRaRERERGpTuIprrgeuAW7MnzCz1wFnAEc458bNrDN3/lDgTOAw4CXAz83sQOdcJuiGi4iIiEjtmbHn0jl3N9BXdvqjwJXOufHcNT2582cAP3DOjTvnngKeAI4KsL0iIiIiUsOq6bms5EDgtWZ2OZAALnbO/SewAbiv6LoduXNT6u7unmMTZm8xv0tERERkMS1Wzunq6pr2/bmGyzDQChwNvBr4oZm9FLAK17r5NDAo3d3di/ZdIiIiIouplnLOXKvFdwC3uawHAB9Ymzu/qei6jcDz82uiiIiIiCwXcw2XPwVOBjCzA4EosBvYBpxpZjEz2x/oAh4IoqEiIiIiUvtmHBY3s5uBk4C1ZrYD+CxwHXBdbnmiJHCuc84BfzCzHwKPAmngQlWKi4iIiKweM4ZL59xZU7z13imuvxy4fD6NEhEREZHlSTv0iIiIiEhgFC5FREREJDAKlyIiIiISGIVLEREREQmMwqWIiIiIBEbhUkREREQCo3ApIiIiIoFRuBQRERGRwChcioiIiEhgFC5FREREJDAKlyIiIiISGIVLEREREQmMwqWIiIiIBEbhUkREREQCo3ApIiIiIoFRuBQRERGRwChcioiIiEhgFC5FREREJDAKlyIiIiISGIVLEREREQmMwqWIiIiIBEbhUkREREQCM2O4NLPrzKzHzB6p8N7FZubMbG3utZnZ18zsCTN7yMxeuRCNFhEREZHaVE3P5fXAaeUnzWwT8Hrg2aLTpwNduceHgW/Mv4kiIiIislzMGC6dc3cDfRXeuhr4a8AVnTsDuNFl3Qe0mNn6QFoqIiIiIjUvPJcPmdnbgJ3Ouf8ys+K3NgDPFb3ekTv3wlT36u7unksT5mQxv0tERERkMS1Wzunq6pr2/VmHSzOrBz4NvKHS2xXOuQrnCmZqYFC6u7sX7btEREREFlMt5Zy59Fy+DNgfyPdabgQeNLOjyPZUbiq6diPw/HwbKSIiIiLLw6yXInLOPeyc63TObXbObSYbKF/pnHsR2Aack6saPxoYdM5NOSQuIiIiIitLNUsR3Qz8B3CQme0ws/Onufx24EngCeA7wAWBtFJEREREloUZh8Wdc2fN8P7moucOuHD+zRIRERGR5Ug79IiIiIhIYBQuRURERCQwCpciIiIiEhiFSxEREREJjMKliIiIiARG4VJEREREAqNwKSIiIiKBUbgUERERkcAoXIqIiIhIYBQuRURERCQwCpciIiIiEhiFSxEREREJjMKliIiIiARG4VJEREREAqNwKSIiIiKBUbgUERERkcAoXIqIiIhIYBQuRURERCQwCpciIiIiEhiFSxEREREJjMKliIiIiARmxnBpZteZWY+ZPVJ07h/M7DEze8jMfmJmLUXvXWpmT5jZ42b2xoVquIiIiIjUnmp6Lq8HTis7dydwuHPuCOC/gUsBzOxQ4EzgsNxn/snMQoG1VkRERERq2ozh0jl3N9BXdu7/OufSuZf3ARtzz88AfuCcG3fOPQU8ARwVYHtFREREpIYFMefyA8C/5p5vAJ4rem9H7pyIiIiIrALh+XzYzD4NpIGt+VMVLnPT3aO7u3s+TZiVxfwuERERkcW0WDmnq6tr2vfnHC7N7FzgLcApzrl8gNwBbCq6bCPw/HwaGJTu7u5F+y4RERGRxVRLOWdOw+JmdhpwCfA259xo0VvbgDPNLGZm+wNdwAPzb6aIiIiILAcz9lya2c3AScBaM9sBfJZsdXgMuNPMAO5zzv25c+4PZvZD4FGyw+UXOucyC9V4EREREaktM4ZL59xZFU5fO831lwOXz6dRIiIiIrI8aYceEREREQmMwqWIiIiIBEbhUkREREQCo3ApIiIiIoFRuBQRERGRwChcioiIiEhgFC5FREREJDAKlyIiIiISGIVLEREREQmMwqWIiIiIBEbhUkREREQCo3ApIiIiIoFRuBQRERGRwChcioiIiEhgFC5FREREJDAKlyIiIiISGIVLEREREQmMwqWIiIiIBEbhUkREREQCo3ApIiIiIoFRuBQRERGRwChcioiIiEhgZgyXZnadmfWY2SNF59rM7E4z684dW3Pnzcy+ZmZPmNlDZvbKhWy8iIiIiNSWanourwdOKzv3KeDfnHNdwL/lXgOcDnTlHh8GvhFMM0VERERkOZgxXDrn7gb6yk6fAdyQe34D8KdF5290WfcBLWa2PqjGioiIiEhtm+ucy3XOuRcAcsfO3PkNwHNF1+3InRMRERGRVSAc8P2swjk33Qe6u7sDbkJtfJeIiIjIYlqsnNPV1TXt+3MNl7vMbL1z7oXcsHdP7vwOYFPRdRuB5+fTwKB0d3cv2neJiIiILKZayjlzHRbfBpybe34u8LOi8+fkqsaPBgbzw+ciIiIisvLN2HNpZjcDJwFrzWwH8FngSuCHZnY+8CzwrtzltwNvAp4ARoHzFqDNIiIiIlKjZgyXzrmzpnjrlArXOuDC+TZKRERERJYn7dAjIiIiIoFRuBQRERGRwChcioiIiEhgFC5FREREJDAKlyIiIiISGIVLEREREQmMwqWIiIiIBEbhUkREREQCo3ApIiIiIoFRuBQRERGRwChcioiIiEhgFC5FREREJDAKlyIiIiISmPBSN0BEREREZikxivX14vX14hqbl7o1JRQuRURERGpJUXC0vl6sryf3vAfrz50fHSlcnjrhTXDC25ewwaUULkVEREQWy/gYtqcHr78X29ObDYt7sqGxECKLguNUXCSCa+vEb+3A32fTIjS8egqXIiIiIkGYNjj24vX1VB8cWzvw2zpxbR3ZEJk7urYO/LYOaGgGs4kPdXcv4A+bHYVLERERkZmMj00ExL7eomHrnmCDY2sHNJYFx2VG4VJERERWt5mCY38vtnd4xtvkg2O2d7EzGyLbs0fX3rkigmM1FC5FRERk5SoEx+JexuyQdWEIW8ExUAqXIiIisjwtRHBsLZ7bOHFUcKyewqWIiIjUnvFEaUFMPiz2FRXHVBMcw7ng2K7guFjmFS7N7BPABwEHPAycB6wHfgC0AQ8C73POJefZThEREVkpZgyOvdjeoRlvMzk4Ti6QcY0tCo6LbM7h0sw2AB8HDnXOjZnZD4EzgTcBVzvnfmBm3wTOB74RSGtFRESktuWDY3kl9VyDY27pHQXHUs5PY15tDkDPt1VhoM7MUkA98AJwMnB27v0bgM+hcCkiIrL8LUZwzBXIuIZm8LxF+FG1w/kZXGoAl+zPPsb7cck+vHgH4X1OAcAfe5GxBy7Aoi3UH3PdEre4sjmHS+fcTjP7MvAsMAb8X+C3wIBzLp27bAewYd6tFBERkYVVEhwrbDm4p8rgGAqX9TJ2TKqsXk3B0TkH/jgWimdfZxKkd94+ESCTfbhkP/54P6QGyc40LOW1HlkIlxZphMwoLhVazJ8xK+bc5B9R1QfNWoFbgf8BDAA/yr3+rHPugNw1m4DbnXMvL/7s4OBg4Uu7a2hFeRERkZXIUuNEhweIDPYRHe4nMpR9RIf6iAz3Ex3qJzy2d8b7+F6IVFMryaY2Uo0tE88L51pJr2kAW/nB0fxxPH8YLzNEKDOUPfrDpCIbSNRvASA6/kfae64hGdvMns6Lcp9Lsn7nX1a8p8PwvQb8UCOZUBO+10Qm1EQ6sp6xNUflLnKYP4rz6pdsSkBXV00Nn/IAABfSSURBVFfheXNz86RGzGdY/FTgKedcL4CZ3QYcC7SYWTjXe7kReL7aBi6k7u7uRfsuERGRRZMcn1joe09P5b2qR2bR41jcy1gYti6a4+h5GBDNPdYs9O9bRM5Pg5/EwvXZ18l+Us/fkRui7ivqbeyHzFjFe4T3eT2xrncB4O+NMtaTpi5amneSoXdAuBGLtmGxViyae0RaMG9uPZK1lHPmEy6fBY42s3qyw+KnAL8Bfgm8k2zF+LnAz+bbSBERkVVpUYJjLjzmguNK45wPqeGSIej8MHSo7ZWE2/8EgPSuuxj/w5WE1r2O+GGXZD+bHiP15A2Vb+xFcqGwbSIcRlvxmg4sXGL1G6k/8aeFIfG86AEfWpgfWyPmM+fyfjP7MdnlhtLA74BvA/8b+IGZfSF37togGioiIrKiFAfH8jmO+UKZaoNj69rJe1Sv8ODo0mOAj4Wzfaf+6E7SL/68UART0svoMhXvYV4EcuHSIk2AB/7E6okWayOy77tLexdzYZLwGmyGYWmzEIRqd27kQpnznMv5KJ5zuVhqqbtYRERWuOT41MUxQQbH1g5cU+uKCY7OT+GSA6XhcDx7DO9zCqHmgwFIPvMjUn+8lsi+7yR6wAcByPQ/ROJ3f135xuGGkt5Fi2UDYqj5MEIth+W+OwOWC4TL0FLlnKDnXIqIiKw+UwbH3PP+Xmx4cMbbTATHid1iXFt2j2rXvnKCY3ZYeig7jJzrZcwM/TfpXXeVVkyP90N66h13vDX7FcKlRZvBi+KKeiStfgOR/d83KUBapAULRWds51znOspkCpciIiJ5k4Jj8TqOCo55zrnscjhFPYvFvY2Rfd+Ft2YTAMnub5LesY1o10eIbHp79vOjO0k/d9vkG5uXDYPF4TD3CLUcXrgsvM+phPd5fcmwtBdrJ7r/exb2h0tVFC5FRGR1qBAciwtkvL6e2QfH/LqNrcUhcvkGR+cncckBLNxQqJjO9D1IuvfXhd7FwjxGf3zK+4TWHlMIl9n5iY3g/ML7XlMXkZedP1EEE8vOZSTSWNWwtK2CpY6WM4VLERFZ/pLjWP/ukt1i5hYcQxNzHFdIcHQuA6kh/AqFLi7ZT6zrz7FoCwDjf7iSTO+viR3+acKdrwUgM/xH0jv/ZfKNvVhR72IbFm0pPPca9i9cFtnvTKKbzyr9aP1Govu9a+F+tCwphUsREaltUwXH4gKZ4YEZbzMpOFbaq3oZBUeX3osb78PiHYWlbtIv/pJM34Nlu78MAv6U9/E3vYNQLlxarAOLrS2prg61HYmFLphULW3huqraOVNFtaw8CpciIrJ0UslJvYyFParnEhxby+Y4LrPg6DLjk3oWixfvjh3+mULhSeL3f4M/9DjxP/kKoeZDAcgMPU76xTsn3zjSVFotHW3LDUW34sU7CpfFDvwoHPjRko+GGg8g1HjAwv1oWXEULkVEZGFMGxxzBTLzDY75oeumFqjhal9/vA+X7MNbs38hHKae+ymZwUdLq6Uzo9Pex6UGsVgbAFa3HksNg58uvB9edwJew+aJHsZYKxZpzq7nKLJIFC5FRGT2Ukmsf3eF9RvnEBxb1lK+zaDf1lnTwdE5B+mRir2LLtmP81PED7+0cH3iNx/Hje+m7tgbsXgnAJmBR8j03lt6YwsXLaUzefeXfJENQPywT01qV6j50EIvpshSUbgUEZFSMwXH/l68of4Zb+M8r+L+1IXg2NaBa26tqeDonMMlXsQlBwtrKgIkn7wBf+SpkqV3cKmpb2Qezl1SqGr21uyLizRCZqLCOrLhLYQ7jisJkoQbNEdRlj2FSxGR1aQQHHMV1Lk9q7N7V88hOLZ24LdPDE9nC2U6ayo4Oj+NSw0UbQs4ULIuo4UbiB18Ue5qn7H/OB/wqT/pXzAv+7/JTN/v8Ye2l944VF+5dzE3HA0Tm9HFt3xxUrtCbVsW5geLLDGFSxGRlWK64JhfEHwFBUfnp3GjO3CZREkv4/j2q/ETPRPrMaamX4LIYhMFLWYhvMaXgXmQGQOvEYDI/u8Bf7x02Z1chbaIlFK4FBFZDlJJbGBPUVgs2qN6rsExv25j2S4ySxkcXSZRNn+xr2Qo2luzqbCXtEsNMvbAn2PRVuqPv7lwj0zf73DjPUV3taI1GEt7F7PP20vaUPfqr09qV7j9VQvye0VWIoVLEZGlNlNw7O/FG6wyOLasLVqCp3aCo8sk8Pc+CzhCTQcVzo0/+uWSoWoyY9Pfp6gX0iItWP0mLNaOc64wVzF60AWlhTGRZu0bLbKIFC5FRBZSPjiW71GdL5SZd3DMhcbWDlxL26IFR+d8SA0Vehj98mrpZD+h1i1EN58JgL/3GRK/uQiv8QDqXn1N9iZelMzu/yhZsBsvOjGHMVY0BJ3vbcxVWgOYF6L+6O9Malt47dEL+ttFZHoKlyIic5VOlc5xnHdwLN1mcMmCY2oYf+RpCMUINR0IgJ/oJfn41yaGqFP9JXtFV2LhhonnsXa8hpfhrdlv4px52UXBww2FIWpC9aqWFlnmFC5FRCopCY5FYbGwGHjP/INjvkBmgYOj85MTFdKTehizi3eH159KZMObAcgMPMz4w58n1P4aQq/4ewDMC5PZ85+lNw43lPUyllZOe3X7FC71YmupO+p/TWpbuOOYBfvdIrI0FC5FZPWZMTj2YkP9mHPT3saZh2ttn9irur2ztMI6P8cxtHD/qfXH9+CPPI1FWwk1vhSAzPAfSXZ/uxAeSY/MfJ9cDyWAxTvxmg/FW7PvxAWRJmJHfK4oQDZjXjTw3yMiy5/CpYisLPng2N+Ll1u3sXwx8GqDo9+2dprg2IFrbgs0ODrnIDM6Ze9i/nVk85mEO18LQKb330n+9z8R3vBmQgd9LHcnH3/gvyZubKGySun8DjBF6zPWrS9cHmo8gLo/+UpJ28xCmssoIlVRuBSR5SPo4Fhpr+q2hQmOAP7o8/h7n8ar31joFUzv+S2pp74/sSajn5z5PmMvFp579RvxWrfg1W0oOreJ+JYvFoaoiTQWdooREVloCpciUhvS6WxY7C/en7q0UKaWgqPzM9ldXwp7S/dP7m1M9hM76CJCrS/P/sQX/g+pZ24hsv85RPc/O3ejFP7QYxM3DtXlehNbSnd+KdoJxivuZWx7JXVtryxpm4XihMrOiYgsFoVLEVl46TQ2ULxX9fyDY3GBTHGhTFA9jpnhP+L2PovXfEihMCX1ws9JP3cr/nh+15fp2wuULObtNbyUUPtrSsNh82HEX/n/TYTIcN282y4ispQULkVkfgrBsWy3mNzWg9bfiw32VRccc8UxQQdHl0nkqqX7ynoZ+0p6HuNbLsdbsyn7s577CekXf0704E9OVD374/gjT+XuahBpLhS4eLGy3V/yvY3xdYV2hNedSHjdiSVts0gjoZbDZv2bRERq1bzCpZm1AN8FDif7V/gPAI8DtwCbgaeBdzvnZl6vQ0Rqz0IFx5ICmaJ1HOfY45gZeAR/dCeh9lfjxdoASD7zI9Iv3IEb74fMaFX3cck+yIVLr/kQQv44FpvYGjDccSxe04HZ8Bhpxjz9/VxEpNx8/8v4VeAO59w7zSwK1AN/A/ybc+5KM/sU8Cngknl+j4gErTg4FoXFiWHrnlkGx44KldWzC47OOUgPV6yW9sf7S3of64+5DguvASD51E34/b8n9orLC+GSzChudGf2uUXK1mIs2/klXzUdW1toS2TDmwvrPuZZtJVQtHUWf8giIqvPnMOlmTUBJwDvB3DOJYGkmZ0BnJS77AbgLhQuRRZXOo0NFu1VHWRwLCmQ6Zx1j2N69/24xC7C+7y+ML9w/L+/Qab333HJAXDpqu7jkgOFcBlqeyVebC0WaSq8H97wZsLrXpetlg43aNcXEZFFYm6G/7lM+UGzLcC3gUeBVwC/BS4CdjrnWoqu63fOlfxVf3BwsPCl3d3dc/p+kVUrkyYyMkh0qJ/IUD/RoT4iw/nn/dnnw4PYDMUmDiPV2EyqqZVkYyupprbs86ZWUo2tJJvaSDU0TR0cXRovM0zIH8bLDBHKDOP5Q7nnuaM/jPnj7NpweeFjHS98gUh6Fz37/A3pSLawpaVvK/V77wPAtzoyoUb8UBO+10Qm1IQfasweC6+b8L0G0PI6IiKLrqurq/C8ubl50t/c5xMuXwXcBxznnLvfzL4KDAEfm024XCzd3d0lfxgiNSmTxgb2TLFXde4x2IfNsKezM8M1t+NyC35XXJKnuR3CpcHROVfo4XPOkdn1C1xygPCmPyucTzx8GZmBhyE1VPXPqj/xp1goDkDyyRtwqSEi+76rUCjjJ3rB+dmh6ZB2fRERma2lyjmVwuV85lzuAHY45+7Pvf4x2fmVu8xsvXPuBTNbD/RMeQeR1aQkOE4s+u3lth60PdUHRz+3V3Vhj+ryyuqi4Jjd9WWsdP3F0T/iBiZXSxNppP413wTAzBh//H9BZpTw+jdApDH3O8ZywdLLbgE4xVqMxbvA4MUK7Y++9NxJv8mLdwTzZywiIktuzuHSOfeimT1nZgc55x4HTiE7RP4ocC5wZe74s0BaKlLLMmlsoG/SbjELGRydy2AWyn4ukyD94i/Af4JI++GF+43958fw9z4L/nh1v6PsuvBL3jjpkujBn8hWSUeaCt8vIiKSN99q8Y8BW3OV4k8C5wEe8EMzOx94FnjXPL9DZGmVB8dcgYzXX7Tl4MAcg2OusjpfYe03tYIbLaqOzj+6cckHcM9M9D56jQdQd+SVuZv7JB//GngxwhvPmBja9pPZwOjFsj2KsQo7v5TsMd1S0uZY10cm/Q71MoqIyHTmFS6dc78HXlXhrVPmc1+RRZMPjv29ePl1G+ccHIsXAJ8IjpmWtbi2dmhdl+1xTA6Q3nUXeJHCUjfOTzP263NwqQGY4bsK35kcmHgRqiP8ktOxSHP287kexfiWy7FQfXZLQVVLi4jIItAKwLJyTRscc+FxLsExv/VgayuZ5jr8hhB+1OEyg2X7TP8x+/yFfsKZE4h1XJy9X2qYZPc3sboNhXBpXjjby+j87HBztHQ9xkq7vxTmQJKdHxk7+KJJbfeK1m0UERFZDAqXsjxl0tk5jGV7VJfsVT2b4JjvZWxdi2trg7b1+G2dZNY4UuN/wOrWEln/eiBb2Tz2wEdhbATGqmuuSw0XnlusnfCGt+LVrSu5pu4138YijZgXmd2fhYiISA1RuJTaM1Nw7O/F+vdUHRz91rX4He1k2hvJtMTxGyL4dYYfzeC8BC41WDSv8TdE9nsn0ZdlZ3b4fb8j9fgNeC1HFMKlRRohPQIWqryXdHRyL2N+sXAAC9cTO+jCSe0t7CwjIiKyjClcyuLyM4XleKZckmdgD+bPHBwzra24tg5oWYff1kmq3SO1Zg/W9FLCG0/FtbST2ftHEr/7K/CfmfhwMveY6t6ZROG51W8ksu87sDX7TZwLxal/7Q9zu75oEW8REZFiCpcSnOmCY75AZobg6AwydUZmbQv+2iYyrfVkmmL49R4u5uOHkviM4TJDkBkleuBbiGx8GwCpF35OcvuXCUVihNZmF+e2SCP4KQjVl63FWLYeY359xkhzdpmdHC/eQfSAD01qZ/E2gyIiIjJB4VKq42cKy/FkC2QqLMkzRXB0gIuBpQFn+M1tjO/fSPIlRtjbRHjNYbjWDsbrXiTRd1PuE4nco0gm98izMC4zsS5jqPEAIpvPxmuc2KHA6vYp2R1GREREFpbCpVQRHHuxgd2TgqML53oZ4+A3GH6HkWltINMUx18Two+DH0nhLAH4xPf9GKHNb4BwhOTTPyD15PVE9n0t7oAzszfs/y/oAyItRdXRldZjzK7TmB2Wnlhex2vYTLRhc0kbzUIQ0kLfIiIii0XhcqXLB8f+yXtUFwpkioKjM/DrwEuA5bJkYl+P5AEesd5GwuyDa+1gdN+9jDVvr/CFaWBk8unwGlxDPYSzldCh1i3w0vcTapnYTcZrPpz6k/4F8xQGRUREliuFy+WsyuCI7+Ni4NdZtqexLvvcbzH8fQy/LkRmTQy/znDRbKJsSLyZUOuh+G0dJMbuINX3CzJvO4/Ipj/NfvWLv4DHnijqWcz3KBbv9lL0CMVKmh5qPphQ88El5xQqRURElj+Fy1rlZ7DB/tJK6qJCGfp2wcgeXMwRGnGYy35stCtEaq2xZneG0ED25NBxDYwdkJ7hC13u4WHRFtJHnohrPQKAUJ+PtR+I13JE4erQuhOpX/c67foiIiIiJRQul8I0wdH6e3CjPbhUPy7qF3oaU3WW7W1sN/yNkIkbRKIAtP7HS/AaXpLtZdznEVKh5+Do8whvOBnXupb0zp/CM7dMsZd08S4wbRBpzM5TLBJq20KobUvJufJrREREREDhMniF4Jgbqu7vhT27YOgFGO4h8vxgdo5jJsPIy0NkWjwafpMilNvpZeCECOP7h6jqH40XxaJtJP7i7/HWbMqe2n0/0fE9WNuRuLr1AET2fSfR/d61QD9YREREZILC5WyUBUfr34kb2IHbuwvGduOnB3HsJRNz2V7G/KMDWGfgOzofGccc+E2tjB+QId2UJNJ4HKGGA3Ftnfje/dj4w1i8vah3sbxiOnskVD9pWDq89jWTmq2haxEREVksCpd5ueDInueh7xls4AXCfZlCocze9c/ihxI0/ypZqKLuOy1Kap0Ha4tvVHm42Ihj0Wb2XvG3sHZfiEQJ9dyLl0mQOf5V+NEWAKKcSnRhf6mIiIjIglk14TIytBO3/RkYfA6GX8CN7cZP9uMywzgbw4+k8ePgYtlePm/c0fGvEwt0j2/J7hKT7mjGi6/DtXdgDS/g+cNYqDHby1jfCY3rsbq1Zb2NLZgXmdSmcOfxi/b7RURERBbDqgmX657/EqNrciXVIaCh/Irc0LEPXiaChetJvOut0LYOv62DiPcUNDSTOOFoLFwPQCT3EBEREZGsVRMuSdYTckk8F8e8xmxvYl0HNOwDLZuw1v2wNZ25amkPyC4HnhdmS+X7ioiIiEjBqgmXO464kq6urpkvFBEREZE585a6ASIiIiKycihcioiIiEhgFC5FREREJDAKlyIiIiISmHmHSzMLmdnvzOxfcq/3N7P7zazbzG4xM60JLiIiIrJKBNFzeRGwvej1l4CrnXNdQD9wfgDfISIiIiLLwLzCpZltBN4MfDf32oCTgR/nLrkB+NP5fIeIiIiILB/zXefyH4G/Bhpzr9uBAedcfv3xHcCG6W7Q3d09zyZUbzG/S0RERGQxLVbOmWnd8DmHSzN7C9DjnPutmZ2UP13hUjfdfRZrYfPu7m4toi4iIiIrUi3lnPn0XB4HvM3M3gTEgSayPZktZhbO9V5uBJ6ffzPnr1b+wEVERESCVks5Z85zLp1zlzrnNjrnNgNnAr9wzr0H+CXwztxl5wI/m3crRURERGRZWIh1Li8BPmlmT5Cdg3ntAnyHiIiIiNQgc27aKZEiIiIiIlXTDj0iIiIiEpgVES7NzJnZ94teh82sN79rkIiIiMhKYWYjM7x/l5m9arHaU25FhEtgL3C4mdXlXr8e2LmE7RERERFZlVZKuAT4V7K7BQGcBdycf8PMjjKzX+f2QP+1mR2UO3+PmW0puu7fzeyIRW21iIiIyCyZ2UnFI7Rmdo2ZvX8Jm1SwksLlD4AzzSwOHAHcX/TeY8AJzrkjgb8Dvpg7/13g/QBmdiAQc849tGgtFhEREVlhVky4zIXCzWR7LW8ve7sZ+JGZPQJcDRyWO/8j4C1mFgE+AFy/KI0VERERWaFWTLjM2QZ8maIh8ZzLgF865w4H3kp2RyGcc6PAncAZwLuB/3/xmioiIiIyZ2lKc1x8qRpSbj7bP9ai64BB59zDRfudQ7bnMl/g8/6yz3wX+GfgHudc34K3UERERGT+ngEONbMY2WB5CnDv0jYpa0X1XDrndjjnvlrhrauAK8zs34FQ2Wd+CwwB31uEJoqIiIjMmZmFgXHn3HPAD4GHgK3A75a0YUVW/Q49ZvYS4C7gYOecv8TNEREREZmSmb0C+I5z7qilbstUVlTP5WyZ2Tlkq8o/rWApIiIitczM/pxsXclnlrot01n1PZciIiIiEpxV3XMpIiIiIsFSuBQRERGRwChcioiIiEhgFC5FREREJDAKlyIiIiISGIVLEREREQnM/wPzU/wQ6wjIYgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,5))\n",
    "plt.plot([\"May\", \"Jul\"], [fl_before, fl_after], label=\"FL\", lw=2)\n",
    "plt.plot([\"May\", \"Jul\"], [poa_before, poa_after], label=\"POA\", lw=2)\n",
    "\n",
    "plt.plot([\"May\", \"Jul\"], [poa_before, poa_before+(fl_after-fl_before)],\n",
    "         label=\"Counterfactual\", lw=2, color=\"C2\", ls=\"-.\")\n",
    "\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See that small difference between the red and the yellow dashed lines? If you really focus you can see the small treatment effect on Porto Alegre. \n",
    "\n",
    "![img](./data/img/diff-in-diff/cant-read.png)\n",
    "\n",
    "\n",
    "Now, what you might be asking yourself is \"how much can I trust this estimator? It is my right to have standard errors reported to me!\". Which makes sense, since estimators without them look silly. To do so, we will use a neat trick that uses regression. Specifically, we will estimate the following model\n",
    "\n",
    "$\n",
    "Y_i = \\beta_0 + \\beta_1 POA_i + \\beta_2 Jul_i + \\beta_3 POA_i*Jul_i + e_i\n",
    "$\n",
    "\n",
    "Noice that \\\\(\\beta_0\\\\) is the baseline of the control. In our case, is the level of deposits in Florianopolis in the month of May. If we turn on the treated city dummy, we get \\\\(\\beta_1\\\\). So \\\\(\\beta_0 + \\beta_1\\\\) is the baseline of Porto Alegre in May, before the intervention, and \\\\(\\beta_1\\\\) is the increase of Porto Alegre baseline on top of Florianopolis. If we turn the POA dummy off and turn the July Dummy on, we get \\\\(\\beta_0 + \\beta_2\\\\), which is the level of Florianópolis in July, after the intervention period. \\\\(\\beta_2\\\\) is then the trend of the control, since we add it on top of the baseline to get the level of the control at the period post intervention. As a recap, \\\\(\\beta_1\\\\) is the increment from going from the treated to the control, \\\\(\\beta_2\\\\) is the increment from going from the period before to the period after the intervention. Finally, if we turn both dummies on, we get \\\\(\\beta_3\\\\). \\\\(\\beta_0 + \\beta_1 + \\beta_2 + \\beta_3\\\\) is the level in Porto Alegre after the intervention. So \\\\(\\beta_3\\\\) is the incremental impact when you go from May to July and from Florianopolis to POA. In other words, it is the Difference in Difference estimator. \n",
    "\n",
    "If you don't believe me, check for yourself. And also notice how we get our standard errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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>  171.6423</td> <td>    2.363</td> <td>   72.625</td> <td> 0.000</td> <td>  167.009</td> <td>  176.276</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>poa</th>       <td> -125.6263</td> <td>    4.484</td> <td>  -28.015</td> <td> 0.000</td> <td> -134.418</td> <td> -116.835</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>jul</th>       <td>   34.5232</td> <td>    3.036</td> <td>   11.372</td> <td> 0.000</td> <td>   28.571</td> <td>   40.475</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>poa:jul</th>   <td>    6.5246</td> <td>    5.729</td> <td>    1.139</td> <td> 0.255</td> <td>   -4.706</td> <td>   17.755</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "smf.ols('deposits ~ poa*jul', data=data).fit().summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non Parallel Trends\n",
    "\n",
    "One obvious problem with Diff-in-Diff is failure to satisfy the parallel trend assumption. If the growth trend from the treated is different from the trend of the control, diff-in-diff will be biased. This is a common problem with non-random data, where the decision to treat a region is based on its potential to respond well to the treatment, or when the treatment is targeted at regions that are not performing very well. Take our marketing example. If we decided to test billboards in Porto Alegre not to check the effect of billboards, but simply because it is performing poorly. Maybe because online marketing is not working there. In this case, It could be that the growth we would see in Porto Alegre without a billboard would be lower than the growth we observe in other cities. This would cause us to underestimate the effect of the billboard there. \n",
    "\n",
    "One way to check if this is happening is to plot the trend using past periods. For example, let's suppose POA had a small decreasing trend but Florianopolis was on a steep ascent. In this case, showing periods from before would reveal those trends and we would know Diff-in-Diff is not a reliable estimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApcAAAErCAYAAACLjKgCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxc9X3v/9d3du2bF9mSsVlk8IJZvIAXwNjOQpqUmybtDdnI1rRJmuampZdwk5T8SlMS2ps0Db8maQIFGhpICL2Q24QGG4zBNgZsCIsNiMUYyZZk2dql2b/3jzMzmpFGliyNRiPp/Xw89NDonDMz3zHH1pvv9jHWWkREREREcsE11Q0QERERkZlD4VJEREREckbhUkRERERyRuFSRERERHJG4VJEREREcsYzFW/a1dWlJeoiIiIi01xFRYUZekw9lyIiIiKSMwqXIiIiIpIzsyZcNjY2TnUTRPJO973MRrrvZTYqpPt+1oRLEREREZl8CpciIiIikjNTslp8JNZaent7icfjOX/tQCBAV1dXzl93tnO5XJSWlmLMsMViIiIiMgsVVLjs7e3F7/fj8/ly/tp+v59AIJDz153twuEwvb29lJWVTXVTREREpAAU1LB4PB6flGApk8fn801KT7OIiIhMTwUVLkVERERkeiuoYfFCUF1dzfLly1M/33333Rw5coRbb72Ve++9dwpbJiIiIlL4FC6HKCoq4oknnsg4duTIkSlqjYiIiMjIrLVECmx2msKliIiIyDTSHY6z61iIR5pDbG8O8oEzi/hwxVS3alBBh8vKf23O6et1frJu1GsGBgbYtGkTAIsXL+buu+/OaRtEREREToe1lhdORtiRCJP7WsNE7eD5p4+HFS4LWbZhcREREZF8OhGM8ejRENubgjx6NETrwODYt8vAJfN8bK3zs60uwAU1Xt54vXMKW5upoMPlWHoaxyoYDObstURERERyKRq3HGgPs705xI6mIAfaI6R1TrKw2MXWugDb6gNcscBPpb9wN/wp6HApIiIiMlMd7YuxoznIjuYQO48G6QwPxkmfCzbU+tm60M/W+gDLKj3TphqewuUYPfbYYxlbFN1xxx2sW7duClskIiIi00koZnmyNZSaO3mwI5px/qwyN1vrA2yrC7Cp1keJt3B7J09F4XKI5ubhi4guu+wyWlpapqA1IiIiMp292R1le3OQ7c0hHj8Woj9tJU6Jx3DZAj9b6/xsrQtwVvnMiGWjfgpjzCLgLqAWiAP/Yq39njGmGrgXWAIcBv7IWtthnD7b7wHvAfqBT1hrD0xO80VEREQKR18kzuMtTu/kjqYgb/TEMs4vr/KwrS7A1roAl8734XdPj6Hu0zGWiBwF/tJae8AYUwbsN8Y8DHwC2GGt/ZYx5ivAV4DrgauAhsTXJcAPEt9FREREZhRrLYc6o+xocnon97aGCKdtal7hM2xZGGBrvZ8tCwMsLHFPXWPzZNRwaa09BhxLPO4xxhwC6oCrgc2Jy+4EduKEy6uBu6y1FnjSGFNpjFmQeB0RERGRaa0zFGfnUWfe5CPNQY72D6ZJA6ye403MnfRz8RwfHtfM6508ldMa3DfGLAEuAvYB85OB0Vp7zBgzL3FZHfB22tOaEseyhsvGxsbU40AggN/vP50mnRZtRzQ5uru7aWtrm+pmyAjS/46JzBa67yWX4hYO9brY2+Fib4ebF3tcxBkMjNVey/qqGBuqYqyrjFHpBeiCLnizK3/tzNd939DQcMrzYw6XxphS4JfA/7DWdp9iOXy2EzbLsWEN7OrqIhAIjLVJpyUYDE7aa8925eXlLFq0aKqbIVk0NjaO+o+AyEyj+15yobU/xiNHQ+xoDvJoc4gTocHeSY+BDfN9bKsLsKXOz8pqL64p3iaokO77MYVLY4wXJ1jeba29P3G4NTncbYxZACS7rpqA9KRRDxzNVYNFREREci0St+xrC/NIc5DtTSGePxnJOH9GqTuxEMfPZQv8lPum5zZB+TCW1eIGuA04ZK39TtqpB4FrgW8lvj+QdvzPjDH34Czk6ZpO8y2rq6tZvnw5sViMpUuX8oMf/IDi4mKam5u57rrreOWVV4jH47zrXe/ipptuwufzpZ57zTXX0N7ezsMPPzyFn0BERETG4q2eKI8k9pzcdSxET2RwoDXghstqnQ3Mt9b5Oad8+mxiPtXGErs3Ah8Dthhjnkt8vQcnVL7DGNMIvCPxM8CvgTeA14AfA5/PfbMnT7K2+N69e/H5fNx+++1Ya/nYxz7G7/3e73HgwAH2799PX18fN910U+p5nZ2dPP/883R1dXH48OGp+wAiIiKS1UDUsr0pyA37Oll3fysX3NfKl/d28p9HgvRELOdVevj8ihLuf2cNb354Ib945xz+dHkpDRVeBcvTMJbV4k+QfR4lwNYs11vgCxNsV0FYv349L730Ert27cLv9/PRj34UALfbzd/93d9xwQUXcMMNN1BcXMyDDz7Iu9/9bubOncv999/PX/zFX0xx60VERGY3ay2NXVGnXndzkN0tIYJp206Wew1XLPSn5k4uKp0Zm5hPtYL+Uyy9dnPuXgvovXPnmK+PRqM8/PDDbNu2jUOHDnHhhRdmnC8vL6e+vp433niDlStX8stf/pLrr7+eefPm8fGPf1zhUkREZAp0h+M8dszZwHzH0RBv92ZuYn5hjTdVEWftPB/eWbZNUD4UdLicCgMDA2zatAlwei4/9rGPcdttt2XtDrfWYoyhra2NN954g/Xr12OMwePxcPDgwYxa5CIiIpJ7cWt54WTEqdfdFOSptjBpFRap8bucMFkf4MqFfuYVzfxNzKdaQYfL0+lpHE0wGGQsGxEl51ymW7ZsGQ8++GDGse7ubpqbmznzzDO566676OrqYtWqVQD09PRw//33K1yKiIhMgvZgjEebk5uYhzgeHNwmyG1g/XwfWxMruy+omfptgmabgg6XheKKK67gG9/4Bj/72c+45ppriMVifO1rX+PDH/4wxcXF/PKXv+S+++5j3bp1ABw+fJj3v//9fO1rX5vilouIiEx/0bhl//Fwau7ks+2RjA2064rdbK13hrqvWOCn0q9tgqaSwuUYGGP46U9/yl/+5V/y93//98Tjcd7xjnfw13/917z11ls0NTWxdu3a1PVLliyhrKyMZ555hjVr1kxhy0VERKan5r4YO5qD7GgOsvNoiK7wYJz0uWBjrT81d/K8Sm0TVEgULodobm7Oery+vp5777132PHFixdz6NChYcd37dqV87aJiIjMVKGYZW9riO1NIR5pDnKwM5px/pxyD1vqnJXdG2t9lHjVO1moFC5FRERkSrzRHWV7k9M7+XhLmP60lTglHsPlC/xsSwx3LylTZJku9F9KRERE8qI3EufxY6FUVZw3ezK3CVpZ7WXrQmdl96XzfPjcGuqejhQuRUREZFJYaznYEWVHc5DtzSH2toaIDC7sptJn2JJY1b2lLsCCYm0TNBMoXIqIiEjOdITi7DzqhMlHmoMc6x9MkwZYO9fL1roA2+oDXFTjxa1NzGcchUsREREZt1jc8tyJCNubg+xoCvFMe5h42j5B84tcTpis87N5oZ/qgHonZzqFSxERETktrf3JbYJCPHI0SEdoME16XbBxvo9t9QG21gVYUaVtgmYbhcshWltbueGGGzhw4AB+v58zzjiDm2++mXPOOScnr//444/j8/m45JJLTvu5n/70pzl06BAf+chH+MIXvjDm53V2dnLffffxmc985rTfM+lzn/sc7373u7n66qvH/RoiIjI9hWOWfW1hHknMnXzhZCTj/OJSN++oD7Clzs9lC/yUaZugWU3hMo21lo9+9KNcc8013H777QA8//zztLW15SxcPvHEE5SWlp5WuIxGo5w4cYJ9+/bx4osvnvZ7dnV1cdttt00oXIqIyOxyuCeaWtW962iI3rRtgorchssW+BLD3QHOKnerd1JSFC7T7Nq1C4/Hw6c+9anUsVWrVmGt5etf/zrbt2/HGMN1113HH/zBH/D4449z6623pjZX/6u/+isuvPBCPvKRj3D++edzzTXX8NBDDxGNRrnjjjvw+/3867/+K263m3vvvZdbbrmFpUuX8uUvf5mmpiYAbr75Zi699FJuvvlmWlpaOHLkCDU1NRw8eJD29nY2bdrELbfcQmNjI3fccQfhcJizzjqLH/3oRxQXF9PW1saXv/xlDh8+DMB3vvMdfvSjH/Hmm2+yadMmrrzySt75zneO2O5vf/vbPPTQQwSDQdatW8c//uM/6h8MEZFZoD8aZ3dLmO1NQR45GqKxK3MT82WVnlS97vXz/QQ8+t0g2RV0uOx75N2ndb2r7ByK1t467PklWx4a0/MPHTrEhRdeOOz4gw8+yAsvvMATTzzBiRMn2LJlCxs2bBj19Wpqati1axc/+clP+P73v8/3v/99PvnJT1JaWsoXv/hFAD7zmc/w+c9/nvXr1/P222/zgQ98gKeeegqA5557joceeoiioiLeeustPvShD/HEE08AcN5553HttdcC8Ld/+7f827/9G3/yJ3/C9ddfz8aNG7n77ruJxWL09vZy4403cujQodRzH3/88RHb/NnPfpbrr78+9fihhx7iqquuGtOfn4iITB/WWl7tijr1upuC7G4NEUrbdrLcZ9i8wM+2+gBbFvqpLy3oyCAFRHfKGDz55JN84AMfwO12M2/ePDZs2MCBAwcoKys75fPe9773AXDhhRfyq1/9Kus1O3fu5OWXX0793NPTQ09PDwBXXXUVRUVFWZ938OBBvvnNb9LV1UVvby9bt24FnN7XH/7whwC43W4qKiro7Owc82fdtWsX//RP/8TAwAAdHR0sW7ZM4VJEZIboCsd57GgotRinqS9zE/OL5nhTK7vXzPXh0TZBMg4FHS7H2uOYq+cvW7aMBx54YNhxa22Wq8Hj8RCPD+7fFQwGM877/X7ACXnRaObwQlI8Hufhhx/OGiJLSkpGbOvnP/957r77bs4//3zuvvvuVK/kWIzU7mAwyHXXXcejjz5KfX09N99887DPJCIi00fcWp4/EWFHYu7kU21hYmm/0uYEXKl63Vcu9DO3SNsEycRpOVeayy+/nHA4zJ133pk6duDAASorK/mP//gPYrEY7e3t7Nmzh9WrV7No0SJefvllQqEQXV1dPPbYY6O+R2lpaapnEmDLli38+Mc/Tv38/PPPj6mtvb291NbWEolE+MUvfpE6fsUVV3DbbbcBEIvF6O7upqysLOM9R2p3MkjW1NTQ29vLgw8+OKa2iIhI4WgPxvj56/18dtdJlt7TwuZfHeemA93sbQ0DsH6+j69fXM7O983l1Q/V8i+XV/NHZxcrWErOFHTPZb4ZY/jpT3/KDTfcwHe/+10CgUBqK6K+vj42bdqEMYa/+Zu/Yf78+QC8//3vZ+PGjZx99tmsWrVq1Pe46qqr+PjHP86vf/1rbrnlFr797W9z3XXXsWHDBmKxGBs2bOC73/3uqK/z1a9+la1bt7Jo0SKWL19Ob28vAN/61rf40pe+xE9/+lNcLhff+c53WLduHZdeeinr169n27Zt3HTTTVnbXVlZybXXXsuGDRs444wzuOiiiybwpykiIvkQjVuePh5mR7Mz3P1ce4T08bb6EjfbEuUVr1jop8KnfiWZXGakId/J1NXVlfVNu7q6qKiomJT3DAaDBAKBSXnt2W4y/7vJxDQ2NtLQ0DDVzRDJq9lw3zf1RnnkaIjtTUF2HgvRHR78tep3w8b5frbWO3Mnl1ZoE/PZYKru+4qKimE3l3ouRUREClwwatnbGkrV6z7UmTmPv6HCw9Y6P1vrAmys9VHsUe+kTJ1Rw6Ux5nbgvUCbtXZl4tiFwA+BABAFPm+tfco4/2v0PeA9QD/wCWvtgclqvIiIyExkreWN7phTr7s5yOPHwgykrcQp9RguX+gsxNlS52dJmfqKpHCM5W68A7gVuCvt2C3A/2et/Y0x5j2JnzcDVwENia9LgB8kvouIiMgp9ETiPH4slJo7ebgnc5ug86u9bKtzhrvXzfXhc2uoWwrTqOHSWrvLGLNk6GGgPPG4AjiaeHw1cJd1JnI+aYypNMYssNYey1F7RUREZgRrLS92RJ163U1BnmwLExncJY5qv7NN0JaFzmKc2mKt5pbpYbz96P8D+C9jzD/gbGeULFdTB7yddl1T4tiI4bKxsXGwMR5n0rHP5xtns05NezbmXjgc5uTJk7S1tU11U2QE6X/HRGaLQr3vuyKwr9PN3g43T3a6aA8Pzo10YTm/LM76qhgbquKcVxrHbZydQHqaoWekFxVJyNd9P9rCofGGy88BX7bW/tIY80fAbcA2IFsf/SmXo6c30FpLb28voVBonM0aWXd3N+Xl5aNfKKfF5XKxZMkSrUQsULNh1azIUIV038filgPtkURFnCD72yPE034rLih2JSriONsEVfm1EEfGp5Du+/GGy2uBLyUe/wL4SeJxE7Ao7bp6BofMR2WMGbWk4ni1tbWxaNGi0S8UERGZgJb+WKq84qNHg3SEBtOk1wWbav3O3Mm6AMurtE2QzDzjDZdHgSuAncAWINkP+yDwZ8aYe3AW8nRpvqWIiMxk4ZjlybYwO5qC7Dga4sWTkYzzS8rcvKMuwNZ6P5tq/ZR61TspM9tYtiL6Gc5K8DnGmCbgRuCPge8ZYzxAEPhs4vJf42xD9BrOVkSfnIQ2i4iITKnDPVF2NAfZ3hTi8WMheqODvZPFHsNltT5nuLs+wFnl2iZIZpexrBa/ZoRTq7Nca4EvTLRRIiIihaQ/GueJY+HUvpOvd2duE7S80pOqiHPpfD9+bRMks5j+d0pERGQIay2vdEXZ3uTMndzTGiKUlicrfIYrFzobmG+tC1BXom2CRJIULkVERIDOUJzHjjkbmO9oCtHcP5gmDXDxHG9iZbef1XN9eFzqnRTJRuFSRERmpbi1/O5EhO1NQR45GuKptjBpFRaZG3Cxtc7PtvoAVy70UxNQ76TIWChciojIrHF8IMYjR0PsSATK9uBgSRyPgQ3zfWyrD7C1zs/51V5c2iZI5LQpXIqIyIwViVuebgvzSHOI7c1BnjuRuU3QolJ3as/Jyxf4KfdpmyCRiVK4FBGRGaUlaNjzSh/bm4M8djREd2RwrDvgdjYx35KYO9lQoU3MRXJN4VJERKaluLUc7onxUkeEg4mvF09GeL27COhMXbe0wpOaO7lhvp8ij8KkyGRSuBQRkYJ3fCDGwY4IL3VEU0Hy5c4o/WmblyeVuC1X1hWxtc7ZKmhxmX7VieST/saJiEjB6I/GebkjmtYb6YTJ42kLb9ItKHaxvMqb9uXB1/4Wy86tz3PLRSRJ4VJERPIuFre80RPlYDJInnTC5Js9MYb3RUKZ17Cs0gmPy6u8LK/2sqLKS5V/+AKcxpOT334RGZnCpYiITBprLW0D8cSQ9uCw9iudEYKx4dd7DDRUeDJ6IpdXeTmj1K2FNyLThMKliIjkRG8kzsudTnh86eTgsPaJUPYh7foS92BPZOKrocKjutwi05zCpYiInJZo3PJ6d3TYApvDPVm6IoFyn2HFkJ7IZZVeKrMMaYvI9KdwKSIiWVlrOdYfT4XHlxI9ka92RQhlyZFelzOkvXLIApu6Eg1pi8wmCpciIkJ3OM6htNXZydXaneFsy2ucyjbLq7ysSBvWbqjw4HUpRIrMdgqXIiKzSCRuaeyKpvVGOo/f7s0+pF3pM4kQOdgTuazKqzKJIjIihUsRkRnIWktzXyzVE5nsjXy1K0oky/oanwvOTWz1syJtWHtBsUtD2iJyWhQuRUSmuc5QnEOdmZuOv9QRoXuEIe0lZe5UeEwOa59d7sGjIW0RyQGFSxGRaSIcs7yaNqSdDJNNfdmHtGv8rtTq7BXVTpg8r9JDqVdD2iIyeRQuRUQKjLWWI72xjJ7Igx0RGruiZCmlTcAN51UOzolMDmvPK9KQtojkn8KliMgU6gjFM8ofHuyIcqgzQk9keIo0wFnJIe3q5CIbD2eVeXBrSFtECoTCpYhIHoRillc6I8MW2Bzrz169Zm7AlbHp+IoqL+dWeijRkLaIFLhRw6Ux5nbgvUCbtXZl2vEvAn8GRIH/tNb+z8TxG4BPAzHgz621/zUZDRcRKUTxxJB2evnDgx0RXuuOEssypF3sMZxX6Rm2wGZukTv/jRcRyYGx9FzeAdwK3JU8YIy5ErgaWGWtDRlj5iWOLwc+BKwAFgLbjTFLrbXZZ5uLiExjJ4KxjPKHBzsiHOqI0pdlYqTLONVr0mtpr6jysqTMjUvzIkVkBhk1XFprdxljlgw5/DngW9baUOKatsTxq4F7EsffNMa8BqwD9uasxSIieTYQdYa0XxqywKZ1IPuQdm2RK6P84fIqL+dWeinyKESKyMw33jmXS4HLjDHfBILAddbap4E64Mm065oSx0bU2Ng4ziacvny+l0ih0H0/djELzUHDa30uXu93vr/W76JpwBBneDAsclnOLolzTrHlnJI45xTHObskTqV3yIUd0NSRn88gDt33Mhvl675vaGg45fnxhksPUAVcCqwFfm6MOQuy/OsL2XfxHWMDc6WxsTFv7yVSKHTfj+z4gLPVz4tpPZEvd0QZyDIx0m3g3HJP5gKbai9nlGpIuxDpvpfZqJDu+/GGyybgfmutBZ4yxsSBOYnji9KuqweOTqyJIiLj1x+N83JHNDGkPTisfTyYfUh7YXH6kLYTJpdWeAloSFtEZEzGGy7/D7AF2GmMWQr4gHbgQeDfjTHfwVnQ0wA8lYuGioicSixueaMnysFkkEys1n6zJ5Z1+KTMazJ6IpNfVX5t9SMiMhFj2YroZ8BmYI4xpgm4EbgduN0Y8yIQBq5N9GK+ZIz5OXAQZ4uiL2iluIjkkrWW1oF4ap/IZE/kK50Rgln+tfEYWFrhYXl15gKbRSVuVa8REZkEY1ktfs0Ipz46wvXfBL45kUaJiAD0RuIcSoTH9GHtk6HsQ9r1Je7UPpHJr4YKDz63QqSISL6oQo+ITLlo3PJadzQxlD04P/Kt3uwDH+U+k6qfnSyBuKzKS4VPQ9oiIlNN4VJE8sZay7H+eEb5w5c6orzaGSGcpTPS63KGtFcMWWBTpyFtEZGCpXApIpOiOxznUNqcyGRvZGc4++5kZ5S6M8ofLq/yck6FB69LIVJEZDpRuBSRCYnELY1d0bTeSOfx2yMMaVf5TUb5w+VVHs6r9FKuIW0RkRlB4VJExsRaS1NfLKP84UsdERq7okSyDGn73XBuhRMeV1R5U6u1a4tcGtIWEZnBFC5FZJjO0OC8yFSY7IzQnWVI2wBnlrmH9UaeVe7BoyFtEZFZR+FSZBYLxSyvJoe0Tw6Gyeb+7EPaNX5XRvnD5VVezqv0UOrVkLaIiDgULkVmAWstR3pjmT2RiSHtaJb1NQE3nFc5uDo7uVp7noa0RURkFAqXIjOMtU5v5J6WMLve9NL8ynEOdUboiWQf0j673J2x6fiKKg9nlnlwa0hbRETGQeFSZJqLW8tLHVH2tITY3RJiT2uY9mByhY0Xp0IrzCtyZZQ/XFHl5dxKD8UeDWmLiEjuKFyKTDPRuOX5ExF2t4TY3Rpmb2uIriELbeYXudhY62cxXWxeupAV1R7mBNxT1GIREZlNFC5FClwoZnm2PczuljB7WkPsaw3TO2Si5KJSNxvm+9hY62fjfD9nlTsVbBobT9Cw0D9FLRcRkdlI4VKkwPRH4zzdFmFPqzPM/czxMMEhi7fPLnezYb6fjbV+NtT6OKNUf5VFRKQw6DeSyBTrDsd5qi2cmi95oD08bFPyZZUeJ0jO97Gh1k9tsYa4RUSkMClciuRZRyjO3tZQapj7dycixNNGuV0GLqjxpoa518/3UaP5kiIiMk0oXIpMsraBGHtawuxODHMf7IhmnPcYWD3Xm+iZ9HPJfB8VqrMtIiLTlMKlSI419UbZ0zo4zN3YlRkm/W5YPSex+KbWx9q5PkpU4UZERGYIhUuRCbDWcrgnluiVDLOnJcRbvZmrb4o9hkvm+VLzJVfP8RHwaINyERGZmRQuRU5DsvpNcr7knpYQR/szV9+Uew3rE0FyY62fC2q8eFXtRkREZgmFS5FTiFvLiycj7Gl1eiUzq984qv2u1OKbDbU+VlZ5VTpRRERmLYVLkTTRuOV3JyLsaQnxRGuYJ09R/WZjrY8N8/2cW+nBZRQmRUREQOFSZrlQzHIgWf2mJcS+tjB9WarfbEwOc6dVvxEREZHhRg2XxpjbgfcCbdbalUPOXQf8PTDXWttunN+43wPeA/QDn7DWHsh9s0XGx6l+E2Z3Ypg7W/Wbc8o9bEj0Sqr6jYiIyOkZy2/NO4BbgbvSDxpjFgHvAI6kHb4KaEh8XQL8IPFdZEp0h+Psa3OC5O6WMM+eGF79ZnmlJ7H4xsf6+ap+IyIiMhGjhktr7S5jzJIsp74L/E/ggbRjVwN3WWst8KQxptIYs8BaeywXjRUZzclgjL2t4dQ+k8+fzF79JjlfUtVvREREcmtc433GmN8Hmq21vxsy96wOeDvt56bEsRHDZWNj43iaMC75fC/Jj/YwPNft5kCXi2e73LzWn7kZudtYzi+Lc1FFnIvLY1xQHic1yh2Gk2/Dyfw3O69038tspPteZqN83fcNDQ2nPH/a4dIYUwx8FXhnttNZjtksx1JGa2CuNDY25u29ZPI09UZT8yVHqn6zZq7TK7mp1seaWV79Rve9zEa672U2KqT7fjw9l2cDZwLJXst64IAxZh1OT+WitGvrgaMTbaTMTtZa3uyJpcoo7m4JcWRI9ZsSj2FdovrNxlo/F6v6jYiIyJQ67XBprX0BmJf82RhzGFiTWC3+IPBnxph7cBbydGm+pYyVtZZXuqLsaUnW5Q5xbGj1G59h/bzkhuWqfiMiIlJoxrIV0c+AzcAcY0wTcKO19rYRLv81zjZEr+FsRfTJHLVTZqBY3PJSRyStlGKYE6HMMFnjd2VsC6TqNyIiIoVtLKvFrxnl/JK0xxb4wsSbJTNRJG55/kSE3S0hdreG2dsaontI9ZvaRPWbDbVO7+TSClW/ERERmU60O7RMmlDMsv/44LZAT2WpfnNGqTs1X3JjrZ8zy1T9RkREZMA20voAAB78SURBVDpTuJSc6YvEeea4U/1md6L6TWiE6jcba/1smO9jkarfiIiIzCj6zS7jlqx+s7vFmS95oD3MkI5Jlld6UsPcG+b7ma/qNyIiIjOawqWM2clgjD2tzuKb3S1hXshS/ebCGu/gApz5PqpV/UZERGRWUbiUEbX2x1KruHe3hDjYmblhucckNixPDHOvm+ejwjd7NywXERERhUtJ83ZvNLX4Zk9LmNe6s1e/2VjrZ+N8Vb8RERGR4RQuZ6lk9ZsnWkLsSWwN9HaW6jeXzPOxIbH4ZvVcH363VnKLiIjIyBQuZ4lk9Ztkr+SI1W/mO72SG2v9rFL1GxERETlNCpczVCxuebEjkpovubf11NVvNtb6WKHqNyIiIjJBCpczRCRu+d2JiDPE3RJib1t4xOo3ya2Bzq3waMNyERERySmFy2kqGLUcaE8svmkNj1j9JrlZuarfiIiISD4oXE4TfZE4Tx8PszsxzL2/fXj1m4YKDxvmDy7AUfUbERERyTeljwLVFY6zL7ktUGuIZ9sjw6vfVHnYON8Z5l4/36fqNyIiIjLlFC4LRLL6TXKY+1TVbzbOd8Kkqt+IiIhIoVG4nCIt/TH2JILknlNUv9lY6wxzXzLPR7mq34iIiEiBU7jMk7d7o+xuSdblDvF6d+aEyYA7WUrR2Wdy7TwfxR6FSREREZleFC4ngbWWN7pj7E4EyT2jVL/ZWOvj4jmqfiMiIiLTn8JlDsSt5ZXOaCpI7mkJ0TKQvfrNpsRqblW/ERERkZlI4XIcktVvdreEU/MmT45Q/Sa5z6Sq34iIiMhsoHA5BpG45bn2SGq+5JOtYbojmfsCLSh2JYKkM8y9VNVvREREZBZSuMwiGLXsb3d6JXcnqt/0D9lkcnGpOzVfcuN8P0tU/UZERERk9HBpjLkdeC/QZq1dmTj298D7gDDwOvBJa21n4twNwKeBGPDn1tr/mqS250yy+s0TiWHuU1W/SQ5z16v6jYiIiMgwY0lIdwC3AnelHXsYuMFaGzXGfBu4AbjeGLMc+BCwAlgIbDfGLLXWDolqU6srHOfJ1uR8yVNUv6n1pzYsV/UbERERkdGNGi6ttbuMMUuGHPtt2o9PAh9MPL4auMdaGwLeNMa8BqwD9uakteMUi1sebXdz+8lOdreEebEje/WbZK/khlo/VX7tMSkiIiJyunIxtvsp4N7E4zqcsJnUlDg2pVwGvv26jxORPgC8Llg715ca5l6n6jciIiIiOTGhcGmM+SoQBe5OHspymc1yLKWxsXEiTRizP1zgIWoNF1fEOL8sTsDd75wYgNa3oDUvrRDJv3z9HRMpJLrvZTbK133f0NBwyvPjDpfGmGtxFvpstdYmA2QTsCjtsnrg6EQamCufpjFv7yVSKBobdd/L7KP7XmajQrrvxzUWbIx5N3A98PvW2v60Uw8CHzLG+I0xZwINwFMTb6aIiIiITAdj2YroZ8BmYI4xpgm4EWd1uB94OLG345PW2j+11r5kjPk5cBBnuPwLhbZSXEREREQmz1hWi1+T5fBtp7j+m8A3J9IoEREREZmetERaRERERHJG4VJEREREckbhUkRERERyRuFSRERERHJG4VJEREREckbhUkRERERyRuFSRERERHJG4VJEREREckbhUkRERERyRuFSRERERHJG4VJEREREckbhUkRERERyRuFSRERERHJG4VJEREREckbhUkRERERyRuFSRERERHJG4VJEREREckbhUkRERERyRuFSRERERHJG4VJEREREckbhUkRERERyRuFSRERERHJm1HBpjLndGNNmjHkx7Vi1MeZhY0xj4ntV4rgxxvyTMeY1Y8zzxpiLJ7PxIiIiIlJYxtJzeQfw7iHHvgLssNY2ADsSPwNcBTQkvj4L/CA3zRQRERGR6WDUcGmt3QWcHHL4auDOxOM7gf+Wdvwu63gSqDTGLMhVY0VERESksI13zuV8a+0xgMT3eYnjdcDbadc1JY6JiIiIyCzgyfHrmSzH7Kme0NjYmOMmFMZ7iRQK3fcyG+m+l9koX/d9Q0PDKc+PN1y2GmMWWGuPJYa92xLHm4BFadfVA0cn0sBcaWxszNt7iRQK3fcyG+m+l9mokO778Q6LPwhcm3h8LfBA2vGPJ1aNXwp0JYfPRURERGTmG7Xn0hjzM2AzMMcY0wTcCHwL+Lkx5tPAEeAPE5f/GngP8BrQD3xyEtosIiIiIgVq1HBprb1mhFNbs1xrgS9MtFEiIiIiMj2pQo+IiIiI5IzCpYiIiIjkjMKliIiIiOSMwqWIiIiI5IzCpYiIiIjkjMKliIiIiOSMwqWIiIiI5IzCpYiIiIjkjMKliIiIiOSMwqWIiIiI5IzCpYiIiIjkjMKliIiIiOSMwqWIiIiI5IxnqhsgIiIiIuPU3YmJhqe6FRkULkVERESmk+5OPPt34XlqJ+5DzxHd/F7Y+L6pblWKwqWIiIhIoRsSKI2NA2DdbgiHprhxmRQuRURERApRdyee/Y/jeepR3C8/h4kPBsroikuIrt1M9OKNUFoOjY1T3NhBCpciIiIihSIZKJ/eifvQs5mBctWQQFmgFC5FREREptIMCJTpFC5FRERE8u1UgfL8dUTXXTmtAmU6hUsRERGRfOjpxLP/CWcOZdZAuZnoxZumZaBMp3ApIiIiMllSgXIn7kMHZmygTDehcGmM+TLwGcACLwCfBBYA9wDVwAHgY9bawtrdU0RERGSyzMJAmW7c4dIYUwf8ObDcWjtgjPk58CHgPcB3rbX3GGN+CHwa+EFOWisiIiJSiE4ZKNcm5lDmPlDaSC82Xlh9eBMdFvcARcaYCFAMHAO2AB9OnL8T+AYKlyIiIjLT9HbheebxvAZKa+PEe14jduIZYif3E+8+hKfuvWDekbP3mKhxh0trbbMx5h+AI8AA8FtgP9BprY0mLmsC6ibcShEREZFCkAyUTz+G++D+4YFy7WaiqzdBaUXO3zp8+F4ib98Pka7Bg8YN0T7w5vztxm0iw+JVwNXAmUAn8AvgqiyX2lO9TmMed5TP53uJFArd9zIb6b6XXHL391L5yrNUHnqGsjdfHiy9aFx0n7WCjmWr6Tr3ImLFpc4TjrUBbRN6T2/4CIH+3zFQfBFRXz0Apd0nKY90EXVXEQosJxRYRiiwFOsqAvJ33zc0NJzy/ESGxbcBb1prjwMYY+4HNgCVxhhPoveyHjg6kQbmSmNjY97eS6RQ6L6X2Uj3veTESD2ULhfRFWudRTmrN+EqraAGqJng28WDbRhPCcZTAkDold8S7fkt1TVz8Z11pfPe4TnYyHsxxfUYYzKeX0j3/UTC5RHgUmNMMc6w+FbgGeBR4IM4K8avBR6YaCNFREREJl1v1+CinKGBcuVgoMzFkLeNhYl1vkDs5DPETuzH9h/Bd95f4F34TgA88y4Dlwd3zZrUc4yvCuOrmvB7T7aJzLncZ4y5D2e7oSjwLPAvwH8C9xhj/jZx7LZcNFREREQk55KB8umduA8ewMRiQFqgXHuFEyjLKif0NtZabH8TsZP7ncU4nS9APDR4gbsYoj2DP1atwl21akLvOVUmtFrcWnsjcOOQw28A6ybyuiIiIiKTJl+BMtpPrON3id7JZ7DB1ozzrtKzcdeswV29BlfFMoxrZtS2mRmfQkRERORU8hQok2y0j/7HPwQ2MnjQW467ejXu6tV4alZPiyHu8VC4FBERkZmptxvPgUQt76GBcsWawTmUE+6h7CP86g+I979NYPU/YozBeEpwlZ0NxuUEypo1uMrOwRh3Lj5ZQVO4FBERkZkjFSgTi3JyHChtPEa85xXivW/hrUvswOguInriGYh0YvuPYEoWAxC4+H9jXDM/TA6lcCkiIiLT2yQHyniondiJ/c7cyZPPQrQXjAvPvMsw3lKMceFf9mVcgfmY4jNSz5uNwRIULkVERGQ6msRAaeNh4p0vETv5DNET+7F9hzPOm6KFuGvWYOMhDM7G6Z45l0z4I80UCpciIiIyPfT1JPahfDRLoFzt1PIeZ6CMB48Ta9/r9FB2PDdkm6AA7qoLcFevwV29Glfxwlx9ohlJ4VJEREQKVzJQPr0T90vPDA+UazcTXX0ZlJ9eoLSxIMRCGJ+zIXqsfR/hV/85dd5VemZiZfcaXJXLMS5f7j7TDKdwKSIiIoUlI1Dux8SiwMQDZVKk+TeEX/1nPPW/j7/hjwGc/SbnXe70TtZcjMs/J2cfZ7ZRuBQREZGpN1KgNOMPlDbSQ6zjWWInnsFVuQrvgm0AzrC2jWLDJ1PXuopqCaz8X7n9TLOUwqWIiIhMjb6ewUU5QwPl8osTcyjHHiitjRPvaXTKK57cT7zrZcCpD+4OdwyGy4oVFG+6JzUkLrmlcCkiIiL5kwyUTz+G+8VnsgTKzcRWX4YtH1v1mnjoJLGTBwa3CYp0DZ40blwVq3DXrMFTs2bwsMsDCpaTRuFSREREJleOAyVA+I1/I9b+JPHe1zOOm8D8VL1ud9UFGE9xTj+KjE7hUkRERHIvh4EyPtBCrONZPAvemSqfGOt8wQmWLl9imyCnZrcprscYM6kfTU5N4VJERERyo68Hz7O7nTmUEwiUNh7O2Pon+OxXsMEWXKVn4S4/FwDfkg9j7X/HXbkS4/ZP3meS06ZwKSIiIuOXg0BprcX2NznzJk88Q6zzRYrW347LXwOAZ97lxAeaAVfqOe7qCyf1Y8n4KVyKiIjI6RktUK69gtiay08dKKN9xDqeS9XstsG2jPPx7ldwzd0AgO+cT03eZ5GcU7gUERGR0U0wUFobJ977Rto2QQfBxgYv8Fbgrr4YT80a3NUXY3xjX9wjhUXhUkRERLLr68Hz7B6nlvfQQLnsosEh74rqU75MrPMlQi/+LTbcMXjQuHBVrHAW4tSswVV2Dsa4Rn4RmTYULkVERGRQMlA+vRP3C0+fdqCMD7QSPfZf4C7Ct/gPAXAV12HDHRj/nER5xTW4qy7EeEvz9rEkfxQuRUREZrvRAuXazcTWZA+U8eBxbKQHd9lZznMiXUQO/zsmMA/vGR/EGIPxVVJ06e2YogXaJmgWULgUERGZjfp78RzY7QTKF5/BRCPA6IHSxsPEO18kemI/sZP7sX2HcVWspGj1PwDgKjsH7xkfwFV5AWABJ0y6ihfm89PJFFK4FBERmS36ezPnUI4hUFprsQNHUwtxYh2/g3ho8DXdAYyvAmvjGOPCGBe+c/44359MCsiEwqUxphL4CbAS539PPgW8AtwLLAEOA39kre0Y4SVERERkMqUC5U7cLz6dGSjPu5DouiuHB8roALHO3zmB8sR+bPBYxku6Ss9KzJ1cjatiOcblzetHksI20Z7L7wEPWWs/aIzxAcXA/wJ2WGu/ZYz5CvAV4PoJvo+IiIiM1YiB0iQCZWJRTqWzSbm1FhI9jwChg7cQa987+Hqe0sSqbqfEYnJzc5Fsxh0ujTHlwOXAJwCstWEgbIy5GticuOxOYCcFEC7n7fkNvmcfxZZXYssqsaUV2PIK53FZJfgDU91EERGR8TvNQJkUfvNuos2/xr/8OtzVFwHgrlmDDXc4q7qr1+Aqb0jV9BYZjbHWju+JxlwI/AtwELgA2A98CWi21lamXddhrc3YCbWrqyv1po2NjeN6/9O17AdfJ3CiZcTzMa+PaHEZseJSosVlRFPf0x+XEi1xjsX8RaAVbyIiMoVcoQEqXv0dVQefoeyNl3AlV3lj6F28lM5lq+k872KipRVg43jDR/AHDzJQcgkxjxMyyzofoKxnOz3l76an4vem8uPINNHQ0JB6XFFRMSwMTSRcrgGeBDZaa/cZY74HdANfPJ1wmS/tv7yLBV6D6enE9HRhujsxvV2JnzsxkchpvZ51u53ez7LKRG9oxZCfK6GswjleVoktLQe31k9JfjU2Nmb8IyAyG8z4+36gL7HK+zHcLzyV0UMZP3cVkXVXpnoo46ETxE4eSCzGOQDRHgB8Sz+Pt/73AYgPtECsH1NyprYJmsam6r7PFi4nknaagCZr7b7Ez/fhzK9sNcYssNYeM8YsANpGfIU86li1njkj/aFbC8GBwaDZ05X5PT2Idie+B/sxXSeh6+SY3t8aA8VlzlB86eBwfCqYlqV/dx7j8+fwT0BERKatUwTK2HkXEF27meiay4mXlxPvOkis/QFir+4n3vt6xsuYQG2iGs7S1DFXUW1eP4rMfOMOl9baFmPM28aYc621rwBbcYbIDwLXAt9KfH8gJy2dTMZAUTG2qBg7b4z7cIVDmN7uIWF0SDDt7kwdp68H09eN6esG3h7TW9hA0fDe0SwhNBVGi0o0VC8iMlOkB8oXn0qNsA0LlKVFRFsfJXbkVmIdz0FsYPA1XH7cVatSVXFM0UL1Tsqkm+g47ReBuxMrxd8APgm4gJ8bYz4NHAH+cILvUZh8fmz1XGz13LFdH49Bb09aCE0G0bTe0Z5ED2mqd3QAExyA9pHniqazHm/mQqXSisEFTFmCKaVl4NIEbRGRgjHQl7YoJ3ugjFx8CTFvMFURh2gf4Vf/GWwMAFOyGHf1Gjw1q3FVrMS4fVP1aWSWmlC4tNY+B6zJcmrrRF53RnK5oTzRAzmW662Fgb7hITRjmD6zd9SEgpjOduhsH1OTrDFQWj4kiGbvFU31jnr1j5SISE4lA+XTO50h7/RAee4FRNc5PZS2soZ48DgDT34K3EUUb/oZxrgxnhK8i/87JjDXWdkdGGOnh8gk0QqTQmUMFJdii0ux8+vH9pxwaGwhNNk72tcNieA6VjZQnBimP1UITYTo0goIaFW9iMgwowTKyLr1hM6uJBp6BRt+jkDl+wEw/jkY3xyMtxQb7sQk9pv0nfXxKfsoIkMpXM4kPj+2Zh62Zt7Yro9FM+aNMnSIvqcrEUTTQmmwHxPsh+NHx/QW1utNC6HJAFqR/Vh5JRSXgcs1gT8EEZECdYpAGT13FaG1Kwgv8hLtO0i8+y54M556ajx0Epe/GmMMRZf8EOPWgk8pXAqXs5nbg62ozij5dUrWQn/v8N7RoSE0fWFTOIQ5eRxOHh/bW7hczhB9aQWUjxBCM3pHy8GjsmMiUqAG+vE8l6jlPSRQRlasIHjxGYTnhIn2vAiRVyH5/+3GhatiZWIT89UYX2qHPwVLKXgKlzJ2xkBJGbakDFs7xqH6UDAtiHalDcuPsO1Tfy+muwO6Owb/kR2FLS4ZvXc0uQVUeSX4i8b/ZyAiMppUoNyJ+4V9mUPeS1fRv+4MgmWvE+9/A3gdEjvaGf/cVJh0V1+E8ZRM3WcQmQCFS5lc/gDWX4udM8Z91KLRwfDZ25XZO9rTOTh0n77Qqb8P098Hrc1jegvr848+RJ++GX5xqeaNisipjRAoo8WG0CWLYcl6WPNBbNUcose2Ez+0HVxe3JXnD24TVLxI2wTJjKBwKYXF48FW1gyrfTuieDwxVD+W3tHOwaH6E61wonVMb+FUYyofFkIZOkSvakwis0uWQGldEPcYWLqK6LrN9NefIHzsl3jqIvir5gDgmbMOc8FNuCvPx7gDU/whRHJPvwFlenO5nO2USsuxC84Y/XprITSQuYI+FURH2Ax/oA/T1QFdHWNuli0pG3klfWLxUvo+pKrGJDJNDPTjeW5vYg7lPohEiJUb+s9xE2qoIVIZxLfwA3iXfQoAV08j7shR3JXnp17CeMvx1Kydqk8gMukULmV2MQYCxc6WSnMXjO05kfDgoqX03tFsc0Z7OqG3G9PXg+nrgZYxVmPyB04dQof0jqoak0geJQPl0ztxP78Pa8OEF7joX+0mtLiMeCCSuLAPgDiD27u5yxpwr7pxChotMnUULkVG4/WdfjWmvp7h+4sO3Qw//VgoiAm1jL0ak9szvFc0I4g634uOn8CU+JyeUa8P6008dnsUTkVOJS1Qup7fR6w0QqjORWiLm8i8gFOLDoAIeMpwV1+cWIxzMS7/GKf1iMxQCpciueZyQ7L3ceHi0a+3FoL92UuAjhRMgwOYzhPQeeKUL33eSG9pjFNtyevDetODp2/wcSKIJh/bREAd9pzE87JfNyTYev3gdivYSmEa0kNpImEscOL9PmLl6VNXXLjKz3XCZM0aXGXnYIxK6YokKVyKTDVjoKgEW1SCnV83tueEQ9mrLw0JoaGuDgLGQCTsDO+HQ873WBTCIed16JnczzeENa7B8OkbKYQOBthhgfVUAdg79Dp/RgDWQisZJpicQ+n0UPYvj9E/30VVLEJs6flE126G6qcxwbdx16x2AmXVRRhv2VS3XKRg6V9akenI58dWz8NWn7oaU2NjIw0NDcNPxGMQiUAkhAkngmckNBhCI2EIpz8OOd+HhtRIGCLpj8Op10u/Lv31TSwG4SCEg5i+SfrzGYF1uYYF1lSPqneE3tdEmB3+nMzwOrzHNjMA41LPVsFIBErXsw8Taz9A4LXw4KmGcmIlYbq/8TXci68AwB95F3hKtE2QyBgpXIrMRi43+N3OPqSJQ/aUT8ihWBSikYzAmhFIh4XZtMfhtJA6WgDOuC7kPDceh1DQ2dw/X583wbrdI0xFGBpsfWO8boQAnAy9accUbIFgP+5nn4CDvyY6cIi+Wkv0XBec68IdOBd7wTuJrrkcT+xVPMaNq+rC1FONt3QKGy4y/Shcikh+uT3Ol79oaoJtKqRm610dGmZPHYCzPmdoz244cT4Wg9gABAemINh6sofUIVMMsk1FGGsATr92sGfX52wXNlWC/bgO/Ib4G78lEj9M/3ywZxog8WU9uCtWEvzcn+AqPRMAD2NcuCciI1K4FJHZIxlsA8X5DbbWQiw2fPpBePjjzF7aLFMWwiMH4PQQnDFlIRaFWBQT7M/Hp8386B7vCGE2M4Rm9sQOnX4wUo/t0IVlfqyNUP3KfxB77iUi3hZi5QYWAolI77JVuOddirvuclyVKzAuX97/TERmOoVLEZHJZgx4PE4FqqLi1OH8BdvoqcNsRiDNDK/DenZHCLPZpiKYSBgTjUA0ghnIzwTb4BkuAlf6CBYDGEzUhce1GNfiLbgWbcYVUM+kyGRTuBQRmcmMAY8XPF5sUUnqcN6C7ZD5tRmLv4ZNL8gegIc9Z+iUhbTX8XaGcHW58JatxHXuezF1l2I051QkrxQuRURkcqTvp5p2eLKDbdNIuySISF5M4UxrEREREZlpFC5FREREJGcULkVEREQkZxQuRURERCRnJhwujTFuY8yzxpj/m/j5TGPMPmNMozHmXmOMNhETERERmSVy0XP5JeBQ2s/fBr5rrW0AOoBP5+A9RERERGQamFC4NMbUA78H/CTxswG2APclLrkT+G8TeQ8RERERmT6MtePfccwYcx9wM1AGXAd8AnjSWntO4vwi4DfW2pXpz+vq6kq9aWNj47jfX0RERETyK30f2YqKCjP0/Lg3UTfGvBdos9buN8ZsTh7Ocukp02u+Nrpt1Ka6MgvpvpfZSPe9zEaFdN+Pu+fSGHMz8DEgCgSAcuA/gHcBtdbaqDFmPfANa+270p+b3nMpIiIiItNTtp7Lcc+5tNbeYK2tt9YuAT4EPGKt/QjwKPDBxGXXAg+M9z1EREREZHqZjH0urwf+whjzGlAD3DYJ7yEiIiIiBWhCC3pERERERNKpQo+IiIiI5MyMC5fGmN6pboNIvhhjrDHm39J+9hhjjicrZonMJLrfRRyjZR1jzE5jzJp8tWeoGRcuRWaZPmClMaYo8fM7gObTeQFjzLi3JBPJswnf7yIy+WZkuDTGlBpjdhhjDhhjXjDGXJ04vsQYc8gY82NjzEvGmN+m/SMlMl39BqdSFsA1wM+SJ4wx64wxe4wxzya+n5s4/gljzC+MMb8Cfpv/JouM23ju98eNMRemXbfbGLMqr60WyTFjzOb0XntjzK3GmE9MYZNSZmS4BILA+621FwNXAv87UZoSoAH4/621K4BO4ANT1EaRXLkH+JAxJgCsAvalnXsZuNxaexHw18DfpZ1bD1xrrd2St5aKTNx47vef4FSQwxizFPBba5/PW4tFZpmZOhxmgL8zxlwOxIE6YH7i3JvW2ucSj/cDS/LfPJHcsdY+b4xZgtOL8+shpyuAO40xDTjVsrxp5x621p7MSyNFcmSc9/svgK8bY/4K+BRwR14aKzJLzdSey48Ac4HV1toLgVacKkIAobTrYszcgC2zy4PAP5A2RJhwE/CotXYl8D4G/x6AM39NZDo6rfvdWtsPPAxcDfwR8O/5a6rIpImSmeMCI12YbzM1WFXg1D2PGGOuBBZPdYNEJtntQJe19gVjzOa04xUMLnj4RL4bJTJJxnO//wT4FfC4euxlhngLWG6M8eMEy63AE1PbJMeM6rlMrHoNAXcDa4wxz+D0Yr48pQ0TmWTW2iZr7feynLoFuNkYsxtw57lZIpNiPPe7tXY/0A38ax6aKDJpklnHWvs28HPgeZzc8+yUNizNjKrQY4y5APixtXbdVLdFREQKhzFmIbATOM9aG5/i5oiM23TIOjOm59IY86c482++NtVtERGRwmGM+TjOqvKvKljKdDZdss6M6rkUERERkak1Y3ouRURERGTqKVyKiIiISM4oXIqIiIhIzihcioiIiEjOKFyKiIiISM4oXIqIiIhIzvw/TYTpMouBtxgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,5))\n",
    "x = [\"Jan\", \"Mar\", \"May\", \"Jul\"]\n",
    "\n",
    "plt.plot(x, [120, 150, fl_before,  fl_after], label=\"FL\", lw=2)\n",
    "plt.plot(x, [60, 50, poa_before, poa_after], label=\"POA\", lw=2)\n",
    "\n",
    "plt.plot([\"May\", \"Jul\"], [poa_before, poa_before+(fl_after-fl_before)], label=\"Counterfactual\", lw=2, color=\"C2\", ls=\"-.\")\n",
    "\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will see how to solve this problem with synthetic control. It will use multiple cities to create a synthetic city that closely follows the trend of the city of interest. But for now, remember that you always need to check if you have parallel trends when applying diff-in-diff. \n",
    "\n",
    "![img](./data/img/diff-in-diff/non-parallel.png)\n",
    "\n",
    "One final issue that it's worth mentioning is that you won't be able to place confidence intervals around your Diff-in-Diff estimator if you only have aggregated thata. Say for instance you don't have data on what each of our customers from Florianópolis or Porto Alegre did. Instead, you only have the average deposits before and after the intervention for both cities. In this case, you will still be able to estimate the causal effect by Diff-in-Diff, but you won't know the variance of it. That's because all the variability in your data got squashed out in the aggregation.\n",
    "\n",
    "## Key Ideas\n",
    "\n",
    "We've explored a technique widely applied when we are estimating causal effects at more macro entities (schools, cities, states, countries...). Difference in Difference takes a treated unit before and after the treatment and compares the trend in the outcome to that of a control unit. Here, we've seen how this could be applied at estimating the effect of a city specific marketing camping.\n",
    "\n",
    "Finally, we looked at how Diff-in-Diff fails if the trend between the treated and control unit is not the same. We also saw how diff-in-diff will be problematic if we only have aggregated data.\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",
    "Other important 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",
    "Finally, I'ld also like to compliment Scott Cunningham and his brilinat work migling Causal Inference and Rap quotes:\n",
    "\n",
    "* [Causal Inference: The Mixtape](https://www.scunning.com/mixtape.html)"
   ]
  },
  {
   "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
}