##############
# Panel Data #
##############
# part 1

rm(list = ls())
ad = "https://www.fsb.miamioh.edu/lij14/411_fatality.txt"
data = read.table(url(ad), header=T)  
head(data)
attach(data)
data[order(year,statename),][1:5,]

plot(beertax[year==1982],y[year==1982],main="1982 cross section",xlab="beertax",ylab="fatality")
abline(lm(y~beertax,subset=(year==1982)))
text(beertax[year==1982],y[year==1982], labels = statename[year==1982], pos = 1,cex = 0.8)

summary(lm(y~beertax,subset=(year==1982)))$coef
summary(lm(y~beertax))$coef

# part 2
plot(beertax[statename=="AL"],y[statename=="AL"],main="AL Time Series",xlab="beertax",ylab="fatality")
abline(lm(y~beertax,subset=(statename=="AL")))
text(beertax[statename=="AL"],y[statename=="AL"], labels = year[statename=="AL"], pos = 1,cex = 0.8)

summary(lm(y~beertax,subset=(statename=="AL")))$coef

# part 3
summary(lm(y~beertax+factor(statename)))$coef

femodel = lm(y~beertax+factor(statename))
poololsmodel = lm(y~beertax)
anova(poololsmodel,femodel)

summary(lm(y~beertax))$r.squared
summary(lm(y~beertax+factor(statename)))$r.squared

# part 4
summary(lm(y~beertax+factor(statename)+factor(year)))$coef
data$youngratio = pop1517/pop
summary(lm(y~beertax+youngratio+mlda+factor(statename)+factor(year),data=data))$coef

# part 5
library(plm)
data_panel = pdata.frame(data, index = c("statename", "year"))
summary(plm(y~beertax, data = data_panel, model = "within"))$coef

# optional advanced codes
mean(y[statename=="AL"])
mean(y[statename=="AR"])
da = data[order(statename,year),]
dy = da$y-rep(tapply(da$y, da$statename, mean),each=7)
dx = da$beertax-rep(tapply(da$beertax, da$statename, mean),each=7)
lm(dy~dx)

# part 6
library(sandwich)  
library(lmtest)
m = plm(y~beertax, data = data_panel, model = "within")
cluster_robust_se = vcovHC(m, type = "HC1", cluster = "group") 
coeftest(m, vcov. = cluster_robust_se)