13-Panel-Data-and-Fixed-Effects.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Controlling What you Cannot See\n",
    "\n",
    "Methods like propensity score, linear regression and matching are very good at controlling for confounding in non-random data, but they rely on a key assumption: conditional unconfoundedness\n",
    "\n",
    "$\n",
    "(Y_0, Y_1) \\perp T | X\n",
    "$\n",
    "\n",
    "To put it in words, they require that all the confounders are known and measured, so that we can condition on them and make the treatment as good as random. One major issue with this is that sometimes we simply can't measure a confounder. For instance, take a classical labor economics problem of figuring out the impact of marriage on men's earnings. It's a well known fact in Economics that married men earn more than single men. What is not as clear is if this relationship is causal or not. It could be that more educated people are both more likely to marry and more likely to have a high earnings job, which would mean that education is a confounder of the effect of marriage on earnings. For this confounder, we could measure the education of the person in the study and run a regression controlling for that. But another confounder could be beauty. It could be that more handsome men are both more likely to get married and more likely to have a high paying job. Unfortunately, beauty is one of those characteristics like intelligence. It's something we can't measure very well.  \n",
    "\n",
    "This puts us in a difficult situation, because if we have unmeasured confounders, we have bias. One way to deal with this is with instrumental variables, like we've seen before. But coming up with good instruments it's no easy task and requires a lot of creativity. Here, let's look at an alternative that takes advantage of time or the temporal structure of data. \n",
    "\n",
    "The idea is to use **panel data**. Panel data is when we have **observations on the same individual for multiple periods of time**. Panel data formats are very common in the industry, where they keep records of customer behavior for the same customer and for multiple time periods. The reason we can leverage panel data is because we can compare the same unit before and after the treatment, and see how they behave with it. Before we dive in the math, let's see how this makes intuitive sense."
   ]
  },
  {
   "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 statsmodels.formula.api as smf\n",
    "import graphviz as gr\n",
    "\n",
    "from linearmodels.panel import PanelOLS\n",
    "\n",
    "\n",
    "%matplotlib inline\n",
    "pd.set_option(\"display.max_columns\", 6)\n",
    "style.use(\"fivethirtyeight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's take a look at the causal graph that we have once we include multiple observations of the same unit across time. Suppose we have a situation where marriage at the first time causes mariage also at the first time. This is also true for times 2 and 3. Also, suppose that beauty is the same across all time periods (a bold statement, but reasonable if time is just a few years) and it causes both marriage and income."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"308pt\" height=\"332pt\"\n",
       " viewBox=\"0.00 0.00 307.94 332.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 328)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-328 303.9355,-328 303.9355,4 -4,4\"/>\n",
       "<!-- Marriage_1 -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>Marriage_1</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"139.4677\" cy=\"-234\" rx=\"51.2156\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"139.4677\" y=\"-229.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Marriage_1</text>\n",
       "</g>\n",
       "<!-- Income_1 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>Income_1</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"62.4677\" cy=\"-162\" rx=\"45.4355\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"62.4677\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Income_1</text>\n",
       "</g>\n",
       "<!-- Marriage_1&#45;&gt;Income_1 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>Marriage_1&#45;&gt;Income_1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M121.2198,-216.937C111.2712,-207.6344 98.8003,-195.9733 87.9105,-185.7906\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"90.0138,-182.9656 80.319,-178.6921 85.2328,-188.0786 90.0138,-182.9656\"/>\n",
       "</g>\n",
       "<!-- Marriage_2 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>Marriage_2</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"177.4677\" cy=\"-162\" rx=\"51.2156\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"177.4677\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Marriage_2</text>\n",
       "</g>\n",
       "<!-- Marriage_1&#45;&gt;Marriage_2 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>Marriage_1&#45;&gt;Marriage_2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M148.861,-216.2022C153.2,-207.981 158.4461,-198.041 163.2546,-188.9301\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"166.4785,-190.3201 168.0508,-179.8425 160.2878,-187.0527 166.4785,-190.3201\"/>\n",
       "</g>\n",
       "<!-- Income_2 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>Income_2</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"45.4677\" cy=\"-90\" rx=\"45.4355\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"45.4677\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Income_2</text>\n",
       "</g>\n",
       "<!-- Marriage_2&#45;&gt;Income_2 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>Marriage_2&#45;&gt;Income_2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M149.4739,-146.7307C129.6083,-135.8949 102.774,-121.258 81.3385,-109.5658\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"82.9733,-106.4708 72.5184,-104.7549 79.6213,-112.6161 82.9733,-106.4708\"/>\n",
       "</g>\n",
       "<!-- Marriage_3 -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>Marriage_3</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"215.4677\" cy=\"-90\" rx=\"51.2156\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"215.4677\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Marriage_3</text>\n",
       "</g>\n",
       "<!-- Marriage_2&#45;&gt;Marriage_3 -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>Marriage_2&#45;&gt;Marriage_3</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M186.861,-144.2022C191.2,-135.981 196.4461,-126.041 201.2546,-116.9301\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"204.4785,-118.3201 206.0508,-107.8425 198.2878,-115.0527 204.4785,-118.3201\"/>\n",
       "</g>\n",
       "<!-- Income_3 -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>Income_3</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"254.4677\" cy=\"-18\" rx=\"45.4355\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"254.4677\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Income_3</text>\n",
       "</g>\n",
       "<!-- Marriage_3&#45;&gt;Income_3 -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>Marriage_3&#45;&gt;Income_3</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M225.1082,-72.2022C229.6082,-63.8944 235.059,-53.8316 240.0363,-44.6427\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"243.1177,-46.3025 244.803,-35.8425 236.9626,-42.9685 243.1177,-46.3025\"/>\n",
       "</g>\n",
       "<!-- Beauty -->\n",
       "<g id=\"node7\" class=\"node\">\n",
       "<title>Beauty</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"178.4677\" cy=\"-306\" rx=\"35.3129\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"178.4677\" y=\"-301.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Beauty</text>\n",
       "</g>\n",
       "<!-- Beauty&#45;&gt;Marriage_1 -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>Beauty&#45;&gt;Marriage_1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M169.0269,-288.5708C164.4957,-280.2055 158.9666,-269.998 153.9215,-260.6839\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"156.9329,-258.8948 149.0925,-251.7689 150.7779,-262.2288 156.9329,-258.8948\"/>\n",
       "</g>\n",
       "<!-- Beauty&#45;&gt;Income_1 -->\n",
       "<g id=\"edge9\" class=\"edge\">\n",
       "<title>Beauty&#45;&gt;Income_1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M147.7259,-297.1258C125.3066,-289.0322 96.0559,-274.7957 79.4677,-252 66.581,-234.2909 62.64,-209.5853 61.7491,-190.5381\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"65.2442,-190.2692 61.5312,-180.3463 58.2458,-190.4189 65.2442,-190.2692\"/>\n",
       "</g>\n",
       "<!-- Beauty&#45;&gt;Marriage_2 -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>Beauty&#45;&gt;Marriage_2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M187.1872,-288.491C191.8437,-278.1547 197.0937,-264.6547 199.4677,-252 202.4179,-236.2743 202.5435,-231.7016 199.4677,-216 197.7054,-207.0032 194.4367,-197.6139 190.948,-189.2358\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"194.0829,-187.6709 186.8286,-179.9491 187.6841,-190.5094 194.0829,-187.6709\"/>\n",
       "</g>\n",
       "<!-- Beauty&#45;&gt;Income_2 -->\n",
       "<g id=\"edge10\" class=\"edge\">\n",
       "<title>Beauty&#45;&gt;Income_2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M148.0237,-296.7162C122.9313,-288.0074 87.3829,-273.1426 61.4677,-252 30.6793,-226.8815 20.8642,-217.7517 8.4677,-180 3.476,-164.7986 3.734,-159.2837 8.4677,-144 11.6451,-133.7413 17.5185,-123.6981 23.6645,-115.1092\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"26.5293,-117.1245 29.8116,-107.051 20.9637,-112.8789 26.5293,-117.1245\"/>\n",
       "</g>\n",
       "<!-- Beauty&#45;&gt;Marriage_3 -->\n",
       "<g id=\"edge8\" class=\"edge\">\n",
       "<title>Beauty&#45;&gt;Marriage_3</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M191.4291,-289.2245C198.9308,-278.9172 208.115,-265.1982 214.4677,-252 229.0372,-221.7308 232.2949,-213.1924 237.4677,-180 239.9315,-164.1908 240.5435,-159.7016 237.4677,-144 235.7054,-135.0032 232.4367,-125.6139 228.948,-117.2358\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"232.0829,-115.6709 224.8286,-107.9491 225.6841,-118.5094 232.0829,-115.6709\"/>\n",
       "</g>\n",
       "<!-- Beauty&#45;&gt;Income_3 -->\n",
       "<g id=\"edge11\" class=\"edge\">\n",
       "<title>Beauty&#45;&gt;Income_3</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M199.1025,-291.1171C211.326,-281.2928 226.1906,-267.4177 235.4677,-252 269.7141,-195.0857 266.4966,-173.8147 275.4677,-108 277.6287,-92.1466 278.4179,-87.7257 275.4677,-72 273.7881,-63.0465 270.6687,-53.67 267.3385,-45.2901\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"270.526,-43.8414 263.4059,-35.9958 264.0794,-46.5692 270.526,-43.8414\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a16cea048>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = gr.Digraph()\n",
    "g.edge(\"Marriage_1\", \"Income_1\")\n",
    "g.edge(\"Marriage_1\", \"Marriage_2\")\n",
    "g.edge(\"Marriage_2\", \"Income_2\")\n",
    "g.edge(\"Marriage_2\", \"Marriage_3\")\n",
    "g.edge(\"Marriage_3\", \"Income_3\")\n",
    "\n",
    "g.edge(\"Beauty\", \"Marriage_1\")\n",
    "g.edge(\"Beauty\", \"Marriage_2\")\n",
    "g.edge(\"Beauty\", \"Marriage_3\")\n",
    "\n",
    "g.edge(\"Beauty\", \"Income_1\")\n",
    "g.edge(\"Beauty\", \"Income_2\")\n",
    "g.edge(\"Beauty\", \"Income_3\")\n",
    "\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that we cannot control beauty, since we can't measure it. But we can still use the panel structure so it is not a problem anymore. The idea is that we can see beauty - and any other attribute that is constant across time - as the defining aspects of a person. And although we can't control them directly, we can control for the individual itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.40.1 (20161225.0304)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"308pt\" height=\"332pt\"\n",
       " viewBox=\"0.00 0.00 307.94 332.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 328)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"#ffffff\" stroke=\"transparent\" points=\"-4,4 -4,-328 303.9355,-328 303.9355,4 -4,4\"/>\n",
       "<!-- Marriage_1 -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>Marriage_1</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"139.4677\" cy=\"-234\" rx=\"51.2156\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"139.4677\" y=\"-229.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Marriage_1</text>\n",
       "</g>\n",
       "<!-- Income_1 -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>Income_1</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"62.4677\" cy=\"-162\" rx=\"45.4355\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"62.4677\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Income_1</text>\n",
       "</g>\n",
       "<!-- Marriage_1&#45;&gt;Income_1 -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>Marriage_1&#45;&gt;Income_1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M121.2198,-216.937C111.2712,-207.6344 98.8003,-195.9733 87.9105,-185.7906\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"90.0138,-182.9656 80.319,-178.6921 85.2328,-188.0786 90.0138,-182.9656\"/>\n",
       "</g>\n",
       "<!-- Marriage_2 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>Marriage_2</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"177.4677\" cy=\"-162\" rx=\"51.2156\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"177.4677\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Marriage_2</text>\n",
       "</g>\n",
       "<!-- Marriage_1&#45;&gt;Marriage_2 -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>Marriage_1&#45;&gt;Marriage_2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M148.861,-216.2022C153.2,-207.981 158.4461,-198.041 163.2546,-188.9301\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"166.4785,-190.3201 168.0508,-179.8425 160.2878,-187.0527 166.4785,-190.3201\"/>\n",
       "</g>\n",
       "<!-- Income_2 -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>Income_2</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"45.4677\" cy=\"-90\" rx=\"45.4355\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"45.4677\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Income_2</text>\n",
       "</g>\n",
       "<!-- Marriage_2&#45;&gt;Income_2 -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>Marriage_2&#45;&gt;Income_2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M149.4739,-146.7307C129.6083,-135.8949 102.774,-121.258 81.3385,-109.5658\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"82.9733,-106.4708 72.5184,-104.7549 79.6213,-112.6161 82.9733,-106.4708\"/>\n",
       "</g>\n",
       "<!-- Marriage_3 -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>Marriage_3</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"215.4677\" cy=\"-90\" rx=\"51.2156\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"215.4677\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Marriage_3</text>\n",
       "</g>\n",
       "<!-- Marriage_2&#45;&gt;Marriage_3 -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>Marriage_2&#45;&gt;Marriage_3</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M186.861,-144.2022C191.2,-135.981 196.4461,-126.041 201.2546,-116.9301\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"204.4785,-118.3201 206.0508,-107.8425 198.2878,-115.0527 204.4785,-118.3201\"/>\n",
       "</g>\n",
       "<!-- Income_3 -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>Income_3</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"254.4677\" cy=\"-18\" rx=\"45.4355\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"254.4677\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Income_3</text>\n",
       "</g>\n",
       "<!-- Marriage_3&#45;&gt;Income_3 -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>Marriage_3&#45;&gt;Income_3</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M225.1082,-72.2022C229.6082,-63.8944 235.059,-53.8316 240.0363,-44.6427\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"243.1177,-46.3025 244.803,-35.8425 236.9626,-42.9685 243.1177,-46.3025\"/>\n",
       "</g>\n",
       "<!-- Person (Beauty, Ingeligence...) -->\n",
       "<g id=\"node7\" class=\"node\">\n",
       "<title>Person (Beauty, Ingeligence...)</title>\n",
       "<ellipse fill=\"none\" stroke=\"#000000\" cx=\"178.4677\" cy=\"-306\" rx=\"119.1136\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"178.4677\" y=\"-301.8\" font-family=\"Times,serif\" font-size=\"14.00\" fill=\"#000000\">Person (Beauty, Ingeligence...)</text>\n",
       "</g>\n",
       "<!-- Person (Beauty, Ingeligence...)&#45;&gt;Marriage_1 -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>Person (Beauty, Ingeligence...)&#45;&gt;Marriage_1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M168.6264,-287.8314C164.205,-279.6688 158.8973,-269.87 154.0292,-260.8827\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"157.0112,-259.0393 149.1708,-251.9134 150.8562,-262.3734 157.0112,-259.0393\"/>\n",
       "</g>\n",
       "<!-- Person (Beauty, Ingeligence...)&#45;&gt;Income_1 -->\n",
       "<g id=\"edge9\" class=\"edge\">\n",
       "<title>Person (Beauty, Ingeligence...)&#45;&gt;Income_1</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M129.1653,-289.5102C110.9306,-281.0369 91.7361,-268.8593 79.4677,-252 66.581,-234.2909 62.64,-209.5853 61.7491,-190.5381\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"65.2442,-190.2692 61.5312,-180.3463 58.2458,-190.4189 65.2442,-190.2692\"/>\n",
       "</g>\n",
       "<!-- Person (Beauty, Ingeligence...)&#45;&gt;Marriage_2 -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>Person (Beauty, Ingeligence...)&#45;&gt;Marriage_2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M187.4059,-288.0042C192.0071,-277.733 197.1308,-264.457 199.4677,-252 202.4179,-236.2743 202.5435,-231.7016 199.4677,-216 197.7054,-207.0032 194.4367,-197.6139 190.948,-189.2358\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"194.0829,-187.6709 186.8286,-179.9491 187.6841,-190.5094 194.0829,-187.6709\"/>\n",
       "</g>\n",
       "<!-- Person (Beauty, Ingeligence...)&#45;&gt;Income_2 -->\n",
       "<g id=\"edge10\" class=\"edge\">\n",
       "<title>Person (Beauty, Ingeligence...)&#45;&gt;Income_2</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M128.9388,-289.5858C106.8146,-280.7032 81.3158,-268.1928 61.4677,-252 30.6793,-226.8815 20.8642,-217.7517 8.4677,-180 3.476,-164.7986 3.734,-159.2837 8.4677,-144 11.6451,-133.7413 17.5185,-123.6981 23.6645,-115.1092\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"26.5293,-117.1245 29.8116,-107.051 20.9637,-112.8789 26.5293,-117.1245\"/>\n",
       "</g>\n",
       "<!-- Person (Beauty, Ingeligence...)&#45;&gt;Marriage_3 -->\n",
       "<g id=\"edge8\" class=\"edge\">\n",
       "<title>Person (Beauty, Ingeligence...)&#45;&gt;Marriage_3</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M192.4781,-287.7759C199.7726,-277.6519 208.4082,-264.5891 214.4677,-252 229.0372,-221.7308 232.2949,-213.1924 237.4677,-180 239.9315,-164.1908 240.5435,-159.7016 237.4677,-144 235.7054,-135.0032 232.4367,-125.6139 228.948,-117.2358\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"232.0829,-115.6709 224.8286,-107.9491 225.6841,-118.5094 232.0829,-115.6709\"/>\n",
       "</g>\n",
       "<!-- Person (Beauty, Ingeligence...)&#45;&gt;Income_3 -->\n",
       "<g id=\"edge11\" class=\"edge\">\n",
       "<title>Person (Beauty, Ingeligence...)&#45;&gt;Income_3</title>\n",
       "<path fill=\"none\" stroke=\"#000000\" d=\"M202.4955,-288.3392C213.9825,-278.7573 227.034,-266.0161 235.4677,-252 269.7141,-195.0857 266.4966,-173.8147 275.4677,-108 277.6287,-92.1466 278.4179,-87.7257 275.4677,-72 273.7881,-63.0465 270.6687,-53.67 267.3385,-45.2901\"/>\n",
       "<polygon fill=\"#000000\" stroke=\"#000000\" points=\"270.526,-43.8414 263.4059,-35.9958 264.0794,-46.5692 270.526,-43.8414\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a16cdceb8>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = gr.Digraph()\n",
    "g.edge(\"Marriage_1\", \"Income_1\")\n",
    "g.edge(\"Marriage_1\", \"Marriage_2\")\n",
    "g.edge(\"Marriage_2\", \"Income_2\")\n",
    "g.edge(\"Marriage_2\", \"Marriage_3\")\n",
    "g.edge(\"Marriage_3\", \"Income_3\")\n",
    "\n",
    "g.edge(\"Person (Beauty, Ingeligence...)\", \"Marriage_1\")\n",
    "g.edge(\"Person (Beauty, Ingeligence...)\", \"Marriage_2\")\n",
    "g.edge(\"Person (Beauty, Ingeligence...)\", \"Marriage_3\")\n",
    "\n",
    "g.edge(\"Person (Beauty, Ingeligence...)\", \"Income_1\")\n",
    "g.edge(\"Person (Beauty, Ingeligence...)\", \"Income_2\")\n",
    "g.edge(\"Person (Beauty, Ingeligence...)\", \"Income_3\")\n",
    "\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Think about it. We don't know the attributes like beauty and intelligence, but we know that person is the same individual. So, we can create a dummy variable indicating that person and add that to a linear model. This is what we mean when we say we can control for the person itself. When estimating the effect of marriage on income with this person dummy in our model, what regression does is it finds the effect of marriage **while keeping the person variable fixed**. Adding this entity dummy is what we call a fixed effect model.\n",
    "\n",
    "\n",
    "## Fixed Effects\n",
    "\n",
    "To make matters more formal, let's first take a look at the data that we have. Following our example, we will try to estimate the effect of marriage on income. Our data contains those 2 variables, `married` and `lwage`, on multiple individuals (`nr`) for multiple years. Note that wage is in log form. In addition to this, we have other controls, like number of hours worked that year, years of education and so on."
   ]
  },
  {
   "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>nr</th>\n",
       "      <th>year</th>\n",
       "      <th>black</th>\n",
       "      <th>...</th>\n",
       "      <th>lwage</th>\n",
       "      <th>expersq</th>\n",
       "      <th>occupation</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>13</td>\n",
       "      <td>1980</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.197540</td>\n",
       "      <td>1</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>13</td>\n",
       "      <td>1981</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.853060</td>\n",
       "      <td>4</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>13</td>\n",
       "      <td>1982</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.344462</td>\n",
       "      <td>9</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>13</td>\n",
       "      <td>1983</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.433213</td>\n",
       "      <td>16</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>13</td>\n",
       "      <td>1984</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.568125</td>\n",
       "      <td>25</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 12 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   nr  year  black  ...     lwage  expersq  occupation\n",
       "0  13  1980      0  ...  1.197540        1           9\n",
       "1  13  1981      0  ...  1.853060        4           9\n",
       "2  13  1982      0  ...  1.344462        9           9\n",
       "3  13  1983      0  ...  1.433213       16           9\n",
       "4  13  1984      0  ...  1.568125       25           5\n",
       "\n",
       "[5 rows x 12 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from linearmodels.datasets import wage_panel\n",
    "data = wage_panel.load()\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generally, the fixed effect model is defined as\n",
    "\n",
    "$\n",
    "y_{it} = \\beta X_{it} + \\gamma U_i + e_{it}\n",
    "$\n",
    "\n",
    "where \\\\(y_{it}\\\\) is the outcome of individual \\\\(i\\\\) at time \\\\(t\\\\), \\\\(X_{it}\\\\) is the vector of variables for individual \\\\(i\\\\) at time \\\\(t\\\\), \\\\(U_i\\\\) is a set of unobservables for individual \\\\(i\\\\) that are unchanging through time, hence the lack of the time subscript, and \\\\(e_{it}\\\\) is the error term. For the education example, \\\\(y_{it}\\\\) is log wages,  \\\\(X_{it}\\\\) are the observable variables that change in time, like marriage and experience and \\\\(U_i\\\\) are the variables that are not observed but constant for each individual, like beauty and intelligence. \n",
    "\n",
    "\n",
    "Now, remember how I've said that using panel data with a fixed effect model is as simple as adding a dummy for the entities. It's true, but in practice, we don't actually do it. Imagine a dataset where we have 1 million customers. If we add one dummy for each of them, we would end up with 1 million columns, which is probably not a good idea. Instead, we use the trick of partitioning the linear regression into 2 separate models. We've seen this before, but now is a good time to recap it. Suppose you have a linear regression model with a set of features \\\\(X_1\\\\) and another set of features \\\\(X_2\\\\).\n",
    "\n",
    "$\n",
    "\\hat{Y} = \\hat{\\beta_1} X_1 + \\hat{\\beta_2} X_2\n",
    "$\n",
    "\n",
    "You can get the exact same \\\\(\\hat{\\beta_1}\\\\) parameter by doing\n",
    "\n",
    "1. regress the the outcome \\\\(y\\\\) on the second set of features \\\\(\\hat{y^*} = \\hat{\\gamma_1} X_2\\\\)\n",
    "2. regress the first set of features on the second \\\\(\\hat{X_1} = \\hat{\\gamma_2} X_2\\\\)\n",
    "3. obtain the residuals \\\\(\\tilde{X}_1 = X_1 - \\hat{X_1}\\\\) and \\\\(\\tilde{y}_1 = y_1 - \\hat{y^*}\\\\)\n",
    "4. regress the residuals of the outcome on the residuals of the features \\\\(\\hat{y} = \\hat{\\beta_1} \\tilde{X}_1\\\\)\n",
    "\n",
    "The parameter from this last regression will be exactly the same as running the regression with all the features. But how exactly does this help us? Well, we can break the estimation of the model with the entity dummies into 2. First, we use the dummies to predict the outcome and the feature. These are steps 1 and 2 above. \n",
    "\n",
    "Now, remember how running a regression on a dummy variable is as simple as estimating the mean for that dummy? If you don't, let's use our data to show how this is true. Let's run a model where we predict wages as a function of the year dummy. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>            <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>       <td>    1.3935</td> <td>    0.022</td> <td>   63.462</td> <td> 0.000</td> <td>    1.350</td> <td>    1.437</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1981]</th> <td>    0.1194</td> <td>    0.031</td> <td>    3.845</td> <td> 0.000</td> <td>    0.059</td> <td>    0.180</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1982]</th> <td>    0.1782</td> <td>    0.031</td> <td>    5.738</td> <td> 0.000</td> <td>    0.117</td> <td>    0.239</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1983]</th> <td>    0.2258</td> <td>    0.031</td> <td>    7.271</td> <td> 0.000</td> <td>    0.165</td> <td>    0.287</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1984]</th> <td>    0.2968</td> <td>    0.031</td> <td>    9.558</td> <td> 0.000</td> <td>    0.236</td> <td>    0.358</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1985]</th> <td>    0.3459</td> <td>    0.031</td> <td>   11.140</td> <td> 0.000</td> <td>    0.285</td> <td>    0.407</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1986]</th> <td>    0.4062</td> <td>    0.031</td> <td>   13.082</td> <td> 0.000</td> <td>    0.345</td> <td>    0.467</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1987]</th> <td>    0.4730</td> <td>    0.031</td> <td>   15.232</td> <td> 0.000</td> <td>    0.412</td> <td>    0.534</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod = smf.ols(\"lwage ~ C(year)\", data=data).fit()\n",
    "mod.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how this model is predicting the average income in 1980 to be 1.3935, in 1981 to be 1.3935 (1.3935+0.1194) and so on. Now, if we compute the average by year, we get the exact same result (add the intercept to the parameters of the other years to get the mean for the year)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "year\n",
       "1980    1.393477\n",
       "1981    1.512867\n",
       "1982    1.571667\n",
       "1983    1.619263\n",
       "1984    1.690295\n",
       "1985    1.739410\n",
       "1986    1.799719\n",
       "1987    1.866479\n",
       "Name: lwage, dtype: float64"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.groupby(\"year\")[\"lwage\"].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means that if we get the average for every person in our panel, we are essentially regressing the individual dummy on the other variables. This motivates the following estimation procedure:\n",
    "\n",
    "1. Create time-demeaned variables by subtracting the mean for the individual:   \n",
    "\\\\(\\ddot{Y}_{it} = Y_{it} -  \\bar{Y}_i\\\\)  \n",
    "\\\\(\\ddot{X}_{it} = X_{it} -  \\bar{X}_i\\\\)\n",
    "2. Run regress \\\\(\\ddot{Y}_{it}\\\\) on \\\\(\\ddot{X}_{it}\\\\)\n",
    "\n",
    "Notice that when we do so, the unobserved \\\\(U_i\\\\) vanishes, since \\\\(\\bar{U_i}=U_i\\\\):\n",
    "\n",
    "\\begin{align}\n",
    "Y_{it} & = \\beta X_{it} + \\gamma U_i + e_{it} \\\\\n",
    "\\bar{Y}_{i} & = \\beta \\bar{X}_{it} + \\gamma U_i + \\bar{e}_{it} \\\\\n",
    "(Y_{it} - \\bar{Y}_{i}) & = (\\beta X_{it} - \\beta \\bar{X}_{it}) + (\\gamma U_i - \\gamma U_i) + (e_{it}-\\bar{e}_{it}) \\\\\n",
    "(Y_{it} - \\bar{Y}_{i}) & = \\beta(X_{it} - \\bar{X}_{it}) + (e_{it}-\\bar{e}_{it}) \\\\\n",
    "\\ddot{Y}_{it} & = \\beta \\ddot{X}_{it} + \\ddot{e}_{it} \\\\\n",
    "\\end{align}\n",
    "\n",
    "To be honest, not only do the unobserved variables that are constant in time vanish. This happens to all the variables that are constant in time. For this reason, you can't include any variables that are constant across time, as they would be a linear combination of the dummy variables and the model wouldn't run. \n",
    "\n",
    "![img](./data/img/fixed-effects/demeaned.png)\n",
    "\n",
    "To check which variables are those, we can group our data by individual and get the sum of the standard deviations. If it is zero, it means the variable isn't changing for the same individual. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "year            1334.971910\n",
       "black              0.000000\n",
       "exper           1334.971910\n",
       "hisp               0.000000\n",
       "hours         203098.215649\n",
       "married          140.372801\n",
       "educ               0.000000\n",
       "union            106.512445\n",
       "lwage            173.929670\n",
       "expersq        17608.242825\n",
       "occupation       739.222281\n",
       "dtype: float64"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.groupby(\"nr\").std().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For our data, we need to remove entinicity dummies, `black` and `hisp`, since they are constant for the individual. Also, we need to remove education. We will also not use occupation, since this is probably mediating the effect of marriage on wage (it could be that single men are able to take more time demanding positions). \n",
    "\n",
    "To run our fixed effect model, first, let's get our mean data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>married</th>\n",
       "      <th>expersq</th>\n",
       "      <th>union</th>\n",
       "      <th>hours</th>\n",
       "      <th>lwage</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>nr</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.000</td>\n",
       "      <td>25.5</td>\n",
       "      <td>0.125</td>\n",
       "      <td>2807.625</td>\n",
       "      <td>1.255652</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>0.000</td>\n",
       "      <td>61.5</td>\n",
       "      <td>0.000</td>\n",
       "      <td>2504.125</td>\n",
       "      <td>1.637786</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>1.000</td>\n",
       "      <td>61.5</td>\n",
       "      <td>0.000</td>\n",
       "      <td>2350.500</td>\n",
       "      <td>2.034387</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>45</th>\n",
       "      <td>0.125</td>\n",
       "      <td>35.5</td>\n",
       "      <td>0.250</td>\n",
       "      <td>2225.875</td>\n",
       "      <td>1.773664</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>110</th>\n",
       "      <td>0.500</td>\n",
       "      <td>77.5</td>\n",
       "      <td>0.125</td>\n",
       "      <td>2108.000</td>\n",
       "      <td>2.055129</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     married  expersq  union     hours     lwage\n",
       "nr                                              \n",
       "13     0.000     25.5  0.125  2807.625  1.255652\n",
       "17     0.000     61.5  0.000  2504.125  1.637786\n",
       "18     1.000     61.5  0.000  2350.500  2.034387\n",
       "45     0.125     35.5  0.250  2225.875  1.773664\n",
       "110    0.500     77.5  0.125  2108.000  2.055129"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y = \"lwage\"\n",
    "T = \"married\"\n",
    "X = [T, \"expersq\", \"union\", \"hours\"]\n",
    "\n",
    "mean_data = data.groupby(\"nr\")[X+[Y]].mean()\n",
    "mean_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To demean the data, we need to set the index of the original data to be the individual identifier, `nr`. Then, we can simply subtract one data frame from the mean data frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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>married</th>\n",
       "      <th>expersq</th>\n",
       "      <th>union</th>\n",
       "      <th>hours</th>\n",
       "      <th>lwage</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>nr</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-24.5</td>\n",
       "      <td>-0.125</td>\n",
       "      <td>-135.625</td>\n",
       "      <td>-0.058112</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-21.5</td>\n",
       "      <td>0.875</td>\n",
       "      <td>-487.625</td>\n",
       "      <td>0.597408</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-16.5</td>\n",
       "      <td>-0.125</td>\n",
       "      <td>132.375</td>\n",
       "      <td>0.088810</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-9.5</td>\n",
       "      <td>-0.125</td>\n",
       "      <td>152.375</td>\n",
       "      <td>0.177561</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.5</td>\n",
       "      <td>-0.125</td>\n",
       "      <td>263.375</td>\n",
       "      <td>0.312473</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    married  expersq  union    hours     lwage\n",
       "nr                                            \n",
       "13      0.0    -24.5 -0.125 -135.625 -0.058112\n",
       "13      0.0    -21.5  0.875 -487.625  0.597408\n",
       "13      0.0    -16.5 -0.125  132.375  0.088810\n",
       "13      0.0     -9.5 -0.125  152.375  0.177561\n",
       "13      0.0     -0.5 -0.125  263.375  0.312473"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "deamed_data = (data\n",
    "               .set_index(\"nr\") # set the index as the person indicator\n",
    "               [X+[Y]]\n",
    "               - mean_data) # subtract the mean data\n",
    "\n",
    "deamed_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can run our fixed effect model on the time-demeaned data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>-8.327e-17</td> <td>    0.005</td> <td>-1.64e-14</td> <td> 1.000</td> <td>   -0.010</td> <td>    0.010</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>married</th>   <td>    0.1147</td> <td>    0.017</td> <td>    6.756</td> <td> 0.000</td> <td>    0.081</td> <td>    0.148</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>expersq</th>   <td>    0.0040</td> <td>    0.000</td> <td>   21.958</td> <td> 0.000</td> <td>    0.004</td> <td>    0.004</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>union</th>     <td>    0.0784</td> <td>    0.018</td> <td>    4.261</td> <td> 0.000</td> <td>    0.042</td> <td>    0.115</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hours</th>     <td> -8.46e-05</td> <td> 1.25e-05</td> <td>   -6.744</td> <td> 0.000</td> <td>   -0.000</td> <td>   -6e-05</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod = smf.ols(f\"{Y} ~ {'+'.join(X)}\", data=deamed_data).fit()\n",
    "mod.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we were to believe that fixed effect eliminates the all omitted variable bias, this model is telling us that marriage increases a man's wage by 11%. This result is very significant. One detail here is that for fixed effect models, the standard errors need to be clustered. So, instead of doing all our estimation by hand (which is only nice for pedagogical reasons), we can use the library `linearmodels` and set the argument `cluster_entity` to True."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Parameter Estimates</caption>\n",
       "<tr>\n",
       "     <td></td>     <th>Parameter</th> <th>Std. Err.</th> <th>T-stat</th>  <th>P-value</th> <th>Lower CI</th>  <th>Upper CI</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>expersq</th>  <td>0.0040</td>    <td>0.0002</td>   <td>16.552</td>  <td>0.0000</td>   <td>0.0035</td>    <td>0.0044</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>union</th>    <td>0.0784</td>    <td>0.0236</td>   <td>3.3225</td>  <td>0.0009</td>   <td>0.0322</td>    <td>0.1247</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>married</th>  <td>0.1147</td>    <td>0.0220</td>   <td>5.2213</td>  <td>0.0000</td>   <td>0.0716</td>    <td>0.1577</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hours</th>   <td>-8.46e-05</td> <td>2.22e-05</td>  <td>-3.8105</td> <td>0.0001</td>   <td>-0.0001</td> <td>-4.107e-05</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from linearmodels.panel import PanelOLS\n",
    "mod = PanelOLS.from_formula(\"lwage ~ expersq+union+married+hours+EntityEffects\",\n",
    "                            data=data.set_index([\"nr\", \"year\"]))\n",
    "\n",
    "result = mod.fit(cov_type='clustered', cluster_entity=True)\n",
    "result.summary.tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how the parameter estimates are identical to the ones we've got with time-demeaned data. The only difference is that the standard errors are a bit larger. Now, compare this to the simple OLS model that doesn't take the time structure of the data into account. For this model, we add back the variables that are constant in time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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>    0.2654</td> <td>    0.065</td> <td>    4.103</td> <td> 0.000</td> <td>    0.139</td> <td>    0.392</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>expersq</th>   <td>    0.0032</td> <td>    0.000</td> <td>   15.750</td> <td> 0.000</td> <td>    0.003</td> <td>    0.004</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>union</th>     <td>    0.1829</td> <td>    0.017</td> <td>   10.598</td> <td> 0.000</td> <td>    0.149</td> <td>    0.217</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>married</th>   <td>    0.1410</td> <td>    0.016</td> <td>    8.931</td> <td> 0.000</td> <td>    0.110</td> <td>    0.172</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hours</th>     <td> -5.32e-05</td> <td> 1.34e-05</td> <td>   -3.978</td> <td> 0.000</td> <td>-7.94e-05</td> <td> -2.7e-05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>black</th>     <td>   -0.1347</td> <td>    0.024</td> <td>   -5.679</td> <td> 0.000</td> <td>   -0.181</td> <td>   -0.088</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hisp</th>      <td>    0.0132</td> <td>    0.021</td> <td>    0.632</td> <td> 0.528</td> <td>   -0.028</td> <td>    0.054</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>educ</th>      <td>    0.1057</td> <td>    0.005</td> <td>   22.550</td> <td> 0.000</td> <td>    0.097</td> <td>    0.115</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod = smf.ols(\"lwage ~ expersq+union+married+hours+black+hisp+educ\", data=data).fit()\n",
    "mod.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This model is saying that marriage increases the man's wage by 14%. A somewhat larger effect than the one we found with the fixed effect model. This suggests some omitted variable bias due to fixed individual factors, like intelligence and beauty not being added to the model. \n",
    "\n",
    "## Visualizing Fixed Effects\n",
    "\n",
    "To expand our intuition about how fixed effect models work, let's diverge a little to another example. Suppose you work for a big tech company and you want to estimate the impact of a billboard marketing campaign on in-app purchase. When you look at data from the past, you see that the marketing department tends to spend more to place billboards on cities where the purchase level is lower. This makes sense right? They wouldn't need to do lots of marketing if sales were skyrocketing. If you run a regression model on this data, it looks like higher cost in marketing leads to less in-app purchase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAEXCAYAAAAJJYvtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3hUZfbA8e8hhNBCqAqClAgoUoSoCFhoKrg2FivoqqC4iG1dy89dUVnErrsqoq4oYndVBGyLKCWCIixGiqgYpEmTABISWkI4vz/uzTAZUu4kM5mb5HyeZ57Mbe89czMzZ9573/u+oqoYY4wxflIt1gEYY4wxoSw5GWOM8R1LTsYYY3zHkpMxxhjfseRkjDHGdyw5GWOM8R1LTqbCEZExIrKqsu2rshCRuSLyUpjbrBWR0dGKyVQ8lpyMr4hILRF5QETSRWSviGwXkf+JyC1Bqz0B9IhVjKUlIt1E5H0RyRCRHPcLebyIHBmy3mQR+aKYchqJyDMiskZE9rvlzRORISXsX93H+YUsm+YuCyupGBMt1WMdgDEhngf6ArcCS4F6QDegZf4KqpoNZMckulISkYHANPdxHrAZOB54EFgsIr1U9VePxU0B6gN/BlYCjYFTgEYetl0PjAA+CoqtGfAHwOv+jYk6qzkZvxkEPK6q01R1jaouVdXJqjo2f4XQU2350yJyqVvj2uPWBOqJyGARWSkiWW6tJSlou8ki8oWI/FVENrrbTRGRxsUFKCJnichXbs1uo4i8IiJFJgYRqQVMBmar6uWqulBV16vqDJxEXAOY4OXgiEh9oDcwWlVnquo6Vf1WVZ9T1Wc9FDEJOEdEmgfNuxaYB6wO2Ve8iDzivsYcEflBRIaGrNNKRGa4x2K9iNxcSMzV3f/RGhHZJyIrROTPXl6vqbosORm/2QwMFJGGYW7XDLgauAg4BzgVeB+4DrgUp2ZwOvD3kO26A32Age46XXC+wAslIv2A6cA77rqDgNbAVBGRIjY7GzgSeCh0garuAp4FzhWRBiW+SqfGmAVcKCJ1PKwf6hcgFRgGICLVcJLTxELWfQinlvUXoBPwBvCGiPR3txVgKk6NrQ9wgftICSnnJWAwTk2vAzAWeFREri1F/KaqUFV72MM3D5yksg7IA5YBLwIXAhK0zhhgVcj0AaBx0LwJbhlNguY9DSwOmp6M82WfFDTvbECBdkXsay7wSEjMLd1tuhbxmu5ylzcoYvlgd/nJQXF9Ucwx+iOwDcgBFruvq5+HY6vAlTjJeg3Oj9OBblk13Nf2krtubWA/MCqkjKk4NUCAM90y2wctbwLsDSqnDXAQOC6knPuAJUHTa3FqgzF/D9rDHw+rORlfUdWvgGNwajmv4tQ4pgAfFlMzAdioqtuCprcAW1Q1I2TeESHb/aCqmUHTX7l/OxSxn5OBv4hIdv4D+MFd1q6IbYqLO2yqOhVojpNYpuBcu5olIp5ODeIkmNo4yeV64FVVzQlZpy1OwvoyZH4q0NF9fjywTVV/DootA+c6WL6TcF7/4pBj9neKPl7GWIMI4z+qegD42n08KSJXAq8DZ+B8ORYmN7SYIuaV9QdZNeBRN55QW4rYJv/LuhPOtZ1QHXFqF56brKvqfmC2+3jYbYb9gIg8rqprS9g2V0ReBe4BeuKcnixy9ZBpCZonhSwPlX+8ewF7SijbmACrOZmK4Ef3b2itJxI6iEi9oOleIfsMtRjoqKqrCnkU1YJwJpAB/C10gbvvm4BPVPX3Ur6G4HibeFz/RZza6Teq+lMhy1fhnNbrHTL/DGCF+3wF0EREAjUgtzFJ+6D1v3X/tizkeP3iMVZTBVnNyfiKiKQCb+MkgQyc00sPATuBOVHYpQKvuTWPhjjXqj5R1fQi1r8PmCki/8I57ZiFc3rqEuAmVd172A5U94jIMOADEXkb+BcFm5LnADeGbFZXRLqGzNuHc0ymAK/gNLXfiVMjexjnOtISTy9adZWbSPYVsXyPiDyDUxvLcMu9BOf631nuarPcGN5wW+nl4NQqD4TsZxIwUUTuAhYAdYATca4HPuolXlP1WHIyfvNf4AqcFl31gK041z2GhVxTipRFwHzgc5x7h2bgXIcplKrOcVvs3Y9ziq4azr1Dn3H4acTg7T4RkV44p9I+xXltm3HuN3pAVX8L2eQU4LuQeSuBE3BOd96Ik7hrueXMBB5U1SJjKCSmHSWscg/O6cancGpkq4ArVXWWu72KyCCcWtiXOA0rHgcSQsq5HrjdLS8Z2IVT6/LS9N1UUaJqp31N1SQik4EWqnpmrGMxxhRk15yMMcb4jiUnY4wxvmOn9YwxxviO1ZyMMcb4jq9b62VmZlq1zhhjKrmkpKTDelGxmpMxxhjfseRkjDHGdyw5xVB6elGdEJjC2PEKjx2v8NjxCk+0j5clJ2OMMb5jyckYY4zv+Lq1njHGqCrZ2dkcPHgwqvupWbMmmZmZJa9ogPCOV7Vq1ahbty7FD8lWkCUnY4yvZWdnk5CQQI0aNaK6n4SEBGrWrBnVfVQm4RyvnJwcsrOzSUxM9Fy+ndYzvrcuK5cRqTsYuSyBEak7WJflueNtUwkcPHgw6onJRFeNGjXCrvlazcn42rqsXAZ9tp01WXlAHN/u2svijBymDWhEq8T4WIdnjIkSqzkZXxuXluUmpkPWZOUxLi0rRhGZqqhhw4acdtpp9OzZk8suu4ydO3fGOqSABx98kLlz55a5nHnz5nHZZZcdNv/mm2/mp58KGyw5uiw5GV/bvCev0PlbiphvTDTUqlWL+fPns2DBAho0aMBLL71U5jLz8iLzHr7nnnvo06dPRMoqzPjx4znuuOOiVn5R7LSe8bVmteMKnd+0iPmm8nu61RMRLe/WdXeEtX737t1ZsWIFAGvWrOGOO+5g27Zt1K5dm6effpr27duzZs0aRowYQV5eHmeeeSbPPfccGzduZN68eTz66KM0bdqU5cuX8/XXXzNmzBjmz5/P/v37GTFiBMOGDWPLli0MGzaMrKws8vLyePLJJznllFO46aabWLJkCSLCFVdcwY033sgNN9zAwIEDufDCC0lNTWX06NHk5eXRrVs3/vnPf5KQkEDnzp0ZMmQIM2bM4MCBA0yePJn27dt7er3nnnsu48aNo1u3bjRv3pyRI0cyY8YMEhISeOeddzjiiCPYtm0bt912Gxs2bADg4YcfpkePHuH9I0JYzcn42uiURNokFkxEbRLjGJ3ivdWPMZGSl5dHamoq55xzDgC33norjz32GKmpqTzwwAPcfvvtANx9992MHDmSOXPm0KxZswJlpKWlMXr0aBYuXMjrr79OvXr1mDNnDnPmzOHVV19l7dq1vPfee/Tv35/58+czf/58OnfuzPLly9m8eTMLFizg66+/5oorrihQ7r59+xg1ahSvvPIKX3/9NQcOHODll18OLG/UqBFffvklw4cPZ/z48aV6/bt37+akk07iq6++okePHrz66quB1ztq1CjmzJnDa6+9xi233FKq8oNZzcn4WqvEeKYNaMS4tCxWb88muVFdRqckWmMIU6727t3Laaedxvr16+natSt9+/YlOzubRYsWcfXVVwfWy8nJAWDRokW8+eabAFx88cXce++9gXVSUlJo3bo1ALNnz2bFihVMnz4dgF27drF69WpSUlK46aabyM3N5dxzz6VLly60bt2atWvXcueddzJgwAD69etXIMb09HRatmxJ27ZtARg6dCgTJ05k1KhRAJx//vkAdO3alY8++qhUx6FGjRoMHDgQgC5duvDVV18BMHfu3ALXpbKyssjKygqr6XgoS07G91olxjOxd0PS07fTrl3LWIdjqqD8a06ZmZlcfvnlTJw4kaFDh5KUlMT8+fPDKqtOnTqB56rKY489Rv/+/Q9b79NPP2XmzJmMHDmSm2++mSFDhjB//nxmzZrFxIkTmTp1KhMmTChQVnESEhIAiIuL48CBA2HFnC8+Pj5wI21wOQcPHuTzzz+nVq1apSq3MJacjDEVSrjXiCIpKSmJRx55hCuuuIJrr72WVq1aMW3aNAYNGoSq8v3339O5c2dOPvlkPvzwQwYPHswHH3xQZHn9+/fn5Zdf5owzziA+Pp5Vq1bRrFkztm/fzlFHHcXVV1/N7t27Wbp0KWeffTbx8fFceOGFtGnTJlAjyte+fXt+/fVXVq9eTXJyMu+88w6nnnpqtA8JAP369WPixImB03nLli2jS5cuZSrTkpMxxoThhBNOoGPHjkyZMoUXX3yR22+/nccff5wDBw4wePBgOnfuzMMPP8z111/Ps88+y9lnn029evUKLeuqq65i/fr19O7dG1WlUaNGvPnmm8yfP5/x48dTvXp16taty/PPP8+mTZu48cYbAzez3n///QXKqlmzJhMmTODqq68ONIgYPnx4WK8tNTWV448/PjA9efJkT9s9+uij3HHHHfTq1Yu8vDx69erFv/71r7D2HUpKqgrGUmUfCTc9PZ127drFOoyAdVm5jEvLYvOePJrVjvPdtR2/HS+/qyzHKzMzk6SkpKjvZ9++fRHrvmjPnj3UqlULEWHKlCm8//77vP322xEp2y/CPV7F/R8LGwnXak4GCO2JwWE9MRhTOkuWLOHOO+9EVUlKSipwbch4Y8nJAMX3xDCxd8MYRWVMxdSrV69ASzZTOnafkwGsJwZjjL9YcjKA9cRgjPEXS04GsJ4YjH9Vq1YtcHOrqZhycnKoVi28dGPXnAxQsCeGLXvyaOrD1nqmaqpbty7Z2dns3bs3qvvZtWtXkU2+zeHCOV75I+GGw5KTCcjvicEYPxGRMnWD49XWrVs5+uijo76fyiLax8tO6xljjPGdcklOIjJJRLaKyPeFLLtDRFREGpdHLMYYY/yvvGpOk4GBoTNF5GjgLGB9OcXhC+uychmRuoORyxIYkbqDdVm5sQ7JGGN8pVyuOanqlyLSupBF/wLuAqaXRxx+ULAnhji+3bXXemIwxpgQMbvmJCIXABtVdWmsYoiF4npiMMYY44hJaz0RqQ3cA5ztdZv09PToBVSOVm9LAA6/sXX19mzS07eXf0AVTGV5H5QXO17hseMVnrIcr5I6JY5VU/JjgDbAUnfgqhZAmoh0V9UthW1QGXpXBkjetINvdx1+v0Zyo7o2kF4JKksv2+XFjld47HiFJ9rHKyan9VR1uaoeoaqtVbU1sAFIKSoxVSbWE4MxxpSs2OQkItVFZLCIvCwii0Vklfv3ZRG5WEQ81bxE5G1gAXCsiGwQkWsjEXxFlN8TwyXJtTgxKY9LkmtZYwhjjAlRZHIRkT/jXBf6EUgFPgaygESgAzAC+KeIPKSqLxS3E1UdUsLy1uGFXbHl98SQnr7dTuUZY0whiqv5tAeKugY0FXhIRJoBt0clMmOMMVVWkclJVUtMOqq6GbgjohEZY4yp8kq8ZuReV+oHdMQ5pZcFrABmq+qB6IZnjDGmKio2OYnICTi9NwiwDMgE6gG3AioiF6rqsqhHaYwxpkopqeb0EvCkqo4PXSAiNwGTgJOiEZgxxpiqq6T7nI4HimqJ9yJOqz1jjDEmokpKTj8CNxSx7M/ucmOMMSaiSjqtdx0wTUTupOA1py5AHjAouuEZY4ypiopNTqq6RETaAX1xTvHVBbKBp4G5qmoDERljjIm4EpuSuwlopvswxhhjoq7Ejl9F5HoR+VpEMkUkz/37tYiMKI8AjTHGVD0l3ef0KHAe8CSwlEPXnLoCfxWRZFX9W9SjLIN1WbmMS8ti8548mtV2ev+2TlaNMcbfSjqtNxzo4nZTFCxNRGbgNJLwbXIqOCS6w4ZEN8YY/yvptJ6UcXlM2ZDoxhhTMZVUc3oZmC0ioaf1TgD+CkyMbnhls3lPXqHztxQx3xhjjD+U1JT8/0RkNTAMp+PX/KbkK4BnVPXf0Q+x9JrVjit0ftMi5htjjPEHL03J/w34OgkVZXRKIoszcgqc2rMh0Y0xxv9KO2TG98Acvw+ZkT8k+ri0LLbsyaOptdYzxpgKodIPmZE/JLoxxpiKw4bMMMYY4zs2ZIYxxhjfsSEzjDHG+I4NmWGMMcZ3vA6Z0YeC9znZkBnGGGOixuuQGZ+7D2OMMSbqShwyoygiUl1EJkUyGGOMMQbKkJyAOODqSAVijDHG5CvpJtzZxSy2DuqMMcZERUnXnE4BHgZCx3MCiAdOi3hExhhjqrySktMS4CdVfT90gYgkAM9FJSpjjDFVWknXnJ4CdhSxLBdnKI0SicgkEdkqIt8HzXtcRH4SkWUiMlVE6nsL2RhjTGVXbHJS1fdUtdDrTqp6UFVf9bifycDAkHmfA51UtQvwM1Ea7n3l9B/5ZOR0nm71BF89No+c3TnR2I2pgtZl5TIidQfn/TeDEak7WJdlt/0ZEykl3ucUCar6pYi0Dpk3M2jyG+DiSO9387ebmHHLJ4HpxRMWsnjCQgCadmtGn7H9ObJL00jv1lQB67JyGfTZ9gJjhS3OyGHagEY2JIsxEVAuycmD4cB/Il3ozrW/F7lsy3ebeef8NwLTp/29N12HpRBXwxohmpKNS8sqkJgA1mTlMS4ty4ZoMSYCRFXLZ0dOzeljVe0UMv8enGE3BmtIMJmZmYHp9PT0sPeZm5XDgtvmk702K6ztmnQ/kg6jOlG3Rd2w92mqhpHLEvh21+E/ZE5MyuOFzvtjEJExFUu7du0Cz5OSkiR0eUxrTiJyNXAe0D80MYUKfiHhaD/zWL648zNWTvfegXrGot/IWPQbANVrVqfP2P4cf0knpNphx69M0tPTS/26qiI/Ha/kTTv4dtfew+c3qku7di1jENHh/HS8KgI7XuGJ9vHynJxEJAk4Fqfz14CiGkx4KG8g8H9Ab1XdU5oyvKieUJ2Bz5zLwGfO5WDeQX5493vm3jeLvJy8kjcGDuw7wBd3fcYXd30GwHF/PJ5T/3YGdY+0WlVVNjolkcUZOQVO7bVJjGN0SmIMozKm8vCUnETkGmACTo/kwYlEgWQP27+N07N5YxHZANyP0zovAfhcRAC+UdWRYcQetmpx1eg0pAudhnQBYMeq7Xw5dg7rUtd6LuOnqT/w09QfAKjbtC59xvYn+ey2uK/BVBGtEuOZNqAR49Ky2LInj6a1ncRkjSGMiQyvNacHgYtV9b+l2YmqDilk9sulKSuSGrZtxKDXnEaCB/YfYMkraXz18Jeet8/eks3H108PTHcdlsIpt/WiZlLNiMdq/KdVYrw1fjAmSrwmp+rAzBLXqsCqJ1TnpJHdOWlkd8BpzTd3zCx+W7LFcxlLXkljyStpADQ6rjF9x/an+SlHRyVeY4ypzLwmp0eB0SLygKoejGZAftG0WzMun34lADnZOSx69hu+fX6R5+23/7SN9y891Dr+lNt6ceKfTya+lp32McaYknhqSi4ivwJNgRxge/AyVY1a06TgpuR+si51DXPum0Xm2p2l2r55jxb0vr8fO+MzrXVQGKw1VXjseIXHjld4Inm8ytKU/MqIRFBJtOrdhmtSrwNgz7bdLHjiK75/e5nn7Td+s4G3znktMN17TD+6/Kkr1aqXZXgtY4ypPMrtJtzS8GvNqSiqys8frWTuvV+wb+e+UpVxzIB2nD66N0ktrR/cUPbLNjx2vMJjxys8Mas5icg9qvqg+3xsUeup6n0Ria4SEBGOveA4jr3gOAB2/ZrJvAfnsuq/3nu3+OWzdH75zFm/Rr0E+o7tz7GDOlhTdWNMlVLcab0WQc+tyVkp1Ds6iXNfuBCAgwcOsvytpcy9d5bn7XN27eezv3zKZ3/5FIDjL+1Er7tOp06TOlGJ1xhj/KLI5KSqNwQ99zRukylaterVOOGqbpxwVTcAtv2UwYz/+4TtS7Z5LuOHd7/nh3edIbHqHZ1E3wf607pvifdAG2NMhVPcab0jVHVrSQWIyJGq+ltkw6r8Gh/XhFOeOJV27dpxYF8uaS8uZsGTX3neftevmUy/5oPAdMr1J9H9lp4kJCZEI1xjjClXxZ3WmyMiqcDrwMLg+5tEpBrQHbgKOAPoVHgRxovqNePpfktPut/SE4CN/9tA6n2zyfihxN8GAWkvLibtxcUAHNH5SPr8oz/NTjwqKvEaY0y0FZecugHXAy8CySKyGsgCEnH600sH/g38JdpBVjXNT27B0P9eBcD+zH0sfHoB3738reftty7/jXcHvxWY7nXX6XS77kSqJ/hl+C5jjCme15twjwY6A/WB34FlqroxyrFVuKbk4Qq3Kaaqsmb2aube+wVZG8Mboypfy9NbccZ9fWnUvnGpto8la+obHjte4bHjFR5f3ISrqr8Cv0YkClNqIkJy/2NI7n8MANm/ZfP1Y/P48f0VnstYP28db5w12SkvTuj7wJl0vLwz1eLsBmBTNa3LymVcWhartyWQvGmH9S7vE3YTbgxF8peHHlR+mvoDc+79gtzduaUqo/35x3Lq33pTr3m9iMQUafbLNjx2vEq2LiuXQZ9tP2xcrmkDGlmCKoEvak7G/6Sa0OGijnS4qCMAv6/5nXnj5rLmi188l/HzRyv5+aOVgel25x3LOc+eZzcAm0prXFpWgcQEsCYrj3FpWTYcSoxZcqqkGrRpwAUv/xGAvNw8lr22hC/HzgmrjPSPV5L+8aFkddXs4TQ4xj6wpvLYvKfwEbG3FDHflJ+wkpPbhPxIVd0cpXhMFMTFx9Ht2hPpdu2JAPy2fAupY+aweXF4bVpe6zcp8PykG0/h1LtOj2icxpS3ZrXjCp3ftIj5pvx4Haa9PvAccDGQC9QRkQuA7qo6OorxmSg4snNTLp3iDE6cuyeHxc8tYtH4b8IqY/GEhSyesDAwfd3iG6xbJVPhjE5JZHFGzmHXnEanJMYwKgPem5K/g9OEfCzwg6o2EJEmwNeqGrUrrtYgovytnP4jM275pExl9Hv4LDoPPSFCER3ix+PlZ3a8vAm01tueTXKjutZaz6NoN4jwmpwygKNUNVdEdqhqQ3d+pqomRSS6Qlhyiq3fV+/gtb6TSl6xGHWOrMvVc4cTX7tGmePx+/HyGzte4bHjFR6/tNbLBBoDgWtNItIyeNpUPg2SG3LrujsA5wbgGbd8ws8f/hRWGbt/y+a5Ds8Epi+Y9EfauPdpGWNMUbwmp5eAKSJyD1BNRHoCDwEvRC0y4ysiwjnjz+Oc8ecBsGnxRt676O2wy/lw+NTA8+antGDwW5faCMDGmMN4TU6PAvuACUA8MAmnX72noxSX8bmjTmoeqFXl5ebx7uC32LosvM7pNy7cwPhj/hmYvvyjKzmyS9OIxmmMqZi8dl+kwFPuw5gC4uLjGPLRnwLT6Z/+zKc3fBh2Oe+c/0bg+fGXduLMxwbYDcDGVFFem5L3Bdaq6hoRaYpTk8oD/q6qW6IZoKl42v2hfaBWtT9rP6/0epH9u/aHVUbwwIoA18y7LqIxGmP8zetpveeAAe7z/PMwB3CG07gg0kGZyiMhMYGRy28OTH836Vu+/Ed4PVUATD79pcDzTkO70P/hsyMSnzHGn7wmp+aqul5EquMkqVZADrApapGZSqnb8BPpNtzpqSJ7SxYvn/LvsMv4/q1lfP/WssC03QBsTOXjNTntEpEjcUa8/UFVs0WkBk7jCGNKpW7TxMDpP4DUMbNZ8kpa2OW8dNLzgeen39OblOtPjkh8xpjY8ZqcxgP/A2pwaOTbU4Hwbnoxphi9x/Sj95h+AGxbmcGbZ78adhnzHkxl3oOpgelRP91KfC37DWVMReO1td6jIjIVyFPV/DEYNgJ2ldpEReNjmxy6AfigMu2q91k/b13Y5Tx33KG7Hc594QLantM+YjEaY6LHc6/kqvpzcdPGRItUE/74xiWB7lJW/fdnPhkZflP14G0SmycybP71SDVrqm6MH3ltSl4PGAP0xunGKPCJVtWWHrafBJwHbFXVTu68hsB/gNbAWuBSVf09rOhNldT2nENN1XP35haoHXmVtTGLZ9o8GZi+bNoVNO3WLGIxGmPKJpym5C1weiV/A7gSuBOY4nH7ycCzwGtB8+4GZqnqIyJytzv9fx7LMwaA+FrxBRpVpE1czLxxc8Mu5z+D3gw8b90vmQtfGRyJ8IwxpeQ1OZ0NdFDV7SKSp6rTRWQx8BHwr5I2VtUvRaR1yOwLgT7u81eBuVhyMmWUMuIkUkacBMCebbuZeOLzJWxxuLWzV/N0qycC09fMu46klvUjFqMxpmRek1M1nJ7JAbLdwQc3A23LsO/AiLqqullEjihDWcYcpnbjOgVqVbP+/jnfv7k07HKCbwA+ceTJnPa33hGJzxhTNK/jOc0CHlLVWSLyNnAQyAZOVNWTPO3IqTl9HHTNaaeq1g9a/ruqNgjeJng8p/T0dC+7McaT7PVZfDl8dpnLOWvqOcQnln2sKmOqmuCxoMoy2GCyu+4v7gi4DwOJwD9U9QcvgRSSnFYCfdxaUzNgrqoeG7yNDTZogkXreKkq71/6HzYt2lCmcs7994W0Heif/6e9v8Jjxys8vhhsUFVXBz3PIDL3N30IXA084v6dHoEyjQmbiHDJe5cHpjcsWM+Uy98Nu5xP/nzoLXxE5yO5dOpQ4uLjIhKjMVWN5/ucRORsoCtQN3i+qt7nYdu3cRo/NBaRDcD9OEnpXRG5FlgPXOI9bGOip0XPloFrVQcPHOTFbhPC7lV96/LfeLbtobZCl34wlGYnHhXROI2pzLze5/QscCkwB9gT7k5UdUgRi/qHW5Yx5ala9WoFelX/8YMVzLztv2GX8+7gtwLP259/LAPHn2djVRlTDK81pyFAV1X9NZrBGON3HQZ3pMPgjoAzVtULncaHXcbPH63k549WBqavmnstDdo0KGYLY6oer8lpO7AzmoEYU9EkJCYUaKq+/K2lzP7b52GX81qflwPPu9/cg553nBaR+PxmXVYu49Ky2Lwnj2a14xidkkirROuU1xSuyOTkttDL9yTwpog8DPwWvF5wYwljqrLOQ0+g89ATANidsbvAUB5eLRr/DYvGfxOYHpE2itqNakcsxlhZl5XLoM+2syYrLzBvcUYO0wY0sgRlClVczWkVoAT1o4fTP14wBaw5kjEh6jQpeAPw14/N438TFoZdzsSU5wLP+z9yNp2GdIlIfOVtXFpWgcQEsCYrj3FpWUzs3TBGUfVDoeQAABpZSURBVBk/KzI5qWq18gzEmMqs112n0+uu0wH4ffUOXus7KewyZt09k1l3zwSg7lGJXDV7eIUZq2rznrxC528pYr4xXlvrNQf2BPcaLiINgFqqakO1GxOGBskND41VpcqMmz8u0EDCi+xNWQV6Y79w8mBa900uZovYala78BMsTYuYb4zXBhHTgOFA8JAWLYCXgFMiHZQxVYWIcM6z53POs+cDsOl/G3nv4rfDLmf6NR8EnrfoeTR/fNNftw2OTklkcUZOgVN7bRKdRhHGFMZr90WZqprkdX6kWPdFJlhVO155OXm8O/gtti7/reSVizHk4z9xROcjIxRV6eW31tuyJ4+mPmytV9XeX2Xli+6LgAwRaauqq/JniEhbnCbmxpgoiKsRx5CP/xSYTv/0Zz69IfwRgN8+7/XA846Xd6b/I2fH5AbgVonx1vjBeOY1OU0CpojIPcBq4BjgAZzTesaYctDuD4dGAN6/az+TTn2RnDC7VVrxznJWvLM8MD3sqxHUaxG1kx/GlJrX5PQIkAs8ARyN0xfey8A/oxSXMaYYCfUSuCGoW6W0lxYz74G5YZfzyqkTA897/d/pnDzKLiEbfyjxmpOIxOHUnK5X1fB+ppWRXXMycOhaxept2SQ3ruu7axV+k7U5i0k9/l2mMqrFV2PEt6OomVQzQlH5n30evYnG57Es4zltBlqqam6ZIgiTJSdTWM8CbRLjrGcBD/LfX3Pvn8XSyd+VqawBT5/LcYM6RCgyf7LPY8mi9XksS3K6C6gPjFHVnFJHECZLTmZE6g7eW733sPmXJNeyi+slKOz9lfFjBm8NfLVM5TZs25Ahn15F9QTPI+5UCPZ5LFm0Po9laa13M9AU+KuIZOB0WwSAqrYsdUTGlMB6FoisJh2aHLoB+KDy4fAPWDtnTVhl7Fi1gwntnwpMX/TOpbToaV8DVUF5fh69JqcrI75nYzywngWiR6oJF06+KDC9fv46pl7xXtjlBI8a3ObMYzh/4iCkmo1VVRmV5+fR02m9WLHTesauOZVeWd5fB/Yd4K0/vMbvv+woUwxXfn4Njdo3LlMZ5cU+jyXz4zWnsUUt8zJMe2lZcjIQ1DpoezbJjay1nleRfH/9OGUFM/8a/gjAwbpddyJn3Ns3IvFEg30evYnG57EsyemVkFlNgd7AVFW9okxRFcOSkwlmxys80Tpee3/fy8SU59CDZft4XrtoJHWPrBuhqMrO3l/h8UX3Rao6LHSeiAzEGb7dGFOF1GpQi1vW3B6YXvTMAhY8+VXY5bzc/YXA8z7/6McJ16REJD5TOZSlLehM4D+RCsQYUzF1v6Un3W/pCUDm+p1MPj38Xs3m3j+buffPBqBmg1oMmz+CGnVrRDROU7F4Hc8pdKCY2sBQ4NeIR2SMqbCSWtYvMFbVF3d+xg/vfR9WGft+38vzHZ8JTJ/34oUcM8BOt1U1XmtOq0Km9wDfAVdHNhxjTGUhIpz1xEDOemIgAFuWbuY/F7wZdjkfXz898Lxpt2Zc/N7lxMXbrQSVnddrTjZkuzGmTJqe0CxQqzp44CBThrzLpkUbwipjy3ebebbtvwLTl04dSrOUoyIap/GHYpOTiBwJ/AvoBKQBt6uqjeFkjCmTatWrccl7lwemV3/xCx9dOzXsct7941uB58cO6sCAp/4Qk7GqTOSVVHN6DqdPvReAi4CngD8Vu4UxxoQp+cxjArWqnN05vNZ3Ert/yw6rjJXTfmTltB8D01enXkv91g0iGqcpPyUlp9OB9qq6U0TexbnOZIwxUVOjTg2uWzQyML3sjSXMueeLsMt5tffLgeen/KUnPW47NSLxmfJRUnKqqao7AVR1m4jUKYeYjDEmoMuVXelyZVcAdmfs5qWTng+7jIVPLWDhUwsC0yPSRlG7Ue2IxWgir6TkFC8iw4D8k7gJIjI8eAVVnRSVyIwxJkSdJnUCp/8A5j+cyrcv/C/sciamPBd4fubjA+h4aeeIxGcip9jui0RkLkHDYxRCVbVfpIPKZ90XmWB2vMJT1Y7XjlXbeb1/aE9r4el+a0+639zDmqp7ENPui1S1T0T2XAwRuQ24DicJLgeGqeq+aO/XGFO5NGzbqMANwJ+O+ohVn/4cVhmLnl7Aoqed03/JZx3D6ff2pX6r+hGP1ZQspkNZikhz4BbgeFXd6za6uByYHMu4jKnIAr1Gb0sgedOOKtmLu4hw7vMXBKY3/m8D71/8TlhlrP78F1Z//gsA8XXi6fvAmRz3x+NtrKpy4odxlqsDtUQkF6dbpE0xjseYCqvgeDtxfLtrL4szcqr8+FfNT24RqFXl5eTxn0FvkrFiq+ftc3fnMvOv/w0MG9Lhko6cetcZ1DnC2ohFS0yTk6puFJEngPXAXmCmqs6MZUzGVGTj0rIKDAQHsCYrj3FpWUzs3TBGUflLXI04hn56VWA6/ZOVfDrqo7DK+PG9Ffz43goAEpsn0ueBM2nTL9luAI6gmI6EKyINgCnAZcBO4D3gfVV9Awo2iEhPT49JjMZUJCOXJfDtrsMv5p+YlMcLnffHIKKKJS8njzXv/8LPk34seeUitL7oGNr96Vji61bdmqoXwY0pSj3YIICItAMuBY7COfX2nqqGd7Xx8DIvAQaq6rXu9FVAD1UdBdZazxRkx6tkI1J38N7qvYfNvyS5ltWcSlDY+2vT4o2k3j+brd//Vqoymxx/BH3G9ueok5tHIkRf8cVggyIyFHgR+ARYB3QG7haRP6vqW8VuXLz1QA8RqY1zWq8/sLgM5RlTpY1OSWRxRk6BU3ttEuMYnZIYw6gqrqNOas6QT5we2/bv2s+i8QtIe9H7V1TGD1t57+K3A9M9bz+VlOtPpnpNP1zu9zevw7SvBq5R1S+D5p0OvK6qrcsUgMg/cE7rHcDpHuk6Vd0PVnMyBdnx8ibQWm97NsmN6lbJ1nqlEe77a83s1cy99wt2bdhVqv216NWS3mP60vjYJqXaPtaiXXPympwygKNUNTdoXjywSVWjdmQtOZlgdrzCY8crPGU5Xru37ubrx+fxw7vhDawYINBnbH86Dz2BatUrxghFvjitB/wTeEhE7lXVfSJSC/iHO98YY6q0OkfU4azHB3LW4wPRg8rKaT8y574vyMnK8VaAwtx7ZzH33lkAtP1De06/pzf1WiRFMWp/85qcRgFNgVtF5HegAU5/e5tF5Ib8lVS1ZeRDNMaYikOqCccNPp7jBh8PwM51O5n/4Fx++Sx0QPGirfr050DvFjUb1KLP2P60P//YKtVU3WtyujKqUZhKJ/+6x+Y9eTSrHWfXPUyVVb9Vfc57cRAAebl5LH9jKaljZnveft/ve5lx88fMuPljADoN7ULPO06r9L2qex2mPTXagZjKo2AvBQ7rpcAYiIuPo+uwFLoOSwEgY8VW5o6ZHdZw9d+/tYzv31oGQFLr+vR94ExandE6GuHGlNem5DWA0cAQDt3n9A7woHXSakJZLwXGeNOk4xGB4epz9+by7QuLCow7VZLMtTuZ9qf3A9MnjerOyTf2oEbdGhGPtbx5Pa33PHAsTiet64BWwN+A5sDwYrYzVdDmPXmFzt9SxHxjDMTXiqfHbacGRuzdsGA9c++bxfaft3suY/Fzi1j83CIAjuzalD5j+9P0hGZRiTfavCanQcAx+aPiAj+IyEJgFZacTIhmtQsfC6dpEfONMYdr0bMlV34+DIB9O/fyzT+/Zumr33ne/rclW/jPBW8Gpk/72xl0HX4icTUqxufQa3LagtNj+M6gebWAzRGPyFR41kuBMZFVs77TYq/P2P6oKr98toq5937B7q27PZcx/+Evmf+w049Cq96t6X1/Pxoc49/T7F6T0+vADBEZD2wAjgZuBF4TkcBIuKrqvQmKqbRaJcYzbUAjxqVlsWVPHk2ttZ4xESMitB3YjrYDnRtgszZn8dXDX7JyuvfOatelruW1fpMAiEuIo8/YM+l4aSdfjVXltYeINR7KUlVNLntIh1gPESaYHa/w2PEKT2U4XgfzDvLDeyuYe98X5O0v3TXeYy/swKl/O4PEZsWf6fBFDxGq2iYiERhjjImaanHV6HR5Zzpd3hmAHau28+UDc1k310v9wrFy+o+BWlidI+rQ54EzOWZA23K/Adi6xjXGmEqqYdtGDHr1IgAO7D/A0lfSAtedvNi9dTef/Hl6YPqEq7vR46+9qFm/VsRjDeX1Pqd6wBigN9AYp+siwLosqkysVwdjKqfAZ7tVG5q91pbRKYkk/LKN1DGz2fKd93ZtS1/9LtBisFFKE5L/kxy11n9ea07PAS2AscAbON0Z3Ykziq2pBKxXB2Mqp6I/2425bNoVAORk5/C/5xayeMJCz+VuT8tgypB3ueT9y6Nyys9r3+xnAxep6nQgz/17GfCniEdkYqK4Xh2MMRWXl892jbo1OPWu07l13R3cuu4OBr1+MfXbNCix7M2LN5JVyvGsSuK15lQNyHSfZ4tIfZx7nNpGJSpT7qxXB2Mqp9J8tlud0Zqr514LwJ7te1jw5Fd8/+bSw9ar16Ieic3rRSbQEF6T01Kc602zgHnABCAb+DkqUZlyZ706GFM5lfWzXbtRbfo/dBb9HzoLVSX9k5Usf2MptTvU5YxRfaJ2b5TX03ojgLXu81uAvUB94KooxGRiYHRKIm0SC75ZrVcHYyq+SH62RYT25x3HRe9cRtuh7anTpE6kwjyM1/ucVgc9zwCui1pEJiasVwdjKqeK+tkO+z4nEVmuqp2jEYyJrVaJ8TakhTGVUEX8bHs9rResVcSjMMYYY4KUJjn5p2dAY4wxlVJpui86J+JRGGOMKaCq99gSVnISkSOATSIS6H08uLGEMcaYsrMeWzye1hORgSKyEefG21VBj/QoxmaMMVWS9dji/ZrTBOABoK6qVgt62B2axhgTYdZji/fTeg2Af6uXkQmNMcaUifXY4r3m9DIwLJqBGGOMcViPLd5rTj2AW0TkbmBL8AJVPSPiURljTBVWUXt1iCSvyekl92GMMaYcVMReHSLJa996r0Y7EGOMMSZfsclJRPqVVICqzi5LAO7YUC8BnQAFhqvqgrKUaYwxpmIrqeb0cgnLFUguYZ2SPA3MUNWLRaQGULuM5RljjKngik1OqtommjsXkXrAGcA17v5ygJxo7tMYY4z/labj10hKBjKAV0TkOxF5SUSiN3qVMcaYCkFieV+tiJwEfAOcqqoLReRpYJeq3guQmZkZCC493XpKMsaYyqJdu3aB50lJSYeNdlGaXskjaQOwQVUXutPvA3cXtmLwC6ks0tPTK+XrihY7XuGx4xUeO17hifbxiulpPVXdAvwqIse6s/oDP8QwJGOMMT4Q65oTwM3Am25LvdVYN0nGGFPlxTw5qeoS4KRYx2GMMcY/Yt1azxhjjDmMJSdjjDG+Y8nJGGOM71hyMsYY4zuWnIwxxviOJSdjjDG+Y8nJGGOM71hyMsYY4zuWnIwxxviOJSdjjDG+Y8nJGGOM78S8bz1jSrIuK5dxaVms3pZA8qYdjE5JpFVifKzDMpWEvb/8yZKT8bV1WbkM+mw7a7LygDi+3bWXxRk5TBvQyL5ATJnZ+8u/7LSe8bVxaVnuF8cha7LyGJeWFaOITGVi7y//suRkfG3znrxC528pYr4x4bD3l39ZcjK+1qx2XKHzmxYx35hw2PvLvyw5GV8bnZJIm8SCXxRtEuMYnZIYo4hMZWLvL/+y5GR8rVViPNMGNOKS5FqcmJTHJcm17GK1iRh7f/mXtdYzvtcqMZ6JvRuSnr6ddu1axjocU8nY+8ufrOZkjDHGdyw5GWOM8R1LTsYYY3zHkpMxxhjfEVWNdQxFyszM9G9wxhhjIiIpKUlC51nNyRhjjO9YcjLGGOM7vj6tZ4wxpmqympMxxhjfseQUAyJytIjMEZEfRWSFiNwa65gqAhGJE5HvROTjWMfidyJSX0TeF5Gf3PdZz1jH5Gcicpv7WfxeRN4WkZqxjslPRGSSiGwVke+D5jUUkc9FJN392yCS+7TkFBsHgNtVtQPQA7hRRI6PcUwVwa3Aj7EOooJ4GpihqscBJ2DHrUgi0hy4BThJVTsBccDlsY3KdyYDA0Pm3Q3MUtV2wCx3OmIsOcWAqm5W1TT3eRbOF0fz2EblbyLSAjgXeCnWsfidiNQDzgBeBlDVHFXdGduofK86UEtEqgO1gU0xjsdXVPVLYEfI7AuBV93nrwKDIrlPS04xJiKtgW7AwthG4ntPAXcBB2MdSAWQDGQAr7inQV8SkTqxDsqvVHUj8ASwHtgMZKrqzNhGVSEcqaqbwfnBDRwRycItOcWQiNQFpgB/UdVdsY7Hr0TkPGCrqn4b61gqiOpACvC8qnYDdhPhUy6ViXut5EKgDXAUUEdEroxtVMaSU4yISDxOYnpTVT+IdTw+dypwgYisBd4B+onIG7ENydc2ABtUNb82/j5OsjKFOxNYo6oZqpoLfAD0inFMFcFvItIMwP27NZKFW3KKARERnOsBP6rqP2Mdj9+p6t9UtYWqtsa5UD1bVe2XbRFUdQvwq4gc687qD/wQw5D8bj3QQ0Rqu5/N/lgDEi8+BK52n18NTI9k4TbYYGycCvwJWC4iS9x5f1fVT2MYk6lcbgbeFJEawGpgWIzj8S1VXSgi7wNpOC1pvwNejG1U/iIibwN9gMYisgG4H3gEeFdErsVJ8JdEdJ/WQ4Qxxhi/sdN6xhhjfMeSkzHGGN+x5GSMMcZ3LDkZY4zxHUtOxhhjfMeSk6mQRKS1iKjbF1pZy7pCRKpUdzUi8rCI/MV9frqIrIx1TOEQkS4i8nWs4zDRY8nJRJWIrBWRHBFpHDJ/iZtcWpdzPIclNVV9U1XPjuI+h4rIYhHJFpHNIvJfETmtjGWuFZEzS7ltE+Aq4N8AqjpPVY8tfqsiy2omIh+KyKbC/p8ikuAOt7BLRLaIyF9Dlvd3h/XY4w4j08rLtqq6DNgpIueXJm7jf5acTHlYAwzJnxCRzkCt0hYWidpSeXG/UJ8CHgKOBFoCz+H05RYr1wCfqureCJR1EJgBXFTE8jFAO6AV0Be4S0QGArg/WD4A7gUaAouB/3jZ1vUm8OcIvAbjR6pqD3tE7QGsBUYD/wua9wRwD6BAa3feuTh35u8CfgXGBK3f2l03/070L4PmVXfXucjdVyecH113A78A24F3gYbueuvd7bLdR0+cL+v5QftTYCSQDvwOTODQDetxwJPANpyke1NwHCGvPcndxyXFHJ8EnOS1yX08BSS4yxoDHwM7cYYrmOe+ttdxksJet/y7gJrAG+7r3Qn8D6fX6ML2ORu4Mmi6D05ffMH/szuAZUAmTsKoWcL/uXrw/zNo/kbg7KDpB4B33OfXA18HLavjvqbjStrWnW7urp8Q6/e5PSL/sJqTKQ/fAPVEpIOIxAGX4XyRBtuNc6qpPk6iukFEQseH6Q10AAYEzxSRYcCjwJmq+j3OwHGD3PWP4lCCAWecI4D6qlpXVRcUEfN5wMk4A/VdGrTPEcA5QFeczlSLG8OmJ07SmFrMOvfgDDjZ1d1Xd5xkDnA7TieuTXBqXX8HVFX/hJNkz3dfw2M4fZslAUcDjXCSa1E1o85ASdeYLsUZXK4N0AUngYfF7e37KGBp0OylQEf3ecfgZaq6G+cHRUcP26LOUBe5QKlOSRp/s+RkysvrOMnnLOAnnF/FAao6V1WXq+pBda4nvI2TXIKNUdXdWvB01F+AO4E+qrrKnfdn4B5V3aCq+3FOD10c5unAR1R1p6quB+bgJA9wvrSfdsv+Had/saI0Arap6oFi1rkCGKuqW1U1A/gHTr+L4HzxNgNaqWquOteGiupvLNfdX1tVzVPVb7XoYVjqA1nFxATwjKpuUtUdwEccev3hqOv+zQyalwkkBi3PpKD85SVtmy8L5/WYSsaSkykvrwNDcX6Bvxa6UEROcS+IZ4hIJs4v/8Yhq/1aSLl3AhNUdUPQvFbAVBHZKSI7cXqYzsOpfXi1Jej5Hg59WR4VEkdhMeXbjtNRZnFJ8ShgXdD0OncewOPAKmCmiKwWkeLGZHod+Ax4x22c8Jg7LEthfufwL/lQRb3+cGS7f+sFzavHocSYHbIseHlJ2+ZLxDmNaSoZS06mXKjqOpxrNH/AuQge6i2cLviPVtUk4AVAQospZLuzgdEiEnxB/lfgHFWtH/So6Z4GKmtPx5uBFkHTRxez7gJgH8Wf+tuEk0zztXTnoapZqnq7qiYD5wN/FZH+7noFXodbs/qHqh6PMxbReTg11cIsA9oXE1NEuDXLzTinK/OdAKxwn68IXuaO1nsMsMLDtojIUUANSj5FaSogS06mPF0L9HOvLYRKBHao6j4R6Y5Ty/JiBc61kQkicoE77wXgwfxmySLSRETyW8dl4DQmSC7la3gXuFVEmotIfeD/ilpRVTOB+9zYBrnjBcWLyDki8pi72ts4ybWJ23rtPtzrcSJynoi0dccY2oVT+8tzt/st+DWISF8R6exe09uFc5ovf91Qn3L4KdNSE5GaOA07ABLc6Xyv4by+BiJyHM41u8nusqlAJxG5yN3mPmCZqv7kYVtwGnLMdk/dmkrGkpMpN6r6i6ouLmLxKGCsiGThfEm9G0a5S3FqChNF5BzgaZxa2Ey3vG+AU9x19wAPAl+5p/16hPkyJgIzcWof3+F80R+giESgzmCSf8Vp5JCBU6u7CZjmrjIOpwn1MmA5zphC49xl7YAvcE5xLQCeU9W57rKHcb64d4rIHUBTnBFvd+Gcxkzl8EYn+V4D/iAipW7OHyK/1SA41xODrwnej9PIYZ0b0+OqOgPAvcZ2Ec7/43ec/9HlXrZ1XYHzQ8RUQjaekzFl4CbDF1S1VYkr+4iIPARsVdWnYh1Labj3yr2oqj1jHYuJDktOxoTBrW30xak9HQlMAb5R1b/ENDBjKhlLTsaEQURq45xiOg7n9NUnwK3FNNs2xpSCJSdjjDG+Yw0ijDHG+I4lJ2OMMb5jyckYY4zvWHIyxhjjO5acjDHG+I4lJ2OMMb7z/8MLB4qQpP8uAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "toy_panel = pd.DataFrame({\n",
    "    \"mkt_costs\":[5,4,3.5,3, 10,9.5,9,8, 4,3,2,1, 8,7,6,4],\n",
    "    \"purchase\":[12,9,7.5,7, 9,7,6.5,5, 15,14.5,14,13, 11,9.5,8,5],\n",
    "    \"city\":[\"C0\",\"C0\",\"C0\",\"C0\", \"C2\",\"C2\",\"C2\",\"C2\", \"C1\",\"C1\",\"C1\",\"C1\", \"C3\",\"C3\",\"C3\",\"C3\"]\n",
    "})\n",
    "\n",
    "m = smf.ols(\"purchase ~ mkt_costs\", data=toy_panel).fit()\n",
    "\n",
    "plt.scatter(toy_panel.mkt_costs, toy_panel.purchase)\n",
    "plt.plot(toy_panel.mkt_costs, m.fittedvalues, c=\"C5\", label=\"Regression Line\")\n",
    "plt.xlabel(\"Marketing Costs (in 1000)\")\n",
    "plt.ylabel(\"In-app Purchase (in 1000)\")\n",
    "plt.title(\"Simple OLS Model\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Knowing a lot about causal inference, you decide to run a fixed effect model, adding the cities indicator as a dummy variable to your model. The fixed effect model controls for city specific characteristics that are constant in time, so if a city is less open to your product, it will capture that. When you run that model, you can finally see that more marketing costs leads to higher in-app purchase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAEXCAYAAAAJJYvtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeXxU9dXH8c/JHgIJ+75DQHZUFEFZAglga7Wbtlrr0qe2PrbWPhVcqtZdpKDVulSrtVpttVZrrbYEwi77vm9h3yEsCWuWmTnPH3eAITNJBsjkTpLzfr3mRXJ/d+49MyQ5872rqCrGGGNMNIlxuwBjjDGmNGtOxhhjoo41J2OMMVHHmpMxxpioY83JGGNM1LHmZIwxJupYczLVhog8ISKbqmA9d4iIpxKWEyci74jIIRFRERnqn/6ciOz3T7vjYtcTjURkqP/1tT6P51TK+25qBmtOJqqIyLv+P2qlH98HJgBXuV0jgIjMKKPO1QGzfQe4BfgG0AKYKyL9gYeBn/in/b2S6tkkIk+EMd8T/jqXhBjrE/A6wm4qxkRCnNsFGBPCV8BNpablq2ohcNyFesryN+D+UtNKAr5OB3ar6tzTE0QkHfCp6udVUF9Z8oBuInKZqi4NmP5TYDvQzp2yjDnLkpOJRsWquq/UozBws544/iMii0Qk3j8tRkSmiMgcEYnzT8vyf39KRHaLyJ9FpNHpFfmX87SIHBCR4yLyEdAgzDpPhajzkH+5M4CngY7+JLJNRN4F3gdiTieUgDq+LyLLRaTQP++LIpISuDIR+ZmIrBWRIn+9nwSsqxPweEDyaV9O3UeBT4C7ApZdByfl/an0zCJylYjM8r+HR0TkbyLStNQ894rILhE5KSKTgLYhlnO5iEz2v895IvJPEbFGaEKy5mSqJXWuu3UH0BIY65/8MHAZcLOqekRkGPA58BHQG/gm0B74TETE/5xfAL8CxvifuxR4vBJK/DbwArANZ/PdFcB9wC8Br39aC3D2tQB/8M/fHbgNyATeOL0wEXkSGAe8DvQCRgHLA9a1zf/808vdWUF9fwRuCWiA3wf24qTWM0SkOTAZ2AVcibOJsifwacA8NwC/A14E+gIfA+NLLac7MBOYB/QDhvnfhxwRSaqgVlMbqao97BE1D+BdwIOz+e70Y7N/7AlgU6n5M/zzP46zSe3bAWMzgOdLzd8WUKCv//tdwLOl5vkE8FRQ5wz/+o6XerwaME+oeu8ovWycxnJ3qWmD/XU2AFKAU8DocurZBDwRxvt7piZgDXCn/+v5OE16qH+9rf3Tn/a/RwkBy+jjn2ew//vZwF9LrWdCqeW8C3xUap5E4CTwzbLeG3vU3oftczLRaAFwe8D3ZR7BparTReQFnD+6b6jqPwOGrwCuEpGfh3hquohsAVoBc0uNzcZJWRX5DPh1qWkFYTzvDBFpgrOP50URmRA45P+3s//rJJwEU5neAu7yHxzRF7gOJxUF6gHMV9Xi0xNUdYWIFPjHZuGkvQ9LPW825+6PuwLoLCKl9xkm4eybM+Yc1pxMNDqlqmEdMi4iscDVOJuIOouIqOrpfTkxOJvC3g/x1H1ArP/rC700/9Fw6yzH6U3r9wHTQ4zvwtkkCRdeZ1new9kk+jvgM1U9eHZr5znKWq+GMc9pMTj/D8+HGDtUwXNNLWTNyVR3TwBdcBrUROBBzv4BXAz0KK+BiMhu/3P/GzD56ohUGoKq7heRnUBXVX0r1DwishYoBEYCq8pYVDFnm2246z7iP6jiVmB4GbOtAe4UkYTT6UlE+gBp/jGAtTjv2esBzyv9Hi7GabKbAz48GFMmOyDCVFsiMgTnIIjbVXUBztFnT4nI6XOhfgPcICK/E5G+ItJJREaJyJ9EJNk/zwvAfSLyQxFJF5H7cQ5GCEeyiDQv9Wha8dOCPAL8QkQeFZGeItJVRL4pIm8CqOpxf51P+I/Y6+I/J+nhgGVsBa4WkbYi0lhEwv3dvgtooqrTyhh/FUgF3vXXdg1OApqtqqcPnngB+J6I3Od/D+8EflhqOc8B3YAPRORKEekgIhki8rKIdAyzVlOLWHMy1ZKINAQ+AF5W1YkAqvop8GfgQxFJU9XpOEeF9cI5Cm0lziasY5w9H+ll4Pf+6cuBAcBTYZZxC84RboGPLef7WlT1fZzzur4OLAQW4STC3QGzPYa/iQGrcfY/XRYw/jhOmtmAcx5T0KHcZay7UFUPljO+HxgBtPbX9aV//d8JmOcznP1LD+C8xz/ASbCBy1kHDATqApNw0tZbQDKQH06tpnYRS9jGGGOijSUnY4wxUceakzHGmKhjzckYY0zUseZkjDEm6kT1eU4FBQV2tIYxxtRwaWlpQWd/W3IyxhgTdaw5GWOMiTrWnFyUm5vrdgnVir1f58fer/Nj79f5ifT7Zc3JGGNM1LHmZIwxJupYczLGGBN1rDkZY4yJOtacjDHGRJ2oPgnXGFVl7cerObzpEM1vauV2OcaYKmLNyUStY3uOMvWhyWyfuQ2AKzvFk56e7m5RxpgqYZv1TNRRVVZ/tIoPRrx7pjEBrHxhOcXHi90rzBhTZaqkOYnIOyJyQERWl5p+r4hsEJE1IvLbqqjFRLeju4/yr9s+ZeqDkyg+dm4jKi4oZv+KfS5VZoypSlW1We9d4FXgL6cniEgGcAPQW1WLRKRpFdViopCqsuajVXz1zIyQ6ahV/9Z0/tkltLk6rLuPG2OquSppTqo6S0Tal5r8v8Dzqlrkn+dAVdRios/RXQVMfWgyO77aHjQWlxzHNQ8Npvdtl7Jp8yYXqjPGuMHNAyK6AINE5FmgEBitqotcrMdUMVVl9Ycrmf3szNBp6arWZI0fRVrb+i5UZ4xxk6hWzS2T/MnpS1Xt6f9+NTANuA+4Avg70FEDCgq8n5NdlLFmObX/JCtfWM6hpXlBY7FJsXS9qzvtvtEBiQm6zYsxpgYIPPI21P2c3ExOu4B/+pvRQhHxAY2B4L9WUCMPIc7Nza2Rr6s8qsqqv65g9nMzKTlREjTe+qo2ZI4fGTIt1cb362LY+3V+7P06P5F+v9xsTv8ChgEzRKQLkAAcdLEeE2FHdxaQ88Akds3dETQWXyeea349hF4/6GNpyRhTNc1JRD4EhgKNRWQX8DjwDvCOf/NeMXC7VtU2RlOl1BeQlk6GSEsD2pD529BpyRhTO1XV0Xo3lzF0a1Ws37inYEc+Ux6YxK55O4PG4lP8aekWS0vGmHPZ5YtMRKhPWfnBcuaMnRUyLbW5ui2Z40aS2ibNheqMMdHOmpOpdAU78pkyZhK75odOS4MeGUrPW3ojYmnJGBOaNSdTadSnrPzLMmY/PwvPKU/QeJtr2pE5bgSprS0tGWPKZ83JVIr87flMeSCb3fN3BY0l1E1g0KND6fH9XpaWjDFhseZkLor6lBXvLWPOuNBpqe2gdgwfN5LUVqkuVGeMqa6sOZkLlr/tCFMemMTuBZaWjDGVy5qTOW9n0tLzs/AUBqeldkPaM2zsCEtLxpgLZs3JnJf8bUfIGTOJPQtDpKV6CQx6NIMe3+tpackYc1GsOZmwqE9Z/uelzP3tV2WmpeHPj6BeS0tLxpiLZ83JVOjI1iNMGZPNnkW7g8YS6iUw+DcZdL8xcmkpbtZ/iZ/+Bd3yDxPXuBnF1/8QX68rIrIuY0x0sOZkyuTz+vxpaTbeohBpaWgHJy21qBexGmIXzSTh728Qc/wosQCH9xPz9jhOPTABbdU+Yus1xrjLmpMJ6ciWw+SMzmbvkj1BYwmpiQz5TQbdvtsj4vuW4qd+Tszxo+dMi8k/SMIXH1B096MRXbcxxj3WnMw5fF4fy99ZytzxodNS+wwnLdVtHrm0FEiKC0NPP3msStZvjHFHuc1JROKA64GvA32A+kA+sAKYCPxLVYP/gplq6chmf1paWkZaejyDbt+JfFoK5GvWitjNa8+Zpgje9F5VVoMxpuqV2ZxE5KfAI8A6YCbwJXAMqAd0A+4CXhSR51T1jSqo1USIz+tj2dtLmPfCnNBpaVhHho/NqrK0FKj45nuI2bmFmF1bEFU0Ng5f5+6UjLyxymsxxlSd8pJTF+BKVd0XYuwz4DkRaQHcH5HKTJU4vOkQOaOz2bdsb9BYYmoiQ54cxiXf6u7aeUua2oBTv3md+Gmfc2LVEpIHZOAdkAmxtkXamJqszN9wVa2w6ajqXmB0pVZkqoTP62PZW4uZ9+IcvEXeoPEOwzsybOwI6jar60J1pSQkUjLqJrZ3upT09HS3qzHGVIEKP3769zsNA3rgbNI7BqwBptn+purpcO4hcsaUnZYy7u1Ol1v7I3WioDEZY2qlig6I6AN8DgiwEigAUoH7ABWRG1R1ZcSrNJXC5/Gx9O3FzC8jLXXuCtd2XkC9RV/i29gQ7xVDKP7e3S5Uaoyp7SpKTm8DL6jqK6UHROTnwDtAv0gUZirXoY0HyRmTzf7lwbsQE9OSGP7dFPps+5CYE8UAxObtJWbq53g7dMV7ZUZVl2uMqeViKhjvDpR1JN4fcY7aM1HM5/Gx6PUFfPj190M2po4jOvPDKXfSK34FMd7ic8ak6BTxX02qqlKNMeaMipLTOuB/gd+HGPupf9xEqUMbD5IzOpv9K4KbUlL9JIY+NZwu11/iHImnvjKWUtZ0Y4yJnIqa04+Bf4nIGM7d59Qb8ALfjGx55kL4PD6WvLmIBS/NxVscvG+p08jOZDyTRUrTlDPTvD0uJzZ3DeI7O7/Gx1PSb2hVlGyMMecotzmp6nIRSQcycDbx1QWOAy8DM1S1JPIlmvNxaONBJt8/kQMr9weNJTVIdtLSN7oGnbdU8o1bidm9jdh1y4g5mo+vfiO8vfvjHXxtVZVujDFnVHgoub8BTfY/TJTyeXwsfmMhC1+eFzotjUon45lMUpqkhHg2EBNL0T2PI3l7idm9HV+7zmiDxhGu2hhjQgvnPKefAHfgnOd0OjmtAf6sqm9FtDoTloMb8si5P5sDq0KnpYynh5N+XXBaCkWbtMDbpEUkyjTGmLBVdJ7TOOA64AWci72e3ufUF/iViHRU1YcrWomIvONfzgFV7VlqbDQwHmiiqgcv6FVUNz4fsSvm0XLeDGKPZ+DtOwAu4PJA3hIvS95YxIKX5+IrCT5wofO1Tlqq07iMtGSMMVGqouT0I6C3/zJFgZaKSDbOQRIVNifgXeBV4C+BE0WkDZAF7Air2pqguIikFx4kdvM6kkuK0KWz8HbuRuGvxkFCYtiLObjen5ZWB6el5IbJZDyTSfrXu1Zm5cYYU2UqOs+poo/zYX3cV9VZwOEQQ78DHgA0nOXUBAn/epfY9cuRkiIApKSI2HXLif/8LxU80+Et8bLw9/P48Lr3Qzam9K934dacO6wxGWOqtYqS05+AaSJSerNeH+BXwAXvcxKR64HdqrrCrSteuyFm8/qgji5A7Oa1VHToY966PHLun0jemgNBY8mNksl42tKSMaZmqOhQ8gdFZAtwJ8EHRPxeVd+8kJWKSB2ce0WNCPc5ubm5F7KqqNOxpIS0ENNPlHjYUsZr9Hl8bP4wl01/3YB6gkNmiyEt6X5vb6gfU2Pep7LU9NdX2ez9Oj/2fp2fi3m/KrrDgKhWzRY1EWkPfKmqPUWkFzAVOOkfbg3sodT9owoKCmrc5r6YFQtIevMZYk6cvc24LyWVwrsfwde7f9D8eWsPkHN/NnlrLS3l5ubaLTPOg71f58fer/NTme9XWlpa0OazC71lxmpg+oXeMkNVVwFNA9axDehXG47W8/XpT/F3fkz8jC/w5h8itn4jSjKuD2pM3mIvi15fwKJX5uPzBB+J1+X6Sxj65DCSG9apqtKNMabKVMktM0TkQ2Ao0FhEdgGPq+qfLrL2assz/AY8w65n0/p1dL6kW9Bh5HlrDjB59EQOrs0Lem5y4zoMeyaTztd2qapyjTGmylXJLTNU9eYKxttXtIwaRwSNiz+nMXmLvSx6bT6LXl1gackYU6tV1JwqumXGuMotp/Y6sHo/OaOzObguOC3VaVKHjGey6DzKtocbY2oHu2WGy7zFXha+Mo/Fry8MmZa63tCNIU8OI7lBsgvVGWOMO+yWGS4qyM3nw5+/z6H1wceB1GlSh2HPZtFppKUlY0ztE+4tM4Zy7nlOdsuMi3A6LS16dQHqCz5a/pJvdWfIExkk1be0ZIypncK9ZUaO/2Eu0v5V+8gZnV1GWkph+NgsOmZ1dqEyY4yJHhU2p7L4z3/6o6r+qBLrqbE8RR4W/n4+i/+wAPWGSEvf7s6Qxy0tGWMMXERzAmKB23GuXG7KsX/lPnLun8ihjYeCxlKapjBs7Ag6ZnZyoTJjjIlOFZ2EO62c4dhKrqXG8RR5WPDyPJa8sTBkWmqV1YbrXriBpLQkF6ozxpjoVVFy6g+MBUrfzwkgHrim0iuqIfat2EvO/dkczg2RlprVZdjYLLxtfdaYjDEmhIqa03Jgvap+UnpARBKB1yNSVSWJWbWQhC8+IOZYAb6UVEpGfhfvFUMiuk5PoYcFL89lyRuLQh6J1+3GHgx+LIOktCS7ArIxxpShoub0EqFvEghQgnMrjagUsy2XpLfHEZPvJJcYIGb/TgpT6uHrfllE1rlv+V5yRpedloY/P4IOwzpGZN3GGFOTVHSe0z/KGfMB71V6RZUk/ssPzjSm02KO5pMw8e8UVnJz8hR6mP/SXJa+GTotdb+xJ4MfG0qibcIzxpiwXMzRelFNTp4IPb3wZMjpF2rfsr1Mvn8iRzYHB8y6zZ201D7D0pIxxpyPGtucfG07w5rFQdO9zdtWyvI9hR7mvziHpW8tDp2WburJ4McySExNrJT1GWNMbVJjm1PxDbcRu2E5Mds2Ij4fiuBr24ni7/3kope9d+keckZnh05LLeo5aWloh4tejzHG1FY1tjmRXIdTv/49cTO+JDZ3Nb72XSgZ/k1IvPD9Pp7CEua9MIdlby8JmZZ6fL8Xgx4ZamnJGGMuUtjNSUTSgK44F389Q1XLO1HXXfEJeLK+jSfr2xe9qL1L9pAzpoy01LIemc+PoN0QS0vGGFMZwmpOInIH8BrOFckDjyhQoEbv7fcUljBvwhyWvr3YebWl9Ly5N9f8eoilJWOMqUThJqdnge+q6sRIFhNt9izaTc4D2eRvORI0VrdlPTLHjaTd4PZVX5gxxtRw4TanOGByJAuJJiWnSpg3fjbL3lkSOi3d4k9L9SwtGWNMJITbnMYBj4rI0/6Tb2us3Yt2MWXMJPK3Bqeleq1TyRw3krbXtHOhMmOMqT3CbU7/BzQHHhCRcy67oKqVc+KQy0pOlTB3/GyWl5GWet3ah2seHkJC3YSqL84YY2qZcJvTrRGtwmW7F+4iZ0w2Bdvyg8YsLRljTNULqzmp6sxIF+KGkpPFTlr689KQaan3D/ty9UODLS0ZY0wVK7M5icgjqvqs/+unyppPVX8TicIibfeCneSMmUTB9uC0lNo6lczxo2gzsEZssTTGmGqnvOTUOuDrNpEupKqUnCxmzrivWPHuspDjvW/zp6UUS0vGGOOWMpuTqv5vwNcXdd8mEXkHuA44oKo9/dPGA98AioHNwJ2qGhxjKtGu+TuZMiabgh0FQWOpbdLIGj+S1gMsLZnw+FT5745CVh4qZmSbZC5vYh9ojKks5W3Wa6qqBypagIg0U9X9Fcz2LvAq8JeAaTnAw6rqEZFxwMPAgxWXfP6KTxQzd9xXrHgvdFrqc/ulDHxwkKUlE7b8Ih/fzTnI6kMlFPrgzbUnGNoykT9nNCRGxO3yjKn2ytusN11EZgLvAwsCz28SkRjgSuA2YDDQs7yVqOosEWlfalrgSb3zge+eV+Vh2rd8LxN//iVHdwanpbS2aWSOH0Xrq2rMVktTRR5ekM/ivJIz3xeUKBN3FvL3TSe5OT3FxcqMqRliyhm7FFgL/BE4JiKrRGSuiKwCjgFvAKuAyrit7I+AiFwaKbFeIicOHA+a3vfOy/jBpNutMZkLsj7fEzSt2Adf7ih0oRpjah5RDXEMdemZRNoAvYD6wBFgparuPq8VOcnpy9P7nAKmPwL0A76tpYopKCg4831ubu75rO4cWz7exPo/rgGgTssUeo/uS8PejS94ecb8z4pEVh6LDZp+bZMSnupaEuIZxphA6enpZ75OS0sL2hYe7nlOO4GdlVeWQ0RuxzlQYnjpxlRa4As5X50e6kTBkiM069OcgWOuIb5OdOxbys3NvajXVdtE0/v1veJjbFpylJPes9OaJsXw6MCWpDe0n6/qyN6v8xPp98u1mw2KyCicAyCGqOrJiua/GDGxMXzno+8RmxD8SdeYC3FP97ocKfTxxfZCjhb7aJocy70969IjShqTMdVdlTQnEfkQGAo0FpFdwOM4R+clAjniHN00X1XvjlQN1phMZRIRHr08jYcvTeW4R0mNF8SO0jOm0lRJc1LVm0NM/lNVrNuYSIqNEdISrCkZU9nKO1oviIjEiEiLSBVjjDHGQJjNSUTqi8jfgEJgk3/a9SLyTCSLM8YYUzuFm5zeAAqAdjiXGwKYB3wvEkUZY4yp3cLd5zQcaKmqJSKiAKqaJyJNI1eaMcaY2irc5FQAnHPWqoi0BfZWekXGGGNqvXCb09vApyKSAcSIyADgPZzNfcYYY0ylCnez3jicgyFeA+KBd4A3gZcjVJcxxphaLNzLFynwkv9hjDHGRFS4h5JniEgH/9fNReQ9EXlHRJpHtjxjjDG1Ubib9V4HRvq/ftH/rwfndhrXV3ZRxhhTVeatzWHl1nkcP3GUpttb8rX+t9Konh2I7LZwm1MrVd0hInE4Ter0+U57IlaZMcZE2OzVE5m+/DOKPUUA5J/M48ixA/z0usdJjE92ubraLdyj9Y6KSDNgCLBWVU/fvS8+MmUZY0zkrdw6/0xjOi2vYC8L1k11qSJzWrjJ6RVgEZAA/NI/7WpgfSSKMsaYqlBcfIqYw/FISQzeZmebVF6BbRRyW7hH640Tkc8Ar6pu9k/eDfw4YpUZY0wEqU9JXNiIBp82xlu/mPy7N0GcEhcTT7d2l7tdXq0X9i0zVHVjed8bY0x1kb89nyljsileEIsAcXlJ1JnRlKIRh+jQohuXtLnU7RJrvbCak4ikAk/g7HNqDJy5gY2qto1IZcYYU8nUp6x4bxlzxs3Cc8pzzlidOU0Y/KMsBgwdToyc192ETAScz6HkrYGngA+AW4ExwKcRqssYYypV/rYj5IyZxJ6Fu4LGEuol0PUn3RmYMdzuaBwlwm1OI4BuqnpIRLyq+rmILAa+AH4XufKMMebiqE9Z/u5S5o77Ck+hJ2i83ZD2DH9+BPtO7LfGFEXCbU4xOFcmBzguIvVxrkjeOSJVGWNMJTiy9QhTxmSzZ9HuoLGEegkMfiyD7jf1RETYl7vfhQpNWcJtTitw9jdNBb7CuQDsccAOijC11slDJzm6q4DmfVq4XYopxef1seLdZcz9bTlpadxI6rWo50J1JhzhNqe7OHsQxC+AsUB94LZIFGVMtMv9zwamPzqFmPhYfphzB4lpSW6XZPyObDlMzphJ7F0cIi2lJjL4saF0v7GnbcKLcuGe57Ql4Os87PwmU0udPHiC6Y9NZdN/z240mPX0DLImjHKxKgNOWlr+zlLmjp+Ntyg4LbXP6MCwsSMsLVUTYZ/nJCIjgL5A3cDpqvqbyi7KmGijquT+ZwMzHpvKqcOnzhlb+4/VdPtOd1oPsLMq3HJk82FyxmSzd0nwlR0SUhMZ8ngG3b7Tw9JSNRLueU6vAjcB04GTEa3ImChzIu8EMx6bwqaJuUFjMfEx9L9vIC36tXKhMuPz+lj+pyXMnTAndFoa1pHhY7Oo29zSUnUTbnK6GeirqjsjWYwx0URV2fjFBmb8ZiqFR04FjTft1YysF0bRuGsTF6ozhzcdImd0NvuW7Q0aS0xNZMgTw7jk290tLVVT4TanQ0B+JAsxJpqcyDvB9EensDk7OC3FJsTS/5cDufynVxATZ1cSqGo+r49lby9h3guz8RZ5g8Y7DO/IsLEjqNusbohnm+qizOYkIh0Dvn0B+KuIjAXOORkg8GCJcpb1DnAdcEBVe/qnNQT+DrQHtgE3qeqR86zfmEqlqmz893onLeUXBo037d2MrAmWltxSblpKS2LoE8Po+q1ulpZqgPKS0yZACbiOHk6DCaRAbBjreRd4FfhLwLSHgKmq+ryIPOT//sEwlmVMRJw4cILpj+awedKmoLHYhFj6/99ALv+JpSU3+Lw+lr61mPkvzgmZljpmdWLYs1mkWFqqMcpsTqpaab+BqjpLRNqXmnwDMNT/9XvADKw5GReoKhs+X8/Mx0OnpWZ9mpM1YRSNujR2oTpzaONBcsZks3/5vqCxxLQkhj45jK7ftLRU04R7tF4r4GTgZjcRaQAkq+qF3pWrmaruBVDVvSLS9AKXY8wFO7H/ONMencKWyaHT0lW/Gshld1lacoPP42PJHxex4Hdz8RaHSEsjOjtpqWmKC9WZSBNVrXgmkUXAj1R1VcC0XsDbqto/rBU5yenLgH1O+apaP2D8iKo2CHxOQUHBmeJyc4N3TBtzoVSVPVN3sfa1VZQcKwkaT7ukPr3HXEq9dqkuVGeObTvKyvHLKNgQfBxWfL14etzbmxYZrSwtVWPp6elnvk5LSwv6jwz3aL0ugY0JQFVXicglF1HbfhFp4U9NLYAD5c0c+EJqitzc3Br5uiKlst6vE/uPM+2RHLbkbA4ai02MZcCvrubSH/er9mmpOv58+Tw+lry5iAUvhU5LnUZ2JuOZyKSl6vh+uSnS71e4zSlPRDqr6pltHyLSGecQ8wv1b+B24Hn/v59fxLKMqZCqsv6ztcx8fBpFR4uCxptf2oKs8aNomN7IherMwQ155IzO5sDK4KuDJzVIZuhTw+nyja6WlmqJcJvTO8CnIvIIsAXoBDwNvB3Ok0XkQ5yDHxqLyC7gcZym9LGI/A+wA7jx/Eo3JnzH9x9n2sOT2To1+MyH2MRYBtx/DZf++HJiYqt3WqqOfB4fi99YyMKX54VOS6PSyXgmk5Qmtm+pNgm3OT0PlAATgDY4zeRPwIvhPFlVby5jaHiY6zfmgqgq6/+5lplPhE5LLS5rSeb4kTTsbGnJDQfX+9PSqtBpKePp4SDzSVAAACAASURBVKRfZ2mpNqqwOYlILE5y+omqjo98ScZUjuP7jjH14Ry2TQuVluIYOPpq+v6PpaWqcuRYHjNX/psThcdo07AzCbMasOiVBfhKfEHzdr7WSUt1Gltaqq0qbE6q6vVfkTz4J8iYKKSqrPtkDTOfmk5xqLR0eUuyxo+iQaeGLlRXO+04sImPZ75OwYlDxO5LYt+/jhG3NzlovuSGyWQ8k0n617u6UKWJJuFu1vsd8KSIPKGqxZEsyJiLcXzfMaY+NJlt07cGjcUmxjHwgWvoe+dllpaq2NSln1Bw9BDJXzWlzqwmiDf4/U//eheGPjXc0pIBwm9O9wLNgV+JSB7OZYsAUFW7iY1xnaqy9h+rmfX0jNBpqV8rssaPpEFHS0tuOLrpOPU/6Bw6LTVKJuNpS0vmXOE2p1sjWoUxF+HY3mNMfXAS22duCxqLS4pj4AOD6HPHpZaWXOAt8bL49QXwUhpxvuCDGtKurstNr9xGnUZ1XKjORLNwb9M+M9KFGHO+VJW1H69m1tPTKT4WvLW55RWtyBw/igYdGoR4tom0vDUHyBmdTd7aA5x7/WjwpXiIufEU33/kZyQlBKcpY8K9tt5TZY3ZbdqNG47tOcrUhyaXnZYeHETfOy5DYuwQ5KrmLfay6LX5LHp1AT5P8HFUMZeX0PiOFL4+7C5rTKZM4W7Wa1Pq++bAEOCzyi3HmPKpKms+WsVXz8yg+HiItHRla7LGj6R+e0tLbshbc4DJoydycG1e0Fhy4zoMeyaTztd2caEyU92Eu1nvztLTRGQUzu3bjakSp/af5F9PfcqOWduCxuKS47j6wcH0uf1SS0su8BZ7WfjqfBa/Fjotdbn+EoY+OYzkhrZvyYQn3OQUymScO9kaE1ETt5/kX28uI/2TRcQXeYLGW/VvTeb4UdRvVz/Es02kHVi9n5zR2RxcF5yW6jSpQ8YzWXQeZRdUrQlUlZKd/8R7YDaNTx2lsKgTiV3uQRIq/3cv3H1OHUtNqgPcAuys9IqMCZC9+ABTHppM99zgG83FJcdxzUOD6X2bpSU3eIu9LHxlHoteW4B6g2+90/WGbgx5chjJDWy/Uk1RsvV9SnZ8Cr4iEgDvgd0UntxDUr+XkJiLyTrBwl1a6TuxnQSW4VxN3JhKp6qs/nAla56cQdvC4PstnbykOfe8dR1pbS0tuWH/qn3kjM7m0PqDQWN1mtRh2LNZdBppaakmUVU8eXPAd+55hL4TW/Hsn0l8i8q9VGq4+5zsBBFTZY7uKmDKg5PZOXt70A9oUXwcU7L6UP+bPa0xucBT5GHh7+ez+A+h09Il3+rOkCcySKpvaanG0RLwnAwx3Yvv+BYq+zre5TYnEWmGc+minsBS4H5VvZh7OBlTJlVl1V9XMPu5mZScCE5LW9s35fPr+3OkYV1GN05wocLabf9Kf1raECotpTB8bBYdszq7UJmpChKTgCQ0RItK7VuMrUNskwGVvr6KktPrQH3gDeA7wEvADyu9ClPrHd1ZQM4Dk9g1d0fQWElCHJOy+rL48s4QI1zaOJ5f9q7nQpW1k6fIw8KX57H4jYWh09K3uzPkcUtLtUF8+5sp3vAKWuzPKBJPbMNLiU3rUenrqqg5DcK5RXu+iHyMs5/JmEqjvoC0dDI4LbUe0IYrnhnB0SMx1N2bz6D2Dbi7e13qVPNbqFcX+1fuI+f+iRzaGLzBJKVpCsPGjqBjZicXKjNuiGtyFZLSFs/2jzmWv4+0DiOIazY0Ivfbqqg5JalqPoCqHhQRu1ywqTQFO/KZ8sAkds0LPugzPiWea349hF639EFihKeA3NwDpKenVn2htZCnyMOCl+expIy01O073Rn8+DCS0pJcqM64KbZOS2K7/ZIdubk0bh65g14qak7xInInZy+MlSgiPwqcQVXfiUhlpsZSn7Lyg+XMGTsrZFpqc3VbMseNJLVNmgvVmX3L95IzOpvDuSHSUrO6DBubRcfhlpZMZFXUnBYAtwV8v5Bz9zkpzl1yjQlLwY58poyZxK75odPSoEeG0vOW3nZbbhd4Cj3Mf2kuS99chPpCpKUbezD4sQxLS6ZKlNucVHVoFdVhajj1KSv/sozZz8/Ccyr4Kg9trmlH5rgRpLa2tOSGfcv2kjN6Ioc3HQ4aS2lWl+HPj6DDsNLn4hsTOZV7Sq8xIeRvz2fKA9nsnr8raCyhbgKDHh1Kj+/3srTkAk+hh/m/m8PSPy4OmZa639iTwY8NJdHSkqli1pxMxKhPWfHeMuaMC52W2g5qx/BxI0ltZQc5uGHv0j3kjM7myObgtFS3uZOW2mdYWjLusOZkIiJ/ez5TxmSze4GlpWjjKSxh/otzWfpWGWnppp4MfiyDxNREF6ozxmHNyVSqM2np+Vl4CoPTUrsh7Rk2doSlJZfsXbKHnDFlpKUW9Zy0NLSDC5UZc66wm5OIpAM3AS2BPcA/VHVjpAoz1U/+tiPkjJnEnoUh0lK9BAY9mkGP7/W0tOQCT2EJ8ybMYenbi51jbEvp8f1eDHpkqKUlEzXCvWXGLcAfgf8A24FewEMi8lNV/VsE6zPVgPqU5X9eytzfflVmWhr+/AjqtbS05IY9i3eTMyab/C1HgsbqtqxH5vMjaDfE0pIJj8/rC3lDycoWbnJ6Bviaqs46PUFEBgHvAxfVnETk/4Af43yeWwXcqaqFF7NMU3WObD3ClDHZ7Fm0O2gsoV4Cg3+TQfcbLS1VJVVl6751rN26lK1/2ci699eFTEs9b+7NNb8eYmnJhO3I5sPkjM6mxeUtaf69VhFdV7jNqR4wr9S0+cBFXc5IRFoBvwC6q+op//X7vg+8ezHLNZHn8/r8aWk23hB3p203tIOTllrYBVqrUmHxSd7PeZEDyw+Q/GkzYg8F36Sxbst6ZI4bSbvB7au+QFMt+bw+lr29hHkvzMFb5GHvsj0M6JkMEbxlV7jN6UXgORF5TFULRSQZeNI/vTJqSBaREpw77O6phGWaCDqyxfn0tHdJ8H9VQmoiQ36TQbfv9rC05IJ/z/wLh94/Tsr8NogGv/89b/GnpXqWlkx4Dm86RM7obPYt23t2osLKCcu4bNRlxCXFR2S94Tane4DmwH0icgRogHO9vb0i8r+nZ1LVtuezclXdLSITgB3AKWCyqk4+n2WYquPz+lj+zlLmjg+dltoP68jwsVnUbW5pyQ27F+1i78MFJOc1DhqLaSTc8Pvv0vaadi5UZqojn9fHsrcWM+/FOXiLvEHjKa1SKDnliVhzEtUQG6NLzyQyJJyFqerM81q5SAPgU+B7QD7wD+ATVf0AoKCg4Exxubm557NoU8mO7zzGyvHLyF8bvFM9rm483e/pSausNpaWXOAt9LDhnXVs+2xLyH1Lp/odotkPGnFNrxuqvjhTLR3f7v99Xx/8+x5fL57uP+tFy+GtL+r3PT397DbBtLS0oAWFe5v282o65yET2KqqeQAi8k9gIPBB6RkDX0hNkZubG/Wvq/S25tI6DO/IsLEjqNusbsRrqQ7vV1XbvXAXOaOzKdieHzTmrV/M8et3kdgjjq8N+znN6kd2B3Z1Zz9fzu/70rcWM7+MtNQxqxPDns0ipVndiL9f4R5KngA8CtzM2fOcPgKevcgj63YAV4lIHZzNesOBxRexPFOJQm5r9ktMTWTIk8O45FvdLS25oORkMXN/O5vl7y4NmZZ04ElOZu6jYeOGDOw+0hqTqdDh3EPkjCnj9z0tiaFPDqPrN7tV2e97uPuc/gB0xTmybjvQDngYaAX8qJznlUtVF4jIJ8BSwINzp90/XujyTOWoaFtzh8xODH/O+fRkqt7uBTvJGTMpZFpKbZ1K5m9H0ax/M9auX03v7pcSE2N3DTZl83l8LH1rEfN/Nzd0WhrR2UlLTav2XrPhNqdvAp1O3xUXWCsiC4BNXERzAlDVx4HHL2YZpvKUm5Zc+PRkzio5WcyccV+x4t1lIcd7/7AvVz88mISUBABSElOtMZlyHdp4kJzR2exfEXzKQVL9JIY+NZwu11/iyu97uM1pH85h3oEf1ZKB4L9gplrL33okZGMK3NYcjiNFPl5eeYzNxzz0bhjPPT3qkhJvfygv1K55O5jywCQKdhQEjaW2SSPztyNpM/C8DpY1tZjP42PJHxex4Hdz8RYHp6VOIzuT8UzVp6VA4Tan94FsEXkF2AW0AX4G/EVEhp2eSVWnVX6Jpip1zOrMJd/qzvrP1gLOp6chTw6n6w3hf3rafdzDtycfYkOBcwDFF9sL+XJHIV+MakxqgjWo81F8opg5z89i5V+Whxzvc/ulDHxw0Jm0ZExFDm08yOT7J3Jg5f6gsaQGyU5a+kZX17eOhNucfur/99elpt/tf4CzW9Zu/lIDDHkigx2zt9PishYX9OnpN4sLzjSm01YcKuGFFcd48gq70224ds510tLRncFpKa2tk5ZaD7C0ZMLj8/hY/MZCFr48L3RaGpVOxjOZpDRxLy0FCvdQcrsqZC2SVD+Zm//zQ1KaplzQp6edx4N/8AHWHCm52NJqheITxcwZO4uV75eRlu64lKsfHER8HUtLJjwHN+SRMzq7zLSU8fRw0q9zPy0Fsvs5mZAu5rylsjbdNUqyTXoV2TlnB1MeyOborqNBY2lt08gcP4rWV7VxoTJTHXlLvCx5YxELXp6LryT4SuKdr3XSUp3G0ZGWAoV7nlMq8AQwBGiMc+ki4PwvWWRqvvt61WXVoRL2F579ZWiTEsODfe2yRmUpPl7M7LEzWfXBipDjfe+8jIEPXGNpyYTt4Po8cu7P5sDq4LSU3DCZjGcySf96VxcqC0+4yel1oDXwFM7VG24FxuBcesiYcwxqkcRrgxrw8upj5Bf5aJIUy6OXpdIxNTLX4KrudszezpQHJ3EsVFpqV5+s8SNp1d/SkgmPt8TLkj8sZMHv54VOS1/rQsbTw6MyLQUKtzmNALqp6iER8arq5yKyGPgC+F3kyjPVVWbrJDJbJ7ldRlQrPl7MV8/NZPVfQ6Qlgb4/upyBY64hPtmauglP3ro8cu6fSN6aA0FjyY2SyXg6utNSoHCbUwxw+pCh4yJSH+ccp84RqcqYGm7H7O1MeSCbY7uPBY2lta9P1oRRtLqitQuVmerIW+Jl8esLWfhK6LSUfl1Xhj41nDqN6rhQ3YUJtzmtwNnfNBX4CngNOA5sjFBdxgWz9xbywsrjHCny0TgphkcvS6VvY9vHUZmKjhUx+7mZrP7byuBBgUt/dDkDLC2Z85C39gA5o7PLTUvt+27Gs+0xTm4pJqZOaxLT70YSovu0jnCb012cPQjiF8BYoD5wWySKMlVv0YFifjzrCPtOnv3UtSH/EP8a2ZhOafaHsjJs/2obUx+cFDIt1e/QgKzxo2h5hV2g1YTHSUsLWPj7+fg8wWmpy/WXMPTJYcQe/YLiTR+B95TzvGO5FJ7YSVK/l5CY6D1gO9zznLYEfJ0H/DhiFRlXvLDy2DmNCWDnCR/PLz/GW0MaulRVzVB0rIjZz85k9YdlpKX/uZwBoy0tmfDlrfGnpbUh0lLjOgx7JpPO13YB4GTujDON6TTf8a149k0jvuWIqij3gpx32xSRVaraKxLFGPfkh7gaMcDBwuBPZCZ822duZcpDkzm+J0Ra6uhPS/0sLZnweIu9LHptPoteXVBuWkpu6OxbUvWC50SoJeE7vjnC1V6cC8l0dp/nGqht3TjmHwi+gkPX+tEb+6NZ0dEivnpmBmv+vip4UOCyu/ox4P6rI3aLa1Pz5K05wOTREzm4Ni9orE6TOmQ8k0XnUefe/E8kFklsgBaVSlgxScQ26h/Jci/ahfzliZ7rW5hK80S/VJYfKmFjwDXxejaI46G+qS5WVT1tm7GVqQ9N5vje4LTUoFNDssaPosXlLV2ozFRH3mIvC1+dz+LXQqelrjd0Y8iTw0hukBzy+fHtvk/xhlfQ4sP+KbHENuhNbMNLI1j1xbuQ5nRtpVdhXNcyJY7srzXmhRXH2HrMQ7cG8dzXq55dRfw8FBUUMuuZGaz9eHXQmMQIl93Vj6t+NdDSkgnbgdX7yRmdzcF1odPSsGez6DSy/FulxzUZgNRphWf7x6jnBLGNriCu5ciouo5eKOfVnESkKbBHRM5cfTzwYAlTvTVMiuXZ/vXdLqNa2jZ9i5OW9h0PGmvQqSFZE0bR4jJLSyY83mIvC1+Zx6LXFqBeDRq/5FvdGfJEBkn1Q6el0mJT2hLbfXRllxlR4V5bbxTwJ6A5527WUyA2AnUZUy0UFRQy6+kZrP1HGWnpJ/246v+uJi7J9t2Z8OxftY+c0dkcWn8waKxOkxSGj82iY1bNv/5BuL8xrwFPA++p6qmKZjamNtg6zUlLJ/YHp6WGnRuSNeFaml/awoXKTHXkKfKw8PfzWfyH0Gmp69ebM/S574Sdlqq7cJtTA+BNVQ1+x4ypZQoLCpn11HTWfbImaExihMvvvoL+9w20tGTCtn/lPnLun8ihjYeCxurUL+Lq29fR/qo4ElOuB6w5BfoTcCfwTgRrMSbqbZ26makP54ROS+mNyJowiuZ9LS2Z8HiKPCx4eR5L3lgYMi2lD9zDVbdsJDHFgxZCya4vSOjwAxcqrXrhNqergF+IyEPAvsABVR1c6VUZE2UKCwqZ9eQ01n26NmhMYoTL//dK+v9igKUlE7Z9K/aSc382h3NDpKUGJVxz2xra9j13v5MWFwTNW1OF+5v0tv9hTK2zZcpmpj08mRMHgs+0b9SlEZkTRtG8j6UlEx5PoYcFL89lyRuLUF9wWup2Yw+uvGEO8SWlDoiISyWu1deqqEr3hXttvfciXYgx0aYw/xQzn5zO+n+GSEuxQr+7r+TK+wYQl2hpyYRn3/K95IwOnZZSmtVl+PMj6DCsI95Tl1O06gn0xA5QDyQ0IK55FrF121d90S4p97dKRIZVtABVnVZ55RgTHbbkbGLqwzmczAuRlro2JmvCKJr1bu5CZaY68hR6mP/SXJa+GTotdb+xJ4MfG0pimnODztjkpiT3ewXvwTn4CvOIazqImKSmVV22qyr6yPenCsYV6FjBPMZUG4X5p5jx+DQ2/Gtd0JjEClfc058r7r3K0pIJ275le5l8/0SObD4cNFa3uZOW2mcE/xmVmFjimtbeXfrl/oapaodIF+C/q+7bQE+cZvcjVZ0X6fUaU9rmyZuY9uvJnMw7GTTW6BJ/WuplacmEx1PoYf6Lc1j61uLQaemmngx+LIPE1EQXqot+0fDx72UgW1W/KyIJQPW5j7CpEU4dOcXMx6ex4fMy0tLP+nPlvQOITbCLoZjw7F26h5zR2aHTUot6TloaGvHP/tWaq81JRFKBwcAdAKpaDBS7WZOpXTZPymXaIzkh01Ljbk3IGj+Kpr2auVCZqY48hSXMe2EOy95eEjIt9fh+LwY9MtTSUhjcTk4dgTzgzyLSB1gC3Keqoe6OZUylOXX4JDMen8bGf68PGouJi6Hfz/pz5c+vsrRkwrZ3yR5yxpSRllrWI/P5EbQbYmkpXOLmFYlEpB8wH7haVReIyMvAUVV9DKCgoOBMcbm5uS5VaWqafV/tYfXLKynOLwoaq9cpld6jLyUt3a7ObsLjLfKy8d11bP1ks7PXvJQ2X2vHJT/pQXxdu1VKoPT0s7f6SEtLC7p/h9vJaRewS1UX+L//BHgo1IyBL6SmyM3NrZGvK1Iu9v06dfgkM34zlY1fbAgai4mL4Yp7r+KKe/rXmLRkP1/n50Lerz2LdpPzQDb5W44EjdVtWY/McSNpN7h9JVUYXSL98+Vqc1LVfSKyU0S6quoGYDgQfMajMRcp978bmf5oDqcOBV9Uv3H3JoyYcC1NetSu80jMhSs5VcK88bNZ9s6SkGmp5y29uebXQ0isZ/uWLpTbyQngXuCv/iP1tuBcYNaYSnHykJOWcr8MnZau/MVV9LunP7HxNSMtmcjbvWgXU8ZMIn9rcFqq16oemb8dRdtr2rlQWc3ienNS1eVAP7frMDVP7n82MP2xKSHTUpMeTcmaMIom3S0tmfCUnCph7vjZLC8jLfW6tQ/XPDyEhLoJVV9cDeR6czKmsp08eMJJS//ZGDQWEx/DlfcOoN89V1paMmHbvXAXOWOyKdiWHzRWr3UqmeNGWlqqZNacTI2y8cv1zHhsKqcOB6elpj2bkTlhFE26NXGhMlMdlZwsdtLSn5daWqpi1pxMjXDy4AmmPzaVTf8NnZb63zeQy+++wtKSCdvuBTvJGTOJgu3BaSm1dSqZvx1Fm6vbulBZ7WDNyVRrqkrulxuY/thUCo+ETktZL4yi8SWWlkx4Sk4WM2fcV6x4d1nI8d4/7MvVDw8mIcXSUiRZczLV1om8E0x/dAqbs4NP0I6Jj6H/Lwdy+U8tLZnw7Zq/kyljsinYEXzH2dTWqWSOH0WbgZaWqoI1J1PtqCobv9jAjN+UkZZ6NyNrwigad7W0ZMJTfKKYNa+sZPvnW0OO976tL1c/ZGmpKllzMtXKiQP+tDQpOC3FJsSeSUsxcTEuVGeqq6+enhGyMaW2SSNr/EhaD7C0VNWsOZlqQVXZPXUX0/4wicL8wqDxZn2akzVhFI26NHahOlPd9f/lANb/ey2eE54z0/rcfikDHxxkackl1pxMtbD8naWsGLskaHpsQixX/Wogl91laclcuLrN69H9nl6sHL+MtLZpZI4fReur2rhdVq1mzclEvW37NrC8wTSoq3D8bANq1rc5WeMtLZmL49k/i5Ld/6FP73zSftaJnj/5Jgn1bX+l26w5mai2bd8GPp75OsdO5ZNwXSqpH7VD45Qr7ruCAfcMtrRkLkrJnokUb/oTeI6TCHTttx3Phq3E93sJibWLtrrJfrNNVJu56guOnXJOgizudpQTGfvJvzuXg5dvscZkLppn90TwHD9nmp7YRsnu/7pUkTnNkpOJaoVF594+/dTQAwAcPRl81r4x50tLNSb/VHzHNlV5LeZc9tHTRLXUOqHvSNu0fssqrsTURJLQIMTEeGIbXV71xZhzWHMyUS3r8ptoWO/c21o0b9CWIX2ud6kiU5PEt/8BktgoYIoQk3YJcU2HuFaTcdhmPRPVGqc158fXPsL05f9iT95OOrTqwpDe3yApoY7bpZkaIK7RZUjfsZRs+4gTBfup1/JK4tt+C4mxS165zZqTiXr16tTn+oF3kJubS3p6utvlmBomNqUtsT0eYGduLo3a289XtLDNesYYY6KONSdjjDFRx5qTMcaYqGPNyRhjTNQRVXW7hjIVFBREb3HGGGMqRVpampSeZsnJGGNM1LHmZIwxJupE9WY9Y4wxtZMlJ2OMMVHHmpMLRKSNiEwXkXUiskZE7nO7pupARGJFZJmIfOl2LdFOROqLyCcist7/czbA7ZqimYj8n/93cbWIfCgiSW7XFE1E5B0ROSAiqwOmNRSRHBHJ9f8b4iq6F86akzs8wP2q2g24CviZiHR3uabq4D5gndtFVBMvA9mqegnQB3vfyiQirYBfAP1UtScQC3zf3aqizrvAqFLTHgKmqmo6MNX/faWx5uQCVd2rqkv9Xx/D+cPRyt2qopuItAa+Drztdi3RTkRSgcHAnwBUtVhV7QZY5YsDkkUkDqgD7HG5nqiiqrOAw6Um3wC85//6PeCblblOa04uE5H2wKXAAncriXovAQ8APrcLqQY6AnnAn/2bQd8WkRS3i4pWqrobmADsAPYCBao62d2qqoVmqroXnA/cQNMK5j8v1pxcJCJ1gU+BX6rqUbfriVYich1wQFWXuF1LNREHXAb8QVUvBU5QyZtcahL/vpIbgA5ASyBFRG51typjzcklIhKP05j+qqr/dLueKHc1cL2IbAM+AoaJyAfulhTVdgG7VPV0Gv8Ep1mZ0DKBraqap6olwD+BgS7XVB3sF5EWAP5/D1Tmwq05uUBEBGd/wDpVfdHteqKdqj6sqq1VtT3OjuppqmqfbMugqvuAnSLS1T9pOLDWxZKi3Q7gKhGp4//dHI4dQBKOfwO3+7++Hfi8MhduNxt0x9XAD4FVIrLcP+3XqvpfF2syNcu9wF9FJAHYAtzpcj1RS1UXiMgnwFKcI2mXAX90t6roIiIfAkOBxiKyC3gceB74WET+B6fB31ip67QrRBhjjIk2tlnPGGNM1LHmZIwxJupYczLGGBN1rDkZY4yJOtacjDHGRB1rTqZaEpH2IqL+a6Fd7LJ+ICK16nI1IjJWRH7p/3qQiGxwu6bzISK9RWSu23WYyLHmZCJKRLaJSLGINC41fbm/ubSv4nqCmpqq/lVVR0RwnbeIyGIROS4ie0Vkoohcc5HL3CYimRf43CbAbcCbAKr6lap2Lf9ZZS6rhYj8W0T2hPr/FJFE/+0WjorIPhH5Vanx4f7bepz030amXTjPVdWVQL6IfONC6jbRz5qTqQpbgZtPfyMivYDkC11YZaSlquL/g/oS8BzQDGgLvI5zLTe33AH8V1VPVcKyfEA28J0yxp8A0oF2QAbwgIiMAvB/YPkn8BjQEFgM/D2c5/r9FfhpJbwGE41U1R72iNgD2AY8CiwKmDYBeARQoL1/2tdxzsw/CuwEngiYv71/3tNnos8KmBbnn+c7/nX1xPnQ9RCwGTgEfAw09M+3w/+84/7HAJw/1rMD1qfA3UAucAR4jbMnrMcCLwAHcZruzwPrKPXa0/zruLGc9ycRp3nt8T9eAhL9Y42BL4F8nNsVfOV/be/jNIVT/uU/ACQBH/hfbz6wCOeq0aHWOQ24NeD7oTjX4gv8PxsNrAQKcBpGUgX/z3GB/58B03cDIwK+fxr4yP/1T4C5AWMp/td0SUXP9X/fyj9/ots/5/ao/IclJ1MV5gOpItJNRGKB7+H8IQ10AmdTU32cRvW/IlL6/jBDgG7AyMCJInInMA7IVNXVODeO+6Z//pacbTDg3OcIoL6q1lXVeWXUfB1wBc6N+m4KWOddwLVAX5yLqZZ3D5sBOE3js3LmeQTnhpN9KHP4lQAABJBJREFU/eu6EqeZA9yPcxHXJjip69eAquoPcZrsN/yv4bc41zZLA9oAjXCaa1nJqBdQ0T6mm3BuLtcB6I3TwM+L/2rfLYEVAZNXAD38X/cIHFPVEzgfKHqE8VzUudVFCXBBmyRNdLPmZKrK+zjNJwtYj/Op+AxVnaGqq1TVp87+hA9xmkugJ1T1hJ67OeqXwBhgqKpu8k/7KfCIqu5S1SKczUPfPc/Ngc+rar6q7gCm4zQPcP5ov+xf9hGc64uVpRFwUFU95czzA+ApVT2gqnnAkzjXXQTnD28LoJ2qlqizb6is642V+NfXWVW9qrpEy74NS33gWDk1AfxeVfeo6mHgC86+/vNR1/9vQcC0AqBewHgB5zo9XtFzTzuG83pMDWPNyVSV94FbcD6B/6X0oIj09+8QzxORApxP/o1LzbYzxHLHAK+p6q6Aae2Az0QkX0Tyca4w7cVJH+HaF/D1Sc7+sWxZqo5QNZ12COdCmeU1xZbA9oDvt/unAYwHNgGTRWSLiJR3T6b3gUnAR/6DE37rvy1LKEcI/iNfWlmv/3wc9/+bGjAtlbON8XipscDxip57Wj2czZimhrHmZKqEqm7H2UfzNZyd4KX9DecS/G1UNQ14A5DSiwnxvBHAoyISuEN+J3CtqtYPeCT5NwNd7JWO9wKtA75vU86884BCyt/0twenmZ7W1j8NVT2mqverakfgG8CvRGS4f75zXsf/t3fvoFFFQRjH/1/hA0HQIihaCL4Iotj5KGMXCSKkESztRIiIYBdREgVttFCCaSSNEBBtFAkiWogpBDGyYGNhJRrwEUEslLGYsxA2JtldNbmG71dlc85dzmXhzs6cWU7JrM5FxA7yLKIeMlP9nQlg+xxr+itKZvmOLFfW7QZq5e/a9LFyWu8WoNbEtUjaACxn/hKl/YccnGwhHQMOlL2FRquBjxHxXdIeMstqRo3cG7km6VD53xAwWG9LltQhqd4dN0k2E2xu8x5GgT5JGyWtAc7MNjEivgD9ZW2Hy3lByyR1S7pUpt0ig2tH6V7rp+zHSeqRtLWcMTRFZn8/y3Xvp9+DpC5Ju8qe3hRZ5qvPbXSfmSXTtklaSTZ2AKwor+tGyPtbK6mT3LO7WcbuADsl9ZZr+oGJiHjdxLWQjRyPSunWlhgHJ1swEfEmIp7PMnwcOC/pK/mQGm3hfV+SmcKwpG7gKpmFjZX3Gwf2lrnfgEHgaSn77WvxNoaBMTL7eEE+6H8wSyCIPEzyFNnkMElmdSeAu2XKANlCPQG8Is8UGihj24CHZInrGXA9Ih6XsYvkg/uzpNPAevLE2ymyjPmEmU0ndSPAQUltt/M3qHcNQu4nTt8TPEs2Obwta7ocEQ8Ayh5bL/l5fCI/oyPNXFscJb+I2BLk85zM/kAJhkMRsWneyRUi6QLwISKuLPZa2lF+K3cjIvYv9lrs33BwMmtByTa6yOxpHXAbGI+Ik4u6MLMlxsHJrAWSVpElpk6yfHUP6JujbdvM2uDgZGZmleOGCDMzqxwHJzMzqxwHJzMzqxwHJzMzqxwHJzMzqxwHJzMzq5xffOethJLCVWUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fe = smf.ols(\"purchase ~ mkt_costs + C(city)\", data=toy_panel).fit()\n",
    "\n",
    "fe_toy = toy_panel.assign(y_hat = fe.fittedvalues)\n",
    "\n",
    "plt.scatter(toy_panel.mkt_costs, toy_panel.purchase, c=toy_panel.city)\n",
    "for city in fe_toy[\"city\"].unique():\n",
    "    plot_df = fe_toy.query(f\"city=='{city}'\")\n",
    "    plt.plot(plot_df.mkt_costs, plot_df.y_hat, c=\"C5\")\n",
    "\n",
    "plt.title(\"Fixed Effect Model\")\n",
    "plt.xlabel(\"Marketing Costs (in 1000)\")\n",
    "plt.ylabel(\"In-app Purchase (in 1000)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take a minute to appreciate the figure above and let it tell you what fixed effect is doing. Notice that fixed effect is fitting **one regression line per city**. Also notice that the lines are parallel. The slope of the line is the effect of marketing costs on in-app purchase. So the **fixed effect is assuming that the causal effect is constants across all entities**, which are cities in this case. This can be a weakness or an advantage, depending on how you see it. It is a weakness if you are interested in finding the causal effect per city. Since the FE model assumes this effect is constant across entities, you won't find any difference in the causal effect. However, if you want to find the overall impact of marketing on in-app purchase, the panel structure of the data is a very useful leverage that fixed effects can explore. \n",
    "\n",
    "## Time Effects\n",
    "\n",
    "Just like we did a fixed effect for the individual level, we could design a fixed effect for the time level. If adding a dummy for each individual controls for fixed individual characteristics, adding a time dummy would control for variables that are fixed across time. One example of such a variable is inflation. Prices and salary tend to go up with time. If the wage marriage proportion also changes with time, we would have time as a confounder. To give a more concrete example, suppose that marriage is increasing with time. Since inflation also makes salary increase with time, some of the positive association we see between marriage and wage would be simply because both are increasing with time. To correct for that, we can add a dummy variable for each time period. In `linear models`, this is as simple as adding `TimeEffects` to our formula and setting the `cluster_time` to ture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Parameter Estimates</caption>\n",
       "<tr>\n",
       "     <td></td>     <th>Parameter</th> <th>Std. Err.</th> <th>T-stat</th>  <th>P-value</th> <th>Lower CI</th>  <th>Upper CI</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>expersq</th>  <td>-0.0062</td>   <td>0.0008</td>   <td>-8.1479</td> <td>0.0000</td>   <td>-0.0077</td>   <td>-0.0047</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>union</th>    <td>0.0727</td>    <td>0.0228</td>   <td>3.1858</td>  <td>0.0015</td>   <td>0.0279</td>    <td>0.1174</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>married</th>  <td>0.0476</td>    <td>0.0177</td>   <td>2.6906</td>  <td>0.0072</td>   <td>0.0129</td>    <td>0.0823</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hours</th>    <td>-0.0001</td>  <td>3.546e-05</td> <td>-3.8258</td> <td>0.0001</td>   <td>-0.0002</td> <td>-6.614e-05</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod = PanelOLS.from_formula(\"lwage ~ expersq+union+married+hours+EntityEffects+TimeEffects\",\n",
    "                            data=data.set_index([\"nr\", \"year\"]))\n",
    "\n",
    "result = mod.fit(cov_type='clustered', cluster_entity=True, cluster_time=True)\n",
    "result.summary.tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this new model, the effect of marriage on wage decreased significantly from `0.1147` to `0.0476`. Still, this result is significant at a 99% level, so man could still expect an increase in earnings from marriage. \n",
    "\n",
    "## When Panel Data Won't Help You\n",
    "\n",
    "Using panel data and fixed effects models is an extremely powerful tool for causal inference. When you don't have random data nor good instruments, the fixed effect is as convincing as it gets for causal inference with non experimental data. Still, it is worth mentioning that it is not a panacea. There are situations where even panel data won't help you.\n",
    "\n",
    "The most obvious one is when you have confounders that are changing in time. Fixed Effects can only eliminate bias from attributes that are constant for each individual. For instance, suppose that you can increase your intelligence level by reading books and eating lots of good fats. This causes you to get a higher paying job and a wife. Fixed effect won't be able to remove this bias due to unmeasured intelligence confounding because, in this example, intelligence is changing in time. \n",
    "\n",
    "![img](./data/img/fixed-effects/time-travel.png)\n",
    "\n",
    "Another less obvious case when fixed effect fails is when you have **reversed causality**. For instance, let's say that it isn't marriage that causes you to earn more. Is earning more that increases your chances of getting married. In this case, it will appear that they have a positive correlation but earnings come first. They would change in time and in the same direction, so fixed effects wouldn't be able to control for that. \n",
    "\n",
    "## Key Ideas\n",
    "\n",
    "Here, we saw how to use panel data, data where we have multiple measurements of the same individuals across multiple time periods. When that is the case, we can use a fixed effect model that controls for the entity, holding all individual, time constant attributes fixed. This is a powerful and very convincing way of controlling for confounding and it is as good as it gets with non random data. \n",
    "\n",
    "Finally, we saw that FE is not a panacea. We understood two situations where it doesn't work: when we have reverse causality and when the unmeasured confounding is changing in time.\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",
    "Another 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'd also like to compliment Scott Cunningham and his brilliant work mingling 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
}