import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from datetime import date
import warnings
sns.set_theme(style="darkgrid", font_scale=1.2)
plt.rcParams['figure.figsize'] = (20.0, 10.0)
warnings.filterwarnings('ignore')
Research Question 1
Data Cleaning
lockdown = pd.read_csv("lockdown_dates.csv")
# Cleaning the lockdown dates dataset.
lockdown_c = lockdown[lockdown['County'].isnull()].drop(columns=['Country','County'])
lockdown_c = lockdown_c[lockdown_c["State"] != "Puerto Rico"]
lockdown_c = lockdown_c.rename(columns={"State": "state", "Date": "lockdown_start", "Type": "lockdown_type"}).replace({'Stay at home': 0, 'Shelter in place': 1})
lockdown_c["lockdown_end"] = ["2020-04-30", "2020-04-24", "2020-05-15", "2021-06-15", "2020-04-27", "2020-05-20", "2020-05-15",
"2020-05-15", "2020-05-04", "2020-04-24", "2020-05-31", "2020-04-30", "2020-05-30", "2020-05-01",
"2020-05-03", "2020-06-29", "2020-05-14", "2020-05-31", "2020-05-15", "2020-05-18", "2020-06-05",
"2020-05-04", "2020-05-11", "2020-05-03", "2020-04-24", "2020-05-15", "2020-06-15", "2020-06-09",
"2020-05-15", "2020-05-15", "2020-05-08", "2020-05-29", "2020-06-19", "2020-05-08", "2020-05-08",
"2020-05-04", "2020-04-30", "2020-05-01", "2020-05-15", "2020-06-10", "2020-05-31", "2020-05-04", "2020-05-26"]
lockdown_c = lockdown_c[["state", "lockdown_start", "lockdown_end", "lockdown_type"]]
lockdown_c["lockdown_start"] = pd.to_datetime(lockdown_c["lockdown_start"], format='%Y-%m-%d')
lockdown_c["lockdown_end"] = pd.to_datetime(lockdown_c["lockdown_end"], format='%Y-%m-%d')
lockdown_c["lockdown_duration"] = (lockdown_c["lockdown_end"] - lockdown_c["lockdown_start"]).dt.days
lockdown_c = lockdown_c.sort_values(by="lockdown_duration", ascending=False)
states_of_interest = ["California", "Kentucky", "Washington", "New Jersey", "New Hampshire", "Georgia", "South Carolina", "Alabama", "Florida", "Texas"]
lockdown_j = lockdown_c[lockdown_c["state"].isin(states_of_interest)]
lockdown_j
stateobject
California10%
Kentucky10%
8 others80%
lockdown_startdatetime64[ns]
2020-03-19 00:00:00 - 2020-04-07 00:00:00
4
California
2020-03-19 00:00:00
44
Kentucky
2020-03-26 00:00:00
74
New Jersey
2020-03-21 00:00:00
73
New Hampshire
2020-03-28 00:00:00
169
Washington
2020-03-23 00:00:00
9
Florida
2020-04-03 00:00:00
126
Texas
2020-04-02 00:00:00
118
South Carolina
2020-04-07 00:00:00
0
Alabama
2020-04-04 00:00:00
21
Georgia
2020-04-03 00:00:00
flights = pd.read_csv("flights.csv")
# Cleaning the flight delays dataset.
airports_of_interest = ["LAX", "CVG", "SEA", "EWR", "MHT", "ATL", "GSP", "BHM", "MCO", "DFW"]
months = ["Jan", "Feb", "Mar", "Apr", "May", "Jun"]
flights_c = flights[["ORIGIN_STATE_NM", "ORIGIN", "DEST_STATE_NM", "DEST", "FL_DATE", "MONTH", "DISTANCE", "CANCELLED", "DEP_DELAY_NEW"]]
flights_c = flights_c.rename(columns={"ORIGIN_STATE_NM": "orig_state", "ORIGIN": "orig_airport", "DEST_STATE_NM": "dest_state",
"DEST": "dest_airport", "FL_DATE": "flight_date", "MONTH": "flight_month",
"DISTANCE": "distance", "CANCELLED": "cancelled", "DEP_DELAY_NEW": "dep_delay"})
flights_c["cancelled"] = flights_c["cancelled"].astype(bool)
flights_c["flight_date"] = pd.to_datetime(flights_c["flight_date"], format='%m/%d/%Y')
flights_c["flight_month"] = flights_c["flight_month"].replace({1: "Jan", 2: "Feb", 3: "Mar", 4: "Apr", 5: "May", 6: "Jun"})
flights_c["flight_month"] = pd.Categorical(flights_c["flight_month"], categories=months, ordered=True)
flights_c
orig_stateobject
orig_airportobject
0
California
ONT
1
California
ONT
2
California
ONT
3
California
ONT
4
California
ONT
5
California
ONT
6
California
ONT
7
California
ONT
8
California
ONT
9
California
ONT
# Merged dataset w/ lockdown dataset focusing on the flights that departed from the states of interest.
flights_orig = flights_c[flights_c["orig_state"].isin(states_of_interest)]
#flights_orig = flights_orig[flights_orig["orig_airport"].isin(airports_of_interest)]
flights_orig_m = flights_orig.merge(lockdown_c, left_on="orig_state", right_on="state", how="left")
flights_orig_m = flights_orig_m.drop(columns=["dest_state", "dest_airport", "state"])
flights_orig_m["lockdown_flight"] = ((flights_orig_m["lockdown_start"] < flights_orig_m["flight_date"]) & (flights_orig_m["lockdown_end"] > flights_orig_m["flight_date"]))
flights_orig_m
orig_stateobject
orig_airportobject
220
Washington
SEA
221
Washington
SEA
222
Washington
SEA
223
Washington
SEA
224
Washington
SEA
225
Washington
SEA
226
Washington
SEA
227
Washington
SEA
228
Washington
SEA
229
Washington
SEA
# Merged dataset w/ lockdown dataset focusing on the flights that arrived (intended to arrive, if cancelled) in the states of interest.
flights_dest = flights_c[flights_c["dest_state"].isin(states_of_interest)]
#flights_dest = flights_dest[flights_dest["dest_airport"].isin(airports_of_interest)]
flights_dest_m = flights_dest.merge(lockdown_c, left_on="dest_state", right_on="state", how="left")
flights_dest_m = flights_dest_m.drop(columns=["orig_state", "orig_airport", "state"])
flights_dest_m["lockdown_flight"] = ((flights_dest_m["lockdown_start"] < flights_dest_m["flight_date"]) & (flights_dest_m["lockdown_end"] > flights_dest_m["flight_date"]))
flights_dest_m
dest_stateobject
dest_airportobject
0
California
SFO
1
California
SFO
2
California
SJC
3
California
SJC
4
California
SJC
5
California
SJC
6
California
SJC
7
California
SMF
8
California
SMF
9
California
SMF
Hypothesis Testing
def avg_difference_in_means(state, binary_col, numerical_col):
state_select = flights_orig_m.loc[flights_orig_m["orig_state"] == state]
observed_avgs = state_select.groupby(binary_col).mean()
observed_diff = abs(observed_avgs[numerical_col][1] - observed_avgs[numerical_col][0])
# Should be True - False
differences = []
repititions = 10000
for i in np.arange(repititions):
shuffled_flights = state_select[[binary_col, numerical_col]]
shuffled_flights[binary_col] = shuffled_flights[binary_col].sample(frac=1).values
shuffled_avgs = shuffled_flights.groupby(binary_col).mean()
differences = np.append(differences, abs(shuffled_avgs[numerical_col][1] - shuffled_avgs[numerical_col][0]))
return np.count_nonzero(differences >= observed_diff) / repititions
cancels_p_values = {"Alabama": avg_difference_in_means("Alabama", "lockdown_flight", "cancelled"),
"California": avg_difference_in_means("California", "lockdown_flight", "cancelled"),
"Florida": avg_difference_in_means("Florida", "lockdown_flight", "cancelled"),
"Georgia": avg_difference_in_means("Georgia", "lockdown_flight", "cancelled"),
"Kentucky": avg_difference_in_means("Kentucky", "lockdown_flight", "cancelled"),
"New Hampshire": avg_difference_in_means("New Hampshire", "lockdown_flight", "cancelled"),
"New Jersey": avg_difference_in_means("New Jersey", "lockdown_flight", "cancelled"),
"South Carolina": avg_difference_in_means("South Carolina", "lockdown_flight", "cancelled"),
"Texas": avg_difference_in_means("Texas", "lockdown_flight", "cancelled"),
"Washington": avg_difference_in_means("Washington", "lockdown_flight", "cancelled")}
cancels_p_values
While we initially had a standard p-value cutoff of 0.05, we can see from these results we would reject the null hypothesis for all ten states because they all had p-values of 0.0. These values imply that the presence of a lockdown order did have a significant effect on cancellation rates of flights. This is corroborated by one of our visualizations from earlier, which indicated that flights that occurred during lockdown had significantly higher cancellation rates than those that didn't for all ten states.
delay_p_values = {"Alabama": avg_difference_in_means("Alabama", "lockdown_flight", "dep_delay"),
"California": avg_difference_in_means("California", "lockdown_flight", "dep_delay"),
"Florida": avg_difference_in_means("Florida", "lockdown_flight", "dep_delay"),
"Georgia": avg_difference_in_means("Georgia", "lockdown_flight", "dep_delay"),
"Kentucky": avg_difference_in_means("Kentucky", "lockdown_flight", "dep_delay"),
"New Hampshire": avg_difference_in_means("New Hampshire", "lockdown_flight", "dep_delay"),
"New Jersey": avg_difference_in_means("New Jersey", "lockdown_flight", "dep_delay"),
"South Carolina": avg_difference_in_means("South Carolina", "lockdown_flight", "dep_delay"),
"Texas": avg_difference_in_means("Texas", "lockdown_flight", "dep_delay"),
"Washington": avg_difference_in_means("Washington", "lockdown_flight", "dep_delay")}
delay_p_values
We also performed similar hypothesis tests under the same procedures for the effect of lockdown orders on flight departure delays. The ten p-values above all fall under the threshold of 0.05 with one exception: New Hampshire. These results imply that lockdown orders also had noticeable impacts on not just cancellation rates, but also the length of departure delays.
Error Rate Control
bonferroni_cancels = dict((k, v) for k, v in cancels_p_values.items() if v <= 0.05 / 10)
bonferroni_cancels
bonferroni_delays = dict((k, v) for k, v in delay_p_values.items() if v <= 0.05 / 10)
bonferroni_delays
For the Bonferroni method, which controls FWER, given that we have 10 tests with a general p-value threshold of 0.05, we would set the cutoff line at 0.05 / 10 = 0.005. However, given that all of our p-values are 0 for the cancellations hypothesis testing, we would still reject the null hypothesis for all cases. There is a similar case for the delays hypothesis testing; we still only fail to reject the null hypothesis for one state, New Hampshire.
def benjamini_hochberg(p_values, alpha):
p_values_sort = sorted(p_values)
kth = 0
for k in np.arange(len(p_values_sort)):
if p_values_sort[k] <= ((k + 1) * alpha / len(p_values)):
kth = k
decisions = (p_values <= p_values_sort[kth])
return decisions
dict(list({k: v for k, v in sorted(cancels_p_values.items(), key=lambda item: item[1])}.items())[0: np.sum(benjamini_hochberg(pd.Series(cancels_p_values), 0.05))])
dict(list({k: v for k, v in sorted(delay_p_values.items(), key=lambda item: item[1])}.items())[0: np.sum(benjamini_hochberg(pd.Series(delay_p_values), 0.05))])
For the Benjamini-Hochberg method, which controls FDR, we see a very similar case to the Bonferroni method, where almost all discoveries are accepted with the exception of New Hampshire for delays.
Research Question 2
Data Cleaning
# California Gas Prices
ca_gas = pd.read_csv("california_gas.csv")
ca_gas.columns = ca_gas.iloc[1].fillna('DROP')
ca_gas = ca_gas.drop([0, 1]).drop("DROP", axis=1)
ca_gas = ca_gas.rename(columns={"Date": "date",
"Weekly California All Grades All Formulations Retail Gasoline Prices (Dollars per Gallon)": "gas_all_grades",
"Weekly California Regular All Formulations Retail Gasoline Prices (Dollars per Gallon)": "regular",
"Weekly California Midgrade All Formulations Retail Gasoline Prices (Dollars per Gallon)": "midgrade",
"Weekly California Premium All Formulations Retail Gasoline Prices (Dollars per Gallon)": "premium",
"Weekly California No 2 Diesel Retail Prices (Dollars per Gallon)": "diesel"})
ca_gas_c = ca_gas[["date", "gas_all_grades", "regular", "midgrade", "premium", "diesel"]]
ca_gas_c["date"] = pd.to_datetime(ca_gas_c["date"], format="%b %d, %Y")
date_range = (ca_gas_c["date"] > "2020-1-1") & (ca_gas_c["date"] <= "2022-1-1")
ca_gas_c = ca_gas_c.loc[date_range]
ca_gas_c["gas_all_grades"] = ca_gas_c["gas_all_grades"].astype(float)
ca_gas_c["regular"] = ca_gas_c["regular"].astype(float)
ca_gas_c["midgrade"] = ca_gas_c["midgrade"].astype(float)
ca_gas_c["premium"] = ca_gas_c["premium"].astype(float)
ca_gas_c["diesel"] = ca_gas_c["diesel"].astype(float)
ca_gas_c["lockdown"] = (("2020-03-19" < ca_gas_c["date"]) & ("2021-06-15" > ca_gas_c["date"]))
ca_gas_c
datedatetime64[ns]
2020-01-06 00:00:00 - 2021-12-27 00:00:00
gas_all_gradesfloat64
2.721 - 4.642
1307
2020-01-06 00:00:00
3.505
1308
2020-01-13 00:00:00
3.497
1309
2020-01-20 00:00:00
3.487
1310
2020-01-27 00:00:00
3.467
1311
2020-02-03 00:00:00
3.455
1312
2020-02-10 00:00:00
3.449
1313
2020-02-17 00:00:00
3.441
1314
2020-02-24 00:00:00
3.443
1315
2020-03-02 00:00:00
3.435
1316
2020-03-09 00:00:00
3.403
cases = pd.read_csv("covid_cases.csv")
cases = cases[["submission_date", "state", "new_case"]]
ca_cases = cases[cases["state"] == "CA"]
ca_cases["submission_date"] = pd.to_datetime(ca_cases["submission_date"], format="%m/%d/%Y")
ca_cases
submission_datedatetime64[ns]
2020-01-22 00:00:00 - 2022-10-18 00:00:00
stateobject
CA100%
15543
2020-06-11 00:00:00
CA
15553
2020-08-18 00:00:00
CA
15562
2021-02-25 00:00:00
CA
15573
2021-02-10 00:00:00
CA
15576
2020-09-06 00:00:00
CA
15639
2020-10-16 00:00:00
CA
15649
2021-02-26 00:00:00
CA
15666
2020-11-10 00:00:00
CA
15698
2021-06-28 00:00:00
CA
15704
2020-08-10 00:00:00
CA
ca_cases["submission_date"] = ca_cases["submission_date"] - pd.to_timedelta(7, unit='d')
week_cases = ca_cases.groupby(["state", pd.Grouper(key="submission_date", freq='W-MON')])["new_case"].sum().reset_index().sort_values("submission_date")
cases_date_range = (week_cases["submission_date"] > "2020-1-1") & (week_cases["submission_date"] <= "2022-1-1")
week_cases = week_cases.loc[cases_date_range]
ca_gas_m = ca_gas_c.merge(week_cases, left_on="date", right_on="submission_date", how="left")
ca_gas_m = ca_gas_m.drop(columns=["state", "submission_date"])
ca_gas_m = ca_gas_m.fillna(0)
ca_gas_m
datedatetime64[ns]
2020-01-06 00:00:00 - 2021-12-27 00:00:00
gas_all_gradesfloat64
2.721 - 4.642
0
2020-01-06 00:00:00
3.505
1
2020-01-13 00:00:00
3.497
2
2020-01-20 00:00:00
3.487
3
2020-01-27 00:00:00
3.467
4
2020-02-03 00:00:00
3.455
5
2020-02-10 00:00:00
3.449
6
2020-02-17 00:00:00
3.441
7
2020-02-24 00:00:00
3.443
8
2020-03-02 00:00:00
3.435
9
2020-03-09 00:00:00
3.403
# Florida Gas Prices
fl_gas = pd.read_csv("florida_gas.csv")
fl_gas.columns = fl_gas.iloc[1].fillna('DROP')
fl_gas = fl_gas.drop([0, 1]).drop("DROP", axis=1)
fl_gas = fl_gas.rename(columns={"Date": "date",
"Weekly Florida All Grades All Formulations Retail Gasoline Prices (Dollars per Gallon)": "gas_all_grades",
"Weekly Florida Regular All Formulations Retail Gasoline Prices (Dollars per Gallon)": "regular",
"Weekly Florida Midgrade All Formulations Retail Gasoline Prices (Dollars per Gallon)": "midgrade",
"Weekly Florida Premium All Formulations Retail Gasoline Prices (Dollars per Gallon)": "premium"})
fl_gas_c = fl_gas[["date", "gas_all_grades", "regular", "midgrade", "premium"]]
fl_gas_c["date"] = pd.to_datetime(fl_gas_c["date"], format="%b %d, %Y")
date_range = (fl_gas_c["date"] > "2020-1-1") & (fl_gas_c["date"] <= "2022-1-1")
fl_gas_c = fl_gas_c.loc[date_range]
fl_gas_c["gas_all_grades"] = fl_gas_c["gas_all_grades"].astype(float)
fl_gas_c["regular"] = fl_gas_c["regular"].astype(float)
fl_gas_c["midgrade"] = fl_gas_c["midgrade"].astype(float)
fl_gas_c["premium"] = fl_gas_c["premium"].astype(float)
fl_gas_c["lockdown"] = (("2020-04-03" < fl_gas_c["date"]) & ("2020-05-04" > fl_gas_c["date"]))
fl_gas_c
datedatetime64[ns]
2020-01-06 00:00:00 - 2021-12-27 00:00:00
gas_all_gradesfloat64
1.817 - 3.396
959
2021-09-27 00:00:00
3.096
960
2021-10-04 00:00:00
3.074
961
2021-10-11 00:00:00
3.22
962
2021-10-18 00:00:00
3.221
963
2021-10-25 00:00:00
3.362
964
2021-11-01 00:00:00
3.348
965
2021-11-08 00:00:00
3.321
966
2021-11-15 00:00:00
3.299
967
2021-11-22 00:00:00
3.396
968
2021-11-29 00:00:00
3.376
cases = pd.read_csv("covid_cases.csv")
cases = cases[["submission_date", "state", "new_case"]]
fl_cases = cases[cases["state"] == "FL"]
fl_cases["submission_date"] = pd.to_datetime(fl_cases["submission_date"], format="%m/%d/%Y")
fl_cases
submission_datedatetime64[ns]
2020-01-22 00:00:00 - 2022-10-18 00:00:00
stateobject
FL100%
78
2022-08-22 00:00:00
FL
93
2021-04-10 00:00:00
FL
102
2020-10-17 00:00:00
FL
112
2021-08-18 00:00:00
FL
133
2021-03-10 00:00:00
FL
191
2022-08-25 00:00:00
FL
230
2021-01-18 00:00:00
FL
272
2021-10-13 00:00:00
FL
316
2020-02-09 00:00:00
FL
334
2021-09-23 00:00:00
FL
fl_cases["submission_date"] = fl_cases["submission_date"] - pd.to_timedelta(7, unit='d')
week_cases = fl_cases.groupby(["state", pd.Grouper(key="submission_date", freq='W-MON')])["new_case"].sum().reset_index().sort_values("submission_date")
cases_date_range = (week_cases["submission_date"] > "2020-1-1") & (week_cases["submission_date"] <= "2022-1-1")
week_cases = week_cases.loc[cases_date_range]
fl_gas_m = fl_gas_c.merge(week_cases, left_on="date", right_on="submission_date", how="left")
fl_gas_m = fl_gas_m.drop(columns=["state", "submission_date"])
fl_gas_m = fl_gas_m.fillna(0)
fl_gas_m
datedatetime64[ns]
2020-01-06 00:00:00 - 2021-12-27 00:00:00
gas_all_gradesfloat64
1.817 - 3.396
0
2020-01-06 00:00:00
2.576
1
2020-01-13 00:00:00
2.52
2
2020-01-20 00:00:00
2.546
3
2020-01-27 00:00:00
2.474
4
2020-02-03 00:00:00
2.401
5
2020-02-10 00:00:00
2.409
6
2020-02-17 00:00:00
2.357
7
2020-02-24 00:00:00
2.444
8
2020-03-02 00:00:00
2.374
9
2020-03-09 00:00:00
2.309
# New York Gas Prices
ny_gas = pd.read_csv("new_york_gas.csv")
ny_gas.columns = ny_gas.iloc[1].fillna('DROP')
ny_gas = ny_gas.drop([0, 1]).drop("DROP", axis=1)
ny_gas = ny_gas.rename(columns={"Date": "date",
"Weekly New York All Grades All Formulations Retail Gasoline Prices (Dollars per Gallon)": "gas_all_grades",
"Weekly New York Regular All Formulations Retail Gasoline Prices (Dollars per Gallon)": "regular",
"Weekly New York Midgrade All Formulations Retail Gasoline Prices (Dollars per Gallon)": "midgrade",
"Weekly New York Premium All Formulations Retail Gasoline Prices (Dollars per Gallon)": "premium"})
ny_gas_c = ny_gas[["date", "gas_all_grades", "regular", "midgrade", "premium"]]
ny_gas_c["date"] = pd.to_datetime(ny_gas_c["date"], format="%b %d, %Y")
date_range = (ny_gas_c["date"] > "2020-1-1") & (ny_gas_c["date"] <= "2022-1-1")
ny_gas_c = ny_gas_c.loc[date_range]
ny_gas_c["gas_all_grades"] = ny_gas_c["gas_all_grades"].astype(float)
ny_gas_c["regular"] = ny_gas_c["regular"].astype(float)
ny_gas_c["midgrade"] = ny_gas_c["midgrade"].astype(float)
ny_gas_c["premium"] = ny_gas_c["premium"].astype(float)
ny_gas_c["lockdown"] = (("2020-03-22" < ny_gas_c["date"]) & ("2020-05-15" > ny_gas_c["date"]))
ny_gas_c
datedatetime64[ns]
2020-01-06 00:00:00 - 2021-12-27 00:00:00
gas_all_gradesfloat64
2.124 - 3.573
1024
2020-01-06 00:00:00
2.722
1025
2020-01-13 00:00:00
2.725
1026
2020-01-20 00:00:00
2.705
1027
2020-01-27 00:00:00
2.673
1028
2020-02-03 00:00:00
2.635
1029
2020-02-10 00:00:00
2.609
1030
2020-02-17 00:00:00
2.587
1031
2020-02-24 00:00:00
2.584
1032
2020-03-02 00:00:00
2.57
1033
2020-03-09 00:00:00
2.547
cases = pd.read_csv("covid_cases.csv")
cases = cases[["submission_date", "state", "new_case"]]
ny_cases = cases[cases["state"] == "NY"]
ny_cases["submission_date"] = pd.to_datetime(ny_cases["submission_date"], format="%m/%d/%Y")
ny_cases
submission_datedatetime64[ns]
2020-01-22 00:00:00 - 2022-10-18 00:00:00
stateobject
NY100%
58450
2022-08-19 00:00:00
NY
8440
2022-08-20 00:00:00
NY
1233
2022-08-21 00:00:00
NY
8544
2022-08-22 00:00:00
NY
52381
2022-08-23 00:00:00
NY
53271
2022-08-24 00:00:00
NY
45058
2022-08-25 00:00:00
NY
56764
2022-08-26 00:00:00
NY
4955
2022-08-27 00:00:00
NY
51429
2022-08-28 00:00:00
NY
ny_cases["submission_date"] = ny_cases["submission_date"] - pd.to_timedelta(7, unit='d')
week_cases = ny_cases.groupby(["state", pd.Grouper(key="submission_date", freq='W-MON')])["new_case"].sum().reset_index().sort_values("submission_date")
cases_date_range = (week_cases["submission_date"] > "2020-1-1") & (week_cases["submission_date"] <= "2022-1-1")
week_cases = week_cases.loc[cases_date_range]
ny_gas_m = ny_gas_c.merge(week_cases, left_on="date", right_on="submission_date", how="left")
ny_gas_m = ny_gas_m.drop(columns=["state", "submission_date"])
ny_gas_m = ny_gas_m.fillna(0)
ny_gas_m
datedatetime64[ns]
2020-01-06 00:00:00 - 2021-12-27 00:00:00
gas_all_gradesfloat64
2.124 - 3.573
70
2021-05-10 00:00:00
3.017
71
2021-05-17 00:00:00
3.061
72
2021-05-24 00:00:00
3.068
73
2021-05-31 00:00:00
3.062
74
2021-06-07 00:00:00
3.088
75
2021-06-14 00:00:00
3.135
76
2021-06-21 00:00:00
3.146
77
2021-06-28 00:00:00
3.156
78
2021-07-05 00:00:00
3.186
79
2021-07-12 00:00:00
3.201
# Texas Gas Prices
tx_gas = pd.read_csv("texas_gas.csv")
tx_gas.columns = tx_gas.iloc[1].fillna('DROP')
tx_gas = tx_gas.drop([0, 1]).drop("DROP", axis=1)
tx_gas = tx_gas.rename(columns={"Date": "date",
"Weekly Texas All Grades All Formulations Retail Gasoline Prices (Dollars per Gallon)": "gas_all_grades",
"Weekly Texas Regular All Formulations Retail Gasoline Prices (Dollars per Gallon)": "regular",
"Weekly Texas Midgrade All Formulations Retail Gasoline Prices (Dollars per Gallon)": "midgrade",
"Weekly Texas Premium All Formulations Retail Gasoline Prices (Dollars per Gallon)": "premium"})
tx_gas_c = tx_gas[["date", "gas_all_grades", "regular", "midgrade", "premium"]]
tx_gas_c["date"] = pd.to_datetime(tx_gas_c["date"], format="%b %d, %Y")
date_range = (tx_gas_c["date"] > "2020-1-1") & (tx_gas_c["date"] <= "2022-1-1")
tx_gas_c = tx_gas_c.loc[date_range]
tx_gas_c["gas_all_grades"] = tx_gas_c["gas_all_grades"].astype(float)
tx_gas_c["regular"] = tx_gas_c["regular"].astype(float)
tx_gas_c["midgrade"] = tx_gas_c["midgrade"].astype(float)
tx_gas_c["premium"] = tx_gas_c["premium"].astype(float)
tx_gas_c["lockdown"] = (("2020-04-02" < tx_gas_c["date"]) & ("2021-05-01" > tx_gas_c["date"]))
tx_gas_c
datedatetime64[ns]
2020-01-06 00:00:00 - 2021-12-27 00:00:00
gas_all_gradesfloat64
1.506 - 3.134
1024
2020-01-06 00:00:00
2.325
1025
2020-01-13 00:00:00
2.324
1026
2020-01-20 00:00:00
2.274
1027
2020-01-27 00:00:00
2.244
1028
2020-02-03 00:00:00
2.18
1029
2020-02-10 00:00:00
2.134
1030
2020-02-17 00:00:00
2.127
1031
2020-02-24 00:00:00
2.204
1032
2020-03-02 00:00:00
2.163
1033
2020-03-09 00:00:00
2.111
cases = pd.read_csv("covid_cases.csv")
cases = cases[["submission_date", "state", "new_case"]]
tx_cases = cases[cases["state"] == "TX"]
tx_cases["submission_date"] = pd.to_datetime(tx_cases["submission_date"], format="%m/%d/%Y")
tx_cases
submission_datedatetime64[ns]
2020-01-22 00:00:00 - 2022-10-18 00:00:00
stateobject
TX100%
12323
2021-09-03 00:00:00
TX
8963
2021-09-04 00:00:00
TX
14305
2021-09-05 00:00:00
TX
59632
2021-09-06 00:00:00
TX
4953
2021-09-07 00:00:00
TX
14616
2021-09-08 00:00:00
TX
47091
2021-09-09 00:00:00
TX
45523
2021-09-10 00:00:00
TX
52431
2021-09-11 00:00:00
TX
54300
2021-09-12 00:00:00
TX
tx_cases["submission_date"] = tx_cases["submission_date"] - pd.to_timedelta(7, unit='d')
week_cases = tx_cases.groupby(["state", pd.Grouper(key="submission_date", freq='W-MON')])["new_case"].sum().reset_index().sort_values("submission_date")
cases_date_range = (week_cases["submission_date"] > "2020-1-1") & (week_cases["submission_date"] <= "2022-1-1")
week_cases = week_cases.loc[cases_date_range]
tx_gas_m = tx_gas_c.merge(week_cases, left_on="date", right_on="submission_date", how="left")
tx_gas_m = tx_gas_m.drop(columns=["state", "submission_date"])
tx_gas_m = tx_gas_m.fillna(0)
tx_gas_m
datedatetime64[ns]
2020-01-06 00:00:00 - 2021-12-27 00:00:00
gas_all_gradesfloat64
1.506 - 3.134
50
2020-12-21 00:00:00
1.984
51
2020-12-28 00:00:00
1.971
52
2021-01-04 00:00:00
1.962
53
2021-01-11 00:00:00
2.057
54
2021-01-18 00:00:00
2.173
55
2021-01-25 00:00:00
2.149
56
2021-02-01 00:00:00
2.154
57
2021-02-08 00:00:00
2.209
58
2021-02-15 00:00:00
2.268
59
2021-02-22 00:00:00
2.403
Causal Inference Modelling
!pip install statsmodels==0.13.5
Requirement already satisfied: statsmodels==0.13.5 in /root/venv/lib/python3.9/site-packages (0.13.5)
Requirement already satisfied: patsy>=0.5.2 in /root/venv/lib/python3.9/site-packages (from statsmodels==0.13.5) (0.5.3)
Requirement already satisfied: scipy>=1.3 in /shared-libs/python3.9/py/lib/python3.9/site-packages (from statsmodels==0.13.5) (1.9.3)
Requirement already satisfied: packaging>=21.3 in /shared-libs/python3.9/py-core/lib/python3.9/site-packages (from statsmodels==0.13.5) (21.3)
Requirement already satisfied: pandas>=0.25 in /shared-libs/python3.9/py/lib/python3.9/site-packages (from statsmodels==0.13.5) (1.2.5)
Requirement already satisfied: numpy>=1.17 in /shared-libs/python3.9/py/lib/python3.9/site-packages (from statsmodels==0.13.5) (1.23.4)
Requirement already satisfied: pyparsing!=3.0.5,>=2.0.2 in /shared-libs/python3.9/py-core/lib/python3.9/site-packages (from packaging>=21.3->statsmodels==0.13.5) (3.0.9)
Requirement already satisfied: python-dateutil>=2.7.3 in /shared-libs/python3.9/py-core/lib/python3.9/site-packages (from pandas>=0.25->statsmodels==0.13.5) (2.8.2)
Requirement already satisfied: pytz>=2017.3 in /shared-libs/python3.9/py/lib/python3.9/site-packages (from pandas>=0.25->statsmodels==0.13.5) (2022.5)
Requirement already satisfied: six in /shared-libs/python3.9/py-core/lib/python3.9/site-packages (from patsy>=0.5.2->statsmodels==0.13.5) (1.16.0)
WARNING: You are using pip version 22.0.4; however, version 22.3.1 is available.
You should consider upgrading via the '/root/venv/bin/python -m pip install --upgrade pip' command.
import statsmodels.api as sm
import seaborn as sns
import itertools
def fit_OLS_model(df, target_variable, explanatory_variables, intercept = False):
target = df[target_variable]
inputs = df[explanatory_variables]
if intercept:
inputs = sm.add_constant(inputs)
fitted_model = sm.OLS(target, inputs).fit()
return(fitted_model)
def mean_squared_error(true_vals, predicted_vals):
return np.mean((true_vals - predicted_vals) ** 2)
ca_gas_m2 = ca_gas_m
ca_gas_m2["lockdown"] = ca_gas_m2["lockdown"].astype(int)
ca_gas_m2
datedatetime64[ns]
2020-01-06 00:00:00 - 2021-12-27 00:00:00
gas_all_gradesfloat64
2.721 - 4.642
0
2020-01-06 00:00:00
3.505
1
2020-01-13 00:00:00
3.497
2
2020-01-20 00:00:00
3.487
3
2020-01-27 00:00:00
3.467
4
2020-02-03 00:00:00
3.455
5
2020-02-10 00:00:00
3.449
6
2020-02-17 00:00:00
3.441
7
2020-02-24 00:00:00
3.443
8
2020-03-02 00:00:00
3.435
9
2020-03-09 00:00:00
3.403
sdo = ca_gas_m2[ca_gas_m2['lockdown']==1]['gas_all_grades'].mean() - ca_gas_m2[ca_gas_m2['lockdown']==0]['gas_all_grades'].mean()
sdo
linear_model = fit_OLS_model(ca_gas_m2, 'gas_all_grades', ['lockdown', 'new_case'])
print(linear_model.summary())
OLS Regression Results
=======================================================================================
Dep. Variable: gas_all_grades R-squared (uncentered): 0.577
Model: OLS Adj. R-squared (uncentered): 0.569
Method: Least Squares F-statistic: 69.61
Date: Tue, 13 Dec 2022 Prob (F-statistic): 8.61e-20
Time: 05:36:56 Log-Likelihood: -237.82
No. Observations: 104 AIC: 479.6
Df Residuals: 102 BIC: 484.9
Df Model: 2
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
lockdown 2.5434 0.353 7.196 0.000 1.842 3.244
new_case 1.342e-05 3.33e-06 4.032 0.000 6.82e-06 2e-05
==============================================================================
Omnibus: 4.560 Durbin-Watson: 0.045
Prob(Omnibus): 0.102 Jarque-Bera (JB): 2.805
Skew: -0.200 Prob(JB): 0.246
Kurtosis: 2.302 Cond. No. 1.26e+05
==============================================================================
Notes:
[1] R² is computed without centering (uncentered) since the model does not contain a constant.
[2] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[3] The condition number is large, 1.26e+05. This might indicate that there are
strong multicollinearity or other numerical problems.
IPW Calculation
import sklearn
from sklearn.linear_model import LogisticRegression as LR
import matplotlib.pyplot as plt
def ipw_calc(state):
lr = LR(penalty='none', max_iter=200, random_state=0)
treatment = state["lockdown"]
outcome = state["gas_all_grades"]
X = state[["new_case"]]
lr.fit(X, treatment)
p_score = pd.DataFrame(data={"new_case": state["new_case"],
"treatment": state["lockdown"],
"outcome": state["gas_all_grades"],
"propensity_score": lr.predict_proba(X)[:, 1]})
p_score = p_score[(p_score["propensity_score"] >= 0.05) & (p_score["propensity_score"] <= 0.95)]
treat_1_df = p_score.loc[p_score["treatment"] == 1]
treat_0_df = p_score.loc[p_score["treatment"] == 0]
ate = ((1 / len(p_score)) * sum(treat_1_df["outcome"] / treat_1_df["propensity_score"])) - ((1 / len(p_score)) * sum(treat_0_df["outcome"] / (1 - treat_0_df["propensity_score"])))
return ate
ate_list = {"California": ipw_calc(ca_gas_m),
"Florida": ipw_calc(fl_gas_m),
"New York": ipw_calc(ny_gas_m),
"Texas": ipw_calc(tx_gas_m)}
ate_list
#Propensity Score Plot for California
plt.hist(ipw_calc(ca_gas_m)[ipw_calc(ca_gas_m)['treatment'] == 1]['propensity_score'], label = "Lockdown(Treatment = 1)",ls='dashed', alpha = 0.5, lw=3, color= 'b')
plt.hist(ipw_calc(ca_gas_m)[ipw_calc(ca_gas_m)['treatment'] == 0]['propensity_score'], label = "No Lockdown(Treatment = 0)",ls='dotted', alpha = 0.5, lw=3, color= 'r')
plt.legend(loc='upper right')
plt.title("Propensity score with Treatment in California")
#Propensity Score Plot for Florida
plt.hist(ipw_calc(fl_gas_m)[ipw_calc(fl_gas_m)['treatment'] == 1]['propensity_score'], label = "Lockdown(Treatment = 1)",ls='dashed', alpha = 0.5, lw=3, color= 'b')
plt.hist(ipw_calc(fl_gas_m)[ipw_calc(fl_gas_m)['treatment'] == 0]['propensity_score'], label = "No Lockdown(Treatment = 0)",ls='dashed', alpha = 0.5, lw=3, color= 'r')
plt.legend(loc='upper right')
plt.title("Propensity score with Treatment in Florida")
#Propensity Score Plot for New York
plt.hist(ipw_calc(ny_gas_m)[ipw_calc(ny_gas_m)['treatment'] == 1]['propensity_score'], label = "Lockdown(Treatment = 1)",ls='dashed', alpha = 0.5, lw=3, color= 'b')
plt.hist(ipw_calc(ny_gas_m)[ipw_calc(ny_gas_m)['treatment'] == 0]['propensity_score'], label = "No Lockdown(Treatment = 0)",ls='dashed', alpha = 0.5, lw=3, color= 'r')
plt.legend(loc='upper right')
plt.title("Propensity score with Treatment in New York")
#Propensity Score Plot for Texas
plt.hist(ipw_calc(tx_gas_m)[ipw_calc(tx_gas_m)['treatment'] == 1]['propensity_score'], label = "Lockdown(Treatment = 1)",ls='dashed', alpha = 0.5, lw=3, color= 'b')
plt.hist(ipw_calc(tx_gas_m)[ipw_calc(tx_gas_m)['treatment'] == 0]['propensity_score'], label = "No Lockdown(Treatment = 0)",ls='dashed', alpha = 0.5, lw=3, color= 'r')
plt.legend(loc='upper right')
plt.title("Propensity score with Treatment in Texas")