rm(list=ls())
library(fst)
library(ggplot2)
library(data.table)
library(stargazer)
library(lfe)
library(tidycensus)
library(R.matlab)
library(rgdal)
library(rgeos)
unempins <- fread("C:/Users/dratnadiwakara2/Documents/OneDrive - Louisiana State University/Raw Data/COVID19/Unemployment Insurance.csv")

fp <- fread("C:/Users/dratnadiwakara2/Documents/OneDrive - Louisiana State University/Raw Data/Crosswalk Files/us-state-ansi-fips.csv")

unempins <- merge(unempins,fp[,c("statefips","statecode")],by.x="state",by.y="statecode")
unempins <- unempins[,c("statefips","basebenefit","weeklywage","state")]
names(unempins) <- c("state","basebenefit","weeklywage","statecode")
unempins[,replacementrate:=(basebenefit+600)/weeklywage]
unempins[,growthrate:=((600)/basebenefit)]

unempins[,state:=ifelse(nchar(as.character(state))==1,paste0("0",state),as.character(state))]
acs1 <- R.matlab::readMat("C:/Users/dratnadiwakara2/Documents/OneDrive - Louisiana State University/Raw Data/COVID19/mat/acs.mat")
median_age <- unlist(acs1$DP05.0018E,use.names = F)
pctwhite <- unlist(acs1$DP05.0037PE, use.names = F)
child <- unlist(acs1$DP02.0013PE, use.names = F)
medianincome <- unlist(acs1$DP03.0062E, use.names = F)
incom1 <- unlist(acs1$DP03.0052PE, use.names = F)
incom2 <- unlist(acs1$DP03.0053PE, use.names = F)
incom3 <- unlist(acs1$DP03.0054PE, use.names = F)
incom4 <- unlist(acs1$DP03.0055PE, use.names = F)
incom5 <- unlist(acs1$DP03.0056PE, use.names = F)

pctless35k <- incom1+incom2+incom3+incom4
pctless50k <- incom1+incom2+incom3+incom4+incom5

pctwithinsurance <- unlist(acs1$DP03.0096PE, use.names = F)

hhsize <- unlist(acs1$DP02.0015E, use.names = F)

acs1 <- data.frame(median_age=median_age,pctwhite=pctwhite,child=child,medianincome=medianincome,pctless35k=pctless35k,pctless50k=pctless50k,pctwithinsurance=pctwithinsurance,hhsize=hhsize)

spatial <- R.matlab::readMat('C:/Users/dratnadiwakara2/Documents/OneDrive - Louisiana State University/Raw Data/COVID19/mat/spatial_data.mat')
fips <- unlist(spatial$FIPS, use.names=FALSE)
pop <- unlist(spatial$Pop, use.names=FALSE)
land <- unlist(spatial$Land,use.names = FALSE)

pcthome <- unlist(spatial$PctHome,use.names = F)
pcthome <- data.table(pcthome)
names(pcthome) <- paste0("week",c("01","02","03","04","05","06","07","08","09",as.character(10:15)))

acs2 <- data.frame(fips=as.numeric(fips),pop=pop,land=land)

acskp <- cbind(acs1,acs2)
acskp <- data.table(acskp)
acskp[,logdensity:=log(pop/land)]

acskp[,fips:=ifelse(nchar(fips)==4,paste0("0",fips),as.character(fips))]
load(file="C:/Users/dratnadiwakara2/Documents/OneDrive - Louisiana State University/Raw Data/Presidential Election/countypres_2000-2016.RData")
election <- x
rm(x)

election <- data.table(election)
election <- election[election$year %in% c(2016)& election$party=="republican"]
election <- election[,c("year","FIPS","candidatevotes","totalvotes")]
election[,rep_pct:=election$candidatevotes/election$totalvotes]

election[,FIPS:=as.character(FIPS)]
election[,FIPS:=ifelse(nchar(FIPS)==4,paste0("0",FIPS),FIPS)]

1 Descriptives

1.1 Unemployment Benefits and Wages

Due to the additional $600 of unemployment benefits provided by CARES Act, in some states unemployed receive more money than they would have typically earned in their jobs.

The graph below shows the replacement rate after CARES (red) and replacement rate before CARES (blue)

1.1.1 Point Plot

test <- unempins
test[,replacementratebeforecares:= basebenefit/weeklywage]
test <- test[,c("statecode","replacementrate","replacementratebeforecares")]
test <- test[order(-test$replacementrate)]

test <- melt(test,"statecode")

test <- test[order(-value)]


ggplot(test) + geom_point(aes(x=statecode,y=value,color=variable),size=4)+theme_minimal()+theme(legend.position = "bottom")+geom_hline(yintercept = 1)

(Source: https://www.nytimes.com/interactive/2020/04/23/business/economy/unemployment-benefits-stimulus-coronavirus.html)

1.1.2 Map

us_states <- readOGR("C:/Users/dratnadiwakara2/Documents/OneDrive - Louisiana State University/Raw Data/Shapefiles/US States","cb_2014_us_state_20m")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\dratnadiwakara2\Documents\OneDrive - Louisiana State University\Raw Data\Shapefiles\US States", layer: "cb_2014_us_state_20m"
## with 52 features
## It has 9 fields
## Integer64 fields read as strings:  ALAND AWATER
us_states <- fortify(us_states,region="GEOID")

us_states <- data.table(us_states)

us_states <- merge(us_states,unempins[,c("replacementrate","state")],by.x="id",by.y="state")


ggplot()+
  geom_polygon(data=us_states[! us_states$id %in% c("02","15","72")], aes(x=long,y=lat,group=group,fill=replacementrate),color="black")+
  scale_fill_gradientn(colors=c("orange","orange3","orangered1","orangered3","orangered4"))+  theme_minimal()+
  theme(axis.title=element_blank(),
        axis.text=element_blank(),
        axis.ticks=element_blank(),
        legend.position = "bottom",panel.grid = element_blank())+ guides(fill=guide_legend(title="Replacement Rate"))

2 Empirical Specification

\[ Outcome_w = \Sigma_w\text{ } \beta_w \times w \times log(ReplacementRate) + County\text{ }Fixed\text{ }Effects + Week\text{ }Fixed\text{ }Effects\]
where \(w\) is a dummy variable that indicates the week.

All the figures below plot the \(\beta_w\) estimates with the corresponding 95% confidence intervals.



3 SG Data (Pct Home)

files  <- list.files(path = "C:/Users/dratnadiwakara2/Documents/OneDrive - Louisiana State University/Raw Data/COVID19/social distancing",pattern = '\\.fst',full.names = TRUE)


sd = lapply(files, read_fst, as.data.table = TRUE) 
sd <- do.call(rbind , sd)

sd[,county:=substr(origin_census_block_group,1,5)]

sd[,dayno:=as.numeric(day-as.Date("2019-12-31"))]
sd[,week:=floor(dayno/7)]

sd <- sd[,.(device_count=sum(device_count,na.rm = T),
            distance_traveled_from_home=mean(distance_traveled_from_home,na.rm=T),
            completely_home_device_count=sum(completely_home_device_count,na.rm = T),
            part_time_work_behavior_devices=sum(part_time_work_behavior_devices,na.rm = T),
            full_time_work_behavior_devices=sum(full_time_work_behavior_devices,na.rm = T)),
         by=.(week,county)]


sd[,state:=substr(county,1,2)]
# sd[,county:=substr(tract,1,5)]

sd[,pcthome:=completely_home_device_count*100/device_count]
sd[,pctparttime:=part_time_work_behavior_devices*100/device_count]
sd[,pctfulltime:=full_time_work_behavior_devices*100/device_count]

ftdata <- sd

sd <- merge(sd,unempins,by="state")

sd <- merge(sd,acskp,by.x="county",by.y="fips")
sd <- merge(sd,election[,c("FIPS","rep_pct")],by.x="county",by.y="FIPS")

sd[,weekno := week]
sd[,week:=as.factor(week)]

sd[,week:= relevel(sd$week,9)]

3.1 Pct Home

r <- list()
r[[1]] <- felm(log(1+pcthome)~log(replacementrate)*week|county+week|0|county,data=sd[sd$medianincome>30000 & sd$medianincome<50000  ])
r[[2]] <- felm(log(1+pcthome)~log(replacementrate)*week|county+week|0|county,data=sd[sd$medianincome>=50000 ])

.coef_plot_2reg(r[[1]],"Income less than 50k",r[[2]],"Income greater than 50k","log(replacementrate):week",8)+ geom_vline(xintercept = 9,color = "grey",linetype="dashed")

4 SG Data (Types of Visits)

core = read_fst(
  "C:/Users/dratnadiwakara2/Documents/OneDrive - Louisiana State University/Raw Data/COVID19/sg/core.fst",
  columns = c("safegraph_place_id","top_category"),  as.data.table = TRUE)
files  <- list.files(path = "C:/Users/dratnadiwakara2/Documents/OneDrive - Louisiana State University/Raw Data/COVID19/sg/weekly/v2",pattern = '\\.fst',full.names = TRUE)

weekly = lapply(files, read_fst, as.data.table = TRUE) 

weekly <- do.call(rbind , weekly)

weekly = merge(weekly,core[,c("safegraph_place_id","top_category")],by="safegraph_place_id")

weekly[,county:=substr(poi_cbg,1,5)]
place_categories <- list()
place_categories[['grocery']] = c("General Merchandise Stores, including Warehouse Clubs and Supercenters","Grocery Stores")
# gasolineother = c("Gasoline Stations","Health and Personal Care Stores")
# place_categories[['gasolineother']] = c("Gasoline Stations","Health and Personal Care Stores")
place_categories[['restaurants']] =c("Restaurants and Other Eating Places")
place_categories[['entertainment1']] = c("Museums, Historical Sites, and Similar Institutions","Gambling Industries","Motion Picture and Video Industries","Spectator Sports")#,"Department Stores"
# place_categories[['entertainment2']] = c("Sporting Goods, Hobby, and Musical Instrument Stores","Amusement Parks and Arcades")
# place_categories[['healthservices']] = c("General Medical and Surgical Hospitals","Offices of Other Health Practitioners")
# place_categories[['airports']] = c("Support Activities for Air Transportation")
category_county_week = weekly[,.(v1=sum(visitsonday1,na.rm=TRUE),
                                 v2=sum(visitsonday2,na.rm=TRUE),
                                 v3=sum(visitsonday3,na.rm=TRUE),
                                 v4=sum(visitsonday4,na.rm=TRUE),
                                 v5=sum(visitsonday5,na.rm=TRUE),
                                 v6=sum(visitsonday6,na.rm=TRUE),
                                 v7=sum(visitsonday7,na.rm=TRUE)),
                              by=.(top_category,date_range_start,county)]

category_county_week[,visits:=v1+v2+v3+v4+v5+v6+v7]

category_county_week[,dayno:=as.numeric(date_range_start-as.Date("2019-12-31"))]
category_county_week[,week:=floor(dayno/7)]

category_county_week <- category_county_week[week>=0]

category_county_week[,state:=substr(county,1,2)]

category_county_week <- merge(category_county_week,unempins,by="state")

category_county_week <- merge(category_county_week,acskp,by.x="county",by.y="fips")
category_county_week <- merge(category_county_week,election[,c("FIPS","rep_pct")],by.x="county",by.y="FIPS")

category_county_week[,weekno := week]
category_county_week[,week:=as.factor(week)]

category_county_week[,week:= relevel(category_county_week$week,9)]

4.1 Difference-in-difference: log(Replacement Rate) x Post CARES

4.1.1 Grocery

i=1 
r <- list()
r[[1]] <- felm(log(1+visits)~log(replacementrate)*week|county+week|0|county,data=category_county_week[ top_category %in% as.vector(place_categories[[i]]) ])


.coef_plot_1reg(r[[1]],"log(replacementrate):week",8)+ geom_vline(xintercept = 9,color = "grey",linetype="dashed")

4.1.2 Restaurants

i=2 
r <- list()
r[[1]] <- felm(log(1+visits)~log(replacementrate)*week|county+week|0|county,data=category_county_week[ top_category %in% as.vector(place_categories[[i]]) ])


.coef_plot_1reg(r[[1]],"log(replacementrate):week",8)+ geom_vline(xintercept = 9,color = "grey",linetype="dashed")

4.1.3 Entertainment

i=3
r <- list()
r[[1]] <- felm(log(1+visits)~log(replacementrate)*week|county+week|0|county,data=category_county_week[ top_category %in% as.vector(place_categories[[i]]) ])


.coef_plot_1reg(r[[1]],"log(replacementrate):week",8)+ geom_vline(xintercept = 9,color = "grey",linetype="dashed")

4.2 Heterogeneity: Split x log(Replacement Rate) x Post CARES

4.2.1 By Pct Less than 35k

4.2.1.1 Grocery

i=1 
r <- list()
r[[1]] <- felm(log(1+visits)~log(pctless35k)*log(replacementrate)*week|county+week|0|county,data=category_county_week[ top_category %in% as.vector(place_categories[[i]]) ])


.coef_plot_1reg(r[[1]],"log(pctless35k):log(replacementrate):week",8)+ geom_vline(xintercept = 9,color = "grey",linetype="dashed")

4.2.1.2 Restaurants

i=2
r <- list()
r[[1]] <- felm(log(1+visits)~log(pctless35k)*log(replacementrate)*week|county+week|0|county,data=category_county_week[ top_category %in% as.vector(place_categories[[i]]) ])


.coef_plot_1reg(r[[1]],"log(pctless35k):log(replacementrate):week",8)+ geom_vline(xintercept = 9,color = "grey",linetype="dashed")

4.2.1.3 Entertainment

i=3
r <- list()
r[[1]] <- felm(log(1+visits)~log(pctless35k)*log(replacementrate)*week|county+week|0|county,data=category_county_week[ top_category %in% as.vector(place_categories[[i]]) ])


.coef_plot_1reg(r[[1]],"log(pctless35k):log(replacementrate):week",8)+ geom_vline(xintercept = 9,color = "grey",linetype="dashed")

4.2.2 By Household Size

4.2.2.1 Grocery

i=1 
r <- list()
r[[1]] <- felm(log(1+visits)~log(hhsize)*log(replacementrate)*week|county+week|0|county,data=category_county_week[ top_category %in% as.vector(place_categories[[i]]) ])


.coef_plot_1reg(r[[1]],"log(hhsize):log(replacementrate):week",8)+ geom_vline(xintercept = 9,color = "grey",linetype="dashed")

4.2.2.2 Restaurants

i=2
r <- list()
r[[1]] <- felm(log(1+visits)~log(hhsize)*log(replacementrate)*week|county+week|0|county,data=category_county_week[ top_category %in% as.vector(place_categories[[i]]) ])


.coef_plot_1reg(r[[1]],"log(hhsize):log(replacementrate):week",8)+ geom_vline(xintercept = 9,color = "grey",linetype="dashed")

4.2.2.3 Entertainment

i=3
r <- list()
r[[1]] <- felm(log(1+visits)~log(hhsize)*log(replacementrate)*week|county+week|0|county,data=category_county_week[ top_category %in% as.vector(place_categories[[i]]) ])


.coef_plot_1reg(r[[1]],"log(hhsize):log(replacementrate):week",8)+ geom_vline(xintercept = 9,color = "grey",linetype="dashed")

4.2.3 By Device Type (android x high replacement x post CARES)

category_county_week = weekly[,.(android=sum(android,na.rm=TRUE),
                                 ios=sum(ios,na.rm=TRUE)),
                              by=.(top_category,date_range_start,county)]


category_county_week[,dayno:=as.numeric(date_range_start-as.Date("2019-12-31"))]
category_county_week[,week:=floor(dayno/7)]

category_county_week <- category_county_week[week>=0]

category_county_week[,state:=substr(county,1,2)]

category_county_week <- merge(category_county_week,unempins,by="state")

category_county_week <- merge(category_county_week,acskp,by.x="county",by.y="fips")
category_county_week <- merge(category_county_week,election[,c("FIPS","rep_pct")],by.x="county",by.y="FIPS")

category_county_week[,weekno := week]
category_county_week[,week:=as.factor(week)]

category_county_week[,week:= relevel(category_county_week$week,9)]


category_county_week <- category_county_week[,c("android","ios","replacementrate","week","county","top_category")]

category_county_week <- melt(category_county_week,id.vars = c("replacementrate","week","county","top_category"))
category_county_week[,top_category:=as.character(top_category)]
category_county_week[,android:=ifelse(as.character(variable)=="ios",0,1)]

4.2.3.1 Grocery

i=1
r <- felm(log(1+value)~android*log(replacementrate)*week|county+week|0|county,data=category_county_week[ top_category %in% as.vector(place_categories[[i]]) ])
print(.coef_plot_1reg(r,"android:log(replacementrate):week",8))

4.2.3.2 Restaurants

i=2
r <- felm(log(1+value)~android*log(replacementrate)*week|county+week|0|county,data=category_county_week[ top_category %in% as.vector(place_categories[[i]]) ])
print(.coef_plot_1reg(r,"android:log(replacementrate):week",8))

4.2.3.3 Entertainment

i=3
r <- felm(log(1+value)~android*log(replacementrate)*week|county+week|0|county,data=category_county_week[ top_category %in% as.vector(place_categories[[i]]) ])
print(.coef_plot_1reg(r,"android:log(replacementrate):week",8))

5 Number of New COVID-19 Cases

CSSE <- fread("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv")

CSSE <- CSSE[CSSE$FIPS>=1000]
CSSE[,c("UID","iso2","iso3","code3","Admin2","Province_State","Country_Region","Lat","Long_","Combined_Key"):=list(NULL)]

CSSE <- melt(data = CSSE,id.vars = c("FIPS"))
names(CSSE) <- c("fips","date","cases")
CSSE[,date:=as.character(CSSE$date)]
CSSE[,date:=as.Date(CSSE$date,"%m/%d/%y")]

CSSE <- CSSE[order(CSSE$fips,CSSE$date)]
CSSE[, lag.cases:=c(NA, cases[-.N]), by=fips]
CSSE[,newcases:=cases-lag.cases]



CSSE_dealth <- fread("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_US.csv")

CSSE_dealth <- CSSE_dealth[CSSE_dealth$FIPS>=1000]
CSSE_dealth[,c("UID","iso2","iso3","code3","Admin2","Province_State","Country_Region","Lat","Long_","Combined_Key"):=list(NULL)]

CSSE_dealth <- melt(data = CSSE_dealth,id.vars = c("FIPS"))
names(CSSE_dealth) <- c("fips","date","deaths")
CSSE_dealth[,date:=as.character(CSSE_dealth$date)]
CSSE_dealth[,date:=as.Date(CSSE_dealth$date,"%m/%d/%y")]
CSSE_dealth <- CSSE_dealth[order(CSSE_dealth$fips,CSSE_dealth$date)]
CSSE_dealth[, lag.deaths:=c(NA, deaths[-.N]), by=fips]
CSSE_dealth[,newdeaths:=deaths-lag.deaths]


CSSE <- merge(CSSE,CSSE_dealth,by=c("fips","date"),all.x=T)


CSSE[,county:=ifelse(nchar(fips)==4,paste0("0",fips),paste0(fips))]
CSSE[,state:=substr(county,1,2)]

CSSE[,dayno:=as.numeric(date-as.Date("2019-12-31"))]
CSSE[,week:=floor(dayno/7)]

CSSE <- merge(CSSE,unempins,by="state")

5.1 County-Level Cases

CSSE_w <- CSSE[,.(newcases=sum(newcases,na.rm = T),
                  newdeaths=sum(newdeaths,na.rm=T),
                  replacementrate=mean(replacementrate,na.rm=T)),
               by=.(county,week)]

CSSE_w[,weekno := week]
CSSE_w[,week:=as.factor(week)]

CSSE_w[,week:= relevel(CSSE_w$week,10)]


CSSE_w <- merge(CSSE_w,acskp,by.x="county",by.y="fips")

CSSE_w[,state:=substr(county,1,2)]
r <- felm(log(1+newcases)~log(replacementrate)*week|county+week|0|county,data=CSSE_w[newcases>=0 & weekno>12 ])
.coef_plot_1reg(r,"log(replacementrate):week",13)

5.2 County-Level Heterogeneity

5.2.1 By Pct Less than 35k

r <- felm(log(1+newcases)~log(pctless35k)*log(replacementrate)*week|county+week|0|county,data=CSSE_w[newcases>=0 & weekno>12 ])
.coef_plot_1reg(r,"log(pctless35k):log(replacementrate):week",13)

5.2.2 By Population Density

r <- felm(log(1+newcases)~logdensity*log(replacementrate)*week|county+week|0|county,data=CSSE_w[newcases>=0 & weekno>12 ])
.coef_plot_1reg(r,"logdensity:log(replacementrate):week",13)

5.3 State-Level Cases

CSSE_w <- CSSE[,.(newcases=sum(newcases,na.rm = T),
                  newdeaths=sum(newdeaths,na.rm=T),
                  replacementrate=mean(replacementrate,na.rm=T)),
               by=.(state,week)]

CSSE_w[,weekno := week]
CSSE_w[,week:=as.factor(week)]

CSSE_w[,week:= relevel(CSSE_w$week,10)]
r <- felm(log(1+newcases)~log(replacementrate)*week|state+week|0|state,data=CSSE_w[newcases>=0 & weekno>10  ])
.coef_plot_1reg(r,"log(replacementrate):week",12)