causal_inference.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "pqsKy--8l_vR"
   },
   "source": [
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/haroldcha/gunviolencedata/blob/master/causal_inference.ipynb)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "CmQVK4ZkCsRM"
   },
   "outputs": [],
   "source": [
    "import numpy as np # linear algebra\n",
    "import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n",
    "import datetime\n",
    "import math\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "8f8hNsBbCsRN"
   },
   "source": [
    "### Using difference-in-difference for causal inference\n",
    "\n",
    "Difference-in-difference is a statistical technique used to study the differential effects in natural experiments. The technique works by comparing a treatment group $T$ and a control group $C$ before and after an intervention. We first look at the slope before and after intervention of the control group. Then, we calculate the slope before and after intervention of the treatment. The difference between these two slopes, $Q$, is our \"normal\" difference between the two groups. That is our difference-in-difference estimate of the intervention effect is\n",
    "\n",
    "$$\n",
    "Q = (T_2 - T_1) - (C_2 - C_1)\n",
    "$$\n",
    "\n",
    "![difference in difference illustration](https://upload.wikimedia.org/wikipedia/commons/d/da/Illustration_of_Difference_in_Differences.png)\n",
    "\n",
    "A major challenge with difference-in-difference is that in application, our calculated $Q$ can include confounding variables that occur before and after the intervention, so our implementation requires a **parallel trend assumption**. That is to say, in order to assume the treatment and control slopes are parallel in the absense of a treatment, we must account for all relevant confounders. This is difficult if we only calculate the naive difference of differences, so we will instead use a regression analysis. \n",
    "\n",
    "Regression models yield the same result when looking at individual-level data. Consider each unit of the natural experiment to have four values:\n",
    "\n",
    "- $y$ is the value of interest\n",
    "- $T$ is a binary value indicating whether they received the treatment\n",
    "- $S$ is whether the unit was measured before or after the intervention\n",
    "- $T \\cdot S$ represents the interaction between time $\\cdot$ treatment\n",
    "\n",
    "See the following OLS model\n",
    "\n",
    "$$\n",
    "y = \\beta_0 + \\beta_1 S + \\beta_2 T + \\beta_3 T \\cdot S \\hspace{0.1cm} [+ \\beta_{i+3} X_{i}] + \\epsilon\n",
    "$$\n",
    "\n",
    "When running our regression analysis, the difference in difference estimate will be the calculated coefficient of the interaction term $T \\cdot S$. This is the same as our naive calculation $Q$ because when calculating the conditional expectation of $y$, the population means are the same as the expected values from the regression.\n",
    "\n",
    "Consider that each optimal coefficient corresponds to a difference in conditional expectations. For example,\n",
    "\n",
    "$$\n",
    "\\hat{\\beta_1} = \\hat{E}(Y | S = 1, T = 0) - \\hat{E}(Y | S = 0, T = 0)\n",
    "$$\n",
    "\n",
    "And in the first expectation term, by calculating the conditional expectation of $y$ given treatment before the intervention, we can obtain over all $y$ in that set the mean in the naive difference-in-difference estimator at the beginning of this section.\n",
    "\n",
    "$$\n",
    "\\hat{E}(Y | S = 1, T = 0) = \\frac{\\hat{E}(y I(S = 1, T = 0))}{\\hat{P}(S = 1, T = 0)}\n",
    "$$\n",
    "\n",
    "$$\n",
    "= \\frac{\\sum_{i=1}^{n}{y_{i,S=1}I(i \\in T = 0)}}{n_{T=0}} = \\bar{y}_{T=0,S=1}\n",
    "$$\n",
    "\n",
    "We can apply these same steps to see that each $\\beta_k$ corresponds to the means of units grouped by treatment and time. So by using a regression, we can easily include each confounder $X_i$ and still get an unbiased estimation of the intervention effects.\n",
    "\n",
    "### Experiment design\n",
    "\n",
    "Using difference-in-difference, we can see if our selected intervention, the 2014 Act 192, caused an increase in injuries + fatalities per incident. Our hypothesis is that if gun advocates can sue Philadelphia municipalities for gun legislation, more violent guns will be available and there will be a higher injury + death rate per gun incident.\n",
    "\n",
    "Our value of interest $y$ is the number of deaths in a given incident. Our treatment value $T$ is a binary value indicating 1 if the incident occurred in South Philadelphia, PA or 0 if in Camden, NJ. We also have a value $I$ to indicate 1 if the violent incident occurred after the intervention, and 0 if it occurred before.\n",
    "\n",
    "Since South Philadelphia and Camden are adjacent, we can exclude confounders such as weather, geography, etc. We'll be controlling for population density, median household income, unemployment rates, and the median age."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "bohxQr6QCsRN"
   },
   "outputs": [],
   "source": [
    "phil_camden = pd.read_csv('https://raw.githubusercontent.com/haroldcha/gunviolencedata/master/philadelphia_camden_gun_violence.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 991
    },
    "id": "euwj8u-xCsRN",
    "outputId": "5bc05bf7-aaca-46cb-e615-977d2a6c81dd"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Unnamed: 0</th>\n",
       "      <th>incident_id</th>\n",
       "      <th>date</th>\n",
       "      <th>state</th>\n",
       "      <th>city_or_county</th>\n",
       "      <th>address</th>\n",
       "      <th>n_killed</th>\n",
       "      <th>n_injured</th>\n",
       "      <th>incident_url</th>\n",
       "      <th>source_url</th>\n",
       "      <th>...</th>\n",
       "      <th>participant_type</th>\n",
       "      <th>sources</th>\n",
       "      <th>state_house_district</th>\n",
       "      <th>state_senate_district</th>\n",
       "      <th>datetime</th>\n",
       "      <th>year</th>\n",
       "      <th>state_pop</th>\n",
       "      <th>city_pop</th>\n",
       "      <th>area (sq. mi)</th>\n",
       "      <th>pop_density</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>288</td>\n",
       "      <td>92423</td>\n",
       "      <td>2014-01-01</td>\n",
       "      <td>Pennsylvania</td>\n",
       "      <td>Philadelphia</td>\n",
       "      <td>34th Street and Haverford Avenue</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/92423</td>\n",
       "      <td>http://www.philly.com/philly/blogs/dncrime/Man...</td>\n",
       "      <td>...</td>\n",
       "      <td>0::Victim||1::Subject-Suspect</td>\n",
       "      <td>http://www.philly.com/philly/blogs/dncrime/Man...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2014-01-01</td>\n",
       "      <td>2014</td>\n",
       "      <td>12788313</td>\n",
       "      <td>1565460.0</td>\n",
       "      <td>46058</td>\n",
       "      <td>277.656715</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>297</td>\n",
       "      <td>92747</td>\n",
       "      <td>2014-01-01</td>\n",
       "      <td>New Jersey</td>\n",
       "      <td>Camden</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/92747</td>\n",
       "      <td>http://www.courierpostonline.com/article/20140...</td>\n",
       "      <td>...</td>\n",
       "      <td>0::Victim||1::Subject-Suspect</td>\n",
       "      <td>http://www.courierpostonline.com/article/20140...</td>\n",
       "      <td>5.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2014-01-01</td>\n",
       "      <td>2014</td>\n",
       "      <td>8864525</td>\n",
       "      <td>75673.0</td>\n",
       "      <td>8722</td>\n",
       "      <td>1016.340862</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>532</td>\n",
       "      <td>92349</td>\n",
       "      <td>2014-01-02</td>\n",
       "      <td>Pennsylvania</td>\n",
       "      <td>Philadelphia</td>\n",
       "      <td>7th Street</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/92349</td>\n",
       "      <td>http://www.philly.com/philly/blogs/dncrime/S-P...</td>\n",
       "      <td>...</td>\n",
       "      <td>0::Victim||1::Victim||2::Subject-Suspect</td>\n",
       "      <td>http://www.philly.com/philly/blogs/dncrime/S-P...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2014-01-02</td>\n",
       "      <td>2014</td>\n",
       "      <td>12788313</td>\n",
       "      <td>1565460.0</td>\n",
       "      <td>46058</td>\n",
       "      <td>277.656715</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>557</td>\n",
       "      <td>220357</td>\n",
       "      <td>2014-01-02</td>\n",
       "      <td>Pennsylvania</td>\n",
       "      <td>Philadelphia</td>\n",
       "      <td>34th and Haverford Avenue</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/220357</td>\n",
       "      <td>http://blog.phillypolice.com/2014/02/wanted-su...</td>\n",
       "      <td>...</td>\n",
       "      <td>0::Victim||1::Subject-Suspect</td>\n",
       "      <td>http://blog.phillypolice.com/2014/02/wanted-su...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2014-01-02</td>\n",
       "      <td>2014</td>\n",
       "      <td>12788313</td>\n",
       "      <td>1565460.0</td>\n",
       "      <td>46058</td>\n",
       "      <td>277.656715</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>604</td>\n",
       "      <td>92727</td>\n",
       "      <td>2014-01-02</td>\n",
       "      <td>New Jersey</td>\n",
       "      <td>Camden</td>\n",
       "      <td>1600 block of Louis Street</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/92727</td>\n",
       "      <td>http://www.nj.com/camden/index.ssf/2014/01/cam...</td>\n",
       "      <td>...</td>\n",
       "      <td>0::Victim</td>\n",
       "      <td>http://www.nj.com/camden/index.ssf/2014/01/cam...</td>\n",
       "      <td>5.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2014-01-02</td>\n",
       "      <td>2014</td>\n",
       "      <td>8864525</td>\n",
       "      <td>75673.0</td>\n",
       "      <td>8722</td>\n",
       "      <td>1016.340862</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>604</th>\n",
       "      <td>235382</td>\n",
       "      <td>1062115</td>\n",
       "      <td>2018-03-04</td>\n",
       "      <td>Pennsylvania</td>\n",
       "      <td>Philadelphia</td>\n",
       "      <td>N 6th St and Spring Garden St</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/106...</td>\n",
       "      <td>http://6abc.com/1-killed-1-injured-in-spring-g...</td>\n",
       "      <td>...</td>\n",
       "      <td>0::Victim||1::Victim</td>\n",
       "      <td>http://philadelphia.cbslocal.com/2018/03/04/po...</td>\n",
       "      <td>181.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2018-03-04</td>\n",
       "      <td>2018</td>\n",
       "      <td>12800922</td>\n",
       "      <td>1583592.0</td>\n",
       "      <td>46058</td>\n",
       "      <td>277.930479</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>605</th>\n",
       "      <td>235559</td>\n",
       "      <td>1063484</td>\n",
       "      <td>2018-03-05</td>\n",
       "      <td>New Jersey</td>\n",
       "      <td>Camden</td>\n",
       "      <td>400 block of Mechanic St</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/106...</td>\n",
       "      <td>http://www.nj.com/camden/index.ssf/2018/03/man...</td>\n",
       "      <td>...</td>\n",
       "      <td>0::Victim</td>\n",
       "      <td>https://www.courierpostonline.com/story/news/c...</td>\n",
       "      <td>5.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2018-03-05</td>\n",
       "      <td>2018</td>\n",
       "      <td>8886025</td>\n",
       "      <td>73725.0</td>\n",
       "      <td>8722</td>\n",
       "      <td>1018.805893</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>606</th>\n",
       "      <td>236316</td>\n",
       "      <td>1066841</td>\n",
       "      <td>2018-03-10</td>\n",
       "      <td>Pennsylvania</td>\n",
       "      <td>Philadelphia</td>\n",
       "      <td>200 block of Lombard St</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/106...</td>\n",
       "      <td>http://6abc.com/parking-lot-attendant-shot-dur...</td>\n",
       "      <td>...</td>\n",
       "      <td>0::Victim</td>\n",
       "      <td>http://6abc.com/parking-lot-attendant-shot-dur...</td>\n",
       "      <td>175.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2018-03-10</td>\n",
       "      <td>2018</td>\n",
       "      <td>12800922</td>\n",
       "      <td>1583592.0</td>\n",
       "      <td>46058</td>\n",
       "      <td>277.930479</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>607</th>\n",
       "      <td>237129</td>\n",
       "      <td>1070055</td>\n",
       "      <td>2018-03-15</td>\n",
       "      <td>Pennsylvania</td>\n",
       "      <td>Philadelphia</td>\n",
       "      <td>700 block of S Columbus Blvd</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/107...</td>\n",
       "      <td>http://6abc.com/2-men-run-from-south-philly-tr...</td>\n",
       "      <td>...</td>\n",
       "      <td>0::Subject-Suspect||1::Subject-Suspect</td>\n",
       "      <td>http://6abc.com/2-men-run-from-south-philly-tr...</td>\n",
       "      <td>175.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2018-03-15</td>\n",
       "      <td>2018</td>\n",
       "      <td>12800922</td>\n",
       "      <td>1583592.0</td>\n",
       "      <td>46058</td>\n",
       "      <td>277.930479</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>608</th>\n",
       "      <td>239385</td>\n",
       "      <td>1080442</td>\n",
       "      <td>2018-03-29</td>\n",
       "      <td>New Jersey</td>\n",
       "      <td>Camden</td>\n",
       "      <td>Kaighn Ave and Haddon Ave</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/108...</td>\n",
       "      <td>https://www.courierpostonline.com/story/news/c...</td>\n",
       "      <td>...</td>\n",
       "      <td>0::Victim||1::Victim||2::Victim||3::Victim</td>\n",
       "      <td>http://philadelphia.cbslocal.com/2018/03/29/ca...</td>\n",
       "      <td>5.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2018-03-29</td>\n",
       "      <td>2018</td>\n",
       "      <td>8886025</td>\n",
       "      <td>73725.0</td>\n",
       "      <td>8722</td>\n",
       "      <td>1018.805893</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>609 rows × 36 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     Unnamed: 0  incident_id        date         state city_or_county  \\\n",
       "0           288        92423  2014-01-01  Pennsylvania   Philadelphia   \n",
       "1           297        92747  2014-01-01    New Jersey         Camden   \n",
       "2           532        92349  2014-01-02  Pennsylvania   Philadelphia   \n",
       "3           557       220357  2014-01-02  Pennsylvania   Philadelphia   \n",
       "4           604        92727  2014-01-02    New Jersey         Camden   \n",
       "..          ...          ...         ...           ...            ...   \n",
       "604      235382      1062115  2018-03-04  Pennsylvania   Philadelphia   \n",
       "605      235559      1063484  2018-03-05    New Jersey         Camden   \n",
       "606      236316      1066841  2018-03-10  Pennsylvania   Philadelphia   \n",
       "607      237129      1070055  2018-03-15  Pennsylvania   Philadelphia   \n",
       "608      239385      1080442  2018-03-29    New Jersey         Camden   \n",
       "\n",
       "                              address  n_killed  n_injured  \\\n",
       "0    34th Street and Haverford Avenue         0          1   \n",
       "1                                 NaN         0          1   \n",
       "2                          7th Street         0          1   \n",
       "3           34th and Haverford Avenue         0          1   \n",
       "4          1600 block of Louis Street         1          0   \n",
       "..                                ...       ...        ...   \n",
       "604     N 6th St and Spring Garden St         1          1   \n",
       "605          400 block of Mechanic St         1          0   \n",
       "606           200 block of Lombard St         0          1   \n",
       "607      700 block of S Columbus Blvd         0          0   \n",
       "608         Kaighn Ave and Haddon Ave         0          4   \n",
       "\n",
       "                                          incident_url  \\\n",
       "0     http://www.gunviolencearchive.org/incident/92423   \n",
       "1     http://www.gunviolencearchive.org/incident/92747   \n",
       "2     http://www.gunviolencearchive.org/incident/92349   \n",
       "3    http://www.gunviolencearchive.org/incident/220357   \n",
       "4     http://www.gunviolencearchive.org/incident/92727   \n",
       "..                                                 ...   \n",
       "604  http://www.gunviolencearchive.org/incident/106...   \n",
       "605  http://www.gunviolencearchive.org/incident/106...   \n",
       "606  http://www.gunviolencearchive.org/incident/106...   \n",
       "607  http://www.gunviolencearchive.org/incident/107...   \n",
       "608  http://www.gunviolencearchive.org/incident/108...   \n",
       "\n",
       "                                            source_url  ...  \\\n",
       "0    http://www.philly.com/philly/blogs/dncrime/Man...  ...   \n",
       "1    http://www.courierpostonline.com/article/20140...  ...   \n",
       "2    http://www.philly.com/philly/blogs/dncrime/S-P...  ...   \n",
       "3    http://blog.phillypolice.com/2014/02/wanted-su...  ...   \n",
       "4    http://www.nj.com/camden/index.ssf/2014/01/cam...  ...   \n",
       "..                                                 ...  ...   \n",
       "604  http://6abc.com/1-killed-1-injured-in-spring-g...  ...   \n",
       "605  http://www.nj.com/camden/index.ssf/2018/03/man...  ...   \n",
       "606  http://6abc.com/parking-lot-attendant-shot-dur...  ...   \n",
       "607  http://6abc.com/2-men-run-from-south-philly-tr...  ...   \n",
       "608  https://www.courierpostonline.com/story/news/c...  ...   \n",
       "\n",
       "                               participant_type  \\\n",
       "0                 0::Victim||1::Subject-Suspect   \n",
       "1                 0::Victim||1::Subject-Suspect   \n",
       "2      0::Victim||1::Victim||2::Subject-Suspect   \n",
       "3                 0::Victim||1::Subject-Suspect   \n",
       "4                                     0::Victim   \n",
       "..                                          ...   \n",
       "604                        0::Victim||1::Victim   \n",
       "605                                   0::Victim   \n",
       "606                                   0::Victim   \n",
       "607      0::Subject-Suspect||1::Subject-Suspect   \n",
       "608  0::Victim||1::Victim||2::Victim||3::Victim   \n",
       "\n",
       "                                               sources state_house_district  \\\n",
       "0    http://www.philly.com/philly/blogs/dncrime/Man...                  NaN   \n",
       "1    http://www.courierpostonline.com/article/20140...                  5.0   \n",
       "2    http://www.philly.com/philly/blogs/dncrime/S-P...                  NaN   \n",
       "3    http://blog.phillypolice.com/2014/02/wanted-su...                  NaN   \n",
       "4    http://www.nj.com/camden/index.ssf/2014/01/cam...                  5.0   \n",
       "..                                                 ...                  ...   \n",
       "604  http://philadelphia.cbslocal.com/2018/03/04/po...                181.0   \n",
       "605  https://www.courierpostonline.com/story/news/c...                  5.0   \n",
       "606  http://6abc.com/parking-lot-attendant-shot-dur...                175.0   \n",
       "607  http://6abc.com/2-men-run-from-south-philly-tr...                175.0   \n",
       "608  http://philadelphia.cbslocal.com/2018/03/29/ca...                  5.0   \n",
       "\n",
       "    state_senate_district    datetime  year state_pop   city_pop  \\\n",
       "0                     NaN  2014-01-01  2014  12788313  1565460.0   \n",
       "1                     5.0  2014-01-01  2014   8864525    75673.0   \n",
       "2                     NaN  2014-01-02  2014  12788313  1565460.0   \n",
       "3                     NaN  2014-01-02  2014  12788313  1565460.0   \n",
       "4                     5.0  2014-01-02  2014   8864525    75673.0   \n",
       "..                    ...         ...   ...       ...        ...   \n",
       "604                   1.0  2018-03-04  2018  12800922  1583592.0   \n",
       "605                   5.0  2018-03-05  2018   8886025    73725.0   \n",
       "606                   1.0  2018-03-10  2018  12800922  1583592.0   \n",
       "607                   1.0  2018-03-15  2018  12800922  1583592.0   \n",
       "608                   5.0  2018-03-29  2018   8886025    73725.0   \n",
       "\n",
       "     area (sq. mi)  pop_density  \n",
       "0            46058   277.656715  \n",
       "1             8722  1016.340862  \n",
       "2            46058   277.656715  \n",
       "3            46058   277.656715  \n",
       "4             8722  1016.340862  \n",
       "..             ...          ...  \n",
       "604          46058   277.930479  \n",
       "605           8722  1018.805893  \n",
       "606          46058   277.930479  \n",
       "607          46058   277.930479  \n",
       "608           8722  1018.805893  \n",
       "\n",
       "[609 rows x 36 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phil_camden"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 991
    },
    "id": "NBAuOI6DCsRO",
    "outputId": "26ba7182-d23d-4770-f513-43df2aecf52b"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Unnamed: 0</th>\n",
       "      <th>incident_id</th>\n",
       "      <th>date</th>\n",
       "      <th>state</th>\n",
       "      <th>city_or_county</th>\n",
       "      <th>address</th>\n",
       "      <th>n_killed</th>\n",
       "      <th>n_injured</th>\n",
       "      <th>incident_url</th>\n",
       "      <th>source_url</th>\n",
       "      <th>...</th>\n",
       "      <th>state_house_district</th>\n",
       "      <th>state_senate_district</th>\n",
       "      <th>datetime</th>\n",
       "      <th>year</th>\n",
       "      <th>state_pop</th>\n",
       "      <th>city_pop</th>\n",
       "      <th>area (sq. mi)</th>\n",
       "      <th>pop_density</th>\n",
       "      <th>median_income</th>\n",
       "      <th>unemployment</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>288</td>\n",
       "      <td>92423</td>\n",
       "      <td>2014-01-01</td>\n",
       "      <td>Pennsylvania</td>\n",
       "      <td>Philadelphia</td>\n",
       "      <td>34th Street and Haverford Avenue</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/92423</td>\n",
       "      <td>http://www.philly.com/philly/blogs/dncrime/Man...</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2014-01-01</td>\n",
       "      <td>2014</td>\n",
       "      <td>12788313</td>\n",
       "      <td>1565460.0</td>\n",
       "      <td>46058</td>\n",
       "      <td>277.656715</td>\n",
       "      <td>45300</td>\n",
       "      <td>0.072</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>297</td>\n",
       "      <td>92747</td>\n",
       "      <td>2014-01-01</td>\n",
       "      <td>New Jersey</td>\n",
       "      <td>Camden</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/92747</td>\n",
       "      <td>http://www.courierpostonline.com/article/20140...</td>\n",
       "      <td>...</td>\n",
       "      <td>5.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2014-01-01</td>\n",
       "      <td>2014</td>\n",
       "      <td>8864525</td>\n",
       "      <td>75673.0</td>\n",
       "      <td>8722</td>\n",
       "      <td>1016.340862</td>\n",
       "      <td>26214</td>\n",
       "      <td>0.093</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>532</td>\n",
       "      <td>92349</td>\n",
       "      <td>2014-01-02</td>\n",
       "      <td>Pennsylvania</td>\n",
       "      <td>Philadelphia</td>\n",
       "      <td>7th Street</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/92349</td>\n",
       "      <td>http://www.philly.com/philly/blogs/dncrime/S-P...</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2014-01-02</td>\n",
       "      <td>2014</td>\n",
       "      <td>12788313</td>\n",
       "      <td>1565460.0</td>\n",
       "      <td>46058</td>\n",
       "      <td>277.656715</td>\n",
       "      <td>45300</td>\n",
       "      <td>0.072</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>557</td>\n",
       "      <td>220357</td>\n",
       "      <td>2014-01-02</td>\n",
       "      <td>Pennsylvania</td>\n",
       "      <td>Philadelphia</td>\n",
       "      <td>34th and Haverford Avenue</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/220357</td>\n",
       "      <td>http://blog.phillypolice.com/2014/02/wanted-su...</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2014-01-02</td>\n",
       "      <td>2014</td>\n",
       "      <td>12788313</td>\n",
       "      <td>1565460.0</td>\n",
       "      <td>46058</td>\n",
       "      <td>277.656715</td>\n",
       "      <td>45300</td>\n",
       "      <td>0.072</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>604</td>\n",
       "      <td>92727</td>\n",
       "      <td>2014-01-02</td>\n",
       "      <td>New Jersey</td>\n",
       "      <td>Camden</td>\n",
       "      <td>1600 block of Louis Street</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/92727</td>\n",
       "      <td>http://www.nj.com/camden/index.ssf/2014/01/cam...</td>\n",
       "      <td>...</td>\n",
       "      <td>5.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2014-01-02</td>\n",
       "      <td>2014</td>\n",
       "      <td>8864525</td>\n",
       "      <td>75673.0</td>\n",
       "      <td>8722</td>\n",
       "      <td>1016.340862</td>\n",
       "      <td>26214</td>\n",
       "      <td>0.093</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>604</th>\n",
       "      <td>235382</td>\n",
       "      <td>1062115</td>\n",
       "      <td>2018-03-04</td>\n",
       "      <td>Pennsylvania</td>\n",
       "      <td>Philadelphia</td>\n",
       "      <td>N 6th St and Spring Garden St</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/106...</td>\n",
       "      <td>http://6abc.com/1-killed-1-injured-in-spring-g...</td>\n",
       "      <td>...</td>\n",
       "      <td>181.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2018-03-04</td>\n",
       "      <td>2018</td>\n",
       "      <td>12800922</td>\n",
       "      <td>1583592.0</td>\n",
       "      <td>46058</td>\n",
       "      <td>277.930479</td>\n",
       "      <td>45300</td>\n",
       "      <td>0.072</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>605</th>\n",
       "      <td>235559</td>\n",
       "      <td>1063484</td>\n",
       "      <td>2018-03-05</td>\n",
       "      <td>New Jersey</td>\n",
       "      <td>Camden</td>\n",
       "      <td>400 block of Mechanic St</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/106...</td>\n",
       "      <td>http://www.nj.com/camden/index.ssf/2018/03/man...</td>\n",
       "      <td>...</td>\n",
       "      <td>5.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2018-03-05</td>\n",
       "      <td>2018</td>\n",
       "      <td>8886025</td>\n",
       "      <td>73725.0</td>\n",
       "      <td>8722</td>\n",
       "      <td>1018.805893</td>\n",
       "      <td>26214</td>\n",
       "      <td>0.093</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>606</th>\n",
       "      <td>236316</td>\n",
       "      <td>1066841</td>\n",
       "      <td>2018-03-10</td>\n",
       "      <td>Pennsylvania</td>\n",
       "      <td>Philadelphia</td>\n",
       "      <td>200 block of Lombard St</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/106...</td>\n",
       "      <td>http://6abc.com/parking-lot-attendant-shot-dur...</td>\n",
       "      <td>...</td>\n",
       "      <td>175.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2018-03-10</td>\n",
       "      <td>2018</td>\n",
       "      <td>12800922</td>\n",
       "      <td>1583592.0</td>\n",
       "      <td>46058</td>\n",
       "      <td>277.930479</td>\n",
       "      <td>45300</td>\n",
       "      <td>0.072</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>607</th>\n",
       "      <td>237129</td>\n",
       "      <td>1070055</td>\n",
       "      <td>2018-03-15</td>\n",
       "      <td>Pennsylvania</td>\n",
       "      <td>Philadelphia</td>\n",
       "      <td>700 block of S Columbus Blvd</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/107...</td>\n",
       "      <td>http://6abc.com/2-men-run-from-south-philly-tr...</td>\n",
       "      <td>...</td>\n",
       "      <td>175.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2018-03-15</td>\n",
       "      <td>2018</td>\n",
       "      <td>12800922</td>\n",
       "      <td>1583592.0</td>\n",
       "      <td>46058</td>\n",
       "      <td>277.930479</td>\n",
       "      <td>45300</td>\n",
       "      <td>0.072</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>608</th>\n",
       "      <td>239385</td>\n",
       "      <td>1080442</td>\n",
       "      <td>2018-03-29</td>\n",
       "      <td>New Jersey</td>\n",
       "      <td>Camden</td>\n",
       "      <td>Kaighn Ave and Haddon Ave</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>http://www.gunviolencearchive.org/incident/108...</td>\n",
       "      <td>https://www.courierpostonline.com/story/news/c...</td>\n",
       "      <td>...</td>\n",
       "      <td>5.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>2018-03-29</td>\n",
       "      <td>2018</td>\n",
       "      <td>8886025</td>\n",
       "      <td>73725.0</td>\n",
       "      <td>8722</td>\n",
       "      <td>1018.805893</td>\n",
       "      <td>26214</td>\n",
       "      <td>0.093</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>609 rows × 38 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     Unnamed: 0  incident_id        date         state city_or_county  \\\n",
       "0           288        92423  2014-01-01  Pennsylvania   Philadelphia   \n",
       "1           297        92747  2014-01-01    New Jersey         Camden   \n",
       "2           532        92349  2014-01-02  Pennsylvania   Philadelphia   \n",
       "3           557       220357  2014-01-02  Pennsylvania   Philadelphia   \n",
       "4           604        92727  2014-01-02    New Jersey         Camden   \n",
       "..          ...          ...         ...           ...            ...   \n",
       "604      235382      1062115  2018-03-04  Pennsylvania   Philadelphia   \n",
       "605      235559      1063484  2018-03-05    New Jersey         Camden   \n",
       "606      236316      1066841  2018-03-10  Pennsylvania   Philadelphia   \n",
       "607      237129      1070055  2018-03-15  Pennsylvania   Philadelphia   \n",
       "608      239385      1080442  2018-03-29    New Jersey         Camden   \n",
       "\n",
       "                              address  n_killed  n_injured  \\\n",
       "0    34th Street and Haverford Avenue         0          1   \n",
       "1                                 NaN         0          1   \n",
       "2                          7th Street         0          1   \n",
       "3           34th and Haverford Avenue         0          1   \n",
       "4          1600 block of Louis Street         1          0   \n",
       "..                                ...       ...        ...   \n",
       "604     N 6th St and Spring Garden St         1          1   \n",
       "605          400 block of Mechanic St         1          0   \n",
       "606           200 block of Lombard St         0          1   \n",
       "607      700 block of S Columbus Blvd         0          0   \n",
       "608         Kaighn Ave and Haddon Ave         0          4   \n",
       "\n",
       "                                          incident_url  \\\n",
       "0     http://www.gunviolencearchive.org/incident/92423   \n",
       "1     http://www.gunviolencearchive.org/incident/92747   \n",
       "2     http://www.gunviolencearchive.org/incident/92349   \n",
       "3    http://www.gunviolencearchive.org/incident/220357   \n",
       "4     http://www.gunviolencearchive.org/incident/92727   \n",
       "..                                                 ...   \n",
       "604  http://www.gunviolencearchive.org/incident/106...   \n",
       "605  http://www.gunviolencearchive.org/incident/106...   \n",
       "606  http://www.gunviolencearchive.org/incident/106...   \n",
       "607  http://www.gunviolencearchive.org/incident/107...   \n",
       "608  http://www.gunviolencearchive.org/incident/108...   \n",
       "\n",
       "                                            source_url  ...  \\\n",
       "0    http://www.philly.com/philly/blogs/dncrime/Man...  ...   \n",
       "1    http://www.courierpostonline.com/article/20140...  ...   \n",
       "2    http://www.philly.com/philly/blogs/dncrime/S-P...  ...   \n",
       "3    http://blog.phillypolice.com/2014/02/wanted-su...  ...   \n",
       "4    http://www.nj.com/camden/index.ssf/2014/01/cam...  ...   \n",
       "..                                                 ...  ...   \n",
       "604  http://6abc.com/1-killed-1-injured-in-spring-g...  ...   \n",
       "605  http://www.nj.com/camden/index.ssf/2018/03/man...  ...   \n",
       "606  http://6abc.com/parking-lot-attendant-shot-dur...  ...   \n",
       "607  http://6abc.com/2-men-run-from-south-philly-tr...  ...   \n",
       "608  https://www.courierpostonline.com/story/news/c...  ...   \n",
       "\n",
       "     state_house_district  state_senate_district    datetime  year state_pop  \\\n",
       "0                     NaN                    NaN  2014-01-01  2014  12788313   \n",
       "1                     5.0                    5.0  2014-01-01  2014   8864525   \n",
       "2                     NaN                    NaN  2014-01-02  2014  12788313   \n",
       "3                     NaN                    NaN  2014-01-02  2014  12788313   \n",
       "4                     5.0                    5.0  2014-01-02  2014   8864525   \n",
       "..                    ...                    ...         ...   ...       ...   \n",
       "604                 181.0                    1.0  2018-03-04  2018  12800922   \n",
       "605                   5.0                    5.0  2018-03-05  2018   8886025   \n",
       "606                 175.0                    1.0  2018-03-10  2018  12800922   \n",
       "607                 175.0                    1.0  2018-03-15  2018  12800922   \n",
       "608                   5.0                    5.0  2018-03-29  2018   8886025   \n",
       "\n",
       "      city_pop area (sq. mi)  pop_density  median_income unemployment  \n",
       "0    1565460.0         46058   277.656715          45300        0.072  \n",
       "1      75673.0          8722  1016.340862          26214        0.093  \n",
       "2    1565460.0         46058   277.656715          45300        0.072  \n",
       "3    1565460.0         46058   277.656715          45300        0.072  \n",
       "4      75673.0          8722  1016.340862          26214        0.093  \n",
       "..         ...           ...          ...            ...          ...  \n",
       "604  1583592.0         46058   277.930479          45300        0.072  \n",
       "605    73725.0          8722  1018.805893          26214        0.093  \n",
       "606  1583592.0         46058   277.930479          45300        0.072  \n",
       "607  1583592.0         46058   277.930479          45300        0.072  \n",
       "608    73725.0          8722  1018.805893          26214        0.093  \n",
       "\n",
       "[609 rows x 38 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phil_camden['median_income'] = 26214\n",
    "phil_camden.loc[phil_camden['city_or_county'] == 'Philadelphia', 'median_income'] = 45300\n",
    "\n",
    "phil_camden['unemployment'] = 0.093\n",
    "phil_camden.loc[phil_camden['city_or_county'] == 'Philadelphia', 'unemployment'] = 0.072\n",
    "\n",
    "phil_camden"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "id": "W7M7RtcSCsRP"
   },
   "outputs": [],
   "source": [
    "def parse_age(age_string):\n",
    "    if type(age_string) != str:\n",
    "        return None\n",
    "    \n",
    "    ages = []\n",
    "    age_split = []\n",
    "    if '||' in age_string:\n",
    "        age_split = age_string.split('||')\n",
    "    elif '|' in age_string:\n",
    "        age_split = age_string.split('|')\n",
    "    else:\n",
    "        age_split = [age_string]\n",
    "    \n",
    "    for age in age_split:\n",
    "        if '::' in age:\n",
    "            ages.append(int(age.split('::')[-1]))\n",
    "        elif ':' in age:\n",
    "            ages.append(int(age.split(':')[-1]))\n",
    "        else:\n",
    "            return None\n",
    "\n",
    "    return np.mean(ages)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "id": "zz3xyHc3CsRP"
   },
   "outputs": [],
   "source": [
    "phil_camden['mean_participant_age'] = phil_camden['participant_age'].apply(parse_age)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 424
    },
    "id": "QJgqQTyECsRQ",
    "outputId": "97bbd04c-7926-48ff-f4ba-de6bb0336c40",
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y</th>\n",
       "      <th>T</th>\n",
       "      <th>I</th>\n",
       "      <th>T*I</th>\n",
       "      <th>median_income</th>\n",
       "      <th>pop_density</th>\n",
       "      <th>unemployment</th>\n",
       "      <th>mean_participant_age</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>45300</td>\n",
       "      <td>277.656715</td>\n",
       "      <td>0.072</td>\n",
       "      <td>18.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>45300</td>\n",
       "      <td>277.656715</td>\n",
       "      <td>0.072</td>\n",
       "      <td>18.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>26214</td>\n",
       "      <td>1016.340862</td>\n",
       "      <td>0.093</td>\n",
       "      <td>39.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>45300</td>\n",
       "      <td>277.656715</td>\n",
       "      <td>0.072</td>\n",
       "      <td>19.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>45300</td>\n",
       "      <td>277.656715</td>\n",
       "      <td>0.072</td>\n",
       "      <td>31.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>602</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>45300</td>\n",
       "      <td>277.930479</td>\n",
       "      <td>0.072</td>\n",
       "      <td>59.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>603</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>45300</td>\n",
       "      <td>277.930479</td>\n",
       "      <td>0.072</td>\n",
       "      <td>17.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>604</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>45300</td>\n",
       "      <td>277.930479</td>\n",
       "      <td>0.072</td>\n",
       "      <td>32.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>605</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>26214</td>\n",
       "      <td>1018.805893</td>\n",
       "      <td>0.093</td>\n",
       "      <td>41.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>606</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>45300</td>\n",
       "      <td>277.930479</td>\n",
       "      <td>0.072</td>\n",
       "      <td>51.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>465 rows × 8 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     y  T  I  T*I  median_income  pop_density  unemployment  \\\n",
       "0    1  1  0    1          45300   277.656715         0.072   \n",
       "3    1  1  0    1          45300   277.656715         0.072   \n",
       "4    1  0  0    0          26214  1016.340862         0.093   \n",
       "5    1  1  0    1          45300   277.656715         0.072   \n",
       "6    1  1  0    1          45300   277.656715         0.072   \n",
       "..  .. .. ..  ...            ...          ...           ...   \n",
       "602  1  1  1    1          45300   277.930479         0.072   \n",
       "603  1  1  1    1          45300   277.930479         0.072   \n",
       "604  2  1  1    1          45300   277.930479         0.072   \n",
       "605  1  0  1    0          26214  1018.805893         0.093   \n",
       "606  1  1  1    1          45300   277.930479         0.072   \n",
       "\n",
       "     mean_participant_age  \n",
       "0                    18.0  \n",
       "3                    18.0  \n",
       "4                    39.0  \n",
       "5                    19.0  \n",
       "6                    31.0  \n",
       "..                    ...  \n",
       "602                  59.0  \n",
       "603                  17.0  \n",
       "604                  32.0  \n",
       "605                  41.0  \n",
       "606                  51.0  \n",
       "\n",
       "[465 rows x 8 columns]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "did_df = phil_camden[['year', 'city_or_county', 'n_killed', 'n_injured', 'median_income', 'pop_density', 'unemployment', 'mean_participant_age']]\n",
    "did_df = did_df.dropna()\n",
    "\n",
    "did_df['T'] = 0\n",
    "did_df.loc[did_df['city_or_county'] == 'Philadelphia', 'T'] = 1\n",
    "\n",
    "did_df['I'] = 0\n",
    "did_df.loc[did_df['year'] > 2014, 'I'] = 1\n",
    "\n",
    "did_df['y'] = did_df['n_killed'] + did_df['n_injured']\n",
    "\n",
    "did_df['T*I'] = did_df['T']\n",
    "\n",
    "did_df = did_df[['y', 'T', 'I', 'T*I', 'median_income', 'pop_density', 'unemployment', 'mean_participant_age']]\n",
    "did_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-BxdHNb-CsRQ"
   },
   "source": [
    "Using OLS, we can find the equivalent value of the difference-in-difference estimate of the policy effect. That is, we can find the mean number of deaths before and after Act 192 in the Philadelphia treatment group, and the mean number of deaths before and after Act 192 in the Camden control group, and calculate the difference between these means. This is because our averages are just expected values from our estimator.\n",
    "\n",
    "The reason we use OLS is to easily include our four confounders as well as calculating errors.\n",
    "\n",
    "To do this, we can estimate the unknown parameters in $y = \\beta_0 + \\beta_1 T + \\beta_2 I + \\beta_3 TI + \\beta_4 X_1 + \\beta_5 X_2 + \\beta_6 X_3 + \\beta_7 X_4$ where $X_i$ is the $i^{th}$ confounder in our dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "zjllI2PGCsRR",
    "outputId": "22769a06-eed6-4126-ffd0-72c1b1ba036d"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "T                       80.069005\n",
       "I                        0.176325\n",
       "T*I                     80.069005\n",
       "median_income           -0.004166\n",
       "pop_density              0.108780\n",
       "unemployment            -7.103068\n",
       "mean_participant_age    -0.001602\n",
       "dtype: float64"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = did_df.loc[:, did_df.columns != 'y']\n",
    "Y = did_df['y']\n",
    "\n",
    "model = sm.OLS(Y, X)\n",
    "results = model.fit()\n",
    "results.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "D7pVwQRrCsRR"
   },
   "source": [
    "We can visualize this result by plotting the mean number of injuries + deaths pre and post intervention, segregated by treatment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 297
    },
    "id": "IMIJg7v8CsRR",
    "outputId": "4145feff-98f0-4d3b-f3fa-c144ed937bda"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAEYCAYAAACkzwftAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAv+klEQVR4nO3de3yV9Znv/c+VMzkBSThIQggEFLBQkQiiqMjB4rRqnaltrfZgdehJZ2tnOrsznU6Z2e0zbseZ5+noWDfbWutux+5narsf2mmnIJ4VrFBUPHOGgAoknJKQrKy1fs8fvztZSVghC1jJnax8369XXmTd6173uhKV9fV3uG5zziEiIiIy1GWFXYCIiIhIOijUiIiISEZQqBEREZGMoFAjIiIiGUGhRkRERDJCTlhvXFFR4WpqasJ6exERERlkNm3adMg5N+ZMXx9aqKmpqWHjxo1hvb2IiIgMMma2+2xer+knERERyQgKNSIiIpIRFGpEREQkI4S2piaZ9vZ26uvraW1tDbuUQa+goICqqipyc3PDLkVERGRQGFShpr6+npKSEmpqajCzsMsZtJxzNDQ0UF9fz+TJk8MuR0REZFAYVNNPra2tlJeXK9D0wcwoLy/XiJaIiEgXgyrUAAo0KdLvSUREpLtBF2pEREREzoRCTYpqamo4dOjQScdXr17N3XffDcDKlSu59957AfjCF77Az3/+8wGtUUREpD9FonGisXjYZfRqUC0UHoquvfZarr322rDLEBEROW2xuONIS4TG5giHmiI0NLfR0BShoamNQ80RGrscO9TUxrHWKI9/5RLmThoddulJKdT0sGvXLpYvX878+fPZvHkz5557Lo8++igA9913H7/61a9ob2/n3//935k+fTqPPPIIGzdu5P777096vXXr1nH//ffzy1/+EoC1a9fygx/8gF/84hcD9jOJiMjw4JzjeFu0M4wcaop0hpSG5oj/agqCS3Mbjc0R4u7k62QZlBXlUV6UT1lRHjMnlFJRnE95UR5jS/IH/gdLUZ+hxsweBj4GHHDOfSjJ8zcB/zV42AR8xTn36tkW9ne/eoM39x8728t0M3NCKd+55vw+z3vnnXf44Q9/yKWXXsoXv/hFHnjgAQAqKir4wx/+wAMPPMC9997LQw891Oe1Fi9ezNe+9jUOHjzImDFj+NGPfsQtt9xy1j+LiIgMD63tsW5h5FBHQOn8s+sIS4RIL9NDpQU5lAfBpKaikLk1o6koyqO82AeX8uK8zuAyqjCP7KyhtyEllZGaR4D7gUd7eX4ncIVz7rCZXQ2sAuanp7xwTJw4kUsvvRSAm2++mX/5l38B4I//+I8BmDt3bsojLWbGZz/7WX7yk59wyy23sH79+s6RHxERGX6isTiNLZHOENIxotLYOc3TfRqoORJLep2C3CzKi/KpKM5jbEkBM8aXUl7sH5cFYaW8yAeV0UW55OdkD/BPOvD6DDXOuWfNrOYUz7/Y5eEGoCoNdaU0otJfem6X7nicn++H3LKzs4lGoylf75ZbbuGaa66hoKCAG264gZwczfqJiGSKeNxxrLU9mOpJPorS9bkjLe1Jr5OTZZ1hpKI4j0llhT6YFOdRHkwFdY6mFOdRmKfPkp7S/Ru5Ffhtb0+a2QpgBUB1dXWa3zp99uzZw/r161mwYAGPPfYYCxcuZPPmzWd8vQkTJjBhwgS++93vsnbt2jRWKiIi6eacoyUS86MmXUZMek71dEwDHW6OEE22MAUYXZjbOWIyfXwp5V1GUTqmfsqL86goyqd0RI56kJ2ltIUaM7sSH2oW9naOc24VfnqKurq65P8GDAIzZszgxz/+MV/60peYNm0aX/nKV7jvvvvO6po33XQTBw8eZObMmWmqUkREUtUWjdEYhJJDTW2J73sJLa3tydelFOfndI6cTCwrZE71qM4FtV1HUcqL8hldmEtOtjqnDCRzru9sEUw//TrZQuHg+dnAL4GrnXPvpvLGdXV1buPGjd2OvfXWW8yYMSOVl/ebXbt28bGPfYzXX389rde9/fbbmTNnDrfeemvarjkYfl8iImGIxR2HW3rs7An+7DrV47cqt3G8NfmSgbycrG4jJh3hpLzHKErHCEtBbuavSwmTmW1yztWd6evPeqTGzKqBXwCfTTXQDDdz586lqKiIf/qnfwq7FBGRQck5x7HWaC87enzPlIYuIyyNLRGS/T+534rs16SUF+fxocqRwWLZxMLZrotpi/M15ZNJUtnS/RiwCKgws3rgO0AugHPuQeBvgXLggeBfjOjZpKyw1dTUpH2UZtOmTWm9nojIUHAiEktM9aTQM6U9lnzmYOSI3M4RkykVxVxUkwgmPUdWRo3IJWsIbkWW9Ehl99ONfTx/G3Bb2ioSEZFBqT0W7xwp6blYtmMUpet25JZetiIX5mV3TvWMLy3g/Aml3bYfd0z1VBTnM7owj7wcrUuR1Gg/mIjIMBWPO46caKex6yhKl+3HPRfTHj2RfCtybrZ1dp4tL85jckVRtzUp5UVdp360FVn6j/7NEhHJEM45miMxvwblFD1TOhq8HW6JEEuyFdkMygrzOkPKjHNKeyym7b4+pbRA61JkcFCoEREZxFrbY0m3H/d2A8JINPlW5JKCnM4QUl1WyJzq0cGalDzKevRMGT1EW+SLKNQMsKeffpp7772XX//612GXIiIhiMbiHG5p774mpSnSYzFtYmSlqS35VuT8nKzO9ScVxXmcN76k186zZUV5w6JFvohCjYjIWXDOcexENKXOs43Nfson2Vbk7I4W+cEC2YllhafsmVKYl60pH5EeFGqSePTRR7n33nsxM2bPns0nP/lJvvvd7xKJRCgvL+enP/0p48aNY+XKlezcuZP33nuPd999l3/+539mw4YN/Pa3v6WyspJf/epX5Obm8p//+Z/ceeedVFRUcOGFF3a+T3NzM3fccQdbtmwhGo2ycuVKrrvuOh555BFWr15NS0sL27dv5/rrr+eee+4J8TciMry0RKInjaL0DC1db0DYW4v8UYW5nYHk3HHFnZ1ne/ZMKS/KY6S2IouctcEban77TXh/S3qvOX4WXH33KU954403+N73vscLL7xARUUFjY2NmBkbNmzAzHjooYe45557Ohvpbd++naeeeoo333yTBQsW8Pjjj3PPPfdw/fXX8x//8R8sX76cP/3TP+XJJ59k6tSpfOpTn+p8r+9973ssXryYhx9+mCNHjjBv3jyWLl0KwCuvvMLmzZvJz8/nvPPO44477mDixInp/X2IDBORaLyzs2zHNE/n3ZCT9Ew50Z58K3JRXnbniEnlqBHMrhzpR1KS9EwZXZRHrlrkiwyowRtqQvLkk0/yiU98goqKCgDKysrYsmULn/rUp3jvvfeIRCJMnjy58/yrr76a3NxcZs2aRSwWY/ny5QDMmjWLXbt28fbbbzN58mSmTZsGwM0338yqVasAWLNmDatXr+bee+8FoLW1lT179gCwZMkSRo4cCcDMmTPZvXu3Qo1IIBZ3HGmJpNR59lBTG8d6a5GfnRWEEh9IaiuKOkNKz54p5UX5jMjTuhSRwWzwhpo+RlT6i3PupHnqO+64g69//etce+21PP3006xcubLzufz8fACysrLIzc3tfG1WVhbRqP+LtLd5b+ccjz/+OOedd1634y+99FLndQGys7M7ryWSiZxzHG+L0titT0rXqZ62bk3fGpsjJJvx8S3y8zp7psycUOqDSbLtyMV5lKhFvkhGGbyhJiRLlizh+uuv56677qK8vJzGxkaOHj1KZWUlAD/+8Y9P63rTp09n586dbN++ndraWh577LHO5z7ykY9w3333cd9992FmbN68mTlz5qT15xEJS2t7rNuUzqFT9ExpaIoQiSXfilxakENFsQ8pNRWFzK0Z3W37cUfn2fKiPEZpK7LIsKZQ08P555/Pt771La644gqys7OZM2cOK1eu5IYbbqCyspKLL76YnTt3pny9goICVq1axUc/+lEqKipYuHBh572lvv3tb3PnnXcye/ZsnHPU1NRoq7cMWtFYnMbOuyKf3Hm2Z8+U5l5a5BfkZnUulh1bUsCM8aWJNSnFeZQVJaZ+RhflaiuyiKTMXLK9hQOgrq7Obdy4sduxt956ixkzZoRSz1Ck35ecjXjccay1/aTOs4d66ZlypCV5i/ycjq3InYtle0z19OiZohb5ItIbM9t0NjfF1t8uIhnCOUdLJHZS59neeqYcbk6+FdkMRo3I7VwsO318aedC2rLivG5TPxVF+ZSO0LoUERkcFGpEBrG2aKzbLp7eeqZ0hJbW9uTrUorzczpHTiaWFTKnelS3GxBWdI6s5DO6MJccbUUWkSFo0IWaZLuP5GRhTRvK2YnFHYdbIt22Hzd2aebW8waEx3vbipyT1W3EZNrYkpN29nR0ni0ryqMgV+tSRCTzDapQU1BQQENDA+Xl5Qo2p+Cco6GhgYKCgrBLGfaccxxrjfayo+fknimNvbTI91uRE4tlZ1eNCnb1dO88WxGElGJtRRYROcmgCjVVVVXU19dz8ODBsEsZ9AoKCqiqqgq7jIx0IhLrcq+evnumtMeSj5qNHJHbOWJSO6aYi2pO7jxbEez2GaUW+SIiZ21QhZrc3Nxu3XpF0qE9Fu8WQk7qmdIUdKYNnmvpZStyYV5257qT8aUFnD+hNGnn2YrifEYX5pGXo3UpIiIDaVCFGpFUxOOOIyfau4+i9OiZ0nUH0NETybci52Zb54hJWVEekyuKuq9J6dIzRVuRRUQGP/0tLaFzztHUFk00cOulZ0rHDQgPt0SI9bIVuawwr3NHz4xzSrstpu16s8Hy4nxKC7QuRUQkkyjUSL9obY+dNGJyqsW0kWjyrcglBTmdIaS6rJA51aM7G7yVFed3Cy2j1SJfRGRYU6iRlERjcQ63tHdfk9JL59mGpghNbcm3IufnZHWuP6kozuO88SWdi2l79kwpK8pTi3wREUmZQs0w5Zzj2Ilot1EU3zMl+WLaIyfak25Fzu5okR8skJ1YVnjSNE/XnimFedma8hERkX6hUJNBWiLRbqMoDc1dpnt6NHhr7KVFPsCowtzOQHLuuGLKp5T32jOltEBbkUVEZHBQqBnEItF4sHg20TOlY7Fsz86zDU0RTrQn34pclJfdOWJSOWoEH64a2eMGhMHISrAuJVct8kVEZAhSqBlAsbjjSEvklJ1nO+/t09TGsd5a5GdndYaQ8qCxW3mXUZTO+/gEj9UiX0REhgOFmrPgnON4WzQIIn13nm1sjpBsxse3yE+MmMycUOqDSbKeKcV5lKhFvoiIyEn6DDVm9jDwMeCAc+5DSZ6fDvwIuBD4lnPu3rRXOYBa22PdwsipOs82NEWIxJJvRS4tyKGi2O/oqakoZG7N6F57pozSVmQREZGzlspIzSPA/cCjvTzfCPwZ8PH0lJRe0Vicxs67Ip/cefZQj2mg5l5a5BfkZnWOnowtKWDG+NLEmpQunWc7goxa5IuIiAysPkONc+5ZM6s5xfMHgANm9tF0Fnam4nHHjf9zQ+cIy5GW5C3yc7KsM4xUFOcxqayw2/bjnj1T1CJfRERkcBvQT2ozWwGsAKiuru6X98jKMorzcxKLZYvyKSvO6zb1U1GUT+kIrUsRERHJJAMaapxzq4BVAHV1dcmbpKTBD79wUX9dWkRERAYpLfwQERGRjKBQIyIiIhkhlS3djwGLgAozqwe+A+QCOOceNLPxwEagFIib2Z3ATOfcsf4qWkRERKSnVHY/3djH8+8DVWmrSEREROQMaPpJREREMoJCjYiIiGQEhRoRERHJCAo1IiIikhEUakRERCQjKNSIiIhIRlCoERERkYygUCMiIiIZQaFGREREMoJCjYiIiGQEhRoRERHJCAo1IiIikhEUakRERCQjKNSIiIhIRlCoERERkYygUCMiIiIZQaFGREREMoJCjYiIiGQEhRoRERHJCAo1IiIikhEUakRERCQjKNSIiIhIRlCoERERkYygUCMiIiIZQaFGREREMoJCjYiIiGQEhRoRERHJCH2GGjN72MwOmNnrvTxvZvYvZrbNzF4zswvTX6aIiIjIqaUyUvMIsPwUz18NTAu+VgA/OPuyRERERE5Pn6HGOfcs0HiKU64DHnXeBmCUmZ2TrgJFREREUpGONTWVwN4uj+uDYyIiIiIDJh2hxpIcc0lPNFthZhvNbOPBgwfT8NYiIiIiXjpCTT0wscvjKmB/shOdc6ucc3XOuboxY8ak4a1FREREvHSEmtXA54JdUBcDR51z76XhuiIiIiIpy+nrBDN7DFgEVJhZPfAdIBfAOfcg8Bvgj4BtQAtwS38VKyIiItKbPkONc+7GPp53wNfSVpGIiIjIGVBHYREREckICjUiIiKSERRqREREJCMo1IiIiEhGUKgRERGRjKBQIyIiIhlBoUZEREQygkKNiIiI9K31GLy5GtpPhF1Jr/psviciIiLDkHNw8G3Yuga2roU96yEehZt+DtOWhV1dUgo1IiIi4rU1wc5nYdtaH2SO7vXHx30ILrkDpi6DifPCrfEUFGpERESGK+egYXswGrMGdr8AsQjkFcOURXD5X/ggM7Iy7EpTolAjIiIynLSfgF3PJ4LM4V3+eMV5MG8FTLsKqhdATl6oZZ4JhRoREZFM17jTTydtW+unl6KtkDMCplyRmFYaPSnsKs+aQo2IiEimibbB7hd9kNm6Bhq2+uNltTD3C36h76SFkFsQapnpplAjIiKSCY7sTSzw3fEMtDdDdj7ULISLbvNBprw27Cr7lUKNiIjIUBRrhz0b/EjMtifgwJv++MhquOBGvzamZiHkFYVb5wBSqBERERkqjr3nA8zWNbDjaWg7Blm5MOkSuOq7PshUnAtmYVcaCoUaERGRwSoWhX0bEzuV3t/ij5dMgPOv9yFmyhWQXxJunYOEQo2IiMhg0nQwMRqz/UloPQKWDdUXw9KVPsiMnTlsR2NORaFGREQkTPE47N+cGI3ZvxlwUDQWpn8Mpi2FKVfCiFFhVzroKdSIiIgMtJZGPwqzdQ1sWwcthwCDqovgym/5nUrjZ0OW7jt9OhRqRERE+ls8Du+/lthyXf8yuDgUlsPUpb75Xe1iKCoPu9IhTaFGRESkP7Qehe1PJTr5Nn3gj0+YA5d/w6+NmTAHsrLDrTODKNSIiIikg3O+V8zWNT7I7NkALgYFI6F2iQ8xU5dA8diwK81YCjUiIiJnqu24797bMa10bJ8/Pn4WLLzTTytVXQTZ+rgdCPoti4iIpMo5OLQ1sVNp94sQb4e8EqhdBIu+6dfIlE4Iu9JhSaFGRETkVCItsOu5RJA5sscfHzMDLv6Kn1aaOB9y8sKtUxRqRERETtKwPbHAd+dzEGuD3EKYsggW3uWnlUZNDLtK6SGlUGNmy4HvA9nAQ865u3s8Pxp4GKgFWoEvOudeT3OtIiIi/aO9FXa/4IPM1jXQuN0fL58GF93q+8ZMuhRy8sOtU06pz1BjZtnAvwLLgHrgZTNb7Zx7s8tpfw284py73symB+cv6Y+CRURE0uLw7sQC353PQnsL5BRAzWUw/8u+k2/ZlLCrlNOQykjNPGCbc24HgJn9DLgO6BpqZgL/AOCce9vMasxsnHPug3QXLCIickaiEdizPuji+wQcfNsfHzUJ5tzs18bULITcEeHWKWcslVBTCezt8rgemN/jnFeBPwaeN7N5wCSgCugWasxsBbACoLq6+gxLFhERSdGx/YkppR1PQ6QJsvP8VNKFn/NBpnyqbg6ZIVIJNcn+Sbsej+8Gvm9mrwBbgM1A9KQXObcKWAVQV1fX8xoiIiJnJxaF+t8nGuB9ECzvLK2CWTf4EDP5csgvDrdO6RephJp6oOsS7ypgf9cTnHPHgFsAzMyAncGXiIhI/zr+gZ9O2rrG35ag7Shk5UD1Alj29z7IjJmu0ZhhIJVQ8zIwzcwmA/uATwOf6XqCmY0CWpxzEeA24Nkg6IiIiKRXPAb7NiWmld57xR8vHg8zr/EhZsoif3sCGVb6DDXOuaiZ3Q78Dr+l+2Hn3Btm9uXg+QeBGcCjZhbDLyC+tR9rFhGR4aa5AbavCxb5roMTjWBZUDUPFn/bB5nxszQaM8yl1KfGOfcb4Dc9jj3Y5fv1wLT0liYiIsNWPO5HYDqmleo3Ag4KK+Dcj/hbEdQuhsKysCuVQUQdhUVEZHA4cdivieno5Nt8EDCovNDfU2naMjhnDmRlhV2pDFIKNSIiEg7n/O6kjp1Ke38PLgYFo/xIzLSrYOoSKKoIu1IZIhRqRERk4LQe8/1iOjr5Hn/PHz/nw3DZ1/09larqICs71DJlaFKoERGR/uMcHHwncYfrPeshHoX8Uqi9MhiNWQol48OuVDKAQo2IiKRXpNnfS6ljWulo0JR+7Pmw4HYfZCbOg+zccOuUjKNQIyIiZ8c5aNgebLdeC7ueh1gEcov8aMzlf+FHY0ZWhV2pZDiFGhEROX3tJ2DXC4lppcNBE/mK82DeCr9TqXoB5OSHW6cMKwo1IiKSmsO7El18dz4H0ROQM8LfS2nB13yQGV0TdpUyjCnUiIhIctE22P1iom/MoXf98dGTYe7nfYiZdCnkjgi3TpGAQo2IiCQcrQ9GY9b6rdftzZCdDzULoe6LfpFveW3YVYokpVAjIjKcxdph70uJnUoH3vTHR1bDhz/tQ8zkyyCvKNw6RVKgUCMiMtwcfz+xNmbH09B2DLJyYNIlcNV3fZCpOFc3h5QhR6FGRCTTxaKwb2MiyLz/mj9ecg6c//FgNOYKKCgNtUyRs6VQIyKSiZoOwvZ1Qe+YddB6BCwbJs6HJd/xQWbc+RqNkYyiUCMikgnicdi/OdEAb98fAAdFY2H6R33zu9orYcTosCsV6TcKNSIiQ1VLI2x/Mthy/QS0HALM3xDyyr/2W67HfxiyssKuVGRAKNSIiAwVzvn1MB07lepfBheHEWV+JGbaVVC7GIrKw65UJBQKNSIig1nrUdj+lJ9S2voENL3vj0+YA5d/A6Yug8oLISs73DpFBgGFGhGRwcQ5OPBWYjRm7waIR6FgpB+FmXaVH5UpHht2pSKDjkKNiEjY2ppg5zOJIHNsnz8+bhZc8mc+yFRdBNn6K1vkVPRfiIjIQHMODm1N7FTa/SLEIpBXArWLYNE3/WhM6YSwKxUZUhRqREQGQqQFdj0fjMasgSO7/fExM2D+l/xozMSLIScv3DpFhjCFGhGR/tK4I9HFd9fzEG2F3ELfvffS/+K3XI+qDrtKkYyhUCMiki7trbD7haBvzFpo2OaPl08N7nC9DKovgdyCcOsUyVAKNSIiZ+PInmA0Zq1f7NveAjkFUHMZzFvh18aU14ZdpciwoFAjInI6ohG/zbpjp9LBt/3xUdVwwU1+bUzNQsgrDLdOkWFIoUZEpC/H9ifWxux4BiLHISsXai6FCz/nG+BVTNPNIUVCllKoMbPlwPeBbOAh59zdPZ4fCfwEqA6uea9z7kdprlVEZGDEolD/+8S00gdb/PHSSpj1J340ZvLlkF8Sbp0i0k2focbMsoF/BZYB9cDLZrbaOfdml9O+BrzpnLvGzMYA75jZT51zkX6pWkQk3ZoO+JtCbl3jbxLZehQsG6oXwNK/80Fm7AyNxogMYqmM1MwDtjnndgCY2c+A64CuocYBJWZmQDHQCETTXKuISPrEY7DvD4kGePs3++PF42DGNT7ETFnkb08gIkNCKqGmEtjb5XE9ML/HOfcDq4H9QAnwKedcvOeFzGwFsAKgulq9GURkgLU0wrZ1QZB5Ak40gmX5WxAs/hsfZMbNgqyssCsVkTOQSqhJNtbqejz+CPAKsBioBdaa2XPOuWPdXuTcKmAVQF1dXc9riIikVzwO77+aWORbvxFwUFjuA8y0Zf4mkYVlYVcqImmQSqipByZ2eVyFH5Hp6hbgbuecA7aZ2U5gOvD7tFQpIpKqE0f8mphtT/gw03wAMKi8MLin0jKYMEejMSIZKJVQ8zIwzcwmA/uATwOf6XHOHmAJ8JyZjQPOA3aks1ARkaScgw/eSPSN2fsSuBgUjIKpS/yITO0SKB4TdqUi0s/6DDXOuaiZ3Q78Dr+l+2Hn3Btm9uXg+QeB/wY8YmZb8NNV/9U5d6gf6xaR4aztOOx4OggyT8DxYPB4/GxYeJcPMpVzIVutuESGk5T+i3fO/Qb4TY9jD3b5fj9wVXpLExEJOAcH30nsVNq9HuLtkF8KtVf6EDN1KZSMD7tSEQmR/jdGRAanSDPsfC4xrXR0jz8+9nxY8FUfZCbOh+zccOsUkUFDoUZEBo+G7UGIWQO7XoBYG+QW+X4xl33d71YaWRV2lSIySCnUiEh42k/48NIxrdQY7C+oOBfm/akPMdULICc/3DpFZEhQqBGRgXV4V+KeSjufhegJyBkBky+Di7/q18aUTQ67ShEZghRqRKR/Rdtgz/pEA7xD7/rjo2v8Ha6nXeXvdp07ItQyRWToU6gRkfQ7Wt9lNOYZiDRBdh7ULIS6L/oGeOW1ujmkiKSVQo2InL1Yu2961xFkDrzhj4+cCLM/6UdjJl8OeUXh1ikiGU2hRkTOzPH3g1sRrIHtT0HbMcjK8Qt7l/03H2TGnKfRGBEZMAo1IpKaeMzfELJjy/X7r/njJefA+R8PRmOugILSUMsUkeFLoUZEetd8CLatC0Zj1sGJw2DZvundkr/1QWbchzQaIyKDgkKNiCTE4/De5sROpX1/ABwUjYFzr/Z9Y2qvhBGjw65UROQkCjUiw11LI2x/MlgfsxZaDgEGVXVw5V/7IDP+w5CVFXalIiKnpFAjMtw4B+9vSdxTqf734OIwogymLvFTSrVLoKg87EpFRE6LQo3IcNB6FHY8HQSZJ6DpfX/8nAvgsr/wQabyQsjKDrNKEZGzolAjkomcgwNvBfdUesJ39I1HIX8kTF2cGI0pGRd2pSIiaaNQI5Ip2pr8vZQ6ppWO1fvj42bBJXf4IFM1D7L1n72IZCb97SYyVDkHDdsSfWN2vwixCOQVw5RFcMVf+ptDjqwMu1IRkQGhUCMylERaYNfzwbTSWn/Ha4Ax02H+l/xozMSLIScv1DJFRMKgUCMy2DXu8It7t66BXc9BtBVyC/29lC65w98ccvSksKsUEQmdQo3IYBNtg90vJBrgNWzzx8tqYe4tvm/MpEshtyDcOkVEBhmFGpHB4MiexB2udz4D7S2QnQ+TL4N5K/zamPLasKsUERnUFGpEwhCNwN4NiSBz8C1/fFQ1XPAZvzam5jLIKwy3ThGRIUShRmSgHHvPL+7duga2Pw2R45CVC5MugTk3+yBTMU03hxQROUMKNSL9JRaF+pcTfWM+2OKPl1bCrD/xIWby5ZBfEm6dIiIZQqFGJJ2aDiRuDLl9nb89gWVD9QJYutIHmbEzNRojItIPFGpEzkY8Bvs3Jxrg7d/sjxePg+nX+J1KUxbBiFFhVikiMiwo1IicrpZG2LbOr4/Z9gS0NIBlQdVFsPhv/GjMuFmQlRV2pSIiw4pCjUhf4nF4/7VE35h9G8HFobDcN76btgxqF0NhWdiViogMaymFGjNbDnwfyAYecs7d3eP5bwA3dbnmDGCMc64xjbWKDJwTR2DHU4kt180HAIMJc+Dyv/SjMRMugKzskAsVEZEOfYYaM8sG/hVYBtQDL5vZaufcmx3nOOf+EfjH4PxrgLsUaGRIcQ4+eCO4p9ITsGcDuBgUjIKpS3yIqV0CxWPCrlRERHqRykjNPGCbc24HgJn9DLgOeLOX828EHktPeSL9qO047HgmseX6+H5/fPxsWHinDzKVdZCtWVoRkaEglb+tK4G9XR7XA/OTnWhmhcBy4PZenl8BrACorq4+rUJFzppzcOjdxE6l3esh3g55JVB7pQ8xU5dC6TlhVyoiImcglVCTrKGG6+Xca4AXept6cs6tAlYB1NXV9XYNkfSJNMPO54JppbX+Hkvge8Us+Kpf6Ft9MWTnhluniIictVRCTT0wscvjKmB/L+d+Gk09Sdgatid2Ku16HmJtkFsEU66AhXf5IDNqYt/XERGRISWVUPMyMM3MJgP78MHlMz1PMrORwBXAzWmtUKQv7a2w+/lEkGnc4Y+XT4OLbvNbriddAjn54dYpIiL9qs9Q45yLmtntwO/wW7ofds69YWZfDp5/MDj1emCNc66536oV6XB4d2KB785nIXoCcgr8vZQu/qpfG1M2OewqRURkAJlz4Sxtqaurcxs3bgzlvWUIikZgz4uJvjGH3vHHR9f4Bb7TroKahZA7ItQyRUTkzJnZJudc3Zm+XntVZfA6us8v7t26FnY8DZEmyM6DSZfC3C/4IFNeq5tDiogIoFAjg0msHfb+PjGtdOANf3zkRJj9yWA05jLILw63ThERGZQUaiRcx9/3HXy3roXtT0HbUcjKgeoFsOzvfZAZM12jMSIi0ieFGhlY8Rjs25RogPfeq/548XiYea0PMVMWQUFpqGWKiMjQo1Aj/a/5EGxb50PM9nVw4jBYFkycD0v+1geZcR/SaIyIiJwVhRpJv3gc3nsl0Tdm3ybAQdEYOPdqmLYUahfDiNFhVyoiIhlEoUbS48Rh2P6kDzLbnoDmg4BB5VxY9Fe+Ad45F0BWVtiViohIhlKokTPjHLy/Jbin0hOw9yVwcT/6MnWpn1KqXQxFFWFXKiIiw4RCjaSu9ZjvF9MRZI6/54+fcwFc9uc+yFTOhazsMKsUEZFhSqFGeuccHHw70Tdmz3qIRyF/JNRe6UPM1KVQMi7sSkVERBRqpIe2Jn8vpY7RmKN7/fFxH4JL7vB3uJ44D7Jzw61TRESkB4Wa4c45aNiW2Km0+wWIRSCv2PeLufwvfJAZWRl2pSIiIqekUDMctZ+AXc8nGuAd3uWPV5wH81b4aaXqBZCTF2qZIiIip0OhZrho3JkYjdn1HERbIWcETLkiMa00elLYVYqIiJwxhZpMFW3zU0lbn/BBpmGrP15WG9zhehlMWgi5BaGWKSIiki4KNZnkyF7YttaPyOx4BtqbITsfahbCRbf5IFNeG3aVIiIi/UKhZiiLtcOeDYkt1wff8sdHVsMFN/q1MTULIa8o3DpFREQGgELNUHPsvS6jMU9D2zHIyoVJl8Ccm3yQqThXN4cUEZFhR6FmsItFYd/GxE6l97f44yUT4PzrfYiZcgXkl4Rbp4iISMgUagajpoO+8d3WNbB9HbQeBcuG6oth6UofZMbO1GiMiIhIFwo1g0E8Bvs3J7Zc798MOCgaC9OvgWlLYcqVMGJU2JWKiIgMWgo1YWlphO1PJm5H0NIAlgWVdXDlt/xOpfGzISsr7EpFRESGBIWagRKPw/uv+dGYbWuh/mVwcSgs9zeFnHYV1C6GwrKwKxURERmSFGr6U+tR2P5UIsg0feCPT7gQLv+GDzIT5kBWdrh1ioiIZACFmnRyDg68megbs2cDuBgUjITaJT7ETF0CxWPDrlRERCTjKNScrbbjvntvx9qYY/v88fGzYOGd/p5KVRdBtn7VIiIi/UmftKfLOTj0bmKn0u4XId4OeSVQuwgWfdOvkSmdEHalIiIiw0pKocbMlgPfB7KBh5xzdyc5ZxHw/wC5wCHn3BVpqzJskRZ/Z+uOBnhH9vjjY2bAxV/x00oT50NOXrh1ioiIDGN9hhozywb+FVgG1AMvm9lq59ybXc4ZBTwALHfO7TGzob9opGF7YjRm1/MQa4PcQpiyCBbe5aeVRk0Mu0oREREJpDJSMw/Y5pzbAWBmPwOuA97scs5ngF845/YAOOcOpLvQftfeCrufh61BJ9/G7f54+TS46FbfN2bSpZCTH26dIiIiklQqoaYS2NvlcT0wv8c55wK5ZvY0UAJ83zn3aM8LmdkKYAVAdXX1mdSbXod3J24OufNZaG+BnAKouQzmf9l38i2bEnaVIiIikoJUQk2yGwy5JNeZCywBRgDrzWyDc+7dbi9ybhWwCqCurq7nNfpfNAJ71ie2XB96xx8fNQnm3OzXxtQshNwRA16aiIiInJ1UQk090HXxSBWwP8k5h5xzzUCzmT0LfBh4l7Ad3ZcYjdnxNESaIDvPTyXN/bwPMuVTdXNIERGRIS6VUPMyMM3MJgP7gE/j19B09f8B95tZDpCHn576v9NZaMriMdj7UmI05oPX/fHSKph1gw8xky+H/OJQyhMREZH+0Weocc5Fzex24Hf4Ld0PO+feMLMvB88/6Jx7y8z+E3gNiOO3fb/en4X3KtYOP/kTiEWgegEs+3sfZMZM12iMiIhIBjPnBn5pC/g1NRs3buyfi+/ZAGNn+NsTiIiIyJBgZpucc3Vn+vrM7ChcfXHYFYiIiMgAywq7ABEREZF0UKgRERGRjKBQIyIiIhlBoUZEREQygkKNiIiIZASFGhEREckICjUiIiKSERRqREREJCMo1IiIiEhGCO02CWZ2ENjdj29RARzqx+uLiIgMN/392TrJOTfmTF8cWqjpb2a28WzuHyEiIiLdDfbPVk0/iYiISEZQqBEREZGMkMmhZlXYBYiIiGSYQf3ZmrFrakRERGR4yeSRGhERERlGFGpEREQkI6Qt1JhZzMxeMbNXzewPZnZJCq/5MzN7y8x+mq46elz/+2a2z8xO+XOa2Sgz++opnn/YzA6Y2es9jn/YzNab2RYz+5WZlQbHl5nZpuD4JjNbnJ6fSEREBMzsejNzZja9y7ExZvaSmW02s8tO9bl2mu81xszazexLKZz7cTOb2ctzlwf5IGpmn+jx3H83s9eDr091Of5TM3snOP6wmeWe6v3TOVJzwjl3gXPuw8BfAf+Qwmu+CvyRc+6mVN7AzHJSLSYIMtcDe4HL+zh9VFBLbx4Blic5/hDwTefcLOCXwDeC44eAa4Ljnwf+V6p1i4iIpOBG4Hng012OLQHeds7NwX/2nVaoMS9ZLrgB2BC8Z18+DiQNNcAe4AvAv/V4348CFwIXAPOBb3QMEgA/BaYDs4ARwG2nevP+mn4qBQ53PDCzb5jZy2b2mpn9XXDsQWAKsNrM7jKzMjP7P8E5G8xsdnDeSjNbZWZrgEeDxPh4cL2XzezSXmq4Engd+AFd/kGY2Tgz+2UwovRqMKJ0N1AbjDT9Y88LOeeeBRqTvMd5wLPB92uBPwnO3+yc2x8cfwMoMLP8lH5zIiIip2BmxcClwK0EocbMLgDuAf7IzF4B/js9Ptd6+SyuCWZMHgD+AExM8pY3An8OVJlZZZc6Phdc61Uz+1/B5+m1wD8G71vb9SLOuV3OudeAeI/rzwSecc5FnXPNwKsEAwnOud+4APB7oOqUvxznXFq+gBjwCvA2cBSYGxy/Cr8FzPAh6tfA5cFzu4CK4Pv7gO8E3y8GXgm+XwlsAkYEj/8NWBh8Xw281Us9DwGfxQesfUBucPx/A3cG32cDI4Ea4PU+fr6TzgFeBK4Lvv86cDzJ6z4BPJGu37O+9KUvfelreH8BNwM/DL5/Ebgw+P4LwP3B990+s3r7LA7OiwMX9/JeE4Gtwff/F/D14PvzgXe6fIaXBX8+Anyij/q7nRPU9gJQiL8Nww7gz3u8Jhcfui471bX7Y/ppOj5hPWpmFhR7FbA5KGg6MC3J6xcSTNM4554Eys1sZPDcaufcieD7pcD9QRJdDZSaWUnXC5lZHvBHwP9xzh0DXgpqAB+YfhC8T8w5d/QsfuYvAl8zs01ACRDpUcf5+LTc5zykiIhIim4EfhZ8/zNSmxY61Wfxbufchl5e92ng/03yXouBnzvnDgE455LNZqTEObcG+A0+oD0GrAeiPU57AHjWOffcqa6V8hqV0yxwvZlVAGPwqfAfnHP/o4+XWbJLBX82dzmWBSzoEnKSWY4fgdnicxWFQAvwHymUnzLn3NsEYcnMzgU+2vGcmVXh19l8zjm3PZ3vKyIiw5OZleMDxYfMzOFnHJyZ/WVfLyXJZ7GZ1dD9M7anG4FxZtax9nWCmU0Lrpe2RnfOue8B3wtq+jdga5cav4PPE30OEPTLmppgNXY20AD8DvhiMAeImVWa2dgkL3sWuCk4ZxFwKBhl6WkNcHuX97ogyTk3Arc552qcczXAZOAqMysE1gFfCV6bHSxGOo4faTndn3Ns8GcW8DfAg8HjUfgA9VfOuRdO97oiIiK9+ATwqHNuUvAZNxHYiZ/t6Krn51qqn8WdzOw8oMg5V9nl8/Qf8KM364BPBiELMyvr5X37FHwWd1xnNjAb/1mPmd0GfAS40TnXcy3OSdIZakYEC4Newa9b+XwwvbMGvw5mvZltAX5O8h94JVBnZq/hF+5+vpf3+bOO88zsTeDLXZ8MgstH6DIq4/zCo+eBa4D/AlwZ1LIJON851wC8EGwZO2mhsJl1DIedZ2b1ZnZr8NSNZvYufh3RfuBHwfHbganAtzt+J339yyMiIpKCG/GzAF09Dnym64Gen2un8Vmcynvd6Jx7Az+y8oyZvQr8c/D8z/C7lzb3XChsZheZWT1+N9X/MLM3gqdygeeCz/RVwM3OuY7ppweBcUHdr5jZ356qYN0mQURERDKCOgqLiIhIRlCoERERkYygUCMiIiIZQaFGREREMoJCjYiIiGQEhRoRERHJCAo1IiIikhH+fyietvTE1vKrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "did_df_treatment = did_df[did_df['T'] == 1]\n",
    "did_df_control = did_df[did_df['T'] == 0]\n",
    "\n",
    "mean_treatment_before = np.mean(did_df_treatment[did_df_treatment['I'] == 0]['y'])\n",
    "mean_treatment_after = np.mean(did_df_treatment[did_df_treatment['I'] == 1]['y'])\n",
    "\n",
    "mean_control_before = np.mean(did_df_control[did_df_control['I'] == 0]['y'])\n",
    "mean_control_after = np.mean(did_df_control[did_df_control['I'] == 1]['y'])\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 4), tight_layout=True)\n",
    "plt.plot([0, 1], [mean_treatment_before, mean_treatment_after], label='philly')\n",
    "ax.set_xticks([0, 1])\n",
    "plt.plot([0, 1], [mean_control_before, mean_control_after], label='camden')\n",
    "ax.set_xticklabels(['Before Act 192', 'After Act 192'])\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4w1NQC6oD3Q0"
   },
   "source": [
    "**Simple Version of Difference-in-Difference Calculation**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "Gi23kJKACsRS",
    "outputId": "c7d3110d-cd70-4b57-eb97-4ca591484f2d"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5126821370750135"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Difference-in-difference (before)\n",
    "mean_treatment_before - mean_control_before"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "dKWxh9S2DQ4n",
    "outputId": "83a2b6e1-7798-4767-97fa-b9b4b0922dc1"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2494172494172493"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Difference-in-difference (after)\n",
    "mean_treatment_after - mean_control_after"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "eAJqs9yKEIoj"
   },
   "source": [
    "**Difference-in-Difference Using Fit OLS Model**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "id": "eXCjX3CgEH19"
   },
   "outputs": [],
   "source": [
    "def fit_OLS_model(df, target_variable, explanatory_variables, intercept = False):\n",
    "    \"\"\"\n",
    "    Fits an OLS model from data.\n",
    "    \n",
    "    Inputs:\n",
    "        df: pandas DataFrame\n",
    "        target_variable: string, name of the target variable\n",
    "        explanatory_variables: list of strings, names of the explanatory variables\n",
    "        intercept: bool, if True add intercept term\n",
    "    Outputs:\n",
    "        fitted_model: model containing OLS regression results\n",
    "    \"\"\"\n",
    "    \n",
    "    target = df[target_variable]\n",
    "    inputs = df[explanatory_variables]\n",
    "    if intercept:\n",
    "        inputs = sm.add_constant(inputs)\n",
    "    \n",
    "    fitted_model = sm.OLS(target, inputs).fit()\n",
    "    return(fitted_model)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "jFm84PG4DZLT",
    "outputId": "b0ae7c72-fb00-4833-c97e-773fa7c24df1"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\yangj\\anaconda3\\lib\\site-packages\\statsmodels\\tsa\\tsatools.py:142: FutureWarning: In a future version of pandas all arguments of concat except for the argument 'objs' will be keyword-only\n",
      "  x = pd.concat(x[::order], 1)\n"
     ]
    }
   ],
   "source": [
    "linear_model = fit_OLS_model(did_df, 'y', ['T', 'I', 'T*I', 'median_income', 'pop_density', 'unemployment', 'mean_participant_age'], intercept=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "R-NyDIbuFbFn",
    "outputId": "e6d022e8-908b-4942-a266-997ad3aae246"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.094\n",
      "Model:                            OLS   Adj. R-squared:                  0.086\n",
      "Method:                 Least Squares   F-statistic:                     11.94\n",
      "Date:                Wed, 04 May 2022   Prob (F-statistic):           3.08e-09\n",
      "Time:                        18:18:16   Log-Likelihood:                -446.12\n",
      "No. Observations:                 465   AIC:                             902.2\n",
      "Df Residuals:                     460   BIC:                             923.0\n",
      "Df Model:                           4                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "========================================================================================\n",
      "                           coef    std err          t      P>|t|      [0.025      0.975]\n",
      "----------------------------------------------------------------------------------------\n",
      "const                  -46.1205     24.137     -1.911      0.057     -93.553       1.312\n",
      "T                       63.3450     33.152      1.911      0.057      -1.803     128.493\n",
      "I                        0.1763      0.065      2.718      0.007       0.049       0.304\n",
      "T*I                     63.3450     33.152      1.911      0.057      -1.803     128.493\n",
      "median_income           -0.0024      0.001     -1.893      0.059      -0.005    9.21e-05\n",
      "pop_density              0.1088      0.057      1.913      0.056      -0.003       0.221\n",
      "unemployment            -5.6195      2.941     -1.911      0.057     -11.399       0.160\n",
      "mean_participant_age    -0.0016      0.003     -0.569      0.570      -0.007       0.004\n",
      "==============================================================================\n",
      "Omnibus:                      134.413   Durbin-Watson:                   1.970\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):              425.322\n",
      "Skew:                           1.334   Prob(JB):                     4.39e-93\n",
      "Kurtosis:                       6.851   Cond. No.                     2.41e+20\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "[2] The smallest eigenvalue is 1.17e-29. This might indicate that there are\n",
      "strong multicollinearity problems or that the design matrix is singular.\n"
     ]
    }
   ],
   "source": [
    "print(linear_model.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qlO1Ff5fHksY"
   },
   "source": [
    "# Analysis\n",
    "\n",
    "Observing the Fit OLS model generated based on our time x treatment interaction $T*I$ against our treatment variable (Act of 192) and confounding variables (median income, population density, unemployment, and mean of participant's age) the constant coefficient demonstrated positive number of 63.345, which signifies positive relationship between our target variable against our confounding variables before and after intervention of Act 192. This also signifies that the general number of injuries and deaths from gun violence incidents have gone up over time.\n",
    "\n",
    "As indicated from the summary of the fit OLS model, the standard error for our interaction variable is fairly high, but it still is within the range of a positive interaction coefficient. The standard error is 33.152 for our $T*I$ term.\n",
    "\n",
    "Since Act 192 allowed individuals to sue cities with strict gun laws, relaxing the gun laws in our treatment region was expected to increase harm of gun violence incidents. So in line with our hypothesis, we believe **Act 192 increased the number of people who were killed or injured in each gun violence incident.**\n",
    "\n",
    "The coefficients of our confounders seemed to have minor effects with the exception of unemployment. Unemployment lead to a decrease in gun violence which is an unexpected result, but considering the limitation of using mean unemployment for our treatment and control groups, it might not be that informative of a term. In a future study, we would recommend researching more granular unemployment statistics to see if the effect of its coefficient still decreases the amount of deaths and injuries from gun violence incidents.\n",
    "\n",
    "Another limitation of our study are the limited number of confounders. While we included four confounding variables, mean participant age was the only one that had individual-level values for each unit in the study. The others were aggregate statistics corresponding to the entirety of South Philadelphia, our treatment, and Camden, our control. For future studies, we would recommend finding other confounders that could explain discrepancies such as rates of policing on a neighborhood level, net worth, educational attainment, and other individual or neighborhood-level variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "QnnWY7BEKGF9"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "Fit_OLS_Model data102final.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}