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

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

## ------------------------------------------------------------------------
cpds <- read.dta("http://uclspp.github.io/PUBLG100/data/CPDS_1960-2013_stata.dta")

## ------------------------------------------------------------------------
cpds <- select(cpds,
               country,
               year,
               rae_ele,
               vturn,
               judrev,
               ud,
               unemp,
               unemp_pmp,
               realgdpgr,
               debt,
               socexp_t_pmp)

## ------------------------------------------------------------------------
hist(cpds$rae_ele, 
     xlab = "Electoral Fractionalization Index (rae_ele)",
     main = "Histogram of Electoral Fractionalization")

## ------------------------------------------------------------------------
range(cpds$vturn, na.rm = TRUE)

## ----echo = FALSE--------------------------------------------------------
# THIS CODE IS HERE FOR DISPLAYING THE CONTENT ON THE COURSE WEBSITE. YOU DO NOT NEED TO WORRY ABOUT WHAT IT DOES
var_ranges <- apply(select(cpds,
                           vturn,
                           ud,
                           unemp,
                           unemp_pmp,
                           realgdpgr,
                           debt,
                           socexp_t_pmp), 2, function(x) {round(range(x, na.rm = TRUE), digits = 2)})
row.names(var_ranges) = c("min", "max")
var_ranges

min_round <- function(x) {
  round(min(x, na.rm = TRUE), digits = 2)
}

max_round <- function(x) {
  round(max(x, na.rm = TRUE), digits = 2)
}

q_round <- function(x, ...) {
  round(quantile(x, na.rm = TRUE, ...), digits = 2)
}

coef_round <- function(model, variable, digits = 3) {
  round(coef(model)[variable], digits = digits)
}

coef_abs <- function(model, variable, digits = 3) {
  round(abs(coef(model)[variable]), digits = digits)
}

coef_pct <- function(model, variable) {
  paste(round(abs(coef(model)[variable] / diff(range(cpds$rae_ele, na.rm = TRUE))) * 100, digits = 2), " %")
}

## ------------------------------------------------------------------------
range(cpds$rae_ele, na.rm = TRUE)
summary(cpds$rae_ele)

## ------------------------------------------------------------------------
country_effects <- plm(rae_ele ~ vturn + judrev + ud + unemp + unemp_pmp + realgdpgr + debt + socexp_t_pmp,
                       data = cpds,
                       index = c("country", "year"), 
                       model = "within", 
                       effect = "individual")  

summary(country_effects)

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

## ------------------------------------------------------------------------
-3.9645e-02 / diff(range(cpds$rae_ele, na.rm = TRUE))

## ------------------------------------------------------------------------
coef(country_effects) / diff(range(cpds$rae_ele, na.rm = TRUE))

## ------------------------------------------------------------------------
time_effects <- plm(rae_ele ~ vturn + judrev + ud + unemp + unemp_pmp + realgdpgr + debt + socexp_t_pmp,
                    data = cpds,
                    index = c("country", "year"), 
                    model = "within", 
                    effect = "time")  

plmtest(time_effects, effect="time")

## ------------------------------------------------------------------------
summary(time_effects)

## ------------------------------------------------------------------------
coef(time_effects) / diff(range(cpds$rae_ele, na.rm = TRUE))

## ------------------------------------------------------------------------
twoway_effects <- plm(rae_ele ~ vturn + judrev + ud + unemp + unemp_pmp + realgdpgr + debt + socexp_t_pmp,
                      data = cpds,
                      index = c("country", "year"), 
                      model = "within", 
                      effect = "twoway")  

summary(twoway_effects)

## ------------------------------------------------------------------------
coef(twoway_effects) / diff(range(cpds$rae_ele, na.rm = TRUE))

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

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

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

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

twoway_effects_hac

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

twoway_effects_scc

## ------------------------------------------------------------------------
dynamic_model <- 
  plm(rae_ele ~ lag(rae_ele) + vturn + judrev + ud + unemp + unemp_pmp + realgdpgr + debt + socexp_t_pmp, 
      data = cpds, 
      index = c("country", "year"), 
      model = "within", 
      effect = "twoways")

## ------------------------------------------------------------------------
pbgtest(dynamic_model)

## ------------------------------------------------------------------------
pcdtest(dynamic_model)

## ------------------------------------------------------------------------
dynamic_model_scc <- coeftest(dynamic_model, 
                              vcov = vcovSCC(dynamic_model, type="HC3", cluster = "group"))

## ------------------------------------------------------------------------
dynamic_model_scc

## ------------------------------------------------------------------------
screenreg(list(twoway_effects, twoway_effects_hac, twoway_effects_scc),
          custom.model.names = c("Twoway", "Twoway (HAC)", "Twoway (SCC)"))

## ------------------------------------------------------------------------
screenreg(list(country_effects, twoway_effects, twoway_effects_hac, twoway_effects_scc, dynamic_model_scc),
          custom.model.names = c("Country Effects", "Twoway Effects", "Twoway (HAC)", "Twoway (SCC)", "Dynamic"))