{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Simple example on using Instrumental Variables method for estimation" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%load_ext autoreload\n", "%autoreload 2" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import patsy as ps\n", "\n", "from statsmodels.sandbox.regression.gmm import IV2SLS\n", "import os, sys\n", "from dowhy import CausalModel" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading the dataset\n", "\n", "We create a fictitious dataset with the goal of estimating the impact of education on future earnings of an individual. The `ability` of the individual is a confounder and being given an `education_voucher` is the instrument." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "n_points = 1000\n", "education_abilty = 1\n", "education_voucher = 2\n", "income_abilty = 2\n", "income_education = 4\n", "\n", "\n", "# confounder\n", "ability = np.random.normal(0, 3, size=n_points)\n", "\n", "# instrument\n", "voucher = np.random.normal(2, 1, size=n_points) \n", "\n", "# treatment\n", "education = np.random.normal(5, 1, size=n_points) + education_abilty * ability +\\\n", " education_voucher * voucher\n", "\n", "# outcome\n", "income = np.random.normal(10, 3, size=n_points) +\\\n", " income_abilty * ability + income_education * education\n", "\n", "# build dataset (exclude confounder `ability` which we assume to be unobserved)\n", "data = np.stack([education, income, voucher]).T\n", "df = pd.DataFrame(data, columns = ['education', 'income', 'voucher'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using DoWhy to estimate the causal effect of education on future income\n", "\n", "We follow the four steps: \n", "1) model the problem using causal graph, \n", "\n", "2) identify if the causal effect can be estimated from the observed variables, \n", "\n", "3) estimate the effect, and \n", "\n", "4) check the robustness of the estimate. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOsAAAD7CAYAAACL3GNOAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOydeVxU9f7/XzMMMwzbDCIM+ypLKirgkhtommhiKmWa5VKZ2uZ27826/m51y0rLbvm1ump9S+3mVqkJikqmrCoooKCyyCIg+zLs6/D+/eHlfB1nUMCZOTNwno/HeSjnfPh8XsyZ1/l8zvucz/vDIyICBweHvpPEZ1sBBwdHz+DMysFhIHBm5eAwEARsC9AlbW1tqK2tRV1dHeRyOdrb29HQ0MAcb2pqQmtrK/OzWCyGiYkJ87OFhQWMjY0hlUohkUhgaWkJY2Njnf4NHH2HiFBWVoaysjLm/Dc2NqKtrQ0ikQimpqYQiUSQSqWQyWSwtbVlW7ISBm/W1tZW5OXlIS8vD2VlZSguLkZZWRlKS0tRUlKC8vJy1NTUoK6uDi0tLRpvXywWQyKRQCqVwtbWFo6OjrC1tYW9vT2zubm5wc3NjTO2jmhvb8eVK1dw7do1pKen4/r168jKykJZWRna29t7XI9QKIS9vT28vb0xfPhwDBs2DAEBARgxYgSMjIy0+Beoh2co0eCCggKkpaUhPT0dOTk5zFZYWIjOzk4Ad43TZRCZTAYHBwfIZDJIpVJYWloyW5e5jIyMIJFImDa6rq5dNDQ0KJ1cuVwOhUKB2tpayOVy1NXVMVtNTQ3Ky8tRVFSE8vJy5qLRdYEQCARwcXGBp6cns40YMQJ+fn6wt7fX0afYf7l69SoiIiIQHR2NhIQENDY2wtLSEkOHDsXw4cPh4+MDBwcHODk5QSaTYdCgQRAIBExv2tLSgubmZrS1taG6uhqlpaW4c+cO7ty5g8zMTKSlpeHmzZtMvZMnT0ZwcDCefvpp+Pj46OJPTNI7sxIRMjIykJCQgJSUFKSlpeHatWuQy+UAAGdnZwwZMkTpS+/p6QkPDw9IpVKW1atSVVWFvLw83Lp1S+kik52djZKSEgDA4MGDMXLkSPj5+cHf3x/jx4+Hl5cXy8r1n7S0NOzbtw9HjhxBbm4uHBwc8MQTTyAoKAhBQUEaN1FnZydu3LiB6OhoxMTE4Pz58ygvL8fQoUMRFhaGpUuXavO8sW/WtrY2XLx4EXFxcUhISMCFCxdQXV0NMzMzjBgxgtn8/Pzg5+enl4bsK1VVVbh69SrS0tKYi1JaWhpaWlpga2uL8ePHY+LEiZg4cSLGjh0LgcDg71oemba2Nhw+fBg7d+5EfHw8PD098cwzzyAsLAxjx44Fj8fTmRaFQoH4+HgcPXoUv/32G4qKijBt2jSsXr0a8+bN0/RQOQnEAjk5ObRr1y5asGABSSQSAkD29vYUGhpKW7ZsodjYWGppaWFDGuu0t7fT5cuX6auvvqIlS5aQi4sLASAzMzOaPn06ffXVV5Sfn8+2TJ3T2tpKe/fuJU9PTzIyMqLQ0FCKioqizs5OtqUREZFCoaCoqChasGABCQQCcnd3p127dlF7e7ummkjUmVlTUlLonXfeIU9PTwJAlpaWNH/+fNq5cyfl5eXpSoZBkpGRQdu3b6ennnqKTE1NCQANGzaMPvjgA7px4wbb8rTO4cOHydnZmUQiEb3xxhtUWFjItqQHkpWVRcuWLSOBQEA+Pj506tQpTVSrXbNmZmbSe++9Rz4+PgSA3Nzc6O2336aYmBhNXnEGFC0tLRQVFUVr1qwhBwcHAkAjRoygjz/+uN/1uLdu3aKQkBDi8Xj00ksv6b1J7yc7O5ueeeYZAkDPPvssFRUVPUp1mjdrR0cHHT9+nKZPn048Ho8cHR1pzZo1FBsbqzdDlv6CQqGg2NhYWrNmDdnZ2RGfz6fp06fT4cOHDf5iuHfvXjI3Nyc/Pz+KjY1lW84jERkZSUOGDCFra2v6/fff+1qN5sxaU1ND//znP8ne3p74fD7Nnj2bIiIiSKFQaKoJjgfQ1tZGv/zyC02bNo14PB65ubnRtm3bqKGhgW1pvaKlpYWWLl1KPB6PNmzYQK2trWxL0ggNDQ308ssvE4/Ho7Vr11JHR0dvq3h0s8rlcvrggw9IKpWSlZUVvfvuu/1uOGZoZGRk0Lp168jc3JxsbW3p888/NwjT1tTUUHBwMEmlUjp58iTbcrTC/v37SSwW09y5c6mpqak3v9p3s7a1tdHnn3/OmPTDDz8kuVze1+o4tEBFRQVt3LiRMe3OnTv1dqRTUVFBfn5+5OTkRGlpaWzL0SoJCQlkbW1NEydOpPr6+p7+Wt/MGhsbS35+fmRiYkLvvfceZ1I9p7y8nNavX08CgYAef/xxSklJYVuSEvX19TR27Fjy8PAwuCBSX7l58ybJZDKaMWNGT4f6vTNrc3Mzvfbaa8Tj8SgkJISys7P7ppSDFa5du0YTJ04kgUBA7777bl/umzSOQqGgWbNmkUwmG3DfpytXrpClpSUtXbq0J8V7btb8/HwaPXo0SSQSOnDgQN8VcrBKZ2cn7dy5k8RiMU2dOpXKyspY1bN161YSCoWUmJjIqg62iIyMJD6fT3v27HlY0Z6ZNSYmhqytrcnPz4+ysrIeXSEH66SkpJCHhwc5OjpSamoqKxqSk5NJKBTS559/zkr7+sJf//pXMjc3p5ycnAcVe7hZ4+LiyNzcnMLCwqixsVFzCg2Qo0ePEgBma25uZlvSI1FdXU1PPPEE2djYsBLUCQoKokmTJmk96FVYWKh03rq2o0ePKpXbtGmTSpmbN29qVRvR3Vcphw4dSgsWLHhQsQeb9dKlS8xrgW1tbZpVaMDMnTu3X5iViKixsZGmTJlCMpmMMjIydNbusWPHiMfj0cWLF3XW5oEDBwgAbdy48YHlgoOD6bvvvtORqrscP36ceDweXbhwobsiid2mdamrq8PChQsxYcIEHDx4kJs43U8xNTVFREQE3NzcsGjRIrS1temk3S+++ALz5s3DuHHjdNKevjNnzhyMHz8eX3zxRbdlujXrunXr0NjYiD179kAoFGpFIId+YGZmhgMHDiAnJwcffPCB1tvLzs5GXFwcVq1apfW2DImVK1fi+PHjqKioUHtcrVnj4+Px448/Yvfu3ZDJZFoVyKEfuLu7Y9u2bfjss8+QkZGh1bYOHjwIBwcHTJ8+XavtGBoLFiyASCTCkSNH1B5Xa9bPPvsMkyZNwrx587QmTC6Xg8fjKW2bN28GAHR0dCjtf/bZZ5nfq6qqwoYNG+Dp6QmhUAgrKyvMmjUL586dY8ps3ryZ+d1JkyYx+0+dOsXsHzx4sIqme+sWiURwcnLC9OnTsWfPHjQ3N6uULy0txcKFCyGVSmFtbY3Q0FDk5OSolKuoqMCaNWvg5uYGoVAIGxsbhIWFITU1lSlz7Ngxpb85MzMTzz33HKytrZl9lZWVffuwe8iKFSvg5eWFf/3rX1ptJzY2Fk888QQreYz0GVNTU4wfPx5xcXHqC9x/FyuXy0koFNK+ffu0e0f9X2bOnEl8Pp9u3bqlcmz8+PG0f/9+5ueSkhJyd3cnmUxG4eHhVFtbS5mZmRQWFkY8Hk8lKGBmZkYTJ05UqTcwMJCsra2V9nXVbWdnR+Hh4VRXV0elpaX00UcfEQD68ssvmbJdAaa5c+dSQkICNTQ00NmzZ8nS0pLGjBmjVG9xcTG5urqSTCajEydOUH19PaWnp1NwcDCZmJhQQkKCUvmuuoODg+ncuXPU2NhIFy9eJCMjI6qoqOj5B9tHtm3bRlZWVlqbtaNQKMjS0pK+/fZbrdT/IPQ5wNTFBx98QJ6enuoOqUaDz5w5QwB09rD8jz/+IAD0+uuvK+2Pi4sjFxcXpS/N8uXLCYDKSxktLS3k4OBAYrGYSktLmf29MWtX3YcOHVIpP3PmTLVmDQ8PVyq3ePFiAqBkqmXLlhEA+vnnn5XKlpSUkEgkosDAQKX9XXWz9SJ7WloaAdDas9eqqioCQGfPntVK/Q/CEMy6f/9+EggE6qaTqkaDc3NzMWjQIJ3lTJ02bRr8/f2xZ88eVFVVMfs///xzrFu3Tinv0NGjRwEAs2fPVqpDJBJh2rRpaG5uxunTp/uko6vuWbNmqRyLjIzEunXrVPaPGTNG6WdHR0cAQHFxMbPv2LFj4PP5CA0NVSprZ2eHYcOG4cqVKygqKlKpe+zYsb3/IzSAj48PeDye2uG8JugayltbW2ul/gfRNexWKBQPLKdQKFgboltbW6OjowO1tbUqx1TM2tTUBLFYrBNhXfzlL39BU1MTvv32WwBAVlYWYmJisGLFCqZMa2sramtrYWJiAgsLC5U6ugJhpaWlvW7/YXV3x71pTAGAz7/7cXalRu2qt7OzExKJROUePTk5GcDd6Oj9mJmZ9frv0ATGxsYQCoVoamrSSv1d9/66/o4BgLm5OYC7jyUfhFwuh6WlpS4kqdClsbGxUeWYilkHDRqEmpoa5gunCxYuXAhnZ2d8/fXXaG1txRdffIFXX31VyTgikQgSiQQtLS2or69XqaOsrAzA3R6rCz6fr/a5YVda057W3Ve6srsLBAK0t7eDiNRuU6dO1Vibj0pdXR1aW1sxaNAgrdRvZWUFAKipqdFK/Q/C29sbAHD9+vVuy7S2tuLWrVuspYLtGl2q+/xVzDp8+HA0NTXh5s2b2lf2XwQCAdauXYvy8nJ88cUXOHjwINasWaNSbv78+QCAEydOKO1vbW3F2bNnIRaLERISwuy3t7fHnTt3lMqWlpaioKCg27pPnjypcszf3x/r16/v/R8GICwsDB0dHYiPj1c5tnXrVri4uKCjo6NPdWuDK1euAAD8/Py0Un/X8Le7Z4naxNPTE76+vrh48aLa0QwAHD58GDY2Nhg+fLiO1d2loqICZmZm6kce99/FdnR0kEwmo48++kjb99JK1NXVkUQiIR6P1+2UofujwXV1dUrR4N27dyuVf/PNNwkA7dixg+rr6+nWrVv03HPPkaOjY7fRYHt7e4qIiKC6ujoqLCyk1157jWQyGd2+fZsp293rhhs3biQASvNFy8rKyNPTkzw8POjkyZMkl8upqqqKdu7cSaampioBLbZfZXzzzTdp6NChWm3Dzc2NPvzwQ6220R2RkZFkbGxMnp6e9Ntvv1FVVRV1dHTQnTt36JtvviFLS0v65ZdfWNFGRPTGG2/QuHHj1B1S/27w3/72N7K3t+9t2olH5m9/+xsBoKtXr3ZbprKyktatW0fu7u5kbGxMEomEQkJC1EYX5XI5rVixguzt7UksFtOkSZMoKSmJAgMDmRe1740M3l+3vb09LVq0iJlpdOHCBZUXvTdt2kREpLJ/9uzZTL1VVVW0YcMG8vDwIGNjY7KxsaEZM2ZQVFQUU0Zd3WqupVqlqqqKLCwstD4LZvHixTRz5kyttvEgrly5Qi+++CK5ubmRSCQioVBITk5OtGDBAoqPj2dNFxFRQEAArVu3Tt0h9WYtLi4miURCf/nLX7SrjEOvWLJkCdnb21NdXZ1W2/nf//1fEovFVFNTo9V2DI28vDzi8/kqjwT/S/ezbn788Ufi8/l07tw5rYnj0B+OHj1KPB6PIiIitN5WfX09mZubs/JihD7z3nvvkUwm626G24OnyM2bN4/s7Ox0OnWKQ/ckJiaSRCKh1atX66zN5cuX02OPPWbw+Y01RUNDA9nb29M777zTXZEHm7WhoYGCgoJIJpPpZBIuh+5JTU2lQYMGUUhIiE7XF8rJySGRSES7du3SWZv6zD//+U+ysLBQegPvPh6eKaK2tpbGjh1LDg4OrN98c2iWyMhIGjRoEM2YMYOV6POaNWtIJpNReXm5ztvWJ3Jzc8nCwoI2b978oGI9y8Ekl8spNDSUjI2Nafv27ZpRyMEaCoWC3n//feLz+bRkyRKdR/27kMvl5ObmRrNmzRqwS6u0t7fThAkTaPjw4Q87Dz3PbtjZ2UmbN28mIyMjmjt3LhUUFDy6Ug6dk5mZSVOnTiWRSKQXAZ64uDgSCAS0ZcsWtqWwwrp168jU1JSuX7/+sKK9T/J97tw58vLyInNzc9q2bRsXIDAQmpub6f333yeRSESjRo2ipKQktiUxbN++nXg8Hu3du5dtKTrls88+Iz6f39PUvn3LyN/c3EwffPABmZiY0NChQ+nQoUN6uyzDQKe9vZ1+/PFH8vDwIAsLC/rXv/6llxfYv//97yQQCJTmL/dnui5QvbitfLSFqbKzs+n5558nPp9Pw4cP50yrR7S3t9MPP/xAnp6eZGxsTCtWrHjU9UG1SmdnJ/31r38lHo/Xr/MId3Z20ttvv008Ho+2bt3am1/VzJKPN27coCVLlpCRkRF5enrSli1bBnyEjy1KSkpoy5Yt5OrqSsbGxrRkyRKDWpbiyy+/JD6fT8uWLevNok0GQVVVFc2dO7evmVg0u5hyRkYGvfXWWySRSMjExISWLFlC0dHRXG+rZTo6Ouj06dP07LPPMu8ev/322wa79GZERAQNHjyYfHx86PLly2zL0Qjnz58nZ2dncnZ2ppiYmL5UofmVz4nuvkyxe/duCggIIADk6OhI69ato4SEhAEbotc0HR0d9Oeff9Lq1avJxsaGANCECRPop59+0unLDdqiqKiIpk6dSkZGRrRmzRqDXamwrKyMli9fTjwej+bOnUuVlZV9rUo7Zr2X9PR0eu+998jX15cAkIuLC7322mt07NixfjfM0TbV1dV0+PBheuWVV8jOzo4A0MiRI+njjz9Wm3DO0Ons7KQffviBbGxsyM7Ojnbs2GEwqyDU19fTp59+SlZWVuTs7KyJaXfaN+u9pKam0v/7f/+PRo8eTXw+n4RCIU2ZMoU+/fRTiouLM5gToSsaGhro3Llz9MEHH9D48ePJyMiIjIyMaMKECbR58+YB8852dXU1vfTSSyQQCMjR0ZG++uorve1pKyoq6OOPP6bBgweTpaUl/eMf/9DUqvO6Neu9lJeX088//0xLly4le3t7AkBCoZDGjx9PGzZsoF9//ZXy8vLYkqdzFAoFZWdn0/79++mtt96iwMBAEggEBICcnZ3plVdeoV9++YWqq6vZlqpzzp07R9bW1jRixAhavXo1mZmZkZmZGb3yyit04cIF1m+tFAoFnT9/nhYvXkwikYgkEglt2rSJqqqqNNlMIo+ISFspKnpDbm4uEhIScOHCBcTHxyM9PR0KhQISiQR+fn4YMWIERowYAT8/PwwZMkRn2Re1QUlJCbKyspCeno5r167h2rVrSE9PR0NDA4yNjTFq1ChMmDAB48ePx4QJE+Ds7My2ZNb4/vvv8cYbb+Dpp5/G3r17YWpqitraWvz000/YtWsX0tPT4erqivnz52P+/PkYP368TtZlamlpQVxcHI4cOYJjx46hpKQEY8eOxapVq7Bw4UJtJLxL0huz3k9DQwPzRb569SquXbuGtLQ0JqGZhYUFPD094enpCQ8PD7i7u8PR0RG2trZwcnKCra0tK2v0tLS0oKysDHfu3EF5eTmKioqQl5eH3Nxc5OTkICcnh8kcKJVKMXLkSOZiNGrUKAwfPpyVzH/6hkKhwKZNm/DZZ5/h7bffxieffMJkj7yX5ORkHDlyBEeOHMHNmzdhamqKxx9/HEFBQQgMDMTw4cPh6uoKHo/XZy2dnZ3Iy8tDWloaLl++jJiYGCQmJqK1tRWjRo3CM888g7CwMAwdOvRR/uSHob9mVQcR4fbt28yXvmvLzc1Ffn6+SsY8GxsbyGQySCQSWFpawtLSEhKJBFKpFJaWlhAIBDAzM2NMbWRkpJSCUi6Xo+vjaWlpQXNzM9rb21FXV4fa2lrI5XLU1dWhrq4OcrkcpaWlqK6uVtIwePBguLm5MReVrguMp6fngO4xH0R9fT1eeOEFREVF4bvvvsOLL77Yo9+7desWYmJiEB0djbi4OOTm5gK4e2H39vaGo6MjHB0dYWdnB4lEArFYDBMTE4jFYjQ1NaG1tRWNjY2oq6tDcXExiouLUVRUhKysLDQ1NYHH48Hb2xuTJ09GUFAQgoOD4eLios2P4l4My6wPo6WlBaWlpSguLkZ5eTnu3LmDsrIy1NbWMqaqq6tDTU0Nk8+3rq6OSfrc3t6OhoYGpr6uHK4CgQBCoRDm5uYQCASwtLSEVCplLgISiQQSiQQymQwODg6QyWRMLy8SiVj5LAyV3NxczJkzBzU1NTh27NgjJTuvra3F9evXkZ6ejuzsbJSWlqKoqAhpaWlMqt3m5ma0tLTA1NQUIpEIZmZmsLS0hL29PRwcHODg4ABvb2/4+flh6NChrOVzBpDEWoDJEMjMzCQAlJyczLaUAUFcXBzZ2trSqFGjlLJJahoXFxf69NNPtVa/luh+MWUOwNXVFXw+H3l5eWxL6fd8//33eOKJJxAUFIT4+HitDS9ra2tRWFiotbzI2oQz6wMQiUSwt7fnzKpFFAoF3nnnHaxcuRLr16/HoUOHYGpqqrX20tPTQUQGaVbBw4sMbNzd3Tmzaol7A0n79u3rcSDpUbh27RokEolBBvc4sz4Ezqza4d5AUnR0tM5WzUtLS4Ofn98jPcphC24Y/BDc3d2Rn5/Ptox+RXx8PMaPHw+hUIiLFy/qdHnLLrMaIpxZH0JXz0r95wkXq+gqkNQd6enpnFn7K+7u7mhubkZ5eTnbUgwaXQeS1FFQUAC5XG6wZuXuWR+Cu7s7ACAvL49ZsJmjd7ARSFJHWloaALC2nOOjwvWsD8HR0RHGxsZckKmP5Obm4vHHH8fly5cRHR3NmlGBu2Z1cXGBVCplTcOjwJn1IRgZGcHZ2Zkzax9gM5CkDkMOLgGcWXsE9/im97AdSFIHZ9YBAGfWnqMPgSR1tLe3IzMz06DNygWYeoC7uzv+/PNPtmXoPfoSSFJHRkYG2traOLP2d9zd3VFYWAiFQgEjIyO25eglbL2R1FPS0tJgbGwMHx8ftqX0GW4Y3APc3d3R3t6OoqIitqXoJfoWSFJHWloafH19Wckeoik4s/aAe5+1ciijj4EkdVy7ds2gh8AAZ9YeIZPJYGZmxpn1HvQ1kNQdhh4JBrh71h7j6urKvdD/X/Q5kKSO2tpaFBUVcWYdKHCPb+5ybyApJiYGY8aMYVvSQ0lLSwMRGexrhl1ww+AewplVNZBkCEYFgJSUFFhZWent/XRP4czaQwa6WQ0lkKSOlJQUBAQEGOSE83vhzNpD3NzcUFxcjJaWFral6BRDCySpIzk5GQEBAWzLeGQ4s/YQd3d3EBEKCgrYlqIz6uvrMX/+fGzfvh379u3Dli1b1GbF12daW1tx48YN+Pv7sy3lkeECTD3Ew8MDwN1nrd7e3iyr0T6GGEhSR1paGtrb27medSDRtezGQLhvNdRAkjqSk5Nhbm4OLy8vtqU8MpxZe8FACDIZciBJHSkpKRg1apTBDd/VYfh/gQ7pz2btD4EkdfSX4BLAmbVX9Fez9odAkjoUCgXS09P7RXAJ4MzaK/qjWe/NkRQTE6P3rw72hhs3bqCpqYnrWQci7u7uqKqqYhZ0NnT6UyBJHcnJyTAxMcFjjz3GthSNwJm1F3RNlesPL/T3t0CSOlJTUzF8+HAYGxuzLUUjcGbtBW5ubuDxeMxQuLGxEenp6QgPD8edO3dYVtcz+msgSR2JiYn9arTAvRTxANra2lBQUIC8vDxmk0qlWLt2LZYvX46amhqmrCG82WRoU9sehY6ODqSkpGDFihVsS9EYnFkfwIYNG/DNN98AAAQCAYyMjNDe3q5kUgAYPHiw3i8h2F/eSOop6enpaG5u1ssUM32FGwY/gHfeeYe53+no6EBrays6OzuVyvD5fIwfP54NeUp8/vnn6OjoUHusvweS1JGYmAhzc3P4+vqyLUVjcGZ9AE5OTli1atUDAxQCgYB1sx4/fhxvv/021q9fr3JsIASS1JGUlISAgID+lY2SOB5IcXExCYVCAtDt9scff7Cmr6mpiZycnIjH4xEA+ve//01ERB0dHbRx40bi8Xi0ceNGUigUrGlkg5EjR9Jf//pXtmVokkTOrD3gzTffJGNjY7VG5fF4VFNTw5q29957jwQCAaPHyMiIjh8/TnPmzCETExP66aefWNPGFo2NjSQQCOjQoUNsS9EkiTwibpXgh1FSUgI3Nze0tbWpHBsyZAiys7NZUHU3aOTr64v29nZmH5/Ph1gsho2NDX755ReMHj2aFW1sEh8fj0mTJiE3N5d5Nt4PSOLuWXuAvb292ntXgUCAiRMnsqQKePPNN1X2dXZ2orW1FTweD0OGDGFBFfskJSXB2toabm5ubEvRKJxZe8i7776rNocPW48GTp06hcjISKVetYuOjg4UFRUhLCys2whxfyYpKQljx441+JxL98OZtYeo6107Ojowbtw4nWtpaWnBypUrHzgzpr29HTExMWojxP2d/vbmUhecWXvB/b2rUChkJXH01q1bUVJSovLMVx1ff/01Dh8+rANV+kFNTQ1ycnI4sw507O3tsXr1aqZ39fPz0/lCRwUFBfj000+7Hd4aGRmBz+fD1NQUzz//PKKiorBgwQKdamSTxMREEBFnVg5g48aNzP/ZCC69/vrrantUgUAAHo+H0aNH49///jfKy8vx008/Yfr06f3u3u1BJCQkYMiQIZDJZGxL0TgD7tFNc3MzKisrUV5ejvLyclRWVqKyshIVFRWor69HY2MjGhoaIJfL0djYiKamJtTV1YGIIJfLAQBNTU1obW2FkZERFAoFU7dQKISZmZlSe1ZWVgDuJlwzMzODmZkZLC0tYWFhwfx/8ODBGDx4MGxsbCCTyWBjY4PBgwdDJBIp1XX69GnMnDmT+dnY2Bjt7e1wc3PDsmXLsGzZsv70qKJPPPnkk3BwcMDevXvZlqJpkvrVi/wKhQKFhYUoKCjA7du3cfv2bRQUFDBbYWEhGhoaNNrevbS1tak8i73/pf/eIJFI4OTkBFdXVzg6OuKXX34BAPB4PJiYmGDx4sV4+eWXMWHChD630Z9QKBRITEzE1q1b2ZaiFQyyZ62vr0dGRgYyMzOZfzMzM5GVlYXW1laNtGFpaQlzc3Om9zMyMoKlpSVzPDc3F/7+/vx/c8QAACAASURBVEoR2Y6ODpUsEnK5nOmVGxsb0djYqLFME6ampvDx8YG3tzd8fX3h6+sLHx8f+Pr6QiwWa6QNQyI1NRX+/v64evUqRowYwbYcTaP/PWttbS2Sk5OZ7cqVK8jOzu5RJLQLsVgMV1dXODk5QSaTMUNOW1tb2NraMsNQS0tLxpwPo6OjAwJB3z++2tpaNDY2oq6uDhUVFaisrERZWRnz/8rKSpSUlKCoqAgFBQVqL0JNTU1ISUlBSkqK0n6BQABfX18EBAQw26hRo2BhYdFnvYbAhQsXYGlpiWHDhrEtRSvoXc+am5uL2NhYxMTEIC4uDtnZ2XiYRD6fD1dXV/j4+MDHxwfu7u5wcXFhNhsbGx2p1x4lJSXMEL+goAC5ubnMaKInE9/5fD58fX0xefJkTJ48GUFBQXo/B7e3LF26FKWlpThz5gzbUrRBEutmvXPnDk6ePIlz584hJibmoelR3N3dERAQgJEjR8LX1xfe3t7w8fGBiYmJjhTrH42NjYxxb968idTUVCQnJ6OoqOiBv+fm5oagoCBMmzYNM2fOhK2trY4Ua4chQ4bgxRdfxAcffMC2FG2ge7N2dnYiMTEREREROHnyJFJTU7vtOZ2dnTFhwgQEBAQgMDAQAQEBTHSV4+GUl5cztw5XrlzBhQsXUFpaqrYsn8/HmDFjEBoaiqeeegr+/v4G9cinvLwcMpkMp06dQkhICNtytIFuzEpEuHjxIg4cOIDDhw+jrKxMbTkfHx9miBYUFARXV1dtSxtwZGVlMbcZsbGx3eZBdnJywqJFi/DCCy9g1KhROlbZe44ePYpnn30W1dXVkEgkbMvRBto1a2ZmJv7zn/9g//79yM3NVTluYWGBGTNm4KmnnsJTTz0FOzs7bUnh6Ibbt28jMjISJ06cwJ9//ommpiaVMkOHDsXixYvx4osv6u0F9O2338apU6dw7do1tqVoiySNTz5XKBR04sQJCgkJYbIX3Ls5OjrS2rVrKSoqilpbWzXdPMcj0NTURCdOnKDXXnuNbG1tVc6dkZERzZ8/n86dO8e2VBUmTpxIq1atYluGNtFcpoimpibasWMHeXt7q5xkKysrWrFiBZ07d27ApRcxVNrb2+nkyZP04osvkrm5uco5HTlyJP3www/U1tbGtlRqbW0lsVhMP/74I9tStMmjm7W9vZ2+++47cnJyUjmhU6ZMoV9//ZVaWlo0IZaDJRobG2nfvn00ZswYlXPs5eVFBw4coM7OTtb0xcXFEQDKzs5mTYMOeDSz/vbbb+Tj46N08kxMTOiVV16h1NRUTYnk0CMuXLhAixYtUslJ5e/vT2fOnGFF08cff0z29vastK1D+mbWsrIyCgsLUzpZIpGI1q1bR2VlZZoWyaGHFBQU0CuvvKKUrA0AvfTSSySXy3WqZcaMGbR48WKdtskCvTfr4cOHycbGRino8NJLL9Ht27e1IZBDz7l58yY9++yzSsFEJycnOnXqlE7ab29vJwsLC9q1a5dO2mORnpu1s7OT3n77baWr6PDhw+ny5cvaFMhhIJw7d47c3d2VUrRu3bpV6+1evHiRAFBGRobW22KZnpm1tbWVFi9ezJwIgUBAGzdu1MvA0YEDB5SG5hy6o7GxkTZu3Eh8Pp85BytWrKD29nattbl161aytbVlNcClIx5u1o6ODpo1axbz4UskElYz0PeUadOmGaRZ6+vraciQITR79my2pfSZX3/9lUxMTJjvzJIlS7TW1lNPPUXPPfec1urXIxIfmtZl06ZNiIyMBAA4ODjg/PnzmDZtWk/fuuBQg7m5OSZNmqT2GBGhs7OzV1MA9Y1nnnkG58+fZ2Y7/fTTT/jyyy813o5CoUBCQgKCg4M1Xrde8iAr//bbb0zgYNCgQZSTk6Oji8ijo889q5mZGU2cOJFtGVonOTmZTE1NmVun8+fPa7T+y5cvEwBKT0/XaL16Svc9a0NDA958800QEfh8Pv7zn//Aw8NDZxcRDsPH398f3333HYC7k/VXrVqlNil5X4mOjoa1tTUee+wxjdWpz3Rr1u+//x4lJSUAgPXr12PWrFk6E1VRUYE1a9bAzc0NQqEQNjY2CAsLQ2pqqkrZjIwMzJs3j0lINnnyZMTFxamU27x5M3g8Hng8ntIQ9NSpU8z+wYMHq/xeVVUVNmzYAE9PT4hEIjg5OWH69OnYs2cPmpubAdz9Ih46dAhPPvkk7OzsIBaL4efnh+3btysNZ7dt2wYej4fGxkbEx8cz7XZlnDh27Bizj8fjoaWlpVstQqEQVlZWmDVrFs6dO8eUub+O/Px8LFy4EFKpFNbW1ggNDUVOTk4vz0jfWbx4MZYsWQLg7sSOX3/9VWN1R0dHIzg4+IHJzvsV3fW5/v7+BIDEYjFVV1frrK8vLi4mV1dXkslkdOLECaqvr6f09HQKDg4mExMTSkhIYMpmZ2eTVColR0dHOnPmDNXX19O1a9doxowZ5ObmpnYY3N0QNDAwkKytrZX2lZSUkLu7O9nZ2VF4eDjV1dVRaWkpffTRRwSAvvzySyIiCg8PJwD0ySefUHV1NVVUVND//M//EJ/PV7vs4MOGwXPnziUA1NzcrKJFJpNReHg41dbWUmZmJoWFhRGPx6PvvvtObR1z586lhIQEamhooKioKBKLxTRmzJhu29YGubm5TIQ4JCREI3UqFAqytram7du3a6Q+A0B9NLimpob5cMPCwnSqaNmyZQSAfv75Z6X9JSUlJBKJKDAwkNm3YMECAkC//vqrUtk7d+6QSCR6ZLMuX76cAKhdOnDmzJlKZp0yZYpKmRdffJGMjY2ptra2Rxq6UGfWLi0HDhxQKtvS0kIODg4kFouptLRUpY7w8HCl8s8++ywBoIqKim7b1waTJk0iAGRmZqaRRzmpqakEgFJSUjSgziBQf8+am5vLDN90vWTgsWPHwOfzERoaqrTfzs4Ow4YNw5UrV5h0JadOnQIAlcwADg4O8Pb2fmQtR48eBQC1twCRkZFYt24dACA0NFRpKNrFyJEj0d7ejuvXr2tMy+zZs5X2i0QiTJs2Dc3NzTh9+rTK792fmb4r71JxcfEja+oNXd+jxsZG5vbqUTh79iysra37YxbDblGbnu/eCcjm5uY6E9Pa2ora2loAeOBs/+zsbNjY2KC+vh4mJiZqNdra2iIrK+uRtZiYmDw0K2BtbS2++OILHD16FEVFRUwy8C7UTejWpJau7PPqUrbc/zl2Lfeh60dD92aMbGxsfOT6zpw5g+nTpw+c+1V0E2AaNGgQ839dXoFFIhGkUikEAgHa29tBRGq3qVOnQiQSwcLCAi0tLWoTd1dXV6ttg8/nq10U+X6DiUQiSCQStLS0PDTP75w5c/DRRx/h1VdfRVZWFjo7O0FEzLNFui8ZR29zGz1MS1eaHH3OtHFvIjx1gbze0Nrairi4ODz55JOPKsugUGtWLy8vZhkIdcM7bdK1pmh8fLzKsa1bt8LFxYVZlKlreNo1HO6isrISmZmZauu3t7dXyaBYWlqqNp3n/PnzAQAnT55UOebv74/169dDoVAgPj4ednZ2WLNmDWxsbBgzdkWL78fU1FTpguHj44Pdu3erLXu/lhMnTijtb21txdmzZyEWi/U2URgRMd8jJyenRzZrXFwcGhsbB97LOd3dzXYFKABQcnKyDu6f71JWVkaenp7k4eFBJ0+eJLlcTlVVVbRz504yNTVVCvbcunWLBg0apBQNvn79OoWEhJCtra3aANObb75JAGjHjh1UX19Pt27doueee44cHR27jQbb29tTREQE1dXVUWFhIb322mskk8mYmUZPPPEEAaDPPvuMKioqqKmpif78809ycXEhABQVFaVU78yZM0kikVBBQQElJCSQQCCgGzduMMd7Eg2uq6tTigbv3r1bqQ11dRARbdy4UeeBmcjISOa79Oqrrz5yfRs3biQfHx8NKDMoun83+OTJk0oZH3SZjqWqqoo2bNhAHh4eZGxsTDY2NjRjxgyVLz0RUWZmJs2bN48sLS2ZxxIRERE0bdo0Rv8rr7zClJfL5bRixQqyt7cnsVhMkyZNoqSkJAoMDGTKb9y4kSlfWVlJ69atI3d3dzI2NiZ7e3tatGgRZWVlMWUqKipo1apV5OzsTMbGxiSTyWj58uX0zjvvMHXeG8XOyMigyZMnk5mZGTk7O9M333xDRERHjx5VycTwwgsvdKtFIpFQSEgInT17lilz4cIFlTo2bdpERKSyXxfvHzc1NZGfnx8zE+fKlSuPXKe/vz+99dZbGlBnUDz4Rf7g4GDmxH788ce6EsXRj1i6dCnzHdLEC/cVFRXE5/Pp+PHjGlBnUDzYrFlZWWRmZkYAiM/n0969e3UljKMfsGXLFsaoLi4uVF5e/sh17t+/nwQCgcqz6wHAg2fdeHl54YcffgCfz0dnZyeWL1/eb5fT49AcRIS//e1veOeddwDcXRjsyJEjGllzKCoqCuPHj+/R4mH9jp5Y+ttvv1WaULx69WqVwAUHBxFRXV0dLVy4kPmuCIVC+uWXXzRWv7OzM3344Ycaq8+A6Hlal6NHj5JYLGZOwtChQ+nSpUvaFMdhYMTExNCQIUOY74iFhYVGczHduHGDANCFCxc0VqcB0buEafHx8UqZ2rvSu9TV1WlLIIcBUFlZSatXr1ZKmubq6krXrl3TaDtfffUVSaVS6ujo0Gi9BkLvsxtWV1fTypUrlR4BWFtb05YtW7ih8QCjoaGBtmzZQlKpVOn7sGDBAq3M1Jo9e7bOJ5boEX1P8n348GGV9VBcXV1p9+7d1NTUpEmRHHpGbW0tffnllySTyZTOv4uLC50+fVorbTY1NZGpqSl9//33WqnfAHi0jPz19fW0ZcsWsrS0VDppEomE1qxZQ/n5+ZoSyqEH3Lp1izZu3EhWVlZK59vKyoq2bNmi1Yv077//Tnw+n0pKSrTWhp6jmYWpSkpK6I033iChUKh0EgUCAT3zzDP0+++/cyvGGShNTU108OBBeuqpp5SeCAAgU1NTeuedd3SSnGDFihU0btw4rbejx2huFTmiu6bdsmULOTo6qrzaJpVKacmSJRQVFcWtJKfnKBQKio2NpZUrV6qMmgCQvb09vf/++zqbwN7Z2UmOjo60efNmnbSnp2jWrF20tLTQjz/+SAEBASonGgA5ODjQihUr6MiRI1RfX68NCRy9pKamhg4ePEhLly5VWh7l3m3SpEl06NAhrSbtVselS5cIAF29elWn7eoZiVpd+RwArl+/jv3792P//v3Iz89XOS4SiRAUFIRZs2YhODgYI0eOhJGRkTYlceBukrfLly8jOjoakZGRiI+PZ6Ye3ouPjw8WL16MF154AZ6eniwoBf7xj39g3759yM/P7/Vc4H5EktbN2gURIT4+HgcPHkRERARu376ttpylpSUmTZqEyZMnIygoCKNHj2ayG3D0naamJiQmJiImJgaxsbG4cOFCtxkbvLy8MGfOHDz//PM6T+ujjlGjRmHSpEn4+uuv2ZbCJroz6/3k5uYiPDwcERERiImJUZu9AQAEAgG8vb0RGBjIbP7+/szkeA5VGhoakJqaihs3buD69eu4cuUKLl++jNbWVrXlBQIBxo0bhzlz5mDOnDkYOnSojhV3T2FhIVxdXREZGam3k+t1BHtmvZeamhrExsYyV/3k5GS1Q7IuBAIBhgwZgsceeww+Pj7w9fVl/v+g3E39jerqamRkZODmzZvIzMxk/n9vwjt1CIVCjBkzBkFBQZg8eTImT56s01xbveHbb7/Fxo0bUVFRARMTE7blsIl+mPV+GhoacOHCBcTGxuLy5ctISUlRmwxMHXZ2dnB3d4eLiwucnZ3h7OwMNzc3ODs7w9HREba2tlpWrxk6OztRUVGBwsJCFBYWoqCgALdv30ZBQQEKCwuRl5eHioqKHtXl5OQEf39/jB49GkFBQRg3bhzEYrGW/wLNMGvWLJiamuK3335jWwrb6KdZ1VFcXIzk5GRmu379OvLy8qBQKHpVj5GREWxsbDB48GAMHjwYdnZ2zM8WFhYwMzODlZUVzMzMYGZmBnNzc0ilUvB4PJiamkIkEjF1icVipat9c3OzUhb9lpYWNDc3o7OzE7W1tairq0NjYyMaGxtRW1uL+vp61NfXo7KyEuXl5SgvL0dFRQUqKytRWVnZ6wyExsbG8PT0xPDhwxEQEAB/f38EBAQYzAXqfhoaGmBjY4N///vfWL58Odty2MZwzKqOtrY2ZGdnIyMjA5mZmbh58yays7NRUFCgkdy0+gifz4ednR3c3Nzg4+PDbI899hg8PDxgbGzMtkSNceTIESxYsADFxcVMutUBjGGb9UG0trYyw8fCwkLk5+fjzp07KC8vR2VlJSoqKlBeXq6SgpRNrK2tmV7exsYGdnZ2cHBwgIuLC1xdXeHi4gInJ6d+ZcgHsXTpUuTk5KjNdDkA6b9m7SltbW3MsLO+vh5NTU2Qy+VoaGhghqw1NTUAVIe5ra2tSgm8zczMlB4zdQ2b+Xw+JBIJzM3NYWZmBlNTU7z//vvIyMgAAPzwww9YsmQJs0AVx91bCDs7O3z44YdYs2YN23L0Ac6sbHHlyhWMGzcOCoUCLi4uuHHjBvc46h66hsAFBQVwdHRkW44+kDRw1h7QMwIDA/HKK68AAAoKCrBlyxaWFekXhw4dQlBQEGfUe+B6Vhaprq6Gj48PKisrIRQKce3aNfj4+LAti3Wamppga2uLbdu2YfXq1WzL0Re4npVNBg0ahI8//hjA3Xtn7t7sLsePH0drayvCwsLYlqJXcD0ry3R2dmLChAm4dOkSgLtLO86bN49lVewyf/58NDU1qV3CcgDD9axsw+fz8c033zAzjdauXauRJRENlbq6Opw6dQoLFy5kW4rewZlVD+CCTf/H0aNHoVAoBvzoQh3cMFhP4IJNd5k5cyaMjY0RHh7OthR9gxsG6wtcsOnugst//PEHli1bxrYUvYQzqx6xYsUKjBs3DgBw5swZHDt2jGVFumXv3r2QSCSYM2cO21L0Es6sesRADzb99NNPeOGFF5RmNnH8H5xZ9YyBGmyKj49HRkYGNxXuAXABJj1kIAabXn31VVy8eBFpaWlsS9FXuACTPjLQgk3Nzc349ddf8fLLL7MtRa/hzKqnDKRg02+//YaGhgY8//zzbEvRa7hhsB4zUKbRTZ8+HRYWFjh69CjbUvQZbhiszwyEYFNeXh7OnTvHBZZ6ANez6jn9Pdj07rvvYu/evbh9+/aASVfTR7ieVd/pz8GmtrY2/Pjjj1i5ciVn1B7AmdUAuD/Y1F/u7Y4cOYLKykouCtxDuGGwgXBvsMnZ2Rk3b940+GDT1KlTIZVK+83FR8tww2BDITAwECtWrABwd/0XQw82ZWRkIDo6GqtWrWJbisHA9awGRH8KNq1fvx6///47bt26BT6f6zN6ANezGhL9JdjU3NyMffv2YeXKlZxRewH3SRkY/SHYdOjQITQ0NOCll15iW4pBwQ2DDRBDDzaNGzcOHh4eOHDgANtSDAluGGyIGHKw6eLFi0hMTMRbb73FthSDg+tZDRRDDTa98MILuHnzJpKTk9mWYmhwPauhcn+wyRB6qpKSEvz6669Yu3Yt21IMEs6sBsyKFSvw+OOPAwCioqL0Pti0c+dOSCQSLidwH+HMasDw+Xx8/fXXD83ZVFtbq7Q0JRu0tbVh9+7dWL16tdJq8Rw9hzOrgfOgYBMRYd++ffDx8UFCQoLONJ0+fRoNDQ1K+w4ePIiqqiqsXLlSZzr6HcRh8FRVVdHgwYMJAAmFQsrIyKDU1FSaOHEiASAA9Mknn+hMj5+fH1lYWNC7775LxcXFREQ0ZswYev7553WmoR+SyJm1n7Br1y7GmM7OziQQCJifAdD8+fN1psXa2poAkEAgIIFAQE8//TQBoISEBJ1p6Ickco9u+gkKhQLe3t7Izc1Ve9zJyQmFhYVa19HZ2QmhUAiFQsHsMzY2Rnt7Ox5//HH8/e9/55J49w3u0U1/4OrVqwgODu7WqABQVFSEkpISrWupqqpSMioAtLe3AwCSkpLw9NNPY+TIkTh8+LBKOY4Hw5nVgJHL5VizZg0CAwMRHx//0PKXL1/WuqaysrJuj3WZMy0tDatWrUJmZqbW9fQnOLMaKAqFAqGhodixY0ePe6ikpCQtqwJKS0sfeJzP50MsFuPMmTMYOnSo1vX0JzizGihGRkY4fvw4pk2b1uPf0UXPWlpa2u20Nz6fD2NjY5w8eRJjxozRupb+BmdWA2bQoEE4depUj+e16qJnLSsrg0AgUNnP4/HA5/Nx9OhRBAcHa11Hf4Qzq4EjEAiwfft27Nq166EZAisrK5Gfn69VPWVlZeDxeCr7+Xw+jhw5glmzZmm1/f4MZ9Z+wsqVK3Hy5ElYWVk9sJy2e9eysjJ0dHQo7ePz+fjPf/7DPbJ5RDiz9iOmT5+OxMREPPbYY92W0bZZi4qKlAJePB4Pu3btwqJFi7Ta7kCAM2s/Y8iQIYiPj8eTTz6p9ri2zVpcXMz8n8fj4ZtvvmHeXeZ4NDiz9kOsrKwQGRmJjRs3qhy7cuUKOjs7tdZ2RUUF8/9t27bhtdde01pbAw3VsB1Hv8DIyAhbtmyBi4sL1q5dy9xH1tfX4/XXX4eRkRHkcjlqa2shl8shl8vR3NyM+vp6pmx7e7vS7BlLS0tmOp5QKISZmRnMzMwglUohkUggkUhQXV0NAJgzZw6GDBmCy5cvw8HBATKZjPldjr7BvRvcD+js7ERubi4yMzORmZmJrKwsZGdnIysrC8XFxVrtSXuKQCCAo6MjvLy84O3tDR8fH3h7e8Pb2xvu7u5qI8gcSiRxZjUwOjo6kJ6ejpSUFCQnJyMlJQVXr15VmT/aU8RiMaRSKaRSKUxNTZV6T5FIBFNTU6ZsQ0MD855vV6/b0NDA9NAtLS190mBpaQl/f3/4+/sjICAAAQEB8PX15XpiZTiz6jvt7e1ITExEdHQ0YmJiEB8f3yNjWlhYwMvLC25ubnBycoKDgwMGDRqEH3/8EU1NTThz5gykUimEQqHGtFZUVKCjowM1NTW4c+cOSkpKmAkEhYWFyM/PR1ZWFpqbmx9al1QqxaRJkzBlyhQEBwfD399/oJuXM6s+cufOHURERCA8PBznzp17YEoWiUTC9Eq+vr7M0NLBwUFteYVCgY8++gibNm1iZZlFIkJhYSGysrKQlZWFGzduIDU19aGjA0tLS0yfPh1z5szB7NmzYWNjo0PVegFnVn0hJycH+/fvx++//47k5GSoOy3GxsYYM2YMJk+ejNGjR8Pf3x8eHh59ut8jIr26T+zs7ERWVhZSUlKQlJSE2NhYpKSkqJ2kYGRkhHHjxmHu3LlYvHgxnJycWFCsczizskltbS1++eUX7N27F/Hx8SoG5fP5GDt2LKZNm4bg4GBMmDDBoDLvPyp1dXWIjY1FdHQ0zp49i5SUFLWf0RNPPIFly5Zh/vz5/fnz4czKBtevX8dXX32Fn3/+WeX+zczMDDNmzEBoaChCQ0Nha2vLkkr9o6ioCBERETh+/DjOnTunEtCysLDA8uXLsXbtWnh6erKkUmskcTmYdMjp06cpJCSEeDyeUn4kgUBATz31FB08eJCam5vZlmkQ1NXV0Z49e2jq1KnE5/OVPk8+n09hYWEUFxfHtkxNwiVM0wUXLlygyZMnK32hAJCXlxd98cUXVFJSwrZEg+b27du0efNmcnZ2VvmMZ82aRVevXmVboibgzKpNbt26Rc8884xKTxocHEzHjh0jhULBtsR+RXt7Ox04cIDGjBmj0tMuW7aMioqK2Jb4KHBm1QYKhYK+/PJLMjU1VfrSPPnkk5SUlMS2vAHB+fPnacKECUqfv0QioR9++IFtaX2FM6umyc7OpkmTJil9Sfz9/enMmTNsSxuQHDlyhHx8fFSGxnfu3GFbWm/hzKpJ/vjjD7KysmK+FGZmZrRjxw5uuMsy7e3t9NFHH5FQKGTOjb29PSUmJrItrTdwZtUUu3btImNjY+bLMGXKFMrJyWFbFsc9XLt2jQICAphzJBaL6fDhw2zL6imcWTXBN998ozTMWrNmDdeb6imtra20fPly5lzxeDzas2cP27J6AmfWR+XAgQPMcz5jY2PavXs325I4esDWrVuVztvp06fZlvQwOLM+CpcuXWLug/h8Pu3fv1/jbdTX19OQIUNo9uzZGq97oLNjxw6mhzU3N6eMjAy2JT2IRC6tSx/p6OjAqlWr0NbWBuBuCpPnn39e4+0QETo7O/ViAnl/480338Tf//53AHfn6q5cuVLtBAq9ge3LhaHyySefMFflBQsWsC2Ho490dnbSlClTmHP53XffsS2pO7glH/tCW1sbHB0dUVlZCUtLS9y4cQOOjo5sy+LoI9nZ2RgxYgRaWlrg7u6OnJwcvZo++F+4JR/7wvHjx1FZWQkAWLdundaMeuzYMfB4PGbrmmVy//78/HwsXLgQUqkU1tbWCA0NRU5Ojkp9VVVV2LBhAzw9PSESieDk5ITp06djz549SrN/7i0nFAphZWWFWbNm4dy5c91qu337NhYuXAgLCwtYW1tjyZIlqKmpQX5+PubMmQMLCwvY29vj1VdfRX19vYq2iooKrFmzBm5ubhAKhbCxsUFYWBhSU1O18Mkq4+XlheXLlwMA8vLyEBsbq/U2+wTbfbshMm/ePCbsr4tnqXPnziUAKjNyuvbPnTuXEhISqKGhgaKiokgsFtOYMWOUypaUlJC7uzvZ2dlReHg41dXVUWlpKX300UcEgL788kulcjKZjMLDw6m2tpYyMzMpLCyMeDyeyjCxS0NYWBhdvnyZGhoaaN++fcybQnPnzqWUlBSqr6+nnTt3EgBav369Uh3FxcXk6upKMpmMTpw4QfX19ZSenk7BwcFkYmKikxXTL126xAyFV65cqfX2+gAXDe4L9vb2BIBGjRqlk/YeZtbwffAFsgAACD1JREFU8HCl/c8++ywBoIqKCmZf17PFQ4cOqdQ/c+ZMxqxd5Q4cOKBUpqWlhRwcHEgsFlNpaamKhhMnTiiVHzZsGAGg6Ohopf3u7u7k4+OjtG/ZsmUEgH7++Wel/SUlJSQSiSgwMFDt56JpHB0ddXpeewkXDe4tHR0dzArivr6+LKu5y/3LJzo7OwNQzo5/9OhRAFC7MFRkZCTWrVunVG727NlKZUQiEaZNm4bm5macPn1apY7Ro0cr/dyVA+r+/Y6Ojkq6gLtDaj6fj9DQUKX9dnZ2GDZsGK5cuYKioiKVNjVN1/ksKCjQelt9gUvy3UvuzU4gFotZVPJ/SCQSpZ+7MhZ2Pe5pbW1FbW0tTExMYGFh0W09Dysnk8kAqF8w2dLSUulnPp8PIyMjpVSmwN38Sfc+hupqU93fcS/Z2dlaz7XUpbUn2RfZgDNrLzEzM4NAIEBHRwcTZNJ3RCIRJBIJamtrUV9f361hH1aurKwMwN0eT5PapFIpGhoa0NzcrHZtV13RtfTHgy4abMINg3sJj8eDt7c3ACA5OdlgXlaYP38+AODkyZMqx/z9/bF+/XqlcidOnFAq09rairNnz0IsFiMkJESj2sLCwtDR0YH4+HiVY1u3boWLi4vKMpKapqWlBenp6QD05/bmfjiz9oFp06YBuJvf9/z58+yK6SGffvop3N3dsX79epw4cQL19fUoKirC66+/jpKSEsasXeXWrVuHiIgI1NfXIysrC4sXL0ZJSQm2b9/ODIc1qc3T0xMvv/wyIiMjUVtbi+rqauzatQsffvghtm3bpvUe9/fff2fyFnedX72D7RCXIXL58mUmzP/cc89prZ2jR4+q5BR64YUX6MKFCyr7N23aRESksv/ed4orKytp3bp15O7uTsbGxmRvb0+LFi2irKwspXbvLyeRSCgkJITOnj3LlOlOQ1JSksr+Tz/9lGJjY1X2v//++0x9VVVVtGHDBvLw8CBjY2OysbGhGTNmUFRUlNY+33uZOnUq8453fn6+TtrsJdwbTH1l1KhRuHr1Kng8Hv744w888cQTbEvi6COHDx/GwoULAdxdkDoqKoplRWrh8gb3lT///BPTp08HEcHLywtXrlx5YKSVQz+pqKjAiBEjUFpaCoFAgKSkJIwaNYptWergXjfsK11Z4IG7jxXmzp2L1tZWllVx9IampibMnTuXeRS1YcMGfTUqAG591keiuroaEyZMQGZmJgBg8eLF2Lt3L6uPHzh6RktLC8LCwhAZGQng7ssb0dHRKs+F9QiuZ30UBg0ahD/++IN5Y2j//v2YPXs25HI5y8o4HkRlZSVCQkIYo3p6eiIiIkKfjQqAe3TzyDg5OSE8PByDBw8GAJw5cwYTJ05EdnY2y8o41JGSkoLRo0cjJiYGwN1XM8+cOaPxx1HagDOrBhg5ciQuXbqEoUOHAgBu3LiBkSNHYuvWrWqXLOTQPR0dHdi6dSsef/xx3L59G8DdiH58fDw8PDxYVtdDWHxu1O+Qy+UUEhKi9Cxx0qRJlJaWxra0Ac2lS5fI399f6bwsWrSImpqa2JbWG7gpcpqms7OTdu3aRebm5kprrSxYsIDy8vLYljegyM/PpyVLliitNSSRSGjXrl1sS+sLnFm1xa1btygoKEjpai4Wi+kvf/kLFRQUsC2vX5OVlUWrVq1SSroOgEJDQw1x2YwuOLNqm+PHj5OXl5fKqmahoaE6yYAwkLh8+TItWbKEjIyMlD7voUOHGlLm/e7gzKoL2tra6JtvvmEyTNy7jR07lnbs2EGVlZVsyzRIiouLadu2bTRixAiVz9bd3Z327dvXX1ZH4MyqS1pbW2nPnj00cuRIlS+WUCikefPm0aFDh6i2tpZtqXpNZWUl7du3j2bOnKnSiwKgxx9/nA4fPkwdHR1sS9Uk3Iv8bPHHH3/g22+/xYkTJ5hE4V0IhUIEBwfj6aefRmhoKNzc3NgRqUdkZ2fj+PHjCA8PR1xcnMojMbFYjPnz5+ONN97AhAkTWFKpVbgX+dlGLpfj8OHD2LdvHxISEtRmhLe3t8ekSZMwffp0TJw4EcOGDWNBqW7Jzc1FXFwc4uPjERUVhby8PLXlAgMDsWTJErz44ouwtrbWsUqdwplVn8jOzsbvv/+OiIgItb1HFy4uLhg9ejQCAgKYzRDewOmOoqIiJCcnM9vly5eZpHT3c++o4+mnn4aLi4uO1bIGZ1Z9paqqCpGRkThz5gyio6MfmnHPwcEBvr6+8Pb2hpeXF7y9veHt7Q1XV1eIRCIdqe6e5uZm5OXlISsrC1lZWcjOzkZ2djauX7/+0FxWQ4YMQVBQEGbOnImQkBCV5GwDBM6shkJ+fj6io6MRHR2NixcvIisrq8evMtra2sLe3h5OTk5wcHCAo6MjrKysIJVKIZFIIJVKIZVKYWFhAaFQCDMzMwB3MxTemzxMLpczw/T6+np0dHSgtrYWtbW1kMvlkMvlTEqWoqIilJSUoLCwECUlJaiqquqRVoFAgKFDh2L8+PEICgrClClTmLSmAxzOrIZKY2MjUlNTmaHjtWvXkJ2drXZpCn1FIpHAy8sLo0aNQkBAAAIDAzFixIj/384d8jAMQkEAviYIEG3BkLShuj9+v2mq6aupxCE6NRLUVNfR3adQ5CEuOcOD1vrq0X4Rw3o3IlJUzWVZsK4rRAQiUuw9PpsxBiEEDMOAaZoQQsgVfZ5neO+/NssNMKz/Zt93bNtW1Nb3OcaIlFLe8nccR/6b2zQNrLX5nrZtoZRC13VFle77Hs45jOMI59wlb7wphpWoEtwUQVQLhpWoEgrA4+ohiOij5wudocsenekDSwAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#Step 1: Model\n", "model=CausalModel(\n", " data = df,\n", " treatment='education',\n", " outcome='income',\n", " common_causes=['U'],\n", " instruments=['voucher']\n", " )\n", "model.view_model()\n", "from IPython.display import Image, display\n", "display(Image(filename=\"causal_model.png\"))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Estimand type: nonparametric-ate\n", "\n", "### Estimand : 1\n", "Estimand name: backdoor\n", "Estimand expression:\n", " d \n", "────────────(Expectation(income))\n", "d[education] \n", "Estimand assumption 1, Unconfoundedness: If U→{education} and U→income then P(income|education,,U) = P(income|education,)\n", "\n", "### Estimand : 2\n", "Estimand name: iv\n", "Estimand expression:\n", "Expectation(Derivative(income, [voucher])*Derivative([education], [voucher])**\n", "(-1))\n", "Estimand assumption 1, As-if-random: If U→→income then ¬(U →→{voucher})\n", "Estimand assumption 2, Exclusion: If we remove {voucher}→{education}, then ¬({voucher}→income)\n", "\n", "### Estimand : 3\n", "Estimand name: frontdoor\n", "No such variable found!\n", "\n" ] } ], "source": [ "# Step 2: Identify\n", "identified_estimand = model.identify_effect(proceed_when_unidentifiable=True)\n", "print(identified_estimand)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "*** Causal Estimate ***\n", "\n", "## Identified estimand\n", "Estimand type: nonparametric-ate\n", "\n", "### Estimand : 1\n", "Estimand name: iv\n", "Estimand expression:\n", "Expectation(Derivative(income, [voucher])*Derivative([education], [voucher])**\n", "(-1))\n", "Estimand assumption 1, As-if-random: If U→→income then ¬(U →→{voucher})\n", "Estimand assumption 2, Exclusion: If we remove {voucher}→{education}, then ¬({voucher}→income)\n", "\n", "## Realized estimand\n", "Realized estimand: Wald Estimator\n", "Realized estimand type: nonparametric-ate\n", "Estimand expression:\n", " \n", "Expectation(Derivative(income, voucher))⋅Expectation(Derivative(education, vou\n", "\n", " -1\n", "cher)) \n", "Estimand assumption 1, As-if-random: If U→→income then ¬(U →→{voucher})\n", "Estimand assumption 2, Exclusion: If we remove {voucher}→{education}, then ¬({voucher}→income)\n", "Estimand assumption 3, treatment_effect_homogeneity: Each unit's treatment ['education'] is affected in the same way by common causes of ['education'] and income\n", "Estimand assumption 4, outcome_effect_homogeneity: Each unit's outcome income is affected in the same way by common causes of ['education'] and income\n", "\n", "Target units: ate\n", "\n", "## Estimate\n", "Mean value: 3.984270555650861\n", "p-value: [0, 0.001]\n", "\n" ] } ], "source": [ "# Step 3: Estimate\n", "#Choose the second estimand: using IV\n", "estimate = model.estimate_effect(identified_estimand,\n", " method_name=\"iv.instrumental_variable\", test_significance=True)\n", "\n", "print(estimate)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have an estimate, indicating that increasing `education` by one unit increases `income` by 4 points. \n", "\n", "Next we check the robustness of the estimate using a Placebo refutation test. In this test, the treatment is replaced by an independent random variable (while preserving the correlation with the instrument), so that the true causal effect should be zero. We check if our estimator also provides the correct answer of zero. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Refute: Use a Placebo Treatment\n", "Estimated effect:3.984270555650861\n", "New effect:-0.004644847508369445\n", "p value:0.43999999999999995\n", "\n" ] } ], "source": [ "# Step 4: Refute\n", "ref = model.refute_estimate(identified_estimand, estimate, method_name=\"placebo_treatment_refuter\", placebo_type=\"permute\") # only permute placebo_type works with IV estimate\n", "print(ref)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The refutation gives confidence that the estimate is not capturing any noise in the data.\n", "\n", "Since this is simulated data, we also know the true causal effect is `4` (see the `income_education` parameter of the data-generating process above)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we show the same estimation by another method to verify the result from DoWhy." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
IV2SLS Regression Results
Dep. Variable: income R-squared: 0.887
Model: IV2SLS Adj. R-squared: 0.887
Method: Two Stage F-statistic: 1370.
Least Squares Prob (F-statistic): 1.69e-189
Date: Tue, 02 Mar 2021
Time: 21:15:13
No. Observations: 1000
Df Residuals: 998
Df Model: 1
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [0.025 0.975]
Intercept 9.8119 0.986 9.953 0.000 7.877 11.746
education 3.9843 0.108 37.020 0.000 3.773 4.195
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 2.309 Durbin-Watson: 1.949
Prob(Omnibus): 0.315 Jarque-Bera (JB): 2.075
Skew: -0.014 Prob(JB): 0.354
Kurtosis: 2.779 Cond. No. 25.2
" ], "text/plain": [ "\n", "\"\"\"\n", " IV2SLS Regression Results \n", "==============================================================================\n", "Dep. Variable: income R-squared: 0.887\n", "Model: IV2SLS Adj. R-squared: 0.887\n", "Method: Two Stage F-statistic: 1370.\n", " Least Squares Prob (F-statistic): 1.69e-189\n", "Date: Tue, 02 Mar 2021 \n", "Time: 21:15:13 \n", "No. Observations: 1000 \n", "Df Residuals: 998 \n", "Df Model: 1 \n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept 9.8119 0.986 9.953 0.000 7.877 11.746\n", "education 3.9843 0.108 37.020 0.000 3.773 4.195\n", "==============================================================================\n", "Omnibus: 2.309 Durbin-Watson: 1.949\n", "Prob(Omnibus): 0.315 Jarque-Bera (JB): 2.075\n", "Skew: -0.014 Prob(JB): 0.354\n", "Kurtosis: 2.779 Cond. No. 25.2\n", "==============================================================================\n", "\"\"\"" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "income_vec, endog = ps.dmatrices(\"income ~ education\", data=df)\n", "exog = ps.dmatrix(\"voucher\", data=df)\n", "\n", "m = IV2SLS(income_vec, endog, exog).fit()\n", "m.summary()" ] } ], "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.8.5" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": false, "sideBar": true, "skip_h1_title": true, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }