
# Load packages
library(MASS)
library(foreign)
library(gplots)

# Set WD
#setwd("C:/Users/spadmin/PhD/Projects/ColRule/Simulation")
rm(list=ls(all=TRUE))

setwd("C:/Users/Julian/Dropbox/TooLittleTooLate/analyses")

## Import data
data <- read.delim("M5b-data.txt", sep = "\t")
data <- na.omit(data)


## Import estimates for naive conflict model (model 5)
beta_m5 <- read.delim("m5b-beta.txt", sep="\t")
vcv_m5 <- read.delim("m5b-vcv.txt", sep="\t")


## Import estimates for bivariate probit model (model 4)
beta_m6 <- read.delim("m6b-beta.txt", sep="\t")
vcv_m6 <- read.delim("m6b-vcv.txt", sep="\t")

vcv_m6 <- vcv_m6[1:13, 1:13]
beta_m6 <- beta_m6[1:13]



### Simulation

k <- 2000 # number of draws
prange <- c(0.025, 0.5, 0.925)

# replace DV with constant
data$onset_do_terr_flag <- 1



### Model 5: naive probit

draw <- mvrnorm(k, data.matrix(beta_m5), data.matrix(vcv_m5))

data0 <- data1 <- data
data0$aut <- 0
data1$aut <- 1

pred0 <- rep(NA, k)
pred1 <- rep(NA, k)
fd.naive <- rep(NA, k)

for (i in 1:k){
	pred0[i] <- mean(pnorm( as.matrix(data0) %*% draw[i,] ))
	pred1[i] <- mean(pnorm( as.matrix(data1) %*% draw[i,] ))
	fd.naive[i] <- pred1[i]-pred0[i]
}

yhat0.naive <- quantile(pred0, probs = prange, na.rm = TRUE)
yhat1.naive <- quantile(pred1, probs = prange, na.rm = TRUE )
fd.naive.quant <- quantile(fd.naive, probs = prange, na.rm = TRUE )




### Model 6: bivariate probit

draw <- mvrnorm(k, data.matrix(beta_m6), data.matrix(vcv_m6))


pred0 <- rep(NA, k)
pred1 <- rep(NA, k)
fd.biv <- rep(NA, k)

for (i in 1:k){
	pred0[i] <- mean(pnorm( as.matrix(data0) %*% draw[i,] ))
	pred1[i] <- mean(pnorm( as.matrix(data1) %*% draw[i,] ))
	fd.biv[i] <- pred1[i]-pred0[i]
}

yhat0.biv <- quantile(pred0, probs = prange )
yhat1.biv <- quantile(pred1, probs = prange )
fd.biv.quant <- quantile(fd.biv, probs = prange )

setwd("C:/Users/Julian/Dropbox/TooLittleTooLate/analyses")
pdf(file="biprobit_new.pdf", onefile=FALSE, paper="special", height=5, width=5*1.668)
#par(mar=c(4.1,4.1,2.1,2.1))
plot(density(fd.biv), main=" ", 
     xlab="Change in predicted probability of conflict onset", xlim=c(-.08,.03), ylim=c(0,100), lwd=3, las=1, bty='n', lty=1, zero.line = FALSE
)
lines(density(fd.naive), lty=2, lwd=3)
#text(-.2, 5.6, "Model 3")
#text(-.58, 6.6, "Model 4")
abline(v=0)
legend("topleft", inset=.05,	 c("Naive", "Instrumented"),  lty=2:1, lwd=3, bty="n")

dev.off()



