{ "cells": [ { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# RESEARCH IN PYTHON: INSTRUMENTAL VARIABLES ESTIMATION\n", "# by J. NATHAN MATIAS March 18, 2015\n", "\n", "# THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\n", "# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\n", "# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\n", "# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\n", "# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\n", "# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN\n", "# THE SOFTWARE." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Instrumental Variables Estimation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This section is taken from [Chapter 10](http://www.ats.ucla.edu/stat/stata/examples/methods_matter/chapter10/default.htm) of [Methods Matter](http://www.ats.ucla.edu/stat/examples/methods_matter/) by Richard Murnane and John Willett. The descriptions are taken from Wikipedia, for copyright reasons.\n", "\n", "In statistics, econometrics, epidemiology and related disciplines, the method of instrumental variables (IV) is used to estimate causal relationships when controlled experiments are not feasible or when a treatment is not successfully delivered to every unit in a randomized experiment.\n", "\n", "In linear models, there are two main requirements for using an IV:\n", "\n", "* The instrument *must* be correlated with the *endogenous explanatory variables*, conditional on the other covariates.\n", "* The instrument *cannot* be correlated with the *error term* in the explanatory equation (conditional on the other covariates), that is, the instrument cannot suffer from the same problem as the original predicting variable.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Example: Predicting Civic Engagement from College Attainment" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Can we use college attainment (COLLEGE) to predict the probability of civic engagement (REGISTER)? College attainment is not randomized, and the arrow of causality may move in the opposite direction, so all we can do with standard regression is to establish a correlation.\n", "\n", "In this example, we use an _Instrumental Variable_ of distance between the student's school and a community college (DISTANCE), to estimate a causal relationship. This is possible only if this variable is related to college attainment and NOT related to the residuals of regressing COLLEGE on REGISTER. \n", "\n", "The python code listed here is roughly parallel to [the code listed in the textbook example](http://www.ats.ucla.edu/stat/stata/examples/methods_matter/chapter10/default.htm) for Methods Matter Chapter 10. If you're curious about how to do a similar example in R, check out \"[A Simple Instrumental Variables Problem](http://www.r-bloggers.com/a-simple-instrumental-variables-problem/)\" by Adam Hyland in R-Bloggers or Ani Katchova's \"[Instrumental Variables in R](https://www.youtube.com/watch?v=OwM3BgWEgUg) video on YouTube." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# THINGS TO IMPORT\n", "# This is a baseline set of libraries I import by default if I'm rushed for time.\n", "\n", "import codecs # load UTF-8 Content\n", "import json # load JSON files\n", "import pandas as pd # Pandas handles dataframes\n", "import numpy as np # Numpy handles lots of basic maths operations\n", "import matplotlib.pyplot as plt # Matplotlib for plotting\n", "import seaborn as sns # Seaborn for beautiful plots\n", "from dateutil import * # I prefer dateutil for parsing dates\n", "import math # transformations\n", "import statsmodels.formula.api as smf # for doing statistical regression\n", "import statsmodels.api as sm # access to the wider statsmodels library, including R datasets\n", "from collections import Counter # Counter is useful for grouping and counting\n", "import scipy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Acquire Dee Dataset from Methods Matter" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import urllib2\n", "import os.path\n", "if(os.path.isfile(\"dee.dta\")!=True):\n", " response = urllib2.urlopen(\"http://www.ats.ucla.edu/stat/stata/examples/methods_matter/chapter10/dee.dta\")\n", " if(response.getcode()==200):\n", " f = open(\"dee.dta\",\"w\")\n", " f.write(response.read())\n", " f.close()\n", "dee_df = pd.read_stata(\"dee.dta\")" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "# Summary Statistics" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "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", " \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", " \n", " \n", "
registercollegedistance
count 9227.000000 9227.000000 9227.000000
mean 0.670857 0.547090 9.735992
std 0.469927 0.497805 8.702286
min 0.000000 0.000000 0.000000
25% 0.000000 0.000000 3.000000
50% 1.000000 1.000000 7.000000
75% 1.000000 1.000000 15.000001
max 1.000000 1.000000 35.000000
\n", "
" ], "text/plain": [ " register college distance\n", "count 9227.000000 9227.000000 9227.000000\n", "mean 0.670857 0.547090 9.735992\n", "std 0.469927 0.497805 8.702286\n", "min 0.000000 0.000000 0.000000\n", "25% 0.000000 0.000000 3.000000\n", "50% 1.000000 1.000000 7.000000\n", "75% 1.000000 1.000000 15.000001\n", "max 1.000000 1.000000 35.000000" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dee_df[['register','college', 'distance']].describe()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "# Cross-Tabulation" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "college 0 1\n", "register \n", "0 1780 1257\n", "1 2399 3791\n", "chi2: 323\n", "p: 0.000\n", "df: 1.000\n", "expected:\n", "[[ 1375.48748239 1661.51251761]\n", " [ 2803.51251761 3386.48748239]]\n" ] } ], "source": [ "print pd.crosstab(dee_df.register, dee_df.college)\n", "chi2 = scipy.stats.chi2_contingency(pd.crosstab(dee_df.register, dee_df.college))\n", "print \"chi2: %(c)d\" % {\"c\":chi2[0]}\n", "print \"p: %(p)0.03f\" % {\"p\":chi2[1]}\n", "print \"df: %(df)0.03f\" % {\"df\":chi2[2]}\n", "print \"expected:\"\n", "print chi2[3]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Correlation Matrix" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": [ "iVBORw0KGgoAAAANSUhEUgAAASYAAADtCAYAAAASnNgZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n", "AAALEgAACxIB0t1+/AAAGdRJREFUeJzt3XmYXFWZx/FvVaeBsEQBQVlkQAZfGRhxQRQDEURFI2Aw\n", "ICMIxgUUFAEdBRlZBhBlccVH8WERlOgoAwrIqjCYCBGVKILgGyIgBlAQIayBdFfNH+cUqRTddatr\n", "u6du/z7Pc5+mlpM6VHe99Z7lvrdUrVYREUlJOe8OiIg0UmASkeQoMIlIchSYRCQ5CkwikhwFJhFJ\n", "zpSMx5cBq/ajIyKyklLeHchTKWMfU5VJ/gZNgN6r1uh9kkwayolIchSYRCQ5CkwikhwFJhFJjgKT\n", "iCRHgUlEkqPAJCLJUWASkeQoMIlIchSYRCQ5CkwikhwFJhFJjgKTiCRHgUlEkqPAJCLJUWASkeQo\n", "MIlIchSYRCQ5CkwikpxCBiYzK5mZ6kqLDKisq6QMHDMbcvfRvPshIu0rVMZkZlPcfdTMymZ2hZnt\n", "nHefRGTiChWY3H3EzMrAT4AngN/m3CURaUMhAlMMRjUzgbcDL3H3x+PjQ7l0TETaMvAXvIzDt5E4\n", "2f0CwrzZHsDpwKXuPqdPXUn+vUqE3ifJNNCBycxK7l6NGdE1wCrAEuDnhL6fDfzA3ffrQ3eSfq8S\n", "ovdJMg3sqlzD6ttZwN+Bi4EfAQ8CXwMOAf6STw9F8nX5sDXNOgDeudyT/JIYyMBkZlsTAtFDMVta\n", "ExgBTgS+BbwQeJm7nxmfX3L3zF+SSJEMTR3cKeSBDEzA9sBxZgawN/Br4FTgVuB24IuEYRwACkoy\n", "GZWGOwtMcVHpm8ArgWeAD7v7n+sefy9wGCEpuJUwQik1a9OqgQqpdatv84ENgZcAjxGGbScAawH7\n", "E96M+dr9LZPZ0NRy5pFhFrCKu78ROAr4Uu0BM5tKGKHs5O47EBaedottVh2rzUQMTGCKq2+VePM+\n", "4AxgAXA58HLgG8AOwEx3/6GCkkx25SmlzCPDdOAqAHe/Cdi27rFlwPbuvizenhLvmw5cOU6blg3E\n", "UM7MVnH3Z+N80ncJwegYYJ14+1bgaeCV7n5/rV1Rh3BZKXZ8zurAz4APurub2SqE4e2/AsuBT7j7\n", "Lf3tef+Y2e6Ev5ER4Fx3P7vh8RcB3wdWA+4HPuDuT5vZbOBIwurhXHf/evy7O4vwBVgFPuruf+zf\n", "/017ujDHNI0wIqkZNbOyu1fiZ+shADM7FFjD3X9mZu8Zr81EXnggMqYYlMrAQkIw3RR4D7A5MAc4\n", "Hti7/sNZ1KAUjZtiA5jZtsA8YDPCBwngQOCp2OZA4Nz+dbe/zGwY+DLwVuBNwEFmtn7D044FLnD3\n", "GcDvgI/Ev7EvALsQ5jEPMbN1gd2BShyyfA74fH/+TzpTGiplHhkeI0yP1KwUYOKpX6cT3q/ZrbRp\n", "VdKByczOj99gEDZNbgRcCHycMOn2deB+dz/B3a+YRFUFmqXYEPZzzQK87r5/q2uzCNjIzKb1vqu5\n", "2BJY7O5L3X058EtgRsNznnsPCUOPt8QP0JbxjIH1gCHgGXf/CfCR+NxNgUd63P+uGBouZx4ZbiCc\n", "SYGZvQH4Q8Pj3wZWBfasG9JltWlJskM5M9uMsAfpajPbhfA//DhhcvtU4G7C7u61gQeg8FlSvXFT\n", "bAB3vxEgrlrW/J4wOfmT+AezHrBGw79TFNOApXW3HydMzo73nCdqj8eTwN9NmLP8KfBU3f3nAXsC\n", "e/Ws5100tErHZ2L9GHirmd0Qb38grsStSTgP9YOEzPy6+Lf21bHatPPCSQYmMxt297vN7HjgM4Sh\n", "2hx338zMDgdeQ0inP+TuD+TX09y0ky6fC2xpZvMJQX4R8M8e9S8XZnYiYQHklcBNdQ+txfOznMcI\n", "wemh+PijtQfc/WIz+zFwHnBA/Im7zzGzI4GbzGxLd3+6N/8n3VEqdzZ4iF/0Bzfcvajuv8eLfI1t\n", "Jiy5wBR3dC+P4/0dCZPbrwJOikOPSwhD0PPc/bra0G0SZUsQAsvuwIUTSJe3A65z90/GOajt3P2Z\n", "Xnay39z9GAgruMDtZrY28CRhGHdaw9NrQ47zgXcA88xsLUKW9NY4r/kkIRvdH9jY3b9AWGSpxCNp\n", "Q6skPVPTVFKBqe6E3DJhj8RnCXNLh7Jir9KIu385Pn8yBiVokmK7+1njtHHgh2Z2NGFZ98A+9DMX\n", "8W/ok8DVhC+xc9z9ATNbBzjL3WcDJwHnm9mBhKxp37gqdwEhSC0HbgEuIKzcnWdmvwCGgcMGIaiX\n", "yoMbmJI7iTcGpV8BdwBvALYA3k3Y3X0G8E13v7affWqRTk5tjd6nPvn923bM/MJ+1TXzk/xdJJUx\n", "RfsBGwA7A+sSvrEuJsz8z4YQvNpZghSZTDqdY8pT7rlenA+oN0o41eTV7n4voZzJAuDHZjYdQEFJ\n", "JFt5ylDmkapcM6Y40V2bU9qNMPdxIfBR4Cozu5awynIQYWVu3bz6KjJoytkbKJOVW2Cq1VOKQek2\n", "4B7g1YS5pT2AIwjLuH8glDHZlHDqiTQws40J79kiwjL4noQh8BJCMD+GsJDwyCBM2vZDs/fM3Zc2\n", "azsoUs6IsuQylKu/mgnhQ3MPcDhhYnQdYLq7H0fYsbtNfM7O7n53Hv0dAFsSVtmOBnYC9iGsYK5H\n", "qH8+k7AK1bjJcDJr9p4VQqlcyjxS1feMKU5c14ZvZxGypJcSNsSdSZhj2svMfh43uv0FuM/d/9bv\n", "vg6QhYQTUtclnGbxCLCFuy82s00Ip/Bs5e4P5tjH1Iz7nuXaqy4qT8l9CrltuWwXiPuPLgNeQZg/\n", "+jRhLulywsbB97n7j7v9uj2mZfDW6H3qkzv3m5m5XWCLuVck+bvIa45pDcLw7e2EuaM5hHPgXkzY\n", "FnDNJN48KdIVLZykm6xcApO7PxF3ID9LGM4d4u6n1x5XUBLpXKc7v9up+xXvW8iKE6TvcvcPTfS1\n", "c1uVc/fHzOw4YCrwLTO7Gviru48qIIl0rgtzTM/V/TKz1xPqfs2qPRjPuTyTUOa6Gu9bDcDdd+7k\n", "hXPN9WLdm88QJmbvqbsck4h0qFQuZx4Z2qn7tQ2wupldbWbXxoA2YbkPQt39cXe/A1YM4USkc13Y\n", "+T1m3a/aDXe/0d2XNLR5EjjN3XclbJSeW9+m5b5PtEEvaQgn0kWlUvbRXDt1vxYBcwHc/U7gYcK5\n", "rxOSVGASke7pQsbUTpncDxBr0JvZhoSsa8LFHFOsLiAiXdCFU1Laqft1DvAdM5tXa9POSffJ1WMa\n", "YHqvWqP3qU/+fuT+mVMjLz7le0n+LpQxiRRUaYBP4lVgEimoUvbkdrIUmEQKapDLnigwiRTVAF+M\n", "QIFJpKDKQ8qYRCQ1CReCy6LAJFJQmmMSkfSUNMckIonRPiYRSY72MYlIcpQxiUh6Opxjaqe0bitt\n", "WjG4s2Mi0lRpaCjzyPBcaV3gKGI5k5pYWncesBmxtG5Wm1YpMIkUVbmUfTTXTmndrDatdb2dRiKS\n", "vi5kTO2U1m3aplWaYxIpqs5PSWmntG47bZ5HGZNIQXXhKintlNZtp83ztJIx6QIBrdN71Rq9T9k6\n", "34TUecbUTmnd57Vp54UzS+v+Zsb27fy7k8rr5i0ooZKxrdL71CfLfnBK5hfAau89MsnfheaYRIpK\n", "ZU9EJDk6JUVEkjM0uB/vwe25iDSnQnEikpyy5phEJDW6GIGIJEcZk4gkRxmTiCRHGZOIJEf7mEQk\n", "NVXt/BaR5OjyTSKSmmqHc0xZ9bvNbHfgGGAEONfdz473LwSWxqfd5e4fmuhrKzCJFFXnk9/P1e82\n", "s9cT6nfPAjCzYeDLhNK5TwE3mNklwOMA7r5zJy88uLmeiDRXLmcfzTWr370lsNjdl7r7cuCXwJuA\n", "bYDVzexqM7s2BrSJd72dRiKSvmp5KPPI0Kx+9zRWDNcgZEovAJ4ETnP3XYGPAnPbqfmtwCRSUNVS\n", "OfPI0Kx+99KGx9YCHgEWAXMB3P1O4GFgg4n2XYFJpKC6kDE1q9/9J2ALM1vbzFYBZgALCKV0vxTb\n", "bEjIrB6YaN81+S1SVJ1vsGxa89vMPglcTUhwznH3B8zsHOA7Zjav1qadq6So5ncXqOb3hOh96pPH\n", "br46s+b3tNfumuTvQhmTSFFpg6WIpKZSGtxTUgY3pErbzGyo4XaS6bx0qFTOPhKljGmSMbMp7j4S\n", "95ZsAyxz9zvy7pd0X2WAy56kGzKlJ2JQGgJuBY4H/mhmx5rZC/PtmXRdqZR9JEqBaZJo2H37KcI+\n", "lE8Bo8BOwMty6Jb0UKU0lHmkSkO5ScDMhtx9NM4llQm7dt8GvB44DRgG9jWz24Dl7p65zCzpa2Fn\n", "d7IUmAquYU7pTOAh4OfATcCbgWWEzGmWuz+bX0+l26oJZ0RZFJgKzMzKdXNKpwHvImRMTwNnE4Zz\n", "WwF7u/uVZlZStlQc1YTnkLIoMBWYu1fi8O0yQjB6G2HC++PA14BPu/vToC0DRZTyHFKWwR2Eyrga\n", "JrpXBdYH/uzutwAfJJSm+DDwYTMbjplVVdlSsVTKQ5lHqpQxFUxtTin+9wsIdXLmAbPM7O/A34Df\n", "AhsB73H3M3LrbB+NVwZ2jOcdDrzY3T9bd99XgD+5+7f70tkuqXZ4SmI7pXWz2rRKgalA4upbbU7p\n", "EmAqoVTFg4QSFofG+w4g/OHsZWZrufvjefW5H8YpA3upuz9Y95zVgHOA1wH/G+9bD/gusAUwcJtQ\n", "uzCUm2hp3UuBHYBVx2ozERrKFYSZrRO3BAwB/0Mo0HURcBhhOHccsD1wPrBPvH1I0YNSNFYZ2BkN\n", "z1kNOA/4PCuqH6xBeJ++xwBWROhCobiJltadEdtcOU6blikwFYCZnQKcbWYvdfdRYHXgUcIk97nA\n", "psAmwD+B+YRvuDe4+8359LjvxisD+xx3f9Tdf9Zw3z3u/us+9K8nurDBsp3Sus3atEyBqRiWEoZm\n", "x5nZOoQ5pEMJWdOvCdsERtz9GXe/DDjM3W/Lrbd9YmYnmtn/EYa10+oeqpWBLbQqpcwjw0RL6z6a\n", "0aZlmmMaYDFT2gTYl/CNdQRwKmEj5erAfxAC1MHu/qvaPqWYVRWeux8DYUEAuN3M1iasSM4g7Osq\n", "tErneccNwO7Ahc1K67Lye1pt0qZlCkwDKtZZHiaM4b9BGLaV488K8F/AicDq7v63ybxPKS4IjFUG\n", "dh3gLHef3dBkrG0TA7eVotp5YGqntO7z2rTzwiqt2wX9Lq1bd+7b6sCngfcT/kAOIQSmY4HLgY+4\n", "+zOxTSq7ulVat0/8z3/N/H3b5i9N8nehjGnA1J37NuTuT5nZycBywjdTLXOaAvy+FpQAEglK0ked\n", "7mPKkwLTAGnYp3ReHN8/DJxCyEQ+CQy7+0Hx+alkSZKDSnVw17YUmAZIXemSGwmrH9cBHwPWBfYk\n", "bJ68ru75CkqT2CBnTIMbUieRuKpUsxlh49+n3P0LhFW3rQkbBI939+sn80S3rFChnHmkShlT4hrq\n", "KR1CCExbAG82syXAeoStAlOBJ0CZUjNmtjGwB+FS1o8RMs0LgCWEzPMY4LPAI/VzdIOoWh3c7ycF\n", "psTVBaXfxbsuA75NOE9pf8CAD9af9yVNbQkcSNhgeRXh9JyXA0cSysLMJJQbPppwjuHAGk04I8qi\n", "wDQY5hA+LNPd/Wkz24/wAfoc8JC7/7Y2fFO2lGkh8H1CdnQVIUBt4e6LzWwTwqrmVkUI9MqYpNce\n", "JmRGOwNXAG8EfgZcUzchrqDUAnd/mJV3ff+h7rHawsGFfe1Uj2hVTnrtRuB64HwzuwbYG9izdmqJ\n", "ApKMZZAD0+D2fBJx94eAgwjf9I8Cu7n75Vp9k2YqLRypUsY0INz9PsIJuoBqdEu2Qc6YFJgGlIZv\n", "kkWT3yKSnNEeBCYzm0rY91XbP/d+d/9Hw3MOJEw9jAAn1U07LCHsHwNY4O5Hj/c6CkwiBdWjU1IO\n", "Bm5x9xPMbB/ClpXDaw+a2UsIZyO8lrDp95dxweZfgJvdfY9WXkSBSaSgKpWeBKbphJPGIewDO6bh\n", "8e2AG2Id8OVmthjYBtgc2MjMriNc4/AId1/EOBSYRAqq0vnlmz5EXTYU/Z0VNb2fVzudUFZ3rFrg\n", "9wMnu/tFZjadMBzcbrzXVmASKahOMyZ3P4dwSavnmNlFrKjpXavzXa+x5netvvodhDkn3P0GM9uw\n", "2WsP7nqiiDTVhYsRjOUGwulQAO8gXEy13q+BHc1s1XjB1S2BPxKqqh4OYGbbAPc2exFlTCIFNdqb\n", "OaZvEc5AmE+40u6+AGZ2BOE6c5eZ2dcJlwkrA0e7+zNm9kXgAjObScic5jR7EdX87oJ+1/wecHqf\n", "+uSKhcsz97rNfM1wkr8LZUwiBdWLfUz9osAkUlDa+S0iyRlN+SzdDApMIgU1yBcjUGASKShlTCKS\n", "nB5tF+gLBSaRgmq+EyhtCkwiBaWhnIgkR9sFRCQ5yphEJDmaYxKR5PQiY2qltG583nqESgRbu/uz\n", "rbarUdkTkYKqVrOPNtRK684AvksorbsSM9sVuAZYfyLt6mVmTK+bt2ACfZ60qg0/pTm9T9k6nrke\n", "He1GN54nq7QuhMvZ7wLcPMF2z8kMTL+6Y2nWUya9ebevyWdmD3HqRb35SyiSz8weApU96YtO55ja\n", "LK2Lu/88tq+/exorSu6O2a6e5phECmq00kpkGv87os3SuuN5jBCcWmqnOSaRgurRHFNWad2utFPG\n", "JFJQPZpjyiytW/fcala78SgwiRRUSyO5CXL3p4H3jHH/V8a472VZ7cajwCRSUJXRzuaY8qTAJFJQ\n", "vciY+kWBSaSgRlvKmNKkwCRSUBWdxCsiqakM8Fm8CkwiBdXa5HeaFJhECkpDORFJTmunpKRJgUmk\n", "oKoKTCKSGm0XEJHkVLUqJyKpUcYkIsnpRcbUQc3vErAEWBSfssDdjx7vdRSYRAqqRxlTrXb3CWa2\n", "D6F290pVLmPN7y+ycs3vzYGb3X2PVl5EheJECqpaqWYebZhOqNlN/PmWMZ5Tq/n9SN19rwU2MrPr\n", "zOxyM3t5sxdRxiRSUKMdXr+pyzW/7wdOdveLzGw6YTi43XivrcAk0gIzK7t7ZbzbKep0H1OXa37/\n", "FhiJ/+4NZrZhsydrKCeSwcymuHvFzNY0s40AUg9KEOaYso42tFvz+1hi9mVm2wD3NnuyMiaRJsys\n", "5O4jZvYy4FLgITNbk/Ahu9Hdk12T79E+pnZrfn8RuMDMZhIypznNXkSBSaQJd6+a2WqED+Rc4EHg\n", "LGAWcBNxeJKiTueYxtJBze+lwO6tvo6GciJjMLOhupsjhPmU7YETgeOAFwNvzaFrLauMVDKPVClj\n", "EmkQ55RGzGx94DWEjYEXA8cDfyBU8H87IUAlS4XiRAqibk5pE2AB4bLWpwOXAMPAbGAnYFd3vzu3\n", "jrag0oOhXL9oKCdSJ84pDQFfAOYD27v7ucC7CEveOwDvcvff5djNllQq1cwjVQpMIqw8p+TutWvY\n", "loGpcfL7w8Cq7r7M3R/Po48TVRmtZB6p0lBOJj0zG3L30Tin9BbgWuAywipchbDTeX3gj/n1cuKq\n", "A1xbV4FJJrW6oLQBYbPgBoRg9ArgAGBvYHVg59TnlBr1YrtAv2goJ5NaDEovIqywXQjsRlh5uwX4\n", "nbvPAvZy91ty7GZbenQSb18oMMmkZGbl+LNECEoHApsAvyDsSn4AuCbOLy3PqZsdGR0dzTxSpaGc\n", "TDp1+5TWJHwGjow/dwWOAL4KvBtY7u7L8utpZ1LOiLIoMMmkEueURsxsM8KmyVHgO8AnCKed/BdQ\n", "cfev5tjNrkh51S2LApNMKnWrbwuAi4DNCEO5p4FD4s+f5tfD7qn0YFWuldK68YTefeLNK2K1y5ZK\n", "8tZojkkmBTN7U/xwAKwNXEeoFlAFFgNnA0e5+6HuvjinbnZVZXQ082hDrbTuDOC7hNK6z4lVGPYl\n", "bEx9A/A2M/v3rHaNFJik8OIE987AXDO7EliDsEfpYMKpJnMJFQN+lFsneyCn0rr3Ek7Xqf3jw8Cy\n", "FtqtREM5Kbx4mskLCaVKbgXucveFZnYAcBJhn9K27v6nPPvZbZ2uurVTWtfdR4B/xi+D04CF7n6n\n", "mU0jnHc4ZrtGCkxSWLXVt3jzF4T5jVcDX4rzIFcQMqbfuPudOXWzZ/IqrRu3WJxLCESHxLsfA6Y1\n", "a1dPgUkKycyG3X15rBJwAPBLd9/PzI4ifFgWAg7s5+6t1q0eKJWRnuxTqpXW/Q1jlNaNmdIlwLXu\n", "fmqr7RopMEkhxaC0MSEAXQYsN7MqYXhxD+FDckpRgxJ0PpQbR9PSusAQMAMYNrN3xDZHjdduPApM\n", "UhhxN/eXgJOBfwCHESZjLwG+CXwMuN7dDzCzi9392dw62we9OIm3xdK6Uxsfj57XbjwKTFIIcQjx\n", "eWBdwvzFxoTg9BTw34QVuL2AZXHuqdBBCWh3O0AStF1AiuIrwEcJGydvi7d/QdjodxXwTmBP4Iy6\n", "CfFCG+STeJUxSSG4++FxI98BwB2EK3IsI5TFXQRsTdj0d2t+veyv+ZfsWMq7D+1SxiQDra5KwJqE\n", "iwa8HNiKcB7cbMLE67XAf06moDToFJhkoMUr5E4DbiZs/tsFeIiwX+l6wiWWlteVy5UBoMAkRTCV\n", "sHt7sbv/nlBbaSPgSWAbd78vz87JxGmOSYrgEeBXwOFmViFcjPIy4HPuviTXnklblDHJwItL/58m\n", "XCzgBOAg4GR3vyvXjknbFJikENz9HuB9hCvnbuvut+XbI+mEhnJSGO7+DPC3vPshnVPGJCLJUWAS\n", "keQoMIlIchSYRCQ5CkwikhwFJhFJjgKTiCRHgUlEkqPAJCLJUWASkeQoMIlIchSYRCQ5CkwikhwF\n", "JhFJjgKTiCRHgUlEkqPAJCLJUWASkeQoMIlIchSYRCQ5pWq1mncfRERWooxJRJKjwCQiyVFgEpHk\n", "KDCJSHIUmEQkOQpMIpKc/we122r+YxVm2AAAAABJRU5ErkJggg==\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sns.corrplot(dee_df[['register','college','distance']])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Linear Regression of REGISTER on COLLEGE" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: register R-squared: 0.035\n", "Model: OLS Adj. R-squared: 0.035\n", "Method: Least Squares F-statistic: 335.9\n", "Date: Wed, 18 Mar 2015 Prob (F-statistic): 1.03e-73\n", "Time: 22:15:12 Log-Likelihood: -5959.0\n", "No. Observations: 9227 AIC: 1.192e+04\n", "Df Residuals: 9225 BIC: 1.194e+04\n", "Df Model: 1 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [95.0% Conf. Int.]\n", "------------------------------------------------------------------------------\n", "Intercept 0.5741 0.007 80.391 0.000 0.560 0.588\n", "college 0.1769 0.010 18.326 0.000 0.158 0.196\n", "==============================================================================\n", "Omnibus: 603.688 Durbin-Watson: 1.955\n", "Prob(Omnibus): 0.000 Jarque-Bera (JB): 1441.712\n", "Skew: -0.689 Prob(JB): 0.00\n", "Kurtosis: 1.640 Cond. No. 2.74\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n" ] } ], "source": [ "result = smf.ols(formula = \"register ~ college\", data = dee_df).fit()\n", "print result.summary()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Two-Stage Least Squares Regression of REGISTER ~ COLLEGE where IV=DISTANCE\n", "###using statsmodels.formula.api.ols" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In two-stage least squares regression, we regress COLLEGE on DISTANCE and use the predictions from that model as the predictors for REGISTER." ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================================================\n", " FIRST STAGE\n", "==============================================================================\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: college R-squared: 0.012\n", "Model: OLS Adj. R-squared: 0.012\n", "Method: Least Squares F-statistic: 115.9\n", "Date: Wed, 18 Mar 2015 Prob (F-statistic): 7.35e-27\n", "Time: 22:15:12 Log-Likelihood: -6598.2\n", "No. Observations: 9227 AIC: 1.320e+04\n", "Df Residuals: 9225 BIC: 1.321e+04\n", "Df Model: 1 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [95.0% Conf. Int.]\n", "------------------------------------------------------------------------------\n", "Intercept 0.6091 0.008 78.812 0.000 0.594 0.624\n", "distance -0.0064 0.001 -10.764 0.000 -0.008 -0.005\n", "==============================================================================\n", "Omnibus: 54.830 Durbin-Watson: 1.703\n", "Prob(Omnibus): 0.000 Jarque-Bera (JB): 1463.752\n", "Skew: -0.190 Prob(JB): 0.00\n", "Kurtosis: 1.086 Cond. No. 19.7\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "\n", "\n", "==============================================================================\n", " SECOND STAGE\n", "==============================================================================\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: register R-squared: 0.001\n", "Model: OLS Adj. R-squared: 0.001\n", "Method: Least Squares F-statistic: 10.35\n", "Date: Wed, 18 Mar 2015 Prob (F-statistic): 0.00130\n", "Time: 22:15:12 Log-Likelihood: -6118.9\n", "No. Observations: 9227 AIC: 1.224e+04\n", "Df Residuals: 9225 BIC: 1.226e+04\n", "Df Model: 1 \n", "Covariance Type: nonrobust \n", "==================================================================================\n", " coef std err t P>|t| [95.0% Conf. Int.]\n", "----------------------------------------------------------------------------------\n", "Intercept 0.5157 0.049 10.632 0.000 0.421 0.611\n", "college_fitted 0.2837 0.088 3.216 0.001 0.111 0.457\n", "==============================================================================\n", "Omnibus: 658.930 Durbin-Watson: 1.938\n", "Prob(Omnibus): 0.000 Jarque-Bera (JB): 1638.898\n", "Skew: -0.726 Prob(JB): 0.00\n", "Kurtosis: 1.532 Cond. No. 23.4\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n" ] } ], "source": [ "print \"==============================================================================\"\n", "print \" FIRST STAGE\"\n", "print \"==============================================================================\"\n", "result = smf.ols(formula = \"college ~ distance\", data = dee_df).fit()\n", "print result.summary()\n", "dee_df['college_fitted'] = result.predict()\n", "\n", "print\n", "print\n", "print \"==============================================================================\"\n", "print \" SECOND STAGE\"\n", "print \"==============================================================================\"\n", "\n", "result = smf.ols(formula = \"register ~ college_fitted\", data=dee_df).fit()\n", "print result.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "^^^^^^ Not sure what's going on with the R2 statistic here (it's 0.001 here, versus 0.022 in the example), although everything else matches what we see from the Stata output in the published example " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Adding Covariates that do not satisfy the requirements of instrumental variables" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the case of the covariate of race/ethicity, we expect that there might be a relationship between race/ethnicity and distance to a community college, as well as a relationship between race/ethnicity and voter registration. \n", "\n", "While race/ethnicity fails the test for instrumental variables, it can still be included as a covariate in a multiple regression model. In such cases, it is essential to include covariates at both stages of a two-stage test." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Correlation Matrix" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": [ "iVBORw0KGgoAAAANSUhEUgAAASYAAADuCAYAAACUCKq3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n", "AAALEgAACxIB0t1+/AAAIABJREFUeJztnXmYXFXR/z8zkxUI+yKCKGsRQUGIyL6ILKLsKAKCoOyC\n", "bC+rsigCKiKKGxJA4CXqjyUoO4iALGERAggC3xAEeRGQANkwBJKZ+f1Rp5NOZ2ZuT6+3Z+rzPP1k\n", "0n1P97e7b9etU6dOVVt3dzdBEAR5or3ZAoIgCEoJwxQEQe4IwxQEQe4IwxQEQe4IwxQEQe4IwxQE\n", "Qe4YkvH4bGB4I4QEQbAAbc0W0EzaMvKYusnfB5RHTZBPXaGpPPKoaVATU7kgCHJHbgyTmbWZWVy1\n", "giDIjDE1BDPrkNTZbB1BEOSDpntMZjZEUqeZtZvZrWa2TbM1BUHQXJpumCTNNbN24I/Au8BjTZYU\n", "BEGTaZphSsaowE7AjsCHJM1Mj3c0RVgQBE2nKekCafo2NwW7l8BjXbsAPwZulHRgozXVgDzqCk3l\n", "kUdNg5qGGyYza5PUnTyiO4FhwKvAXen1LgV+L2m/RmmqEXnUFZrKI4+aBjUNXZUrWX0bC/wHGA9c\n", "A7wJ/Aw4EvhXI3UFQZAvGmaYzGxd3BBNSd7SYsBc4Gzg18CSwGqSLk7Ht0mK8ppBMAhpZPB7E+AJ\n", "M3sV2Ah4FPgK8AHwLLAH8H7h4DBKQTB4qbthKlp9ux/4MPAhYAY+bfseMArYHzhY0v2R/R0EQV2n\n", "coXVt/TffwM/BzYAbgG+APwCuASYLemdWhqlZBB/BXwS98QOlvRiyTGLAH8Gvi5JZjYMD76vAcwB\n", "viXpqRpo2Rk4HZ+6Xi7p0pLHlwV+B4wAXgMOkvSeme0JnIwHZ8dJuihNg8cCa6X7D5f0jyq09fk5\n", "9aXdzJYHHge2lTSpnq/Z2xgzWwO4AugCngG+mRZXfgZsBsxMT7sLMAv4CbAhvuhyhqTby9Bb0fdX\n", "9PglwNuSTjWzA4GvpYdGAusBK0iakfnBDSLq5jGZ2bCUEtBhZuOAwpe7P/AK8DQe5B4p6Z3CuBpO\n", "4XYDhknaFDgFuKBE3xjgPmBV/AcOcAgwK405BLi8WhFmNhT/MWwHbAUcmn7QxZwBXC1pS+AJ4LD0\n", "QzwP2BafBh9pZsvgn2OXpM2B7wDnVCmx18+pL+3psd8A/23Qa+4GDO9hzE+A09Jn1wbsmu7fANhe\n", "0jbpNhM/94akz243YHSW0Eq/v6LxhwHrks4xSVcUNOHJxEeHUVqYuhkmSR+kH9dE3DP7GPBlYHXg\n", "QOAs4EvFV8oax5U2A25Pz/sIMKbk8WH4yami+z5eNGYSsJKZLV6ljtHAZEnTJc0BHgC27E0rcBvw\n", "OUldwOj0g1oO6ADel/RH5p/4HwOmVqmvr8+pL+3n44sWrzfoNTfDP5vSMRtIui/9fRvwueR5rwmM\n", "NbMHzOyg9Pj2wL/N7Gbc6/xTGVor+v4AzGxTPJ76G0rSEdKFcZ1S7ytwam6YzOzKNAUBd59XAq4F\n", "jgKOAS4CXpP0PUm31rGqwOJ4LKtAZ3G2uaQJkl4tGfMk8MX0PjbGDcKiNdAxvej/M/Gk0t6Oebfw\n", "eNpDuAd+Fb4Hn4oU7r8C/yx/VwN9vX1OPWpP05Epku5M9/f3++v3a/YypqPktQuf3aL4Z7MfvqPg\n", "SDP7BLAssLqkLwI/BH5bptZ+f39m9iHckzqKnj+f0/CLc9ADNTVMZrYqPj27w8y2BR7Ev8j9gR/h\n", "U49hwFKFMZK667QCNwMPrBdoT15IX1wOzDCz+3FvahLwTt9DesbMzjaze/CrcrHXNYqFvZwZRceM\n", "AqYVHpA0Hjfuw4EDiu4/EI8zjTWzkZVoLHrt3j6n6SWPFbQdBGyX3t/6wJVmtkKdX7OnMZ14bKn0\n", "2FnARZJmS3oXuBuP5byNxzdJXtZavQmswff3JdwQ3orHCfc1swPScy8JrCXpr729/mCnZobJzIZK\n", "egm/ChyFnwDbSloVv9pvgF/dT5dUifvfXx7E9+AVvJ+/lzFmI+BuSVsA1wGvS3o/Y0yPSDo9xRFW\n", "ANYws6VScH1L4KHetAKfB+4zs1Fm9tcUq+vGYzmdZra/mZ2ajn0P/2FmGdy+6Otzeh5Ys0T7BElb\n", "Sdo6vb8ngQMk/aeer9nHmCfMbKv09+fxuOFawAOpYsVQYHM8SP9A0XOsRx+JvNV+f5J+LmlMeo4f\n", "AL+TdFU6ZkvgL2V8ToOWmqzKpYzuOckd3wK4Cr+Sfj/FaP6EG8ErJN1dmLrVOVfpBvyq/mD6/0Fm\n", "tg+wmKSxvYwR8P/M7DS83vkh1YpICwDHA3fgn8Flkl43s6WBsZL2BL6Pex2HAFOAfdOq3NW4kZoD\n", "PAVcja/8XGFmfwWGAsdUajwTfX5OPWmv4rUqfk0zW2hM+vcE3GschufDXZdW5a7CDcgc/Lx7zsxe\n", "BH5tZgXDcniW0Eq/vx6eqvhcXwt4sYdjgkTVe+Vs/obcdjyL+1Q8tvQInqv0WXzV5PJ0fLVGKa/7\n", "mvKoKzSVRx41DWpqsok3GaWHgeeAjfEVkT3w7O6fA7+SVCvXNa8nUR51habyyKOmQU2tEiz3A1YE\n", "tgGWwacc44Hdk6uLmZUTfA6CIKgs+G1mpQatE99q8ilJr+DlTB4CbjCzzQDCKAVBUC799phSoLsQ\n", "U/oiHiS+Fg8k3m5mf8FXQQ7FV+iWqZ3cIAgGA/0yTMkoFZLhngFeBj6Fx5Z2AY7D8zj+jpcx+Ri+\n", "9aTpmNnKuMZJeN7J7viU81XceJ6OB+6nVrnKVSuNT0ma0gwdPWhZ4POSNL2vsYOBVjifWpmyDVPJ\n", "6tvpuFE6FrgXWBrYTNKZKVP5QDypcpuU25QHRuPL/1Px7QN748u2J+NbFXbCp6Sn4UXrmq3xf/Cl\n", "52bR1+c16A0TrXE+tSxlGaYUuC4YpbG4l/QRPCXgYvwL2MvM7pI03sz+Bfxb0hv1El4BE/EEz2Xw\n", "E2kqsKakyWa2Cp4Uuo6kZp5ExRqb7Wn2+nk1VVV+yP35dMtQy0zJ+cIc5XI1sux0gZR/dBOwNh4/\n", "OhGPJd2C73j/qqQb6qq2RFPOyKOu0FQeedRUNbcusnamYdpp1vO5fN/9iTEtik/fdsRjRwfi07UV\n", "8LSAOxuU0R0EQRl0jKyuA5pl183aB9+YPxf38I/EDXyfddDKoWzDJOndtFXjA3w6d6SkHxeJDKMU\n", "BDmifUjVztC8ullm9hm8BtZuAGnj+NnAupJmm9nv8FX6oaS6WaVj+qW9PwenglZn4lUnf21mH0ul\n", "J+pZJSAIggroGNmeecugr7pZs4FNJM1O/x+S7uutbla/6HeCZSpcdhIe2HtZ89sxBUGQI9qGtmXe\n", "Mui1blZyRKYAmNnRwKKS/tzXmP5Q0ZaUZJyeS6KizVIQ5JCO4VVXNeqzplkyOD/Ca+TvWc6Ycqla\n", "eRilIMgn7R1tmbcMsmqa/QYvYLh70ZSukjpoC5GZLvD41ptW8rx1Y8N7J0A+l3bzuOQcmsojj5qq\n", "ZsKYT2c6DZs+9rde33da0CqssIHXwNoQb1b7WLrdVzTkp8CNpWPUjw46BcrJYwqCoPFUbSgf/sxG\n", "mb/fjR95NJcGOTPGFB5T2eTxqhuayiOPmqqmY2gjG23Xlro2vAyCoHm0ZceQcksYpiAYoLQPqS7z\n", "u5mEYQqCAUoZq265pXUnoXWkkM1e9P/W/YaDQUv7kI7MW14Jj6mEkrpT6wGzJT3XbF1B0F/a2lv3\n", "ehqGqYRklDrwxLDJwM5mdhbe2XVan4ODIEe08qpc6yqvMSX7eU7AO8KegBfB2xpYrQmygqBi2trb\n", "M295JTwmFqhl3oYb6+l4edTPAOfjpRz2NbNngDmxDSdoBdqH5NfwZDHoDVNJTOlivM72XXjZ4M/i\n", "pRxOAHaT9EHzlAZB/8hzcDuLQW2YimqZd+Ce0a64x/QecCk+nVsH+JKk26KSQtBKtHK6wKA2TJK6\n", "imqZv4dP387CC8n/DDhR0nsQKQNB61Gtx5RVWjcdswjwZ+DrkpTum8j8Tjr/lPSN/r72oDRMJe3K\n", "hwPLA3dLesrMvo7vmj4YeNfMLgY6o5Nw0GrUIF2g19K6AGY2Bg9/fJi04d/MRgBI2qaaFx50hqkQ\n", "U0p/LwHMxEs37GZm/wHewA3TSsCXJf28Dhp2xnvzzQUul3RpL8cdC6wg6dSi+y4Enpf0m0ZoMbNl\n", "8TZFI4DX8DIWBS+yp6vlcnhNnnUriclVogdYAvh90WHr4/3dfgtchhcymwN8S9JT/dSTVZB/Ib0p\n", "NDAW7zPXDRwu6R9mtj5wEb7S+z5wQD3bO9UgxrRAad1kiIoZhhuq/y26bz1gETO7A7cvp6USu/2i\n", "dcP2FVDU3rzDzG4GxuNXgSn4j+lo4CfA5cCfgGFmNqrXJ6xMw9D0GtsBWwGHmtnyJceMMLNxeNeJ\n", "wpVoOTO7DW+VVZM4VzlagDOAqyVtCTwBHJbGjsEN+qpFGncA7sQ90IbpkfSGpG3SVfo04HHcMBwK\n", "zJK0Kd6c8vIKZM3zGoBT8PMlS+/OQJekzYHvAOekIT8Fjko6x+PGs260tbdl3jLos0yupAmSXi0Z\n", "81/gfEk7AIcD4yoprTtoDJOZLZ1SAjqAPwBvA9fj7WeG400WNgGuxLuqnol3gplZYymjgcmSpkua\n", "AzwAbFlyzAjgCvyELpw9iyZN/0vtSnSUo2XeVRMvMv+59HfhaqmiYzuBbfHmj43WU4gDXgQckRYp\n", "Ps78K/4kYCUzW7yfmvoqyN+jXkl/JBlwvNVZ4fP4iqRCRceheFyzbrR1dGTeMqikTO4kYByApBfw\n", "39mK/dU+KAyTmf0QuNTMPpKaJywCTMOD3JfjJ88qwDvA/cAsYGNJj9dBzuIs2GJ7Jj4VmYekaamw\n", "e/F9L0t6tNFaSo55t/B4T1dLSXdJeqcZehI7A8+kHwTAk3hLoUKZ1+VwA99fTb15Db3qTRfBK3BD\n", "+bt03xtJy6bAN4EL+6mlX9TAY6qkTO5BJK/SzD6Mf0av91f7oDBM+MnzSeBMM1sajyEdjVvzR/E0\n", "gbmS3pd0E3CMpGdqKcDMzjaze/ApYvFVexSVexiN0DKj6JhRuEHPq5798NZiBS4HZpjZ/bh3Nwm/\n", "+PSHvryG6SWPLaBX0oF4nGls6sOGme0N/BrYSdLb/dTSL2qwifcGYLaZPYgbm+PMbB8zO6SPMZcB\n", "i5vZffjM5KBKFo4GdPA7eUqrAPviV7Pj8K4OF+Ne01dwA3WEpIcLeUr1aEkl6fSkaQjwrJkthc/H\n", "t8RzqBpGP7UUrppXAp9nwRrPedMzRtJDRf/fCF9tPT7FxDaS9H4/5T2Ie2LX9uA1PA+sWarXzPYH\n", "VpZ0Hj5d6wK6zOyreNxra0l1vxhVu+UkTYePKLl7ofrdxStwaWFp/6pemAHsMZnZMHwePwb4Rbr9\n", "jBSkxA3UOsB2kv7QqDyl9MUdD9wBTAAuk/S6mS1tZtf3MKSnQHdNgt9lavk+8BUzewDfovOLMp66\n", "In3V6EmrgdNLnxI4xswm4N93X1f63ujVa0hxpYX0AtcB65vZX/H41DH4quDP8EL+483snrQ5vG60\n", "D2nPvOWVAdklpWjv2yLAicDX8JPnSDyudAZwC76i834aU21Wdx7rRoem8sijpqp589sHZp7Py59z\n", "RS7f94CbyhXtfeuQNMvMzsWvVgfhV9ej8Pf9ZLFbH1tNgoFGe/aqW24ZUIapOE8JuCLN/d8Gfohf\n", "FY8Hhko6NB0fe9+CAUsUissJRaVLJuCrKXfjy7LLALsDI9N9hePDKAUDlrYWri6Q3+hXP0irOQVW\n", "xXNVTkirIkcD6+JJi2dJujc25AaDgba29sxbXml5j6mkntKRuGFaE/ismb2KJ9XNxL2ldyE/npKZ\n", "rQzsgi/BzsC9uquBV3Ev73TgVGBqBcvcA4ZW+JxKND4laUozdBTTyh5TyxumIqP0RLrrJuA3+B6m\n", "/QHDN5nWbbNkFYzGl7Cn4svKe+MJeSfjJVh2wrd5nAbkUX+jaIXPqVjj/+D7L5tKGVtOckvLG6bE\n", "gfiJuZmk98xsP/xk/Q4wRdJjhelbXrylxER8u8Iy+A9uKrCmpMlmtgq+grhOTo1qI2mFz6lY49NN\n", "1DGPVg5+D4g8JjPbFT8pviTpVjP7JR4/O6ooIF5vo5THXJjQVB551FQ1M392Qub5PuqYC3L5vgeK\n", "xzQBuBe40szuBL4E7F7YWpIzLykIGkOOu6Bk0brKi0iBxkPxfVXTgC9KuiVW34LBTA3KnjSNgeIx\n", "Ienf+H4oIGp0BwFVxpiyqnemYxaoYlrOmLKkV6U8x6QqATGFCwYtNfCYeq3eCT1XMc0aUy4D1jAF\n", "wWCnBp14+6reCT1XMc0aUxZhmIJgoNLRkX3rm0pqfvc5plwGTIwpCIIS2qoOs1ZS87uSMQuRaZjm\n", "/Pq2/j5n3Tl/fFfuYkcn7hHOZ5Av2jqq9jv6qt5ZyzELkal849GlteCbT46NQO4MJqGpXPKmqfpV\n", "5eozv28AtkvVOwEOMrN9gMUkjS13TCUvnJn5/fBzpdVKm8vGo5fg/PH5a4qbjGXeUhTymNEcmhrE\n", "7P/3o0xjO2Lvk3L5vnPregS9UxpMrCS4GAwC2tqzbzklgt8tRlGZl8WAJST9u5LgYjAIyHFmdxZh\n", "mFqIVAp4rpmtBtwITEkG6lhgQiSUBguQY48oizBMLYSkbjMbgTdMHIfXHhqLJ7k9Asxtorwgb7Sw\n", "x9S6JnUQkZorFJiL54lsApwNnAmsAGzXBGlBnqk+wbJphMeUc4piSssDG+DlZMcDZ+E5Im3AjriB\n", "CoL5VJ9g2TTCMOWYopjSKsBDeKfZHwN/wrsM7wlsDewg6aWmCQ3ySfUJlk0jpnI5JsWUOoDzgPuB\n", "TSRdDuwKPAZsDuwq6Yk+niYYrLS1Zd9yShimHFIcUypU4cS/q5Ep+H0wMFzSbEkzm6ExaAHaO7Jv\n", "OaV1fb0BSuom3JliSp8D/oJ3fhkHdOGbJJcH/tE8lUFL0MKldcMw5Ygio7QiXoBrRdwYrQ0cgNcy\n", "XwTYJmJKQRbdOfaIsgjDlCOSUVoWX2G7FrgT+D7wFLCVpHFmNkzSB83UGbQIVSZYZpXJNbOd8Waj\n", "c4HLJV2a7p+IL9QA/FPSN/r72mGYcoCZtUvqSnXKz8QbJ/4e+CveM288cKeZrYGfILV+/R5PsKLH\n", "l8XbY40AXgMOSv379gFOBGYD10q6sBlagCXwz6vA+sDJki5JYz4D/EDSNo3UlT6jPfHGnN3AOEkX\n", "FY2pWldf1MBjmlcmN2m9IN2HmQ3Fm8qOAWYBD5rZn/Cu11T7nlp3EjpASHlKXYW9b/hJfAmwKXAc\n", "8E9gD7yZ5+xabzspOsG2A7YCDk3xrWLOAK6WtCXe8fgwM1sGOBf4LF5OdVcz+1QztEh6Q9I26cdw\n", "GvA4nhGPmZ2U/h7eaF3J4zgP2BZPiD3SzJaula5Mql+V66tM7mhgsqTpkuYAD+CfzXrAImZ2h5n9\n", "JRm0fhOGqYmkmNJcM1sVTwe4C/cAvgXcDXwb+JakyZL+VScZPZ1gW5YcM+8EBW7Dg/KrAU9JmpaM\n", "5cM9jGuUFmBeZ5yLgCOKDPhk3LBXszZeka60uXp0WjldDugACtPwWujqk+6OjsxbBn2VyV2c+dM1\n", "cE9pCeC/wPmSdgAOB8ZVUv0iDFMTKVp9ewhv2vkmPpXbHzgSnxrcXGcZvZ1gvR3zbnr8BWAdM1s+\n", "tfDZFg/MN0NLgZ2BZyS9ULhD0niq30NYsa70He+Be1H34NOeWunqm+rLnvRVJnd6yWOj8Nbtk/AV\n", "ZNL38Da+iNMvwjA1ATPbysxGpv8uhXtHN+JxiMnApcApko6WNLlOGs42s3vwLPLFix4qnGDFzCg6\n", "ZhQwTdI0fKp5PW5AJwJvNUNL0WP74dPgmlArXckIrYRP2w6olb4suts7Mm8ZPAjsBNBDmdzngTXN\n", "bCkzG4Z7kA/hHv8FacyH8c/k9f5qj+B3g0nTjW2AY5Jx+jZ+hTmC+VtNngauqacOSacnPUOAZ81s\n", "KdwN3xLvaFxM4QS9Evg8cF8aN0bSFmY2HA/U/7AZWooeGyPpoUo01EOXmY3CPd7tJH1gZv8FOmkQ\n", "3dWXPemztK6ZHQ/cgTs4l0l63cwuA35rZoXv5aC6NCMIakvaZrIkvrrxNL6cOtHMDsBTAxbBf2DP\n", "N0jP3F5OsKWBsZL2TLquNLNDgCnAvmlcp5k9jv/YLpb0z2ZoATCz5VhwulVKxYsGVXxG75nZ1biR\n", "moOnfVxdK12ZVLkql+J0R5TcPano8ZspCTVImouHIqoian7XiKya34UqAenv3YG9gE/h7u9xwO7A\n", "HOBvxTGSKsljLevQ1CBmPH5HptFbfMMdcvm+w2NqAGY2VNKcVCXgAOABSfuZ2Sl4kHsi3s10vxS7\n", "CYKq6W6LzO+gD5JRWhk3QDcBc8ysG49TvIzHJn4YRimoJd3trfvzbl3lOSflblyAJyG+BRwDvIIH\n", "uH8FfBO4V9IBZjY+tpkEtaY7x2VNsgjDVAfSyts5wDL4svHKuHGaBXwXDyjuBcxOsafcGKXk2e2C\n", "Bzln4LGvq/HKmcvg2zJOBaZKqvn2mFbR1ArEJt6glAuBr+HbI57BV99+jP+YvgV8Af9xbVIIiOeI\n", "0fheval4JvPewFr4Vpnt8WlnJ/7e3hzEmnJPDdIFmkbrKs8xko7F40kHAM/hGclH4dsSJuG5MJtI\n", "erppIntnIp4w+ShuBM4DXkyJns/j7+NVSY00AHnUlHu62joyb3kl0gVqxIl7tBf2vhU25P4S+CLw\n", "Dl4Gdzc8vnQK8H9FlSnrSR6XwUNTg5jy7KOZ6QLLfXyjXL7v8JhqSDJKi+O72/+D7x+bgucr3Yvv\n", "Tp/TIKMUDHJa2WMKw1R7RuLZ25MlPYnHRlbCp2/rSfp3M8UF9cPMljaz65uto0B3W3vmLa9E8Lv2\n", "TMVLgBxrZl14M8qbgO9IerWpyoK6IukdvKVWLsizR5RFGKYakzZrnognT34P32byhWr3kQVBf6k2\n", "j6mS0rpZY8olv75cCyPpZeCreOfcMZKeaa6iYDBSgxjTvNK6+KLNBYUH+qjquRveWmyhMf0hPKY6\n", "kRL93mi2jmDw0l39QuMCpXXNrMfSugBmVqjquQlewbOnMWUTHlMQDFBq4DFVUlq3rzFlk+kxjZ71\n", "SH+fs85szxc/2e+CeA1gJS69q++ksEZz8OdymaISNIgarLr1t7TutIwxZZNpmJbYcPv+PmfdGb3G\n", "Ss2W0CM5NQS5MpaJ0JRN1SdTV/WG6UF818K1fZXWZcGqnt19jCmbTMM0/fE7K3neurHEhtvz3OT8\n", "pQKNXmMlLr0rX+d2MpR5s5Z5zLLOo6aq6e6u+i1VUlp3oTGVvHAEv4NggNJFU0rr9jSm34RhCoIB\n", "SlcLr22FYQqCAUoN0gWaRhimIBigdHWHxxQEQc5oZY+pdU1q0DRS6eAg53R1t2fe8kp+lQW5JLWi\n", "6k5dZoMc00Vb5i2vhGEKysLM2sxss9SK6sN4d9lPmFnr1tYY4ITHFAwGvgPcb2Y74vuhlgReiWqc\n", "+aWbtsxbXongd1Au1+GdSW7CS1ksCtxjZtcA7+L1eGY1UV9QQp49oixaV3nQEFLdHSQ9h3tNlwMn\n", "4f3cHgJ2xYuFLdMsjUHPdHe3Zd7ySnhMQa+kri9zzOyjwLHAk8A4vCrnwcB1kr5pZqMkzWym1mBh\n", "OnNseLIIwxT0iqROM1sReAR4FnhN0pVmNgyvZX6dmX0Mn8rVnJ5Kt/Zy3LHACpJOTf/fBzgRmA1c\n", "K+nCeuows2XxvncjgNeAgyS9Z2a74004u9O4i9Pxp+I78IcCv5B0ZTX6eiPPMaQswjAFWWwHvAJs\n", "lwzVqXgI4Higq16eUlHp1jF4a/UHzezG4qaWZjYCuAz4NB4Dw8yWAc7FW2ZNx+Ng90p6ol46gDOA\n", "qyVdZWYnA4cBP03jPoWXBXnWzH6f/r+JpE3NbFF8WlwX6hFjMrOReEfp5fDicF+T9FbJMYcAh+KG\n", "/PuSbkm5b68yfxPwQ5JO6+11IsYULICZDSn8m06mucB6wNbpkM3wdlSv1rkV1bzSrZLmAIXSrcWM\n", "AK4AzmF+2ZLVgackTUs73R/uYVytdcwrQYuXlf1c+nsOvnq5SJG+HYCnzeyP+ELCjVVo65Ou7uxb\n", "BRyBf75bAlfhccd5mNmHgKOBTfH3el4y7qsDj0vaJt16NUoQHlNQhJkNkTTXzD6Ce0T/B/wC2Be4\n", "zcweA9YBNk8/+nrSW+nWeUiaBvzZzL5WdPcLwDqpMP67eNPR8fXUUXLMu0WPX4A3P/0vcL2k6Wna\n", "9xG8S/NquGFauwp9vVKnVbnNgB+mv2/Hp7jFbAQ8mIz4HDObjF/YVgdWMrO7gfeA4yRNohfCMAXz\n", "SEbpo8BE3O1eCVgX2APYG/8B3iVJ9dJgZmcDm+Ptf4rrOo/Ce/b1iaSpZnYccD3wNv5e3up7VNU6\n", "ZuCfzZT0+LRk3I8CPopPAa82s72SluckzQUmmdlsM1u2dDpUC6ot9Gxm38AXPYr5D/NrevdkpEfR\n", "syF/DThX0vVmthk+Hdyot9cOwzTISdO1Y4Bf4tO2I4E/pf//HE8HWAn4kqQZvT1PrZB0etI1BI/L\n", "lJZu7ZM0boykLcxsOPBX5l/h66XjQWAn4Erg88B9+DSzE3g/tY5/E5/WPYB/3j9JGfSL4ga05nRW\n", "6TFJugyP4c0jdRoubEcq1PkuprTmd8GQP4efX0h6ML33XokYU7ATHqS9AI+DXA9MxoO3Z+KlU5cD\n", "lm2kqORRFEq3TmB+6dbe2nB3F43rNLPHgfuBS6ppNlqmju8DX0ktjD6Dr7S9gBuqCWZ2P+41XCHp\n", "FuAJM3sUn8YdWa9pcVd3W+atAgpGGOYb4WIeBbYws+FmtgQeo/sHvkBwLICZrYcvqPRKW3ff/l53\n", "1Pwuj1as+Z3ylDrN7CDcQxqLBzM/iidSvoHHQXaoYaA7j/W186ipam5/8oPME3LH9Yf1632nVbkr\n", "gRXxTrv7SnozTZ8nS7rJzA7GV+XagXMk3ZCM1NXAYrjn9M2IMQU9kozS8ngw8n/wQPdM/MR7FvgQ\n", "8NU6r76SLsiEAAAOHklEQVQFdaKrq/a2VtJ7wJd7uP/Cor8vBS4teXw6nrtVFmGYBiFpanEN7iU9\n", "hl/F/oUHNk/G4wJHA3MkzW6WzqA68lzWJIuIMQ0yzGwRfEn9p/hy71X4Uu7f8KnbSfgVcdEwSq1N\n", "V1db5i2vhMc0yJA0y8yOAd4ETgBeKlpifwlf8l670JM+D5jZysAueNbwDGB3PF7xKr55+HTgVGCq\n", "pPcHq6ZSKkygzAVhmAYhkmaa2XlAB3BCShl4CQ8AT8iTUUqMBg7Bl51vx3Oq1sKnndvjq0Sd+L60\n", "N3t5jsGgaQEqXHXLBWGYBikpC/l7+I/nJDxN4NOSXm6qsJ6ZiG+SXQY3AlOBNSVNNrNV8ETGdUr2\n", "rw1GTQvQyh5TpAvUiFZMFwAwsyXxQPc19czoLiKPS/N51FQ11zyUbZq+vEl7Lt93BL8HOWm/2TkN\n", "MkpBA6lTgmVDiKlcgKSuZmsIak9XC3+rYZiCYIDSyjGmMExBMEDJc55SFmGYgmCAEh5TEAS5ox4x\n", "pnJK66bjlsMrEawr6YNyxxWIVbkgGKA0o7QugJntANwJLN+fccVkekxvL7laPzTXnyWAia/3WWOq\n", "KYxeA5ZcLG/LIB3scvjzuXLob7y4LlVkgx7orE+P5KzSuuBJu9viZYX7M24emYZptdXXyDqk4ey3\n", "RT6Dentt3NFsCQuRU0OQK2OZyJumqk/yaqdyFZbWRdJdaXzx3cV10XscV0ymYfrni5OzDmkoq62+\n", "BuPuz9s55Mbyuofrc4mqlL027mCXw59vtowFSIYyb1eWAZn5XW3wu8LSur1RqIte1riIMQUDmrRB\n", "eVDS1dWdeauArNK6NRkXq3LBgKXQjqrZOppFnTK/fw1cmYoNvo+39qK4tG7Rsd1Z43ojDFMwIDGz\n", "oZLmmFkH3ixgFvDztDdwUNBZB8NUTmndovtWK/q7x3G9EVO5YECSjFI7Xpnzs3jrpJmphfigoLs7\n", "+5ZXwjAFAwozOyU1mwRv1T0c2FbSCcB+lARzBzKdnd2Zt7wShikYMJjZlni527dSnal/4+23v5EO\n", "WRkYYWaLNUliQwmPKQiajJm1S7oP75R7HfB7vMvtz/Cut/fhTRcvkfRu85Q2jvCYgqCJpNW3LpjX\n", "OfcqvO72wXi7853xrrfbSxpvZm2DIY2glT2mWJULWprkKc1Nq28/Al7H40g7ArcCSwHHSro1Hd8G\n", "UK+23Hkizx5RFmGYgpalqMV5G97AcxM8u3gLPK60M26czsXjTYPCIBXIqOefa2IqF7QsySi14zGl\n", "WXig+0i8gec4fBPpcpKeGAxTt1IixhQEDcTMij39JYFPAPsAG0q6Cg94G7CKpLebIDEX1KnsSUMI\n", "wxS0FIVtJmbWbmY7AkOBrYFHgF+b2aaSxgKbS3q8HjElM9vZzB41swlmdnAfxx2bGosW33ehmR1W\n", "Ky190dXZnXnLKxFjCloGM2srCnRPxKdufwdOAfYC7gCuMbN1SDGlOmgYCvwEGINPHx80sxuLG1ua\n", "2Qg8AP9pfJpZqOh4FbAm8Fw9tJVS4SbdXBCGKWgZiryeK4B/AsfjHYTPAH6M717/SHGL8zoEu0fj\n", "m1WnA5jZA8zPnSowImm8Ezee4DlVZ+I76xsS76qHR1RFad024FVgUjrkIUmn9fY6MZULck/ykAp/\n", "jwSeBu4FDgN+AHwcb3M+RNJDdc5TKi54Bj0UPZM0TdKfS+57WdKjddLUI51dXZm3Cqi0tO7qwOOS\n", "tkm3Xo0ShMcU5JyilIB2YFPgbeAiYH88vvQy8BRwhqSXoT4pAWZ2NrA58Ek8nlVgFDC11q9XC+o0\n", "lau0tO6GwEpmdjfwHnCcpEk9jAXCMAU5pjjQDTyM73X7ED51uwn4DZ7hvZ+kCfVMnpR0ekET8KyZ\n", "LQX8F5/GnV/r16sF1U7lalxa9zXgXEnXm9lm+HRwo95eOwxTkFuKjNJpeGziK8AeeIb3X/BA8jBJ\n", "zzUqoztpOh4PtLcDl0l63cyWBsZK2rNkSE96GhKV7qoywbLGpXUfA+am533QzPrsKBKGKcgdaZtJ\n", "IQByAvA94M/AO8BP8fIl60m6OB3f0G0mkm4Gbi657x1gz5L7ruxh7Hfrq24+XXPrUsKyUCL3b/Sv\n", "tO4Z+Pd3vpmtB7zS18FhmIJcUTR9a8Ozts83s+WBrwG/ACYAawDzuiwMpm0m/aFO2QKVltb9AXC1\n", "me2Ee04H9vUiYZiC3FCSp3Qv8KSZXSXpRDObC5wMrAh8XdK9hcB4MzXnma461NatorTudHzvYlmE\n", "YQpygZktL+nN5CndgwdLJwJHmtlM4BhgMeCLwLpmdk9fLaaDSLAMgqows5vxDO7TcI9oGTw7+iw8\n", "hvFV4DpJR6fEvb3waV3TMLOV8WqZk/BVqt3xlaZXcf2nA6cCUyW93wyN9fCYGkUYpiAPnAVMMrNv\n", "44l54/H4xDjgbtxLmg4g6StmtqKkKU3SWmA0cAiew3Q7sDewFj7d3B4PEHfixvbNXp6jrnS3sMcU\n", "md9B0yiqEvA4/qM+G/+x/xaPMa0J/AE4StJTaZ8akl5vvNqFmAj8DngUN0znAS9KmowH5o8CXi3e\n", "Q9doOju7Mm95JTymoCmY2V6SrkuB7j/iGd1n41tL2vH8me8B35X0TIo95aZ5ZSqnUpxY+feix+5O\n", "f17bUFEltHKhuDBMQcMxs+3wKgA/ADqAkcD6+H6qK/A9cCOBb0iaXRgXaQH9o7M+eUwNIQxT0Azu\n", "AY7Dc1tewHOUngZuwWMz1wPjwihVRyt7TBFjChpO6mTyKzyr2/C8pLn4KtdLeCvvWwdjOdxa0jm3\n", "M/OWV8IwBU1B0hxgLB5TOsTMfp6W1beRdF9KtmzdS34O6O7uzrzllTBMQdNIxumXeM7PN81s7Ywh\n", "QT/omtuVecsrEWMKmkragnIhcG2hnlK6P7+X8xahszO/U7UswjAFTSfFl16GefvlwijVgHokWJZT\n", "Wjdt6N07/fdWSd8rtyRvgZjKBbkijFLt6OzszLxVQJ+ldc1sNbziwCaSNga2N7NPZI0rJdNjWm31\n", "NSoRX1f22yKfizV7bdyRfVCDufHiXIZt8mh8mqqpu7ubtrYFzuuqT/I6bUnJKq37CrBD0QVmKDC7\n", "jHELkGmY7n76vTL1NobPfmIk7zz512bLWIil19+KR54vt5hfY/jM2ktyxT35sgEHbtPGj67PV+zj\n", "pD07oEGdSwoU6k4V/l9ilGpCV5UxpkpK66b39E5K9TgfmCjpBTMrbuLQY0neYiLGFAQNJtWRmmtm\n", "q+DVOBcDLgFeSyuVNaFaj6nS0rqpr97luCE6Mt09A+8w0+u4YiLGFAQNJnV9WR7fvLw4/uO9Ad/I\n", "TKpzXjV1ijEVSutCD6V1k6f0J+BJSUcUTen6HFdKeExB0CBKVhw3A/4B3AYcDPwLb3X+j6J651XR\n", "XVnfuCz6LK2L733cEhhqZp9PY07pbVxvhGEKggZQVMt8RNoD+CKwCV7e5RjgY8BOqZTwzFq8ZrUx\n", "pp4os7TuyF6GLzSuN8IwBUGdKYopfRS41MxeBK4BLgC+hef8jAa2rZVRgiitGwRBH6SY0tJ4VYUH\n", "gUck3Z3iTNvisaUHJL1Uy9ftyvEm3SzCMAVBnSiJKX0Uj62cmxp07g98WdLOLNhyvGZ0d+d3L1wW\n", "YZiCoA4UxZRGAsOAt/AW598ws98AnwaWNLNFJM2qh4b7btg8n5nIZRCGKQhqTFFM6WPAxfjv7FDg\n", "cLzUyy54guGO9TJKrU4YpiCoMSmm9GG8UcF9eGD7YWADYG1gJbxRwf81T2W+iQTLIKgRqbECqZvL\n", "2sADwI+Af+L7xV4B1pP0UBilvgnDFAQ1IMWUCp7Sw3hM6Va8rvkv8Qadt+L5S0EGMZULgioxs/YU\n", "U/owXtLjSUl/N7PJwD7ApcBwYIykfzVTa6sQHlMQVEEKdHelPKXPA58FhqXp3BzgStxYbR1GqXzC\n", "MAVBhSSjVNiQ+zTwX7y5wn7ACZLmSLoK+LakfzRTa6sRU7kgqJAio3QP3jL8JtxLGgqca2ZzJF0Q\n", "VTn7T3hMQdBPCqtvieHAIsDGwCclfQBcCJyIVw4IKiA8piDoB0UZ3cvhzTrfwsuVXAf8yswOk/Qo\n", "vkE3qJDwmIKgHxRVCZgInAncjBfa3wtYAbjQzIY3UeKAIDymIOgHqULjiXhO0h+B8cBq+G9pY2Bo\n", "6igcVEF4TEGQgZkNS/8WqgW8B3wOTwM4CejC9729IikSKGtAGKYg6AMzW1bSB2a2EnComW2A5yat\n", "im8zWR/PX7q7iTIHHDGVC4JeMLOrgS4z+xHwB7yL7HLAdviG3G/iNa63kPRc04QOQMJjCoLemQZs\n", "CPwUX/rfBJ++3QmsJekQ4GBJf2+exIFJGKYgKMHMNkhB7hOAO4B1caP0Dt4A8gZgrJktQT67Crc8\n", "MZULgiLSHrdDgU/gQe7D8M6xX8XbWp8DHAKMlDS9t+cJqiM8piBIpL1vc/CWSp/BPaXXJZ2Jx5gO\n", "AE4Fpkt6rWlCBwFhmIKABeopLQV8Eu8muxRwiZktC/wWL19yca0aUga9E1O5YNDT1dVVyOheFbgd\n", "34j7Bl6Bcl/gC3gnky/Vsu9b0DvhMQWDnvb2dlI3k0vxfKStgb8Di+KF3sYDJ4dRahxhmILAaQeW\n", "ByZLegv4Lr737QVJ35D0VFPVDTJiKhcEzvt4sbdDzOwDPJFyCDClqaoGKeExBQFeNQD4Nm6cvgVs\n", "D+wa3UyaQ3hMQZCQ9JKZ7YcXfkPStCZLGrSEYQqCIlIFyg+arWOw09bdHRn1QRDki4gxBUGQO8Iw\n", "BUGQO8IwBUGQO8IwBUGQO8IwBUGQO8IwBUGQO8IwBUGQO/4/XMeBCWj26M0AAAAASUVORK5CYII=\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sns.corrplot(dee_df[['register','college','distance', 'black','hispanic','otherrace']])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Two-Stage Least Squares Regression of REGISTER ~ COLLEGE + BLACK + HISPANIC + OTHERRACE where IV=DISTANCE" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================================================\n", " FIRST STAGE\n", "==============================================================================\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: college R-squared: 0.022\n", "Model: OLS Adj. R-squared: 0.021\n", "Method: Least Squares F-statistic: 51.24\n", "Date: Wed, 18 Mar 2015 Prob (F-statistic): 9.75e-43\n", "Time: 22:15:20 Log-Likelihood: -6554.4\n", "No. Observations: 9227 AIC: 1.312e+04\n", "Df Residuals: 9222 BIC: 1.315e+04\n", "Df Model: 4 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [95.0% Conf. Int.]\n", "------------------------------------------------------------------------------\n", "Intercept 0.6431 0.009 71.039 0.000 0.625 0.661\n", "distance -0.0069 0.001 -11.636 0.000 -0.008 -0.006\n", "black -0.0577 0.016 -3.613 0.000 -0.089 -0.026\n", "hispanic -0.1162 0.013 -8.766 0.000 -0.142 -0.090\n", "otherrace 0.0337 0.024 1.404 0.160 -0.013 0.081\n", "==============================================================================\n", "Omnibus: 53.726 Durbin-Watson: 1.715\n", "Prob(Omnibus): 0.000 Jarque-Bera (JB): 1408.886\n", "Skew: -0.188 Prob(JB): 1.16e-306\n", "Kurtosis: 1.123 Cond. No. 62.3\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "\n", "\n", "==============================================================================\n", " SECOND STAGE\n", "==============================================================================\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: register R-squared: 0.005\n", "Model: OLS Adj. R-squared: 0.004\n", "Method: Least Squares F-statistic: 10.53\n", "Date: Wed, 18 Mar 2015 Prob (F-statistic): 1.65e-08\n", "Time: 22:15:20 Log-Likelihood: -6103.0\n", "No. Observations: 9227 AIC: 1.222e+04\n", "Df Residuals: 9222 BIC: 1.225e+04\n", "Df Model: 4 \n", "Covariance Type: nonrobust \n", "==================================================================================\n", " coef std err t P>|t| [95.0% Conf. Int.]\n", "----------------------------------------------------------------------------------\n", "Intercept 0.5266 0.047 11.200 0.000 0.434 0.619\n", "college_fitted 0.2489 0.082 3.041 0.002 0.088 0.409\n", "black 0.0617 0.015 4.007 0.000 0.032 0.092\n", "hispanic 0.0283 0.015 1.879 0.060 -0.001 0.058\n", "otherrace -0.1067 0.023 -4.604 0.000 -0.152 -0.061\n", "==============================================================================\n", "Omnibus: 654.082 Durbin-Watson: 1.939\n", "Prob(Omnibus): 0.000 Jarque-Bera (JB): 1619.291\n", "Skew: -0.723 Prob(JB): 0.00\n", "Kurtosis: 1.543 Cond. No. 22.6\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n" ] } ], "source": [ "print \"==============================================================================\"\n", "print \" FIRST STAGE\"\n", "print \"==============================================================================\"\n", "result = smf.ols(formula = \"college ~ distance + black + hispanic + otherrace\", data = dee_df).fit()\n", "print result.summary()\n", "dee_df['college_fitted'] = result.predict()\n", "\n", "print\n", "print\n", "print \"==============================================================================\"\n", "print \" SECOND STAGE\"\n", "print \"==============================================================================\"\n", "\n", "result = smf.ols(formula = \"register ~ college_fitted + black + hispanic + otherrace\", data=dee_df).fit()\n", "print result.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Interactions Between the Endogenous Question Predictor and Exogenous Covariates in the Second Stage Model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case, we explore whether interactions between college and race/ethnicity are significant predictors of voter registration. Here, it's important to meet the \"rank condition\": that \"for every endogenous predictor included in the second stage, there must be at least one instrument included in the first stage.\"\n", "\n", "To do this, we need to create a series of stage-one instruments, one for the main effect, and one for each interaction. In " ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================================================\n", " FIRST STAGE\n", "==============================================================================\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: college R-squared: 0.022\n", "Model: OLS Adj. R-squared: 0.021\n", "Method: Least Squares F-statistic: 29.57\n", "Date: Wed, 18 Mar 2015 Prob (F-statistic): 1.12e-40\n", "Time: 22:15:20 Log-Likelihood: -6553.3\n", "No. Observations: 9227 AIC: 1.312e+04\n", "Df Residuals: 9219 BIC: 1.318e+04\n", "Df Model: 7 \n", "Covariance Type: nonrobust \n", "======================================================================================\n", " coef std err t P>|t| [95.0% Conf. Int.]\n", "--------------------------------------------------------------------------------------\n", "Intercept 0.6452 0.010 64.147 0.000 0.625 0.665\n", "distance -0.0071 0.001 -9.830 0.000 -0.009 -0.006\n", "black -0.0651 0.023 -2.858 0.004 -0.110 -0.020\n", "hispanic -0.1276 0.019 -6.588 0.000 -0.166 -0.090\n", "otherrace 0.0590 0.035 1.681 0.093 -0.010 0.128\n", "distance:black 0.0009 0.002 0.442 0.658 -0.003 0.005\n", "distance:hispanic 0.0013 0.002 0.818 0.413 -0.002 0.004\n", "distance:otherrace -0.0030 0.003 -1.018 0.309 -0.009 0.003\n", "==============================================================================\n", "Omnibus: 53.396 Durbin-Watson: 1.715\n", "Prob(Omnibus): 0.000 Jarque-Bera (JB): 1407.217\n", "Skew: -0.188 Prob(JB): 2.67e-306\n", "Kurtosis: 1.124 Cond. No. 93.0\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: register R-squared: 0.005\n", "Model: OLS Adj. R-squared: 0.004\n", "Method: Least Squares F-statistic: 6.691\n", "Date: Wed, 18 Mar 2015 Prob (F-statistic): 6.29e-08\n", "Time: 22:15:20 Log-Likelihood: -6100.6\n", "No. Observations: 9227 AIC: 1.222e+04\n", "Df Residuals: 9219 BIC: 1.227e+04\n", "Df Model: 7 \n", "Covariance Type: nonrobust \n", "=====================================================================================\n", " coef std err t P>|t| [95.0% Conf. Int.]\n", "-------------------------------------------------------------------------------------\n", "Intercept 0.4640 0.056 8.360 0.000 0.355 0.573\n", "college_fitted 0.3587 0.097 3.703 0.000 0.169 0.549\n", "black 0.2780 0.163 1.702 0.089 -0.042 0.598\n", "hispanic 0.1765 0.121 1.456 0.145 -0.061 0.414\n", "otherrace 0.1742 0.176 0.988 0.323 -0.171 0.520\n", "collegeXblack -0.3986 0.303 -1.315 0.189 -0.993 0.196\n", "collegeXhispanic -0.2929 0.249 -1.177 0.239 -0.781 0.195\n", "collegeXotherrace -0.4634 0.285 -1.623 0.105 -1.023 0.096\n", "==============================================================================\n", "Omnibus: 653.574 Durbin-Watson: 1.940\n", "Prob(Omnibus): 0.000 Jarque-Bera (JB): 1616.409\n", "Skew: -0.723 Prob(JB): 0.00\n", "Kurtosis: 1.545 Cond. No. 88.9\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n" ] } ], "source": [ "print \"==============================================================================\"\n", "print \" FIRST STAGE\"\n", "print \"==============================================================================\"\n", "# generate the stage one main effect instrument\n", "result = smf.ols(formula = \"college ~ distance + black + hispanic + otherrace +\" +\n", " \"distance:black + distance:hispanic + distance:otherrace\", data = dee_df).fit()\n", "dee_df['college_fitted'] = result.predict()\n", "print result.summary()\n", "\n", "# generate the stage one interaction instrument for distance:black\n", "# note that we have DROPPED the irrelevant terms. \n", "# The full form for each interaction, which gives the exact same result, is:\n", "# result = smf.ols(formula = \"college:black ~ distance + black + hispanic + otherrace +\" +\n", "# \"distance:black + distance:hispanic + distance:otherrace\", data = dee_df).fit()\n", "\n", "result = smf.ols(formula = \"college:black ~ distance + black + distance:black\", data = dee_df).fit()\n", "dee_df['collegeXblack'] = result.predict()\n", "\n", "\n", "# generate the stage one interaction instrument for distance:hispanic\n", "result = smf.ols(formula = \"college:hispanic ~ distance + hispanic + distance:hispanic\", data = dee_df).fit()\n", "dee_df['collegeXhispanic'] = result.predict()\n", "\n", "# generate the stage one interaction instrument for distance:hispanic\n", "result = smf.ols(formula = \"college:otherrace ~ distance + otherrace + distance:otherrace\", data = dee_df).fit()\n", "dee_df['collegeXotherrace'] = result.predict()\n", "\n", "# generate the final model, that includes these interactions as predictors\n", "result = smf.ols(formula = \"register ~ college_fitted + black + hispanic + otherrace +\" +\n", " \"collegeXblack + collegeXhispanic + collegeXotherrace\", data = dee_df).fit()\n", "print result.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "^^^ in this particular case, we find no significant interactions and fall back on our previous model, which simply included race/ethnicity as a covariate" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Binomial Regression: Logistic Model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this example, we use a logistic model with a two-stage least-squares regression. NOTE: This is not attempted in the textbook example, so I cannot be completely certain about this, unlike the above results." ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================================================\n", " FIRST STAGE\n", "==============================================================================\n", " Generalized Linear Model Regression Results \n", "==============================================================================\n", "Dep. Variable: college No. Observations: 9227\n", "Model: GLM Df Residuals: 9222\n", "Model Family: Binomial Df Model: 4\n", "Link Function: logit Scale: 1.0\n", "Method: IRLS Log-Likelihood: -6253.9\n", "Date: Wed, 18 Mar 2015 Deviance: 12508.\n", "Time: 22:15:20 Pearson chi2: 9.23e+03\n", "No. Iterations: 6 \n", "==============================================================================\n", " coef std err z P>|z| [95.0% Conf. Int.]\n", "------------------------------------------------------------------------------\n", "Intercept 0.5840 0.038 15.447 0.000 0.510 0.658\n", "distance -0.0283 0.002 -11.456 0.000 -0.033 -0.023\n", "black -0.2372 0.066 -3.621 0.000 -0.366 -0.109\n", "hispanic -0.4744 0.055 -8.699 0.000 -0.581 -0.368\n", "otherrace 0.1426 0.101 1.414 0.157 -0.055 0.340\n", "==============================================================================\n", "\n", "\n", "==============================================================================\n", " SECOND STAGE\n", "==============================================================================\n", " Generalized Linear Model Regression Results \n", "==============================================================================\n", "Dep. Variable: register No. Observations: 9227\n", "Model: GLM Df Residuals: 9222\n", "Model Family: Binomial Df Model: 4\n", "Link Function: logit Scale: 1.0\n", "Method: IRLS Log-Likelihood: -5825.1\n", "Date: Wed, 18 Mar 2015 Deviance: 11650.\n", "Time: 22:15:20 Pearson chi2: 9.23e+03\n", "No. Iterations: 6 \n", "==================================================================================\n", " coef std err z P>|z| [95.0% Conf. Int.]\n", "----------------------------------------------------------------------------------\n", "Intercept 0.0572 0.211 0.272 0.786 -0.356 0.470\n", "college_fitted 1.1312 0.367 3.079 0.002 0.411 1.851\n", "black 0.2899 0.073 3.995 0.000 0.148 0.432\n", "hispanic 0.1284 0.068 1.886 0.059 -0.005 0.262\n", "otherrace -0.4589 0.101 -4.566 0.000 -0.656 -0.262\n", "==================================================================================\n" ] } ], "source": [ "print \"==============================================================================\"\n", "print \" FIRST STAGE\"\n", "print \"==============================================================================\"\n", "result = smf.glm(formula = \"college ~ distance + black + hispanic + otherrace\", \n", " data=dee_df,\n", " family=sm.families.Binomial()).fit()\n", "print result.summary()\n", "dee_df['college_fitted'] = result.predict()\n", "\n", "print\n", "print\n", "print \"==============================================================================\"\n", "print \" SECOND STAGE\"\n", "print \"==============================================================================\"#\n", "result = smf.glm(formula = \"register ~ college_fitted + black + hispanic + otherrace\",\n", " data=dee_df,\n", " family=sm.families.Binomial()).fit()\n", "print result.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Binomial Regression: Probit Model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this example, we use a probit model with a two-stage least-squares regression. NOTE: This is not attempted in the textbook example, so I cannot be completely certain about this, unlike the above results." ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==============================================================================\n", " FIRST STAGE\n", "==============================================================================\n", "Optimization terminated successfully.\n", " Current function value: 0.677791\n", " Iterations 4\n", " Probit Regression Results \n", "==============================================================================\n", "Dep. Variable: college No. Observations: 9227\n", "Model: Probit Df Residuals: 9222\n", "Method: MLE Df Model: 4\n", "Date: Wed, 18 Mar 2015 Pseudo R-squ.: 0.01585\n", "Time: 22:15:20 Log-Likelihood: -6254.0\n", "converged: True LL-Null: -6354.7\n", " LLR p-value: 1.862e-42\n", "==============================================================================\n", " coef std err z P>|z| [95.0% Conf. Int.]\n", "------------------------------------------------------------------------------\n", "Intercept 0.3640 0.023 15.566 0.000 0.318 0.410\n", "distance -0.0176 0.002 -11.504 0.000 -0.021 -0.015\n", "black -0.1474 0.041 -3.607 0.000 -0.227 -0.067\n", "hispanic -0.2959 0.034 -8.709 0.000 -0.363 -0.229\n", "otherrace 0.0888 0.062 1.423 0.155 -0.033 0.211\n", "==============================================================================\n", "\n", "\n", "==============================================================================\n", " SECOND STAGE\n", "==============================================================================\n", "Optimization terminated successfully.\n", " Current function value: 0.631312\n", " Iterations 4\n", " Probit Regression Results \n", "==============================================================================\n", "Dep. Variable: register No. Observations: 9227\n", "Model: Probit Df Residuals: 9222\n", "Method: MLE Df Model: 4\n", "Date: Wed, 18 Mar 2015 Pseudo R-squ.: 0.003563\n", "Time: 22:15:20 Log-Likelihood: -5825.1\n", "converged: True LL-Null: -5845.9\n", " LLR p-value: 1.962e-08\n", "==================================================================================\n", " coef std err z P>|z| [95.0% Conf. Int.]\n", "----------------------------------------------------------------------------------\n", "Intercept 0.0427 0.129 0.330 0.742 -0.211 0.296\n", "college_fitted 0.6903 0.226 3.060 0.002 0.248 1.133\n", "black 0.1752 0.044 4.021 0.000 0.090 0.261\n", "hispanic 0.0782 0.042 1.879 0.060 -0.003 0.160\n", "otherrace -0.2834 0.063 -4.533 0.000 -0.406 -0.161\n", "==================================================================================\n" ] } ], "source": [ "import patsy\n", "print \"==============================================================================\"\n", "print \" FIRST STAGE\"\n", "print \"==============================================================================\"\n", "a,b = patsy.dmatrices(\"college ~ distance + black + hispanic + otherrace\",\n", " dee_df,return_type=\"dataframe\")\n", "result = sm.Probit(a,b).fit()\n", "print result.summary()\n", "dee_df['college_fitted'] = result.predict()\n", "\n", "\n", "print\n", "print\n", "print \"==============================================================================\"\n", "print \" SECOND STAGE\"\n", "print \"==============================================================================\"#\n", "\n", "a,b = patsy.dmatrices(\"register ~ college_fitted + black + hispanic + otherrace\",\n", " dee_df,return_type=\"dataframe\")\n", "result = sm.Probit(a,b).fit()\n", "\n", "print result.summary()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.6" } }, "nbformat": 4, "nbformat_minor": 0 }