rm(list=ls())
library(stargazer)
library(data.table)
library(ggplot2)
library(lfe)
library(rgdal)
library(rgeos)
library(gridExtra)
library(scales)
library(zoo)
library(foreign)
library(fst)
library(stringr)
library(ggrepel)
library(lubridate)
library(tidycensus)
library(arrow)
library(dplyr)
library(purrr)
library(jsonlite)
library(quantmod)
library(tidyverse)
library(R.matlab)
# library(directlabels)
fips <- fread("fips.csv")
names(fips) <- "fips"

states <- floor(fips/1000)

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

statefips <- statefips[,c("statefips","statecode")]

states <- merge(states,statefips,by.x="fips",by.y="statefips")
statenames <- c('AL','AK','AZ','AR','CA','CO','CT','DE','DC','FL','GA','ID','IL','IN','IA','KS','KY','LA','ME','MD','MA','MI','MN','MS','MO','MT','NE','NV','NH','NJ','NM','NY','NC','ND','OH','OK','OR','PA','RI','SC','SD','TN','TX','UT','VT','VA','WA','WV','WI','WY')

bmat <- fread("bmaterra.csv")
names(bmat) <- paste0("week",c("01","02","03","04","05","06","07","08","09","10","11","12"))

statefe <- bmat[1:50,]

statefe <- cbind(statenames,statefe)
names(statefe) <- c("State", paste0("week",c("01","02","03","04","05","06","07","08","09","10","11","12")))

statefe <- melt(statefe,id.vars = "State")

exovars <- bmat[51:55,]

exovars[,var:=c("pop","child","hs","white","inc")]

exovars <- melt(exovars,id.var="var")
ggplot(exovars,aes(x=variable,y=value,color=var,group=var))+geom_line()+theme_minimal()+theme(legend.position="bottom",legend.title = element_blank())+ xlab("")+ylab("")+labs(title="Exogenous Variables (Dep. Var: Pct Home)",subtitle="SDEM Estimates")

# ggplot(statefe,aes(x=variable,y=value)) + geom_text(aes(label=State),hju st=0, vjust=0,size=2)+theme_minimal()

ggplot(statefe,aes(x=variable,y=value))+geom_boxplot()+ geom_text(aes(label=State),hjust=1.5, vjust=0,size=1.5,check_overlap = TRUE)+theme_minimal()+ xlab("")+ylab("")+labs(title="State Fixed Effects (Dep. Var: Pct Home)",subtitle="SDEM Estimates")

fips <- fread("fips.csv")
names(fips) <- "fips"

xsub <- fread("xsub.csv")
names(xsub) <- c("pop","child","hs","white","inc")


pcthome <- fread("pcthome.csv")
pcthome <- log(pcthome)
pcthome <- pcthome - pcthome$V1
names(pcthome) <- paste0("pcthome_week",c("00","01","02","03","04","05","06","07","08","09","10","11","12"))


df <- cbind(pcthome,fips,xsub)
df[,state:=floor(fips/1000)]

df <- merge(df,statefips,by.x="state",by.y="statefips")
statefe <- NULL
coefs <- NULL

for(w in c("01","02","03","04","05","06","07","08","09","10","11","12") ) {
  reg <- felm(as.formula(paste0("pcthome_week",w,"~pop+child+hs+white+inc|statecode")),data=df)
  
      
  # t <- summary(reg)
  # rsq <-c(rsq,t$adj.r.squared)
  # print(t$adj.r.squared)
  
  temp <- data.frame(getfe(reg)[,c("effect","idx")])
  temp$effect <- as.numeric(temp$effect)
  temp$idx <- as.character(temp$idx)
  names(temp) <- c(paste0("pcthome_week",w),"state")
  
  temp2 <- data.frame(reg$coefficients)
  temp2['var'] <- row.names(temp2)
  row.names(temp2) <- NULL
  
  if(is.null(statefe)) {
    statefe = temp
    coefs <- temp2
  } else{
    statefe <- merge(statefe,temp,by="state")
    coefs <- merge(coefs,temp2,by="var")
  }
  

}

statefe <- melt(statefe,id.vars = "state")
coefs <- melt(coefs,id.vars = "var")
# ggplot(statefe,aes(x=variable,y=value)) + geom_text(aes(label=state),hjust=0, vjust=0,size=2)+theme_minimal()



ggplot(coefs,aes(x=substr(variable,13,14),y=value,color=var,group=var))+geom_line()+theme_minimal()+theme(legend.position="bottom",legend.title = element_blank())+ xlab("Week")+ylab("")+labs(title="Exogenous Variables (Dep. Var: Pct Home)",subtitle="OLS Estimates")

ggplot(statefe,aes(x=substr(variable,13,14),y=value))+geom_boxplot()+ geom_text(aes(label=state),hjust=1.5, vjust=0,size=1.5,check_overlap = TRUE)+theme_minimal()+theme(legend.position="bottom",legend.title = element_blank())+ xlab("Week")+ylab("")+labs(title="State Fixed Effects (Dep. Var: Pct Home)",subtitle="OLS Estimates")

# .printtable(reg)
zip_county <- fread("C:/Users/dratnadiwakara2/Documents/OneDrive - Louisiana State University/Raw Data/Crosswalk Files/ZIP_COUNTY_092016.csv")

zip_county<- zip_county[order(zip_county$ZIP,-zip_county$RES_RATIO)]

zip_county <- zip_county[!duplicated(zip_county$ZIP),c("ZIP","COUNTY")]
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")
core = read_fst(
  "C:/Users/dratnadiwakara2/Documents/OneDrive - Louisiana State University/Raw Data/COVID19/sg/core.fst",
  columns = c("safegraph_place_id","location_name","top_category","naics_code","latitude","longitude","region","postal_code","iso_country_code"),
  as.data.table = TRUE)
files  <- list.files(path = "C:/Users/dratnadiwakara2/Documents/OneDrive - Louisiana State University/Raw Data/COVID19/sg/weekly",pattern = '\\.fst',full.names = TRUE)


sg = lapply(files, read_fst, columns = c("safegraph_place_id","location_name","postal_code","naics_code","date_range_start","raw_visit_counts","raw_visitor_counts",paste0("visitsonday",1:7),"android","ios"),
  as.data.table = TRUE) 

sg <- do.call(rbind , sg)

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

sg = merge(sg,zip_county,by.x="postal_code",by.y="ZIP",all.x=TRUE)

sg[,week:=as.numeric(format(date_range_start,format="%Y%m%d"))]

sg[,place_category:=NA]


sg[,category:=ifelse(top_category %in% place_categories[['grocery']],"grocery",
                         ifelse(top_category %in% place_categories[['restaurants']],"restaurants",
                               ifelse(top_category %in% place_categories[['entertainment1']],"entertainment",NA )))]
wkstart <- c("3_01","3_08","3_15","3_22","3_29","4_05","4_12","4_19","4_26")

category_county_week = sg[,.(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=.(category,week,COUNTY)]

category_county_week <- category_county_week[!is.na(category) & !is.na(COUNTY)]

category_county_week[,visits:=v1+v2+v3+v4+v5+v6+v7]
category_county_week[,c("v1","v2","v3","v4","v5","v6","v7"):=list(NULL)]

grocery <- dcast(category_county_week[category=="grocery",c("COUNTY","week","visits")],COUNTY~week,value.var = "visits")
grocery[,2:ncol(grocery)] <- log(1+grocery[,2:ncol(grocery)] )
grocery[,2:ncol(grocery)] <- grocery[,2:ncol(grocery)] - grocery$`20200301`
names(grocery) <- c("fips",paste0("grocery_week",wkstart))

df <- merge(df,grocery,by="fips",all.x=T)



restaurants <- dcast(category_county_week[category=="restaurants",c("COUNTY","week","visits")],COUNTY~week,value.var = "visits")
restaurants[,2:ncol(restaurants)] <- log(1+restaurants[,2:ncol(restaurants)] )
restaurants[,2:ncol(restaurants)] <- restaurants[,2:ncol(restaurants)] - restaurants$`20200301`
names(restaurants) <- c("fips",paste0("restaur_week",wkstart))

df <- merge(df,restaurants,by="fips",all.x=T)


entertainment <- dcast(category_county_week[category=="entertainment",c("COUNTY","week","visits")],COUNTY~week,value.var = "visits")
entertainment[,2:ncol(entertainment)] <- log(1+entertainment[,2:ncol(entertainment)] )
entertainment[,2:ncol(entertainment)] <- entertainment[,2:ncol(entertainment)] - entertainment$`20200301`
names(entertainment) <- c("fips",paste0("enterta_week",wkstart))

df <- merge(df,entertainment,by="fips",all.x=T)
cats <- c("grocery","restaur","enterta") 
titles <- c("Grocery Visits","Restaurant Visits","Entertainment visits")

  i=1
for(cat in cats) {
  statefe <- NULL
  coefs <- NULL
  rsq <- NULL
  gr1<- list()
  gr2<-list()
  
  

  print(titles[i])
  for(w in wkstart[2:length(wkstart)] ) {
    reg <- felm(as.formula(paste0(cat,"_week",w,"~pop+child+hs+white+inc|statecode")),data=df)
    
    t <- summary(reg)
    rsq <-c(rsq,t$adj.r.squared)
    
    temp <- data.frame(getfe(reg)[,c("effect","idx")])
    temp$effect <- as.numeric(temp$effect)
    temp$idx <- as.character(temp$idx)
    names(temp) <- c(paste0("grocery_week",w),"state")
    
    temp2 <- data.frame(reg$coefficients)
    temp2['var'] <- row.names(temp2)
    row.names(temp2) <- NULL
    
    if(is.null(statefe)) {
      statefe = temp
      coefs <- temp2
    } else{
      statefe <- merge(statefe,temp,by="state")
      coefs <- merge(coefs,temp2,by="var")
    }
    
  
  }
  
  statefe <- melt(statefe,id.vars = "state")
  coefs <- melt(coefs,id.vars = "var")
  # ggplot(statefe,aes(x=variable,y=value)) + geom_text(aes(label=state),hjust=0, vjust=0,size=2)+theme_minimal()
  
  
  
  print(gr1[[i]] <- ggplot(coefs,aes(x=substr(variable,13,16),y=value,color=var,group=var))+geom_line(size=1)+geom_point(size=2)+theme_minimal()+theme(legend.position="bottom",legend.title = element_blank())+ xlab("Week")+ylab("")+labs(title=paste0("Exogenous Variables (Dep. Var: ",titles[i],")"),subtitle="OLS Estimates"))
  
   print(gr2[[i]] <- ggplot(statefe,aes(x=substr(variable,13,16),y=value))+geom_boxplot()+ geom_text(aes(label=state),hjust=1.5, vjust=0,size=2,check_overlap = TRUE)+theme_minimal()+theme(legend.position="bottom",legend.title = element_blank())+ xlab("Week")+ylab("")+labs(title=paste0("State Fixed Effects (Dep. Var: ",titles[i],")"),subtitle="OLS Estimates"))
  
  i=i+1
}
## [1] "Grocery Visits"

## [1] "Restaurant Visits"

## [1] "Entertainment visits"