# JW June 30, 2014 
library(mgcv)
library(foreign)
library(arm)
library(MASS)

rm(list=ls(all=TRUE))

setwd("C:/Users/Julian/Dropbox/TooLittleTooLate/data")
source("../analyses/mod.summary.R")

ra <- read.dta("autonomy_RR.dta")


RA <- subset(ra, isrelevant == 1 & status_monopoly == 0 & status_dominant==0 
	& !(pys_family==0 & family_warhist>0) & terr_relevant==1)

f1 <- as.formula(onset_do_terr_flag ~ egip +status_autonomy +warhist1 +b +b2 +excl_groups_count 
	+federal +ln_rgdppc_lag +ln_pop +c_incidence_flagl +pys_family 
	+inc_spline1 +inc_spline2 +inc_spline3)

f2 <- as.formula(onset_do_terr_flag ~ egip +status_autonomy + egipXwarhist + autonomyXwarhist +warhist1 +b +b2 +excl_groups_count 
	+federal +ln_rgdppc_lag +ln_pop +c_incidence_flagl +pys_family 
	+inc_spline1 +inc_spline2 +inc_spline3)

f3 <- as.formula(onset_do_terr_flag ~ egip_aut +egip_naut +status_autonomy +egip_autXwarhist +egip_nautXwarhist +autonomyXwarhist +warhist1 +b +b2 +excl_groups_count 
	+federal +ln_rgdppc_lag +ln_pop +c_incidence_flagl +pys_family 
	+inc_spline1 +inc_spline2 +inc_spline3)

f4 <- as.formula(onset_do_terr_flag ~ egip_aut +egip_naut +status_autonomy +egip_autXwarhist +egip_nautXwarhist +autonomyXwarhist + cent_aut2 +other_downgraded2+warhist1 +b +b2 +excl_groups_count 
	+federal +ln_rgdppc_lag +ln_pop +c_incidence_flagl +pys_family 
	+inc_spline1 +inc_spline2 +inc_spline3)


m1 <- glm(f1, data= RA,
	family=binomial(link="logit") )
m2 <- glm(f2, data= RA,
	family=binomial(link="logit") )
m3 <- glm(f3, data= RA,
	family=binomial(link="logit") )
m4 <- glm(f4, data= RA,
	family=binomial(link="logit") )

mod.summary(m1, white.adjust="cluster", cluster=ra$cowcode)
mod.summary(m2, white.adjust="cluster", cluster=ra$cowcode)
mod.summary(m3, white.adjust="cluster", cluster=ra$cowcode)
mod.summary(m4, white.adjust="cluster", cluster=ra$cowcode)


#draw simulations
k <- 1000

su <- mod.summary(m3, white.adjust="cluster", cluster=ra$cowcode)
pars <- m3$coefficients
vcov <- su$vcov  # clustered standard errors
m3.sim <- mvrnorm(k, pars, vcov)


# Relevant data
d3 <- as.data.frame(model.matrix(m3))

###########
# PRE WAR
###########
d3.pre <- subset(d3, egip_aut==0 & egip_naut==0 & status_autonomy==0 & warhist1==0) 


# (0) BASELINE: no autonomy or inclusion
d3.0 <- as.matrix(d3.pre)
pred.0.pre <- array(NA, dim=c(dim(m3.sim)[1] ,nrow(d3.0)))
for (i in 1:nrow(d3.0)){
  pred.0.pre[,i] <- invlogit(
    m3.sim %*% d3.0[i,]     
	)
}

# (1) AUTONOMY
d3.1 <- d3.pre
d3.1$status_autonomy <- 1
d3.1 <- as.matrix(d3.1)

pred.1.pre <- array(NA, dim=c(dim(m3.sim)[1] ,nrow(d3.1)))
for (i in 1:nrow(d3.1)){
  pred.1.pre[,i] <- invlogit(
    m3.sim %*% d3.1[i,]     
	)
}

# (2) EGIP w/o Autonomy
d3.2 <- d3.pre
d3.2$egip_naut <- 1
d3.2 <- as.matrix(d3.2)

pred.2.pre <- array(NA, dim=c(dim(m3.sim)[1] ,nrow(d3.2)))
for (i in 1:nrow(d3.2)){
  pred.2.pre[,i] <- invlogit(
    m3.sim %*% d3.2[i,]     
	)
}

# (3) EGIP w Autonomy
d3.3 <- d3.pre
d3.3$egip_aut <- 1
d3.3 <- as.matrix(d3.3)

pred.3.pre <- array(NA, dim=c(dim(m3.sim)[1] ,nrow(d3.3)))
for (i in 1:nrow(d3.3)){
  pred.3.pre[,i] <- invlogit(
    m3.sim %*% d3.3[i,]     
	)
}


############
# POST WAR
############
d3.post <- subset(d3, egip_aut==0 & egip_naut==0 & status_autonomy==0 & warhist1==1) 

# (0) BASELINE: no autonomy or inclusion
d3.0 <- as.matrix(d3.post)
pred.0.post <- array(NA, dim=c(dim(m3.sim)[1] ,nrow(d3.0)))
for (i in 1:nrow(d3.0)){
  pred.0.post[,i] <- invlogit(
    m3.sim %*% d3.0[i,]     
	)
}

# (1) AUTONOMY
d3.1 <- d3.post
d3.1$status_autonomy <- 1
d3.1$autonomyXwarhist <- 1
d3.1 <- as.matrix(d3.1)

pred.1.post <- array(NA, dim=c(dim(m3.sim)[1] ,nrow(d3.1)))
for (i in 1:nrow(d3.1)){
  pred.1.post[,i] <- invlogit(
    m3.sim %*% d3.1[i,]     
	)
}

# (2) EGIP w/o Autonomy
d3.2 <- d3.post
d3.2$egip_naut <- 1
d3.2$egip_nautXwarhist <- 1
d3.2 <- as.matrix(d3.2)

pred.2.post <- array(NA, dim=c(dim(m3.sim)[1] ,nrow(d3.2)))
for (i in 1:nrow(d3.2)){
  pred.2.post[,i] <- invlogit(
    m3.sim %*% d3.2[i,]     
	)
}

# (3) EGIP w Autonomy
d3.3 <- d3.post
d3.3$egip_aut <- 1
d3.3$egip_autXwarhist <- 1
d3.3 <- as.matrix(d3.3)

pred.3.post <- array(NA, dim=c(dim(m3.sim)[1] ,nrow(d3.3)))
for (i in 1:nrow(d3.3)){
  pred.3.post[,i] <- invlogit(
    m3.sim %*% d3.3[i,]     
	)
}


resmp <- array(NA, dim=c(k,8))
for (i in 1:k){
  resmp[i,1] <- mean(pred.0.pre[i,],na.rm=T)
  resmp[i,2] <- mean(pred.1.pre[i,],na.rm=T)
  resmp[i,3] <- mean(pred.2.pre[i,],na.rm=T)
  resmp[i,4] <- mean(pred.3.pre[i,],na.rm=T)
  resmp[i,5] <- mean(pred.0.post[i,],na.rm=T)
  resmp[i,6] <- mean(pred.1.post[i,],na.rm=T)
  resmp[i,7] <- mean(pred.2.post[i,],na.rm=T)
  resmp[i,8] <- mean(pred.3.post[i,],na.rm=T)	
}


pdf("../R&R/Latex/model3_plot_cl2.pdf", width=5*1.668, height=5, onefile=TRUE)
#par(mfrow=c(2,1)) 
plot(density(resmp[,2]-resmp[,1]), zero.line = FALSE, main="", #Pre Conflict", 
     xlab="Change in probability of conflict onset", xlim=c(-0.05,0.02), ylim=c(0,500), lwd=3, las=1, bty='n', lty=3,
)
lines(density(resmp[,3]-resmp[,1]), lty=2, lwd=3)
lines(density(resmp[,4]-resmp[,1]), lty=1, lwd=3)
abline(v=0)
legend("topleft", #inset=.02, # title="Access to Power",
  	 #c("Autonomy", "Included", "Included & Autonomy"),
  	 c("Regional autonomy", "Inclusion at the center", "Inclusion and autonomy"),  
     lty=3:1, lwd=3, bty="n")

plot(density(resmp[,6]-resmp[,5]), zero.line = FALSE, main="", #Post Conflict", 
     xlab="Change in probability of conflict onset", xlim=c(-0.05,0.02), ylim=c(0,80), lwd=3, las=1, bty='n', lty=3,
)
lines(density(resmp[,7]-resmp[,5]), lty=2, lwd=3)
lines(density(resmp[,8]-resmp[,5]), lty=1, lwd=3)
abline(v=0)
legend("topleft", #inset=.02, # title="Access to Power",
  	 #c("Autonomy", "Included", "Included & Autonomy"),  
  	 c("Regional autonomy", "Inclusion at the center", "Inclusion and autonomy"),  
  	 lty=3:1, lwd=3, bty="n")
# 
# plot(density( (resmp[,2]-resmp[,1])- (resmp[,6]-resmp[,5]) ), main="Change from Pre to Post Conflict", 
#      xlab="change in probability of conflict onset", xlim=c(-0.04,0.04), ylim=c(0,100), lwd=3, las=1, bty='l', lty=3,
# )
# lines(density( (resmp[,3]-resmp[,1]) -(resmp[,7]-resmp[,5]) ), lty=2, lwd=3)
# lines(density( (resmp[,4]-resmp[,1]) -(resmp[,8]-resmp[,5]) ), lty=1, lwd=3)
# abline(v=0)
# legend("topleft", inset=.02, title="Access to Power",
#   	 c("Autonomy", "Included", "Included + Autonomy"),  lty=3:1, lwd=3, bty="n")
 dev.off()





