Title IX Health Behaviors Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Importing the NHIS data into a pandas dataframe**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "        year  serial  strata  psu         nhishid  hhweight     region  \\\n",
      "481367  2004   42087    5058    1  00002004071562      2491  Northeast   \n",
      "481368  2004   42088    5058    2  00002004071567      2487  Northeast   \n",
      "481369  2004   42088    5058    2  00002004071567      2487  Northeast   \n",
      "481370  2004   42088    5058    2  00002004071567      2487  Northeast   \n",
      "481371  2004   42088    5058    2  00002004071567      2487  Northeast   \n",
      "\n",
      "        pernum           nhispid     hhx  ...  \\\n",
      "481367       2  0020040715620102  071562  ...   \n",
      "481368       1  0020040715670101  071567  ...   \n",
      "481369       2  0020040715670102  071567  ...   \n",
      "481370       3  0020040715670103  071567  ...   \n",
      "481371       4  0020040715670104  071567  ...   \n",
      "\n",
      "                                            usborn racenew  height  \\\n",
      "481367  No, born outside U.S. and U.S. territories   White     NIU   \n",
      "481368                           Yes, born in U.S.   White     NIU   \n",
      "481369                           Yes, born in U.S.   White      70   \n",
      "481370                           Yes, born in U.S.   White     NIU   \n",
      "481371                           Yes, born in U.S.   White     NIU   \n",
      "\n",
      "                                     alcstat1  smokestatus2 mod10dmin  \\\n",
      "481367                                    NIU           NIU         0   \n",
      "481368                                    NIU           NIU         0   \n",
      "481369  Current drinker (1+ drinks past year)  Never smoked         0   \n",
      "481370                                    NIU           NIU         0   \n",
      "481371                                    NIU           NIU         0   \n",
      "\n",
      "       mod10fwk vig10dmin vig10fwk  strongfwk  \n",
      "481367        0         0        0          0  \n",
      "481368        0         0        0          0  \n",
      "481369       95         0       95         95  \n",
      "481370        0         0        0          0  \n",
      "481371        0         0        0          0  \n",
      "\n",
      "[5 rows x 31 columns]\n"
     ]
    }
   ],
   "source": [
    "# The data was downloaded from IPUMS and then formatted into a STATA data file\n",
    "# Importing necessary libraries for the entire project\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import sqlite3\n",
    "from matplotlib import style\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import statsmodels.formula.api as smf\n",
    "\n",
    "style.use(\"seaborn-colorblind\")\n",
    "\n",
    "# The rdd library allows us to calculate the optimal bandwidth to select our sample\n",
    "from rdd import rdd\n",
    "\n",
    "# Defining the filepath that contains the NHIS data, reading it in, and putting into a pandas dataframe\n",
    "\n",
    "dtafile = r'C:\\Users\\devin\\Documents\\School\\ECON6060FinalProject\\NHIS2000-2004.dta'\n",
    "df = pd.read_stata(dtafile)\n",
    "\n",
    "# Outputting a sample of the dataframe to make sure it loaded in correctly\n",
    "\n",
    "print(df.tail())\n",
    "\n",
    "# Putting the data into a SQLite database\n",
    "\n",
    "con = sqlite3.connect(r'C:\\Users\\devin\\Documents\\School\\ECON6060FinalProject\\ProjectData.sqlite')\n",
    "df.to_sql(\"initial_data\", con, if_exists=\"replace\")\n",
    "con.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cleaning the data and creating necessary variables**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEICAYAAABfz4NwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZa0lEQVR4nO3de7SddX3n8fdHoggoKqCnGUKFLuMFWJZIiri0ncPClqgoOAOdUJGoTKOIa3Qt6Ah2rdFOzVp1NYjDWKixKJGqmOVlQCutFDheplwMwshNaiqMpEQQuZjgQAn5zh/P78jmcC77nORkJznv11p77Wd/n+f37N/zy8n+nOeyn5OqQpKkZwy6A5KkHYOBIEkCDARJUmMgSJIAA0GS1BgIkiTAQJB2WkkqyUva9EVJPtqmfzfJHYPtnXZGBoJ2SElGkjyYZPcB9uHA9qE7r49l35HkiSSb2uPOJJ9N8tLt0ddeVfXdqnrZ9n5f7fwMBO1wkhwI/C5QwFsG25tpuaaqngM8D3g98P+AG5IcOthuSf0xELQjOgW4FrgIWDZaTLJvkq8n+WWS7yf5aJLv9cx/eZIrkjyQ5I4kfzjVGyV5U5Ib2zrvTvKRntnfac8Ptd/6X9NP56vqiar6l6p6L/Bt4NfrTHJkkn9K8lCS/5NkuGfeO5PcnmRjkp8kefeYvv5Jkg1J7knyrkm2aTjJ+p7XdyU5M8kPkzyc5EtJnt0z/9gkN7U+/VOSV/azndr1GAjaEZ0CfL49jkky1Op/BTwC/AZdUPSGxV7AFcAXgBcBJwHnJzlkivd6pL3f84E3AaclOb7N+732/Pyqek5VXTODbfkq3d4OSfYH/g74KLAPcCbwlSQvbMveBxwL7A28Ezg3yata2yVt+d8HFtLtgUzHHwJLgIOAVwLvaOt9FfAZ4N3AvsCngMsGeahOg2MgaIeS5HXAi4E1VXUD8C/AHyXZDfiPwIer6ldVdRuwuqfpscBdVfXZqtpcVT8AvgKcMNn7VdVIVd1cVVuq6ofAF4F/vw036R66D3+Ak4FvVtU32/tdAawF3tj68ndtz6Kq6tvAt2hhQveB/tmquqWqHqFnr6NP51XVPVX1APB14LBW/2PgU1V1XduzWQ08Bhw5o63VTs1A0I5mGfCtqrq/vf5Cq70QmAfc3bNs7/SLgVe3wx4PJXkIeBvd3sSEkrw6ydVJfp7kYeA9wH7bZlMA2B94oKePJ47p4+uA+a0vb0hybTvk9RBdUIz25d/x1O39v9Psx896pn8FPKenT2eM6dMB7f00x0x59YS0vSTZg+434d2SjH6A7U53OGcI2AwsAP65zTugp/ndwLer6ven+bZfAD4JvKGqHk3yCZ78EN4WtwJ+K/Ddnj5eXFV/PHahdojmK3SHry6tqseT/C8gbZENPHV7f3Mb9G20TyuqasU2Wp92Yu4haEdyPPAEcDDdIY3DgFfQfaCeQnc8/iNJ9kzy8lYb9Q3gpUnenuSZ7fE7SV4xxXs+F3ighcERwB/1zPs5sAX4relsRJLdkhyU5H8Cw8CftVl/C7w5yTFtmWe3E8ALgGfRhd/Pgc1J3gD8Qc9q1wDvSHJwkj2BD0+nT5P4NPCetqeUJHu1E+3P3Ubr107EQNCOZBndcfKfVtXPRh90v8G/DXgf3SWdPwMupjve/xhAVW2k+wBdSnfc/mfAx+g+ZCfzXuC/J9kI/De6D17aOn8FrAD+dzucMtVx9dck2QT8EhihOzn8O1V1c1vf3cBxwIfoPvjvBv4EeEbr/39p7/8gXTBd1tOXy4FPAFcB69rzVquqtXTnET7Z3ncd7YSz5p74B3K0s0ryMeA3qmrZlAtLmpJ7CNpptO8ZvLId2jgCOBX42qD7Je0qDATtTJ5Ldx7hEbpDK+cAl07VKMmtefKWEr2Pt03nzZP89QTr+esZbY20g/GQkSQJ6GMPoV0JcX37mv2tSf6s1fdJd5uAH7fnF/S0OTvJunS3Dzimp354kpvbvPOSpNV3b1+nX5fkunT3spEkbUdT7iG0D+29qmpTkmcC3wPeD/wHusv1/iLJWcALquqDSQ6mu/rjCLovt/wj8NKqeiLJ9a3ttcA36b49eXmS9wKvrKr3JFkKvLWq/tNk/dpvv/3qwAMP3IpNn9wjjzzCXnvtNWvr39k5PhNzbCbn+ExutsfnhhtuuL+qXjjuzKrq+wHsCfwAeDVwBzC/1ecDd7Tps4Gze9r8A/CatsyPeuon0X1l/tfLtOl5wP20sJrocfjhh9dsuvrqq2d1/Ts7x2dijs3kHJ/Jzfb4AGtrgs/Vvk4qty/R3ER3860rquo6YKiqNrRQ2UB3QzHovqrf+xX79a22f5seW39Km6raDDxMd6MtSdJ20tetK6rqCeCwJM8HvpbJ7++ecWo1SX2yNk9dcbIcWA4wNDTEyMjIJN3YOps2bZrV9e/sHJ+JOTaTc3wmN8jxmda9jKrqoSQjdLfRvTfJ/KrakGQ+3d4DdL/5995zZQHdN0fXt+mx9d4269P9darn8eQNwXrffxWwCmDx4sU1PDw8ne5Py8jICLO5/p2d4zMxx2Zyjs/kBjk+/Vxl9MK2ZzB687HXAz+i+1r96DdEl/Hk9eCXAUvblUMH0d27/fp2WGljuj8QEtpNvHrajK7rBOCqdqxLkrSd9LOHMB9Yne5+9M+gu0/9N5JcA6xJcirwU+BEgKq6Ncka4Da6u1Oe3g45AZxG91ew9gAubw+AC4GLk6yj2zNYui02TpLUvykDobo/GrJonPovgKMnaLOC7qZgY+trgaedf6iqR2mBIkkaDG9dIUkCDARJUmMgSJIA/4SmdmE54+sDed+r3+wfG9POyT0ESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElq5uQX0/r5wtLKRVs4akBfbJoNdc6bB90FSTs49xAkSYCBIElqDARJEjBHzyFIs+mG9Q8P7PyT54q0NdxDkCQBBoIkqTEQJEmAgSBJajyprFk3qL9cJml63EOQJAEGgiSpMRAkSYCBIElqDARJEtBHICQ5IMnVSW5PcmuS97f6R5L8a5Kb2uONPW3OTrIuyR1JjumpH57k5jbvvCRp9d2TfKnVr0ty4CxsqyRpEv3sIWwGzqiqVwBHAqcnObjNO7eqDmuPbwK0eUuBQ4AlwPlJdmvLXwAsBxa2x5JWPxV4sKpeApwLfGzrN02SNB1TBkJVbaiqH7TpjcDtwP6TNDkOuKSqHquqO4F1wBFJ5gN7V9U1VVXA54Dje9qsbtNfBo4e3XuQJG0f0/piWjuUswi4Dngt8L4kpwBr6fYiHqQLi2t7mq1vtcfb9Ng67flugKranORhYF/g/jHvv5xuD4OhoSFGRkam0/1fW7loy5TLLNizv+V2FjMdq4ls2rSp73XuSuPYj0H+7Gzrf+fZMJ2fnblokOPTdyAkeQ7wFeADVfXLJBcAfw5Uez4HeBcw3m/2NUmdKeY9WahaBawCWLx4cQ0PD/fb/afo59bEKxdt4cwbd51z7nXy8DZd38jICP2O/670p0j7McifnW397zwbpvOzMxcNcnz6+qlN8ky6MPh8VX0VoKruraonqmoL8GngiLb4euCAnuYLgHtafcE49ae0STIPeB7wwEw2SJI0M/1cZRTgQuD2qvp4T31+z2JvBW5p05cBS9uVQwfRnTy+vqo2ABuTHNnWeQpwaU+bZW36BOCqdp5BkrSd9HPI6LXA24Gbk9zUah8CTkpyGN2hnbuAdwNU1a1J1gC30V2hdHpVPdHanQZcBOwBXN4e0AXOxUnW0e0ZLN2ajZIkTd+UgVBV32P8Y/zfnKTNCmDFOPW1wKHj1B8FTpyqL5Kk2bPrnDWVJG0V/x7CHLGt/ybBykVb5tzVQ9Kuzj0ESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRLQRyAkOSDJ1UluT3Jrkve3+j5Jrkjy4/b8gp42ZydZl+SOJMf01A9PcnObd16StPruSb7U6tclOXAWtlWSNIl+9hA2A2dU1SuAI4HTkxwMnAVcWVULgSvba9q8pcAhwBLg/CS7tXVdACwHFrbHklY/FXiwql4CnAt8bBtsmyRpGuZNtUBVbQA2tOmNSW4H9geOA4bbYquBEeCDrX5JVT0G3JlkHXBEkruAvavqGoAknwOOBy5vbT7S1vVl4JNJUlW11VsozSE54+sDed86580DeV9tW1MGQq92KGcRcB0w1MKCqtqQ5EVtsf2Ba3uarW+1x9v02Ppom7vbujYneRjYF7h/zPsvp9vDYGhoiJGRkel0/9dWLtoy5TIL9uxvubnK8ZnYXByb6fxf3LRp04z/784FgxyfvgMhyXOArwAfqKpftsP/4y46Tq0mqU/W5qmFqlXAKoDFixfX8PDwFL0e31F9/Ba1ctEWzrzRc+4TcXwmNhfHpk4e7nvZkZERZvp/dy4Y5Pj09VOb5Jl0YfD5qvpqK9+bZH6bPx+4r9XXAwf0NF8A3NPqC8apP6VNknnA84AHprsxkqSZ6+cqowAXArdX1cd7Zl0GLGvTy4BLe+pL25VDB9GdPL6+HV7amOTIts5TxrQZXdcJwFWeP5Ck7aufQ0avBd4O3Jzkplb7EPAXwJokpwI/BU4EqKpbk6wBbqO7Qun0qnqitTsNuAjYg+5k8uWtfiFwcTsB/QDdVUqSpO2on6uMvsf4x/gBjp6gzQpgxTj1tcCh49QfpQWKJGkw5taZL0nShAwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSmikDIclnktyX5Jae2keS/GuSm9rjjT3zzk6yLskdSY7pqR+e5OY277wkafXdk3yp1a9LcuA23kZJUh/62UO4CFgyTv3cqjqsPb4JkORgYClwSGtzfpLd2vIXAMuBhe0xus5TgQer6iXAucDHZrgtkqStMGUgVNV3gAf6XN9xwCVV9VhV3QmsA45IMh/Yu6quqaoCPgcc39NmdZv+MnD06N6DJGn7mbcVbd+X5BRgLXBGVT0I7A9c27PM+lZ7vE2PrdOe7waoqs1JHgb2Be4f+4ZJltPtZTA0NMTIyMiMOr5y0ZYpl1mwZ3/LzVWOz8Tm4thM5//ipk2bZvx/dy4Y5PjMNBAuAP4cqPZ8DvAuYLzf7GuSOlPMe2qxahWwCmDx4sU1PDw8rU6POuqMr0+5zMpFWzjzRs+5T8TxmdhcHJs6ebjvZUdGRpjp/925YJDjM6Of2qq6t6qeqKotwKeBI9qs9cABPYsuAO5p9QXj1J/SJsk84Hn0f4hKkrSNzCgQ2jmBUW8FRq9AugxY2q4cOoju5PH1VbUB2JjkyHZ+4BTg0p42y9r0CcBV7TyDJGk7mvKQUZIvAsPAfknWAx8GhpMcRndo5y7g3QBVdWuSNcBtwGbg9Kp6oq3qNLorlvYALm8PgAuBi5Oso9szWLoNtkuSNE1TBkJVnTRO+cJJll8BrBinvhY4dJz6o8CJU/VDkjS75taZL0nShAwESRJgIEiSGgNBkgQYCJKkZmtuXSFJAKSPb/+PWrloS193C+hXnfPmbbauuc49BEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSp8fbXknZq07n19ra0K9522z0ESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpmTIQknwmyX1Jbump7ZPkiiQ/bs8v6Jl3dpJ1Se5IckxP/fAkN7d55yVJq++e5Eutfl2SA7fxNkqS+tDPHsJFwJIxtbOAK6tqIXBle02Sg4GlwCGtzflJdmttLgCWAwvbY3SdpwIPVtVLgHOBj810YyRJMzdlIFTVd4AHxpSPA1a36dXA8T31S6rqsaq6E1gHHJFkPrB3VV1TVQV8bkyb0XV9GTh6dO9BkrT9zPSbykNVtQGgqjYkeVGr7w9c27Pc+lZ7vE2PrY+2ubuta3OSh4F9gfvHvmmS5XR7GQwNDTEyMjKjzq9ctGXKZRbs2d9yc5XjMzHHZnK7yvjM9PNnKps2bZq1dU9lW9+6Yrzf7GuS+mRtnl6sWgWsAli8eHENDw/PoItwVB9fdV+5aAtn3ug594k4PhNzbCa3y4zPjRtnZbUrF23hzKsnX/ds3TZjpv8q97bDQLTn+1p9PXBAz3ILgHtafcE49ae0STIPeB5PP0QlSZplMw2Ey4BlbXoZcGlPfWm7cuggupPH17fDSxuTHNnOD5wyps3ouk4ArmrnGSRJ29GUh4ySfBEYBvZLsh74MPAXwJokpwI/BU4EqKpbk6wBbgM2A6dX1RNtVafRXbG0B3B5ewBcCFycZB3dnsHSbbJlkqRpmTIQquqkCWYdPcHyK4AV49TXAoeOU3+UFiiSpMHZBc7sSJK2BQNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkZqsCIcldSW5OclOSta22T5Irkvy4Pb+gZ/mzk6xLckeSY3rqh7f1rEtyXpJsTb8kSdO3LfYQjqqqw6pqcXt9FnBlVS0ErmyvSXIwsBQ4BFgCnJ9kt9bmAmA5sLA9lmyDfkmSpmE2DhkdB6xu06uB43vql1TVY1V1J7AOOCLJfGDvqrqmqgr4XE8bSdJ2Mm8r2xfwrSQFfKqqVgFDVbUBoKo2JHlRW3Z/4Nqetutb7fE2Pbb+NEmW0+1JMDQ0xMjIyIw6vXLRlimXWbBnf8vNVY7PxBybyTk+k+tnfGb62TeVrQ2E11bVPe1D/4okP5pk2fHOC9Qk9acXu8BZBbB48eIaHh6eZnc7R53x9SmXWbloC2fe6Dn3iTg+E3NsJuf4TK6f8amTh2flvbfqX6Wq7mnP9wFfA44A7m2HgWjP97XF1wMH9DRfANzT6gvGqUuStqMZB0KSvZI8d3Qa+APgFuAyYFlbbBlwaZu+DFiaZPckB9GdPL6+HV7amOTIdnXRKT1tJEnbydYcMhoCvtauEJ0HfKGq/j7J94E1SU4FfgqcCFBVtyZZA9wGbAZOr6on2rpOAy4C9gAubw9J0nY040Coqp8Avz1O/RfA0RO0WQGsGKe+Fjh0pn2RJG09z+xIkgADQZLUGAiSJMBAkCQ1BoIkCTAQJEmNgSBJAgwESVJjIEiSAANBktQYCJIkwECQJDUGgiQJMBAkSY2BIEkCDARJUmMgSJIAA0GS1BgIkiTAQJAkNQaCJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAiSJGAHCoQkS5LckWRdkrMG3R9Jmmt2iEBIshvwV8AbgIOBk5IcPNheSdLcskMEAnAEsK6qflJV/wZcAhw34D5J0pySqhp0H0hyArCkqv5ze/124NVV9b4xyy0HlreXLwPumMVu7QfcP4vr39k5PhNzbCbn+ExutsfnxVX1wvFmzJvFN52OjFN7WlJV1Spg1ex3B5KsrarF2+O9dkaOz8Qcm8k5PpMb5PjsKIeM1gMH9LxeANwzoL5I0py0owTC94GFSQ5K8ixgKXDZgPskSXPKDnHIqKo2J3kf8A/AbsBnqurWAXdruxya2ok5PhNzbCbn+ExuYOOzQ5xUliQN3o5yyEiSNGAGgiQJMBCeJslfJvlRkh8m+VqS5/fMO7vdWuOOJMcMsJsD4y1GnirJAUmuTnJ7kluTvL/V90lyRZIft+cXDLqvg5JktyQ3JvlGe+3YNEmen+TL7TPn9iSvGeT4GAhPdwVwaFW9Evhn4GyAdiuNpcAhwBLg/HbLjTnDW4yMazNwRlW9AjgSOL2NyVnAlVW1ELiyvZ6r3g/c3vPasXnS/wD+vqpeDvw23TgNbHwMhDGq6ltVtbm9vJbuOxHQ3Urjkqp6rKruBNbR3XJjLvEWI2NU1Yaq+kGb3kj3H3p/unFZ3RZbDRw/kA4OWJIFwJuAv+kpOzZAkr2B3wMuBKiqf6uqhxjg+BgIk3sXcHmb3h+4u2fe+labSxyDSSQ5EFgEXAcMVdUG6EIDeNEAuzZInwD+K7Clp+bYdH4L+Dnw2XZI7W+S7MUAx2dOBkKSf0xyyziP43qW+VO6wwGfHy2Ns6q5ds2uYzCBJM8BvgJ8oKp+Oej+7AiSHAvcV1U3DLovO6h5wKuAC6pqEfAIAz58tkN8MW17q6rXTzY/yTLgWODoevKLGt5ewzEYV5Jn0oXB56vqq618b5L5VbUhyXzgvsH1cGBeC7wlyRuBZwN7J/lbHJtR64H1VXVde/1lukAY2PjMyT2EySRZAnwQeEtV/apn1mXA0iS7JzkIWAhcP4g+DpC3GBkjSeiOAd9eVR/vmXUZsKxNLwMu3d59G7SqOruqFlTVgXQ/K1dV1ck4NgBU1c+Au5O8rJWOBm5jgOPjN5XHSLIO2B34RStdW1XvafP+lO68wma6QwOXj7+WXVf7be8TPHmLkRWD7dFgJXkd8F3gZp48Tv4huvMIa4DfBH4KnFhVDwykkzuAJMPAmVV1bJJ9cWwASHIY3Qn3ZwE/Ad5J94v6QMbHQJAkAR4ykiQ1BoIkCTAQJEmNgSBJAgwESVJjIEiSAANBktT8f5Newcb1QmLEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Querying the database and limiting the sample to just women who have a birth year that isn't missing\n",
    "\n",
    "con = sqlite3.connect(r'C:\\Users\\devin\\Documents\\School\\ECON6060FinalProject\\ProjectData.sqlite')\n",
    "restricted = pd.read_sql_query(\"\"\"SELECT * \n",
    "                               FROM initial_data \n",
    "                               WHERE sex = 'Female' AND birthyr < 2001\"\"\", con)\n",
    "con.close()\n",
    "\n",
    "# Creating dummy variables for region of residence at time of survey\n",
    "restricted['Northeast'] = np.where(restricted['region'] == 'Northeast', 1, 0)\n",
    "\n",
    "restricted['Midwest'] = np.where(restricted['region'] == 'North Central/Midwest', 1, 0)\n",
    "\n",
    "restricted['South'] = np.where(restricted['region'] == 'South', 1, 0)\n",
    "\n",
    "restricted['West'] = np.where(restricted['region'] == 'West', 1, 0)\n",
    "\n",
    "# Checking if the dummies are mutually exclusive and exhaustive\n",
    "\n",
    "assert (restricted['Northeast'] + restricted['Midwest'] + restricted['South'] + restricted['West']).eq(1).all() == 1\n",
    "\n",
    "# Creating the age variable which will be used as the running variable for the regression discontinuity analysis\n",
    "# It is calculated to be the age of the individual in July of 1978\n",
    "\n",
    "restricted.loc[restricted['birthmo'] == 'January', 'birthmo_num'] = 1\n",
    "restricted.loc[restricted['birthmo'] == 'February', 'birthmo_num'] = 2\n",
    "restricted.loc[restricted['birthmo'] == 'March', 'birthmo_num'] = 3\n",
    "restricted.loc[restricted['birthmo'] == 'April', 'birthmo_num'] = 4\n",
    "restricted.loc[restricted['birthmo'] == 'May', 'birthmo_num'] = 5\n",
    "restricted.loc[restricted['birthmo'] == 'June', 'birthmo_num'] = 6\n",
    "restricted.loc[restricted['birthmo'] == 'July', 'birthmo_num'] = 7\n",
    "restricted.loc[restricted['birthmo'] == 'August', 'birthmo_num'] = 8\n",
    "restricted.loc[restricted['birthmo'] == 'September', 'birthmo_num'] = 9\n",
    "restricted.loc[restricted['birthmo'] == 'October', 'birthmo_num'] = 10\n",
    "restricted.loc[restricted['birthmo'] == 'November', 'birthmo_num'] = 11\n",
    "restricted.loc[restricted['birthmo'] == 'December', 'birthmo_num'] = 12\n",
    "\n",
    "restricted['Age_at_Deadline'] = (1978 + (7/12)) - (restricted['birthyr'] + (restricted['birthmo_num'] / 12))\n",
    "\n",
    "# Restricting the sample to only those with a birth month\n",
    "    \n",
    "restricted = restricted[restricted.birthmo != 'Unknown-refused']\n",
    "restricted = restricted[restricted.birthmo != 'Unknown-not ascertained']\n",
    "restricted = restricted[restricted.birthmo != \"Unknown-don't know\"]\n",
    "\n",
    "hist = restricted.hist('Age_at_Deadline')\n",
    "\n",
    "# Creating dummy variables for race and ethnicity\n",
    "\n",
    "restricted['Hispanic'] = np.where(restricted['hispyn'] == 'Yes, of Hispanic ethnicity', 1, 0)\n",
    "\n",
    "restricted['White'] = np.where(((restricted['racenew'] == 'White') & (restricted['Hispanic'] == 0)), 1, 0)\n",
    "\n",
    "restricted['Black'] = np.where(((restricted['racenew'] == 'Black/African American') & (restricted['Hispanic'] == 0)), 1, 0)\n",
    "\n",
    "restricted['Other'] = np.where(((restricted['Hispanic'] == 0) & (restricted['White'] == 0) & (restricted['Black'] == 0)), 1, 0)\n",
    "\n",
    "# Checking if the dummies are mutually exclusive and exhaustive\n",
    "\n",
    "assert (restricted['Hispanic'] + restricted['White'] + restricted['Black'] + restricted['Other']).eq(1).all() == 1\n",
    "\n",
    "# Creating a dummy variable for if they were born in the US or not\n",
    "\n",
    "restricted['US_born_dummy'] = np.where(restricted['usborn'] == 'Yes, born in U.S.', 1 , 0)\n",
    "\n",
    "# Creating a dummy variable and a continuous variable for height\n",
    "\n",
    "restricted['Exceptional_Height'] = np.where(restricted['height'] == \\\n",
    "                                            'Exceptionally short or tall (2019-forward: or sex is unknown)', 1, 0)\n",
    "\n",
    "# Creating the lifestyle variables\n",
    "    \n",
    "    # Drinking\n",
    "    \n",
    "restricted['Drinker'] = np.where(restricted['alcstat1'] == 'Current drinker (1+ drinks past year)', 1, 0)\n",
    "\n",
    "    # Smoking\n",
    "    \n",
    "restricted['Smoker'] = np.where(((restricted['smokestatus2'] == 'Current every day smoker') | \\\n",
    "                                 (restricted['smokestatus2'] == 'Current some day smoker')), 1, 0)\n",
    "\n",
    "    # Exercise\n",
    "    \n",
    "        # Values of 997, 998, and 999 indicated unknown values, so these observations are dropped\n",
    "    \n",
    "restricted = restricted[restricted.mod10dmin < 997]\n",
    "\n",
    "restricted = restricted[restricted.vig10dmin < 997]\n",
    "    \n",
    "restricted['moderate_exercise'] = restricted['mod10dmin'] * restricted['mod10fwk']\n",
    "\n",
    "restricted['vigorous_exercise'] = restricted['vig10dmin'] * restricted['vig10fwk']\n",
    "\n",
    "        # Values of 94 (less than once per week), 95 (never), and 96 (unable to do this activity) are treated as zeroes\n",
    "        # I'm creating both a dummy variable (which takes a value of 1 if the respondent does any strength-increasing\n",
    "        # exercises, and 0 if not) and a variable which contains the number of times the respondent does\n",
    "        # strength-increasing exercises per week\n",
    "\n",
    "restricted['strong_dummy'] = np.where(((restricted['strongfwk'] < 94) & (restricted['strongfwk'] > 0)), 1, 0)\n",
    "\n",
    "restricted['strong_freq'] = restricted['strongfwk']\n",
    "restricted.loc[restricted['strong_freq'] > 93, 'strong_freq'] = 0\n",
    "\n",
    "# Rescaling the weights due to pooling 5 years of data\n",
    "\n",
    "restricted['perweight'] = restricted['perweight'] / 5\n",
    "\n",
    "# Putting the dataframe into the SQLite database\n",
    "\n",
    "con = sqlite3.connect(r'C:\\Users\\devin\\Documents\\School\\ECON6060FinalProject\\ProjectData.sqlite')\n",
    "restricted.to_sql(\"restricted\", con, if_exists=\"replace\")\n",
    "con.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Calculating the difference in means for \"baseline\" characteristics and running the regression discontinuity**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal bandwidth :  2.0677778526186157\n",
      "                            WLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                Drinker   R-squared:                       0.147\n",
      "Model:                            WLS   Adj. R-squared:                  0.090\n",
      "Method:                 Least Squares   F-statistic:                     2.230\n",
      "Date:                Sun, 09 May 2021   Prob (F-statistic):             0.0977\n",
      "Time:                        22:07:02   Log-Likelihood:                 102.80\n",
      "No. Observations:                  49   AIC:                            -197.6\n",
      "Df Residuals:                      45   BIC:                            -190.0\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:                  HC1                                         \n",
      "===========================================================================================\n",
      "                              coef    std err          z      P>|z|      [0.025      0.975]\n",
      "-------------------------------------------------------------------------------------------\n",
      "Intercept                   0.3087      0.011     28.285      0.000       0.287       0.330\n",
      "Age_at_Deadline            -0.0059      0.010     -0.587      0.557      -0.026       0.014\n",
      "treated                     0.0322      0.018      1.778      0.075      -0.003       0.068\n",
      "Age_at_Deadline:treated     0.0242      0.015      1.577      0.115      -0.006       0.054\n",
      "==============================================================================\n",
      "Omnibus:                        1.516   Durbin-Watson:                   2.198\n",
      "Prob(Omnibus):                  0.469   Jarque-Bera (JB):                1.335\n",
      "Skew:                          -0.250   Prob(JB):                        0.513\n",
      "Kurtosis:                       2.365   Cond. No.                         7.56\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors are heteroscedasticity robust (HC1)\n",
      "                            WLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                 Smoker   R-squared:                       0.062\n",
      "Model:                            WLS   Adj. R-squared:                 -0.001\n",
      "Method:                 Least Squares   F-statistic:                     1.690\n",
      "Date:                Sun, 09 May 2021   Prob (F-statistic):              0.183\n",
      "Time:                        22:07:02   Log-Likelihood:                 122.37\n",
      "No. Observations:                  49   AIC:                            -236.7\n",
      "Df Residuals:                      45   BIC:                            -229.2\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:                  HC1                                         \n",
      "===========================================================================================\n",
      "                              coef    std err          z      P>|z|      [0.025      0.975]\n",
      "-------------------------------------------------------------------------------------------\n",
      "Intercept                   0.1266      0.008     15.105      0.000       0.110       0.143\n",
      "Age_at_Deadline             0.0045      0.008      0.589      0.556      -0.011       0.020\n",
      "treated                     0.0194      0.011      1.788      0.074      -0.002       0.041\n",
      "Age_at_Deadline:treated     0.0063      0.009      0.680      0.496      -0.012       0.024\n",
      "==============================================================================\n",
      "Omnibus:                        0.863   Durbin-Watson:                   2.259\n",
      "Prob(Omnibus):                  0.650   Jarque-Bera (JB):                0.895\n",
      "Skew:                          -0.173   Prob(JB):                        0.639\n",
      "Kurtosis:                       2.435   Cond. No.                         7.56\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors are heteroscedasticity robust (HC1)\n",
      "                            WLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:      moderate_exercise   R-squared:                       0.053\n",
      "Model:                            WLS   Adj. R-squared:                 -0.010\n",
      "Method:                 Least Squares   F-statistic:                    0.7299\n",
      "Date:                Sun, 09 May 2021   Prob (F-statistic):              0.540\n",
      "Time:                        22:07:02   Log-Likelihood:                -255.23\n",
      "No. Observations:                  49   AIC:                             518.5\n",
      "Df Residuals:                      45   BIC:                             526.0\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:                  HC1                                         \n",
      "===========================================================================================\n",
      "                              coef    std err          z      P>|z|      [0.025      0.975]\n",
      "-------------------------------------------------------------------------------------------\n",
      "Intercept                  73.1605     15.813      4.626      0.000      42.167     104.154\n",
      "Age_at_Deadline            21.9961     19.312      1.139      0.255     -15.854      59.847\n",
      "treated                    25.7833     20.182      1.278      0.201     -13.773      65.339\n",
      "Age_at_Deadline:treated   -13.0774     21.520     -0.608      0.543     -55.255      29.100\n",
      "==============================================================================\n",
      "Omnibus:                       12.399   Durbin-Watson:                   2.486\n",
      "Prob(Omnibus):                  0.002   Jarque-Bera (JB):               12.777\n",
      "Skew:                           1.067   Prob(JB):                      0.00168\n",
      "Kurtosis:                       4.307   Cond. No.                         7.56\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors are heteroscedasticity robust (HC1)\n",
      "                            WLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:      vigorous_exercise   R-squared:                       0.028\n",
      "Model:                            WLS   Adj. R-squared:                 -0.036\n",
      "Method:                 Least Squares   F-statistic:                    0.6623\n",
      "Date:                Sun, 09 May 2021   Prob (F-statistic):              0.580\n",
      "Time:                        22:07:03   Log-Likelihood:                -260.34\n",
      "No. Observations:                  49   AIC:                             528.7\n",
      "Df Residuals:                      45   BIC:                             536.3\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:                  HC1                                         \n",
      "===========================================================================================\n",
      "                              coef    std err          z      P>|z|      [0.025      0.975]\n",
      "-------------------------------------------------------------------------------------------\n",
      "Intercept                 102.7744     29.982      3.428      0.001      44.011     161.538\n",
      "Age_at_Deadline            -3.3660     21.029     -0.160      0.873     -44.582      37.850\n",
      "treated                   -12.2790     33.846     -0.363      0.717     -78.617      54.059\n",
      "Age_at_Deadline:treated     9.8720     24.039      0.411      0.681     -37.244      56.987\n",
      "==============================================================================\n",
      "Omnibus:                       28.907   Durbin-Watson:                   2.125\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               62.103\n",
      "Skew:                           1.689   Prob(JB):                     3.27e-14\n",
      "Kurtosis:                       7.359   Cond. No.                         7.56\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors are heteroscedasticity robust (HC1)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            WLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:           strong_dummy   R-squared:                       0.006\n",
      "Model:                            WLS   Adj. R-squared:                 -0.060\n",
      "Method:                 Least Squares   F-statistic:                   0.09127\n",
      "Date:                Sun, 09 May 2021   Prob (F-statistic):              0.964\n",
      "Time:                        22:07:03   Log-Likelihood:                 116.49\n",
      "No. Observations:                  49   AIC:                            -225.0\n",
      "Df Residuals:                      45   BIC:                            -217.4\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:                  HC1                                         \n",
      "===========================================================================================\n",
      "                              coef    std err          z      P>|z|      [0.025      0.975]\n",
      "-------------------------------------------------------------------------------------------\n",
      "Intercept                   0.1106      0.010     10.549      0.000       0.090       0.131\n",
      "Age_at_Deadline            -0.0033      0.008     -0.386      0.699      -0.020       0.013\n",
      "treated                     0.0003      0.014      0.022      0.983      -0.028       0.029\n",
      "Age_at_Deadline:treated     0.0054      0.011      0.496      0.620      -0.016       0.027\n",
      "==============================================================================\n",
      "Omnibus:                        0.161   Durbin-Watson:                   1.748\n",
      "Prob(Omnibus):                  0.923   Jarque-Bera (JB):                0.287\n",
      "Skew:                           0.123   Prob(JB):                        0.866\n",
      "Kurtosis:                       2.716   Cond. No.                         7.56\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors are heteroscedasticity robust (HC1)\n",
      "                            WLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:            strong_freq   R-squared:                       0.047\n",
      "Model:                            WLS   Adj. R-squared:                 -0.017\n",
      "Method:                 Least Squares   F-statistic:                    0.8689\n",
      "Date:                Sun, 09 May 2021   Prob (F-statistic):              0.464\n",
      "Time:                        22:07:03   Log-Likelihood:                 50.287\n",
      "No. Observations:                  49   AIC:                            -92.57\n",
      "Df Residuals:                      45   BIC:                            -85.01\n",
      "Df Model:                           3                                         \n",
      "Covariance Type:                  HC1                                         \n",
      "===========================================================================================\n",
      "                              coef    std err          z      P>|z|      [0.025      0.975]\n",
      "-------------------------------------------------------------------------------------------\n",
      "Intercept                   0.4127      0.047      8.764      0.000       0.320       0.505\n",
      "Age_at_Deadline            -0.0441      0.033     -1.318      0.188      -0.110       0.021\n",
      "treated                    -0.0440      0.061     -0.724      0.469      -0.163       0.075\n",
      "Age_at_Deadline:treated     0.0388      0.044      0.887      0.375      -0.047       0.125\n",
      "==============================================================================\n",
      "Omnibus:                        2.041   Durbin-Watson:                   2.060\n",
      "Prob(Omnibus):                  0.360   Jarque-Bera (JB):                1.170\n",
      "Skew:                           0.304   Prob(JB):                        0.557\n",
      "Kurtosis:                       3.451   Cond. No.                         7.56\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors are heteroscedasticity robust (HC1)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAATD0lEQVR4nO3df7DldV3H8efLBTFYEgjYNnZzaQQTythc0Sati1piSZCFYaZrUZRiZSNO4DRRjTuDtdgvdZQC2REVt8AWK6eIWNESjVUTkcg1DXZYQUOQxaJg3/1xvguH5Z6999x7zp57P/f5mLlzzvn++ry/nz3ndT7nc35sqgpJUlueMOkCJEmjZ7hLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJcWgCSV5Knd9cuTvLm7/rwkt022Oi1GhrvGLsnWJF9PctAEa1jTBegBs9j21UkeTrKr+/tSkncnOX5/1Nqvqj5aVU/b3+1q8TPcNVZJ1gDPAwr4iclWM5SPV9Vy4MnAC4H/BrYl+Z7JliXNjuGucXsVcCNwObB+z8Ik35bkQ0m+keRfkrw5ycf61n93kmuT3JPktiQvm6mhJD+e5NPdMe9I8jt9q2/oLu/tRuM/MJviq+rhqvpiVb0W+AjwyDGTPCfJPye5N8m/JpnqW/fzSW5Ncn+S/0jyy3vV+sYkO5PcmeQX9nFOU0l29N3+cpLzknw2yX1JPpDkSX3rX5LkM11N/5zkGbM5T7XHcNe4vQp4b/f3oiQruuVvBx4Avp1e6PcH/yHAtcD7gKOBlwPvSHLiDG090LV3GPDjwGuSnNGt+6Hu8rCqWl5VH5/DuVxN71UISY4B/gZ4M3AEcB5wVZKjum3vBl4CfCvw88AfJvn+bt9Tu+1/BDiO3iuDYbwMOBU4FngG8OruuN8PXAb8MvBtwLuAayY5HabJMdw1NkmeCzwF2FxV24AvAj+bZBnwU8CFVfXNqvo8sKlv15cAX66qd1fVQ1X1KeAq4Kf31V5Vba2qm6tqd1V9Fng/8MMjPKU76QU5wM8Bf1tVf9u1dy1wE/BjXS1/0434q6o+Avw93RMDvXB+d1V9rqoeoO/VwCz9SVXdWVX3AB8CTuqW/xLwrqr6RPeKYxPwIPCcOZ2tFjXDXeO0Hvj7qvpad/t93bKjgAOAO/q27b/+FODZ3dTCvUnuBV5Bb5Q/UJJnJ7k+yVeT3Af8CnDkaE4FgGOAe/pqPHOvGp8LrOxqeXGSG7tppXvphf6eWr6Dx57vfw5Zx1f6rn8TWN5X0xv2qml1156WmBk/OSDNRZJvoTdCXZZkTxgdRG/KZAXwELAK+Pdu3eq+3e8APlJVPzJks+8D3ga8uKr+J8kf8WigjuLnT38S+Ghfje+pql/ae6NuGuQqelNEW6rq/5L8FZBuk5089ny/cwS17alpQ1VtGNHxtIg5cte4nAE8DJxAb9rgJODp9MLxVfTmr38nycFJvrtbtsdfA8cneWWSA7u/ZyV5+gxtHgrc0wX7ycDP9q37KrAb+K5hTiLJsiTHJvlTYAr43W7VFcBpSV7UbfOk7s3PVcAT6T2RfRV4KMmLgR/tO+xm4NVJTkhyMHDhMDXtw58Bv9K9gkmSQ7o3mQ8d0fG1iBjuGpf19OaVb6+qr+z5ozeyfgXwOnofM/wK8B568+MPAlTV/fTC8Cx689xfAd5CLzD35bXA7yW5H/hteiFKd8xvAhuAf+qmLGaah/6BJLuAbwBb6b0x+qyqurk73h3A6cCb6IX4HcAbgSd09f9a1/7X6T3JXNNXy4eBPwL+EdjeXc5bVd1Eb979bV272+nebNXSE/+zDi0ESd4CfHtVrZ9xY0kzcuSuieg+x/6MbvrgZOBs4IOTrktqheGuSTmU3rz7A/SmLy4Gtsy0U5Jb8ujPAvT/vWKYxpO8c8Bx3jmns5EWGKdlJKlBjtwlqUEL4nPuRx55ZK1Zs2bO+z/wwAMccsghoyuocfbXcOyv4dhfw5lPf23btu1rVXXUdOsWRLivWbOGm266ac77b926lampqdEV1Dj7azj213Dsr+HMp7+SDPx2s9MyktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUoAXxDVVpockbPvTI9Y1rd3NKd7suPm1SJUlDceQuSQ0y3CWpQU7LqHn9Uyz9nGJRyxy5S1KDDHdJapDhLkkNcs5dkoa0GN7HceQuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGuRHISXNaDF89E+PZbg3yAeiJMNd0pxNN5DYuHY3U/u/lAVv0KDr+tMOHUt7hvuYOHqWNEkzhnuSJwE3AAd12/9lVV2Y5AjgA8Aa4MvAy6rq690+FwBnAw8Dv1ZVfzeW6jvbdtz3yP+U06+FIB30JAFtnJ+k8ZjNyP1B4PlVtSvJgcDHknwYeClwXVVdlOR84HzgN5OcAJwFnAh8B/APSY6vqofHdA4DjWr0bMBq1CZ9n/KVZftm/Chk9ezqbh7Y/RVwOrCpW74JOKO7fjpwZVU9WFVfArYDJ4+yaEnSvqWqZt4oWQZsA54KvL2qfjPJvVV1WN82X6+qw5O8Dbixqq7oll8KfLiq/nKvY54DnAOwYsWKZ1555ZVzPom77rmPHd+c/fbPXPXkoY6/bcd9Qx9r0D4Loe1du3axfPnyoeqYhHH34b6O07/PqoN55P41bNvD1jTKNubS/rD3qemsOhhWHDH+c5ik+d6n+j3tsGVzfjyecsop26pq3XTrZvWGajelclKSw4APJvmefWye6Q4xzTEvAS4BWLduXU1NTc2mlGldfMUWzvv07L+PVT83XFvTzefPdKxB+wzafvDL9MHnNde2t27dynz6e38Ztg9HeZz+fTau3f3I/WvYtoetaZRtzKX9Ye9T09m4djc/89Lpj9OK+d6n+l1/2vKxPB6H+rRMVd2bZCtwKnBXkpVVtTPJSuDubrMdwOq+3VYBd46iWGmU9jXvLS12s/m0zFHA/3XB/i3AC4G3ANcA64GLusst3S7XAO9L8lZ6b6geB3xyDLWPnA92Sa2Yzch9JbCpm3d/ArC5qv46yceBzUnOBm4HzgSoqluSbAY+DzwEnDuJT8pI0lI2Y7hX1WeBtdMs/y/gBQP22QBsmHd12qdhX2ns2X7j2t2Pmf9bbB9/82N80sz8hqrmxICVFjbDXRqBVp7sfN+pHUsy3L0DS2qd/1mHJDVoSY7c9ViT/p0TLS3jnsLy/txjuO9nTgktLf57a1IMd2mRGPaJYimNUvV4zrlLUoMMd0lqkNMyaobz29KjDHftk4G5/y3VPm/li2ALheGu/cIHrrR/Ge7SEJbqqLoVS2mQYbhLWtCWUiCPkuGukXJkq5bM9We1FwLDXRO1kB4MWrpavB8a7pIWpRYDeZT8EpMkNchwl6QGGe6S1CDn3OfJeT8tVN43lzZH7pLUIMNdkhrktIykkXNKaPIcuUtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNWjGcE+yOsn1SW5NckuSX++WH5Hk2iRf6C4P79vngiTbk9yW5EXjPAFJ0uPNZuT+EPCGqno68Bzg3CQnAOcD11XVccB13W26dWcBJwKnAu9IsmwcxUuSpjdjuFfVzqr6VHf9fuBW4BjgdGBTt9km4Izu+unAlVX1YFV9CdgOnDziuiVJ+3DAMBsnWQOsBT4BrKiqndB7AkhydLfZMcCNfbvt6JbtfaxzgHMAVqxYwdatW4et/RGrDoaNa3fPef+lxv4ajv01HPtrOLt27ZpX/g0y63BPshy4Cnh9VX0jycBNp1lWj1tQdQlwCcC6detqampqtqU8zsVXbOG8T/ve8GxtXLvb/hqC/TUc+2s415+2nPnk3yCz+hdIciC9YH9vVV3dLb4rycpu/Urg7m75DmB13+6rgDtHU64kaTZm82mZAJcCt1bVW/tWXQOs766vB7b0LT8ryUFJjgWOAz45upIlSTOZzbTMDwKvBG5O8plu2ZuAi4DNSc4GbgfOBKiqW5JsBj5P75M251bVw6MuXJI02IzhXlUfY/p5dIAXDNhnA7BhHnVJkubBdz0kqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkho0Y7gnuSzJ3Uk+17fsiCTXJvlCd3l437oLkmxPcluSF42rcEnSYLMZuV8OnLrXsvOB66rqOOC67jZJTgDOAk7s9nlHkmUjq1aSNCszhntV3QDcs9fi04FN3fVNwBl9y6+sqger6kvAduDk0ZQqSZqtA+a434qq2glQVTuTHN0tPwa4sW+7Hd2yx0lyDnAOwIoVK9i6descS4FVB8PGtbvnvP9SY38Nx/4ajv01nF27ds0r/waZa7gPkmmW1XQbVtUlwCUA69atq6mpqTk3evEVWzjv0743PFsb1+62v4Zgfw3H/hrO9actZz75N8hc/wXuSrISoLu8u1u+A1jdt90q4M65lydJmou5hvs1wPru+npgS9/ys5IclORY4Djgk/MrUZI0rBmnZZK8H5gCjkyyA7gQuAjYnORs4HbgTICquiXJZuDzwEPAuVX18JhqlyQNMGO4V9XLB6x6wYDtNwAb5lOUJGl+fNdDkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBYwv3JKcmuS3J9iTnj6sdSdLjjSXckywD3g68GDgBeHmSE8bRliTp8cY1cj8Z2F5V/1FV/wtcCZw+prYkSXtJVY3+oMlPA6dW1S92t18JPLuqXte3zTnAOd3NpwG3zaPJI4GvzWP/pcb+Go79NRz7azjz6a+nVNVR0604YO717FOmWfaYZ5GqugS4ZCSNJTdV1bpRHGspsL+GY38Nx/4azrj6a1zTMjuA1X23VwF3jqktSdJexhXu/wIcl+TYJE8EzgKuGVNbkqS9jGVapqoeSvI64O+AZcBlVXXLONrqjGR6Zwmxv4Zjfw3H/hrOWPprLG+oSpImy2+oSlKDDHdJatCiCvcklyW5O8nn9lr+q91PHdyS5PcnVd9CM11/JflAks90f19O8pkJlrigDOivk5Lc2PXXTUlOnmSNC82APvu+JB9PcnOSDyX51knWuFAkWZ3k+iS3dln1693yI5Jcm+QL3eXho2hvUYU7cDlwav+CJKfQ+/brM6rqRGDjBOpaqC5nr/6qqp+pqpOq6iTgKuDqCdS1UF3OXv0F/D7wu11//XZ3W4+6nMf32Z8D51fV9wIfBN64v4taoB4C3lBVTweeA5zb/SzL+cB1VXUccF13e94WVbhX1Q3APXstfg1wUVU92G1z934vbIEa0F8AJAnwMuD9+7WoBWxAfxWwZ+T5ZPy+xmMM6LOnATd0168Ffmq/FrVAVdXOqvpUd/1+4FbgGHqD003dZpuAM0bR3qIK9wGOB56X5BNJPpLkWZMuaJF4HnBXVX1h0oUscK8H/iDJHfReFV4w2XIWhc8BP9FdP5PHfqFRQJI1wFrgE8CKqtoJvScA4OhRtNFCuB8AHE7vZc4bgc3dqFT79nIctc/Ga4DfqKrVwG8Al064nsXgF+hNOWwDDgX+d8L1LChJltObEn19VX1jXO20EO47gKur55PAbno/xKMBkhwAvBT4wKRrWQTW8+j7En9B7xdPtQ9V9W9V9aNV9Ux6A4gvTrqmhSLJgfSC/b1Vted+dVeSld36lcBIppZbCPe/Ap4PkOR44In4i3QzeSHwb1W1Y9KFLAJ3Aj/cXX8+4DTWDJIc3V0+Afgt4J2TrWhh6GYULgVuraq39q26ht4ggu5yy0jaW0zfUE3yfmCK3sj8LuBC4D3AZcBJ9F7+nVdV/zihEheU6fqrqi5NcjlwY1X5oOsz4P51G/DH9Kb//gd4bVVtm1SNC82APlsOnNttcjVwQS2moBmTJM8FPgrcTG+GAeBN9ObdNwPfCdwOnFlV034QYqj27HNJak8L0zKSpL0Y7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalB/w+fMmxpFjEa3AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Now that the variables are ready, it's time to run the regression discontinuity analysis\n",
    "\n",
    "# Creating the function which will do the regression discontinuity analysis and create most of the graph\n",
    "\n",
    "def regression_discontinuity(variable):\n",
    "    \n",
    "    # Estimating the RDD model\n",
    "\n",
    "    model = smf.wls(f\"{variable}~Age_at_Deadline * treated\", final_rdd).fit(cov_type = \"HC1\")\n",
    "    \n",
    "    # Printing the results\n",
    "\n",
    "    print(model.summary())\n",
    "    \n",
    "    # Storing the predicted results so they can be plotted\n",
    "\n",
    "    final_rdd['predictions'] = model.fittedvalues\n",
    "    \n",
    "    # Creating two separate dataframes so that two lines can be plotted\n",
    "\n",
    "    final_rdd_treat = final_rdd[final_rdd['Age_at_Deadline'] < 0]\n",
    "    final_rdd_untreat = final_rdd[final_rdd['Age_at_Deadline'] >= 0]\n",
    "    \n",
    "    # Plotting the data points, a dividing line for treated and untreated groups, and the two lines\n",
    "\n",
    "    graph = final_rdd.plot.scatter(x=\"Age_at_Deadline\", y=f\"{variable}\")\n",
    "    plt.plot(final_rdd_treat['Age_at_Deadline'], final_rdd_treat['predictions'], color = \"C2\", label = \"Linear Fit\")\n",
    "    plt.plot(final_rdd_untreat['Age_at_Deadline'], final_rdd_untreat['predictions'], color = \"C2\")\n",
    "    plt.xlabel(\"Centralized Age at Deadline\", fontsize = \"xx-large\")\n",
    "    plt.rc('xtick', labelsize = 16)\n",
    "    plt.rc('ytick', labelsize = 16)\n",
    "    plt.rc('legend', fontsize = 18)\n",
    "    plt.legend()\n",
    "    plt.axvline(x = 0, ymin = 0, color = \"C1\")\n",
    "    plt.gcf().set_size_inches(16, 12)\n",
    "\n",
    "# Creating weighted averages of the necessary variables by Age_at_Deadline\n",
    "\n",
    "con = sqlite3.connect(r'C:\\Users\\devin\\Documents\\School\\ECON6060FinalProject\\ProjectData.sqlite')\n",
    "final_rdd = pd.read_sql_query(\"\"\"SELECT SUM(perweight*Drinker) / SUM(perweight) AS Drinker, \n",
    "                              SUM(perweight*Smoker) / SUM(perweight) AS Smoker,\n",
    "                              SUM(perweight*moderate_exercise) / SUM(perweight) AS moderate_exercise,\n",
    "                              SUM(perweight*vigorous_exercise) / SUM(perweight) AS vigorous_exercise,\n",
    "                              SUM(perweight*strong_dummy) / SUM(perweight) AS strong_dummy,\n",
    "                              SUM(perweight*strong_freq) / SUM(perweight) AS strong_freq, COUNT(*) AS N, Age_at_Deadline\n",
    "                              FROM restricted\n",
    "                              GROUP BY Age_at_Deadline\"\"\", con)\n",
    "con.close()\n",
    "\n",
    "# Calculating the optimal bandwidth\n",
    "\n",
    "optimal_bandwidth = rdd.optimal_bandwidth(final_rdd['strong_freq'], final_rdd['Age_at_Deadline'], 18)\n",
    "print('Optimal bandwidth : ', optimal_bandwidth)\n",
    "\n",
    "# Restricting the sample to those in the optimal bandwidth for \n",
    "\n",
    "final_rdd = final_rdd[final_rdd['Age_at_Deadline'] < (18 + optimal_bandwidth)]\n",
    "final_rdd = final_rdd[final_rdd['Age_at_Deadline'] > (18 - optimal_bandwidth)]\n",
    "\n",
    "# Creating a histogram which shows how many individual observations fall into each age \"bin\"\n",
    "\n",
    "final_means = restricted[restricted['Age_at_Deadline'] < (18 + optimal_bandwidth)]\n",
    "final_means = final_means[final_means['Age_at_Deadline'] > (18 - optimal_bandwidth)]\n",
    "\n",
    "hist2 = final_means.hist('Age_at_Deadline', bins = 49)\n",
    "\n",
    "# Creating a couple dummy variables for treatment/non-treatment status\n",
    "    \n",
    "final_rdd['treated'] = np.where(final_rdd['Age_at_Deadline'] < 18, 1, 0)\n",
    "final_rdd['untreated'] = np.where(final_rdd['Age_at_Deadline'] >= 18, 1 , 0)\n",
    "\n",
    "assert (final_rdd['treated'] + final_rdd['untreated']).eq(1).all() == 1\n",
    "\n",
    "# Normalizing the Age_at_Deadline variable\n",
    "\n",
    "final_rdd['Age_at_Deadline'] = final_rdd['Age_at_Deadline'] - 18\n",
    "\n",
    "# Running the regression discontinuity analysis using the function defined previously, calling it for each variable of interest\n",
    "\n",
    "regression_discontinuity(\"Drinker\")\n",
    "plt.title(\"Prevalence of Women who Currently Drink\", fontsize = \"xx-large\")\n",
    "plt.ylabel(\"Percentage who are Current Drinkers\", fontsize = \"xx-large\")\n",
    "plt.savefig(\"Drinker_rdd.png\")\n",
    "plt.close()\n",
    "\n",
    "regression_discontinuity(\"Smoker\")\n",
    "plt.title(\"Prevalence of Women who Currently Smoke\", fontsize = \"xx-large\")\n",
    "plt.ylabel(\"Percentage who are Current Smokers\", fontsize = \"xx-large\")\n",
    "plt.savefig(\"Smoker_rdd.png\")\n",
    "plt.close()\n",
    "\n",
    "regression_discontinuity(\"moderate_exercise\")\n",
    "plt.title(\"Amount of Moderate Exercise by Women\", fontsize = \"xx-large\")\n",
    "plt.ylabel(\"Minutes of Moderate Exercise per Week\", fontsize = \"xx-large\")\n",
    "plt.savefig(\"moderate_exercise_rdd.png\")\n",
    "plt.close()\n",
    "\n",
    "regression_discontinuity(\"vigorous_exercise\")\n",
    "plt.title(\"Amount of Vigorous Exercise by Women\", fontsize = \"xx-large\")\n",
    "plt.ylabel(\"Minutes of Vigorous Exercise per Week\", fontsize = \"xx-large\")\n",
    "plt.savefig(\"vigorous_exercise_rdd.png\")\n",
    "plt.close()\n",
    "\n",
    "regression_discontinuity(\"strong_dummy\")\n",
    "plt.title(\"Prevalence of Women Who do Strength Improving Activities\", fontsize = \"xx-large\")\n",
    "plt.ylabel(\"Percentage who do Strength Improving Activities at Least Once a Week\", fontsize = \"xx-large\")\n",
    "plt.savefig(\"strong_dummy_rdd.png\")\n",
    "plt.close()\n",
    "\n",
    "regression_discontinuity(\"strong_freq\")\n",
    "plt.title(\"Frequency of Strength Improving Activities by Women\", fontsize = \"xx-large\")\n",
    "plt.ylabel(\"Number of Times per Week\", fontsize = \"xx-large\")\n",
    "plt.savefig(\"strong_freq_rdd.png\")\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}