#########################
###  Participant Models
#########################
rm(list=ls())
data1 <- read.csv("main participant data.csv", header=T)

###  Main Results, Table 2
#  First Stage
first <- glm(onset ~ minpow + rival + contig_land, data=data1, family="binomial"(link="logit"), maxit=200)
data1$resid <- first$residuals

#  Main Paper, Table 2, Model 3
fit3 <- glm(esc_s2 ~ jme_new2 + atop_sy2 + jointdemo + cinc + unga + dy_trade + escalate_l + resid + factor(year), data=data1, family="binomial"(link="logit"), maxit=200)
summary(fit3)


#  Main Paper, Table 2, Model 4
fit4 <- glm(esc_s2 ~ jnally2 + jally_new2 + atop_sy2 + jointdemo + cinc + unga + dy_trade + escalate_l + resid + factor(year), data=data1, family="binomial"(link="logit"), maxit=200)
summary(fit4)


#  Heckman Selection Models
#  Main Paper, Table 3, Participants Model
library(sampleSelection)
heckvan4 <- heckit(onset~minpow + rival + contig_land, 
                   esc_s2 ~ jnally2 + jally_new2 + atop_sy2 + jointdemo + cinc + unga + dy_trade + escalate_l + factor(year),  data = data1)
summary(heckvan4)


#  Figure 2
library(ggplot2)
probabilities <- read.csv("pred_prob_part.csv", header=T)

ggplot(data=probabilities, aes(x=ally_jme, y=fit, group=1)) +
  geom_bar(stat='identity', fill="gray75") +
  geom_errorbar(aes(ymin = cilo, ymax = cihigh), width=0.5) + 
  xlab("Ally JME") + ylab("Predicted Probability of Participant Escalation (%)")

#Figure 3
data1$year73 <- 0
data1$year74 <- 0
data1$year75 <- 0
data1$year76 <- 0
data1$year77 <- 0
data1$year78 <- 0
data1$year79 <- 0
data1$year80 <- 0
data1$year81 <- 0
data1$year82 <- 0
data1$year83 <- 0
data1$year84 <- 0
data1$year85 <- 0
data1$year86 <- 0
data1$year87 <- 0
data1$year88 <- 0
data1$year89 <- 0
data1$year90 <- 0
data1$year91 <- 0
data1$year92 <- 0
data1$year93 <- 0
data1$year94 <- 0
data1$year95 <- 0
data1$year96 <- 0
data1$year97 <- 0
data1$year98 <- 0
data1$year99 <- 0
data1$year00 <- 0
data1$year01 <- 0
data1$year02 <- 0
data1$year03 <- 0
data1$year73[data1$year==1973] <- 1
data1$year74[data1$year==1974] <- 1
data1$year75[data1$year==1975] <- 1
data1$year76[data1$year==1976] <- 1
data1$year77[data1$year==1977] <- 1
data1$year78[data1$year==1978] <- 1
data1$year79[data1$year==1979] <- 1
data1$year80[data1$year==1980] <- 1
data1$year81[data1$year==1981] <- 1
data1$year82[data1$year==1982] <- 1
data1$year83[data1$year==1983] <- 1
data1$year84[data1$year==1984] <- 1
data1$year85[data1$year==1985] <- 1
data1$year86[data1$year==1986] <- 1
data1$year87[data1$year==1987] <- 1
data1$year88[data1$year==1988] <- 1
data1$year89[data1$year==1989] <- 1
data1$year90[data1$year==1990] <- 1
data1$year91[data1$year==1991] <- 1
data1$year92[data1$year==1992] <- 1
data1$year93[data1$year==1993] <- 1
data1$year94[data1$year==1994] <- 1
data1$year95[data1$year==1995] <- 1
data1$year96[data1$year==1996] <- 1
data1$year97[data1$year==1997] <- 1
data1$year98[data1$year==1998] <- 1
data1$year99[data1$year==1999] <- 1
data1$year00[data1$year==2000] <- 1
data1$year01[data1$year==2001] <- 1
data1$year02[data1$year==2002] <- 1
data1$year03[data1$year==2003] <- 1

heckvan4_prob <- heckit(onset~minpow + rival + contig_land, 
                        esc_s2 ~ jnally2 + jally_new2 + atop_sy2 + jointdemo + cinc + unga + dy_trade + escalate_l  + year74 + year75 + year76 + year77 + year78 + year79 + year80 + year81 + year82 + year83 + year84 + year85 + year86 + year87 + year88 + year89 + year90 + year91 + year92 + year93 + year94 + year95 + year96 + year97 + year98 + year99 + year00 + year01 + year02 + year03,  data = data1)
vcov( heckvan4_prob, part = "outcome" )
ndata <- data.frame(jally_new2=c(0,1), jnally2=0, minpow=mean(data1$minpow), rival=mean(data1$rival), contig_land=mean(data1$contig_land), atop_sy2=mean(data1$atop_sy2), jointdemo=mean(data1$jointdemo), cinc=mean(data1$cinc), unga=mean(data1$unga), dy_trade=mean(data1$dy_trade), escalate_l=mean(data1$escalate_l), 
                    year74=mean(data1$year74), year75=mean(data1$year75), year76=mean(data1$year76), year77=mean(data1$year77), year78=mean(data1$year78), year79=mean(data1$year79), year80=mean(data1$year80), year81=mean(data1$year81), 
                    year82=mean(data1$year82), year83=mean(data1$year83), year84=mean(data1$year84), year85=mean(data1$year86), year86=mean(data1$year86), year87=mean(data1$year87), year88=mean(data1$year88), year89=mean(data1$year89), year90=mean(data1$year90),
                    year91=mean(data1$year91), year92=mean(data1$year92), year93=mean(data1$year93), year94=mean(data1$year94), year95=mean(data1$year95), year96=mean(data1$year96), year97=mean(data1$year97), year98=mean(data1$year98), 
                    year99=mean(data1$year99), year00=mean(data1$year00), year01=mean(data1$year01), year02=mean(data1$year02), year03=mean(data1$year03))
#probabilities <- data.frame(predict(heckvan4, ndata, type="response", se.fit=TRUE))
probabilities2 <- predict(heckvan4_prob, ndata, type="conditional", part="outcome",se.fit=TRUE)
probabilities <- data.frame(predict(heckvan4_prob, ndata, type="conditional", part="outcome",se.fit=TRUE))
probabilities$ally_jme <- as.factor(c(0,1))
probabilities$fit <- probabilities$E.yo.ys.1 * 100
probabilities$se.fit <- .0001069191 * 100
probabilities$cihigh <- probabilities$fit + 1.96*probabilities$se.fit
probabilities$cilo <- probabilities$fit - 1.96*probabilities$se.fit 

ggplot(data=probabilities, aes(x=ally_jme, y=fit, group=1)) +
  geom_bar(stat='identity', fill="gray75") +
  geom_errorbar(aes(ymin = cilo, ymax = cihigh), width=0.5) + 
  xlab("Ally JME") + ylab("Predicted Probability of Participant Escalation (%)")


#  Figure 4
library(ggplot2)
library(scales)

bin <- read.csv("jcr2_part newvar 2.csv", header=T)
bin <- subset(bin, atop != 0)
bin$lo <- bin$mean - 1.96*bin$sd
bin$hi <- bin$mean + 1.96*bin$sd

test <- ggplot(data=bin, aes(x=atop, y=mean)) +
  geom_line()

test <- test + 
  geom_ribbon(aes(ymin = lo, ymax = hi), linetype=2, alpha=0.1) +
  xlab("Number of Alliances") + 
  ylab("Marginal Effect on Participant Escalation") + scale_y_continuous(labels = percent)
test


################################
###  Appendix Robustness Checks

#  Table A1:  Regional JMEs
regdata <- read.csv("regional participant data.csv", header=T)
fitrp <- glm(esc_s2 ~ rnajme + rjally + atop_sy2 + jointdemo + cinc + unga + dy_trade + escalate_l + resid, data=regdata, family="binomial"(link="logit"), maxit=200)
summary(fitrp)


#  Table A2:  Combining Reciprocity and Escalation
fitc <- glm(escpar ~ jnally2 + jally_new2 + atop_sy2 + jointdemo + cinc + unga + dy_trade + escalate_l + resid + factor(year), data=data1, family="binomial"(link="logit"), maxit=200)
summary(fitc)


#  Table A3:  Multinomial Logit  |  Warning!  This takes a while.
library(nnet)

#  Setting Descalate as baseline
data1$shift_new <- ifelse(data1$shift_s2 == 0, 5, data1$shift_s2)
data1$shift_new <- ifelse(data1$shift_new == 1, 0, data1$shift_new)
data1$shift_new <- ifelse(data1$shift_new == 5, 1, data1$shift_new)

rob2t <- multinom(shift_new ~ jnally2 + jally_new2 + atop_sy2 + jointdemo + cinc + unga + dy_trade + escalate_l + resid + factor(year), data=data1, maxit=300)
summary(rob2t)


#  Table A4:  Rare Events
library(Zelig)
fit2 <- zelig(esc_s2 ~ jme_new2 + jally_new2 + atop_sy2 + jointdemo + cinc + unga + dy_trade + escalate_l + resid + factor(year), data=data1, model="relogit", maxit=200)
summary(fit2)


#  Table A5:  K-Adic Correction
kdata <- read.csv("participant kadic data.csv", header=T)
kadic2 <- glm(esc_s2 ~ jnally2 + jally_new2 + atop + jointdemo + cinc + rival + majpow + contig_land + contig_sea, data=kdata, family="binomial"(link="logit"))
summary(kadic2)


#  Figure A1:  Leading JMEs
part <- read.csv("jcr2_lag_part_bin.csv", header=T)
part$hi <- part$coef + (part$error * 1.96)
part$lo <- part$coef - (part$error * 1.96)

pd <- position_dodge(0.2)

ggplot(part, aes(x=month, y=coef, color=Variable)) +
  geom_errorbar(aes(ymin=lo, ymax=hi), width=.1, position=pd) +
  geom_line(position=pd) +
  geom_point(position=pd) +
  xlab("Months until MID") +
  ylab("Partner")


#  Figure A2:  System Effects
xiao1 <- read.csv("sys effect part.csv", header=T)
ggplot(data=xiao1, aes(x=year, y=xiao, group=1)) +
  geom_line() +
  geom_point() +
  xlab("Year") + ylab("Total Predicted Reduction in Participant Conflict Escalation (%)") + scale_x_discrete(breaks=seq(1975,2005,5))