##Attacks Data

#install the following packasges. Then comment them out.
install.packages("arm") #this installs a package for summarizing regression results. after you do it once, you don't have to do it again. 
install.packages("stargazer") #this installs stargazer
install.packages("sandwich")
install.packages("lmtest")
install.packages("multiwayvcov")
install.packages("latest")
install.packages("AER")
install.packages("systemfit")

#load the packages
library(stargazer)  
library(arm)
library(xlsx)
library(sandwich)
library(foreign)
library(lmtest)
library(multiwayvcov)
library(AER)
library(systemfit)

rm(list=ls())  #this clears your data window
setwd("~/Desktop/THESIS") #this sets your working directory. You need to set it to 
#wherever you want to store your data. 


mydat <- read.csv("3.27attack_final.csv", header=TRUE)#load the packages

names(mydat) # this lists the variables you have in the dataframe. 

hist(mydat$lji,
     main="Histogram of Latent Judicial Independence",
     xlab = "Latent Judicial Independence",
     ylab = "Frequency",
     border = "black",
     col = "blue",
     las=1)

hist(mydat$civil,
     main="Histogram of Civil Society Index ",
     xlab = "Civil Society Index",
     ylab = "Frequency",
     border = "black",
     col = "blue",
     las=1)

#This is the model without any interactions or covariates
fit.1  <- lm(lji ~ safe + civil, data=mydat)
m1.vcovCL<-cluster.vcov(fit.1, mydat$cowcode)
robust_se.1 <-sqrt(diag(m1.vcovCL))
summary(fit.1)

fit.2 <- lm(lji ~ safe + civil + safeXcivil, data=mydat)
m2.vcovCL<-cluster.vcov(fit.2, mydat$cowcode)
robust_se2.2 <-sqrt(diag(m2.vcovCL))
summary(fit.2)

fit.3 <- lm(lji ~ safe + civil + safeXcivil + democracy + lngdp, data=mydat)
m3.vcovCL<-cluster.vcov(fit.3, mydat$cowcode)
robust_se2.3 <-sqrt(diag(m3.vcovCL))
summary(fit.3)

fit.4 <- lm(lji ~ safe + civil + safeXcivil + democracy + lngdp + as.factor(cowcode), data=mydat)
m4.vcovCL<-cluster.vcov(fit.4, mydat$cowcode)
robust_se2.4 <-sqrt(diag(m4.vcovCL))
summary(fit.4)


fit.5 <- lm(lji ~ safe_2 + civil_2 + safeXcivil_2 + democracy_2 + lngdp_2 + as.factor(cowcode), data=mydat)
m5.vcovCL<-cluster.vcov(fit.5, mydat$cowcode)
robust_se2.5 <-sqrt(diag(m5.vcovCL))
summary(fit.5)

fit.6 <- lm(lji ~ safe_1 + civil_1 + safeXcivil_1 + democracy_1 + lngdp_1 + as.factor(cowcode), data=mydat)
m6.vcovCL<-cluster.vcov(fit.6, mydat$cowcode)
robust_se2.6 <-sqrt(diag(m6.vcovCL))
summary(fit.6)

fit.7 <- lm(lji ~ safe_5 + civil_5 + safeXcivil_5 + democracy_5 + lngdp_5 + as.factor(cowcode), data=mydat)
m7.vcovCL<-cluster.vcov(fit.7, mydat$cowcode)
robust_se2.7 <-sqrt(diag(m7.vcovCL))
summary(fit.7)

#Here is a table with cluster corrected standard errors. Compare it to the next table. 
stargazer(fit.1, fit.2,fit.3,fit.4, type="text",
         se = list(robust_se.1, robust_se2.2, robust_se2.3, robust_se2.4),
          omit.stat = "f",
          add.lines = list(c("State fixed effects?", "No", "No", "No", "Yes")),
          covariate.labels = c("Judicial Safety Index", "Civil Society Index", "Judicial Safety X Civil Society",
                               "Demoracy Indicator", "Ln(GDP/cap)","Constant"),
          dep.var.caption = "Latent Judicial Independence",
          dep.var.labels = "",
          omit="cowcode",
          title = "Results with Standard Errors Clustered on State")

#Here is a table with uncorrected standard errors. See how much smaller they are. The assumption that errors are uncorrelated within states 
#is very very likely to be wrong. So this table suggests overly optimistic results. 
stargazer(fit.1, fit.2,fit.3,fit.4, type="text",
          #          se = list(robust_se.1, robust_se2.2, robust_se2.3, robust_se2.4),
          omit.stat = "f",
          add.lines = list(c("State fixed effects?", "No", "No", "No", "Yes")),
          covariate.labels = c("Judicial Safety Index", "Civil Society Index", "Judicial Safety X Civil Society",
                               "Demoracy Indicator", "Ln(GDP/cap)","Constant"),
          dep.var.caption = "Latent Judicial Independence",
          dep.var.labels = "",
          omit="cowcode",
          title = "Results with OlS standard errors")


#Here is another table with cluster corrected standard errors with a TWO-year lag. Compare it to the next table. 
stargazer(fit.1, fit.2,fit.3,fit.4,fit.5, type="text",
          se = list(robust_se.1, robust_se2.2, robust_se2.3, robust_se2.4,robust_se2.5),
          omit.stat = "f",
          add.lines = list(c("State fixed effects?", "No", "No", "No", "Yes","Yes")),
          covariate.labels = c("Judicial Safety Index", "Civil Society Index", "Judicial Safety X Civil Society",
                               "Judicial Safety Index (2 yr lag)", "Civil Society Index  (2 yr lag)", 
                               "Judicial Safety X Civil Society  (2 yr lag)",
                               "Demoracy Indicator", "Ln(GDP/cap)","Constant"),
          dep.var.caption = "Latent Judicial Independence",
          dep.var.labels = "",
          omit="cowcode",
          title = "Results with Standard Errors Clustered on State")

#Table with cluster corrected standard errors with ONE year lag
stargazer(fit.1, fit.2,fit.3,fit.4,fit.6, type="text",
          se = list(robust_se.1, robust_se2.2, robust_se2.3, robust_se2.4,robust_se2.6),
          omit.stat = "f",
          add.lines = list(c("State fixed effects?", "No", "No", "No", "Yes","Yes")),
          covariate.labels = c("Judicial Safety Index", "Civil Society Index", "Judicial Safety X Civil Society",
                               "Judicial Safety Index (1 yr lag)", "Civil Society Index  (1 yr lag)", 
                               "Judicial Safety X Civil Society  (1 yr lag)",
                               "Demoracy Indicator", "Ln(GDP/cap)","Constant"),
          dep.var.caption = "Latent Judicial Independence",
          dep.var.labels = "",
          omit="cowcode",
          title = "Results with Standard Errors Clustered on State")

#Table with cluster corrected standard errors with FIVE year lag
stargazer(fit.1, fit.2,fit.3,fit.4,fit.7, type="text",
          se = list(robust_se.1, robust_se2.2, robust_se2.3, robust_se2.4,robust_se2.7),
          omit.stat = "f",
          add.lines = list(c("State fixed effects?", "No", "No", "No", "Yes","Yes")),
          covariate.labels = c("Judicial Safety Index", "Civil Society Index", "Judicial Safety X Civil Society",
                               "Judicial Safety Index (5 yr lag)", "Civil Society Index  (5 yr lag)", 
                               "Judicial Safety X Civil Society  (5 yr lag)",
                               "Demoracy Indicator", "Ln(GDP/cap)","Constant"),
          dep.var.caption = "Latent Judicial Independence",
          dep.var.labels = "",
          omit="cowcode",
          title = "Results with Standard Errors Clustered on State")


#Here is the Instrumental variables regression where we need to instrument for safe and safeXcivil, and where we include fixed effects.
#You'll see that we don't have strong results. 

iv1 <- ivreg(lji ~ safe + civil +safeXcivil+as.factor(cowcode)| break_down + breakXcivil + civil +as.factor(cowcode), data = mydat)
summary(iv1)

iv2 <- ivreg(lji ~ safe + civil +safeXcivil+as.factor(cowcode)| newpower + newXcivil + civil +as.factor(cowcode), data = mydat)
summary(iv2)


#Now, let's subset the data to see which countries had large attacks on the judiciary. Remember a big attack is a year in which the safe index
#dropped by 1. 
big_attack_mydat <- subset(mydat,bigattack==1,select=c(country_name, year,bigattack,break_down,civil)) 
big_attack_mydat


#Now, let's subset the data to see which countries had VERY large attacks on the judiciary. Remember a GIANT attack is a year in which the safe index
#dropped by 1.5. 
giant_attack_mydat <- subset(mydat,giantattack==1,select=c(country_name, year,bigattack,giantattack,break_down,civil)) 
giant_attack_mydat

#Finally, let's see which countries had giant attack with a robust civil society -- say something at least .7 on the scale.
big_attack_robustcivil_mydat <- subset(mydat,bigattack==1 & civil>.7,select=c(country_name, year,bigattack,giantattack,break_down,civil)) 
big_attack_robustcivil_mydat

#GIANT attack w/robust civil
GARC <- subset(mydat,giantattack==1 & civil>.7,select=c(country_name, year,bigattack,giantattack, break_down,civil)) 
#large attack w/robust civil
BARC <- subset(mydat,bigattack==1 & civil>.7,select=c(country_name, year,bigattack,giantattack, break_down,civil)) 


#number of total attacks
table(mydat$giantattack) #total giant attacks
table(mydat$bigattack) #total big attacks