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
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
# 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
# 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
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
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
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
# 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
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
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
# 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
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
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
# 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
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
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
Causal Inference Modelling
!pip install statsmodels==0.13.5
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
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())
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")