###########################################################################
# 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("CPDS_1960-2013_stata.dta")

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

## ------------------------------------------------------------------------
summary(cpds$rae_ele)

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

## ----echo = FALSE--------------------------------------------------------
# THIS CODE IS HERE FOR DISPLAYING THE CONTENT ON THE COURSE WEBSITE. YOU DO NOT NEED TO WORRY ABOUT WHAT IT DOES
coef_round <- function(model, variable) {
  sprintf("`%.3f`", coefficients(model)[variable])
}

coef_abs <- function(model, variable, digits = 3) {
  sprintf("`%.3f`", abs(coefficients(model)[variable]))
}

coef_pct <- function(model, variable) {
  sprintf("`%.2f%%`", abs(coefficients(model)[variable] / diff(range(cpds$rae_ele, na.rm = TRUE)) * 100))
}

## ------------------------------------------------------------------------
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))

## ------------------------------------------------------------------------
coefficients(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)

## ------------------------------------------------------------------------
coefficients(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)

## ------------------------------------------------------------------------
coefficients(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_hac <- coeftest(twoway_effects, 
                               vcov = vcovHC(twoway_effects, 
                                             method = "arellano", 
                                             type = "HC3"))

twoway_hac

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

twoway_scc

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

## ------------------------------------------------------------------------
screenreg(list(country_effects, twoway_effects, twoway_hac, twoway_scc),
          custom.model.names = c("Country Effects", 
                                 "Twoway Effects", 
                                 "Twoway (HAC)", 
                                 "Twoway (SCC)"))

## ------------------------------------------------------------------------
twoway_scc

## ------------------------------------------------------------------------
twoway_scc[,1] / diff(range(cpds$rae_ele, 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
coef_round <- function(model, variable) {
  sprintf("`%.3f`", model[,1][variable])
}

coef_abs <- function(model, variable, digits = 3) {
  sprintf("`%.3f`", abs(model[,1][variable]))
}

coef_pct <- function(model, variable) {
  sprintf("`%.2f%%`", abs( model[,1][variable] / diff(range(cpds$rae_ele, na.rm = TRUE)) * 100))
}