###########################################################################
# PUBLG100: Introduction to Quantitative Methods
#
# Week 8 Seminar: Panel Data, Time-Series Cross-Section Models
#
#

# Set your working directory
#
# CAUTION: Make sure the directory you specify here matches the working directory on your computer.
#          We're using N:/PUBLG100 only for illustration purposes and it would only work if you're 
#          using a UCL dekstop. If you're using your own laptop, then replace N:/PUBLG100 with the 
#          appropriate directory (or folder)
setwd("N:/PUBLG100")


# Verify that your working directory is set correctly
getwd()


## ----eval = FALSE--------------------------------------------------------
## install.packages("plm")

## ----message = FALSE-----------------------------------------------------
library(plm)
library(lmtest)
library(texreg)

## ------------------------------------------------------------------------
rm(list = ls())

## ------------------------------------------------------------------------
guns <- read.csv("guns.csv")

## ------------------------------------------------------------------------
guns$shall <- factor(guns$shall, levels = c(0, 1), labels =c("NO", "YES"))

## ------------------------------------------------------------------------
state_effects <- plm(mur ~ shall + incarc_rate + pm1029, 
                     data = guns, 
                     index = c("stateid", "year"), 
                     model = "within", 
                     effect = "individual")

summary(state_effects)

## ------------------------------------------------------------------------
plmtest(state_effects, effect="individual")

## ------------------------------------------------------------------------
time_effects <- plm(mur ~ shall + incarc_rate + pm1029, 
                    data = guns, 
                    index = c("stateid", "year"), 
                    model = "within", 
                    effect = "time")

summary(time_effects)

## ------------------------------------------------------------------------
plmtest(time_effects, effect="time")

## ------------------------------------------------------------------------
twoway_effects <- plm(mur ~ shall + incarc_rate + pm1029, 
                      data = guns, 
                      index = c("stateid", "year"), 
                      model = "within", 
                      effect = "twoways")

summary(twoway_effects)

## ------------------------------------------------------------------------
screenreg(list(state_effects, time_effects, twoway_effects), 
          custom.model.names = c("State Fixed Effects", 
                                 "Time Fixed Effects", 
                                 "Twoway Fixed Effects"))

## ------------------------------------------------------------------------
pbgtest(twoway_effects)

## ------------------------------------------------------------------------
twoway_effects_hac <- coeftest(twoway_effects, 
                               vcov = vcovHC(twoway_effects, 
                                             method = "arellano", 
                                             type = "HC3"))

screenreg(list(twoway_effects, twoway_effects_hac),
          custom.model.names = c("Twoway Fixed Effects", 
                                 "Twoway Fixed Effects (HAC)"))

## ------------------------------------------------------------------------
pcdtest(twoway_effects)

## ------------------------------------------------------------------------
twoway_effects_pcse <- coeftest(twoway_effects, 
                                vcov = vcovBK(twoway_effects, 
                                              type="HC3", 
                                              cluster = "group")) 

twoway_effects_pcse

## ------------------------------------------------------------------------
twoway_effects_scc <- coeftest(twoway_effects,
                               vcov = vcovSCC(twoway_effects, 
                                              type="HC3", 
                                              cluster = "group"))

twoway_effects_scc

## ------------------------------------------------------------------------
screenreg(list(state_effects, 
               time_effects, 
               twoway_effects, 
               twoway_effects_pcse, 
               twoway_effects_scc), 
          custom.model.names = c("State Effects", 
                                 "Time Effects", 
                                 "Twoway Fixed Effects", 
                                 "PCSE", 
                                 "SCC"))