#Overall setup
rm(list = ls())                     # delete all objects in memory
library(foreign)                    # library for foreign (non-R) data files
#library(arm)

setwd("") # Define your working directory

qra=read.dta("koenigmaeder_strategiccompliance_ajps.dta")
attach(qra)


#Model 1
#u11:znumber_of_provisions,zdelegationratio,zefficient_score
#u13:const
#u14:zmsdis_sum,zmss_sum,zcouconf_sum,typodir,zinterest,numissues
#u23:0
#u24:const,zcomdis_sum,zcomss_sum,zcouconf_sum,typodir,zinterest,dummytime

depvar=cbind(i,c,d)         # dummies for dependent variable
indvar=cbind(
  znumber_of_provisions,zdelegationratio,zefficient_score,zmsdis_sum,zmss_sum,
  zcouconf_sum,typodir,zinterest,numissues,zcomgain_sum,
  zcomss_sum,zcouconf_sum,typodir,zinterest,dummytime
) 

# log-likelihood procedure (NOT the negative of it!)

llik=function(b,X,Y) {
  i=as.matrix(Y[,1])              # i, c, d are dep var dummies
  c=as.matrix(Y[,2])
  d=as.matrix(Y[,3])
  r=nrow(Y)
  o=rep(1,r)
  x11=as.matrix(cbind(X[,1],X[,2],X[,3]))     
  x13=as.matrix(cbind(o))                           
  x14=as.matrix(cbind(X[,4],X[,5],X[,6],X[,7],X[,8],X[,9]))               
  x24=as.matrix(cbind(o,X[,10],X[,11],X[,12],X[,13],X[,14],X[,15]))        
  b11=as.matrix(b[1:3])
  b13=as.matrix(b[4])                                          
  b14=as.matrix(b[5:10])                                           
  b24=as.matrix(b[11:17])                                       
  u11=x11 %*% b11
  u13=x13 %*% b13                                               
  u14=x14 %*% b14
  u24=x24 %*% b24
  p4=pnorm(u24/sqrt(2))              
  p3=1-p4                            
  p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1))  
  p1=1-p2                            
  lnpi=log(p1)                
  lnpc=log(p2*p3)
  lnpd=log(p2*p4)
  llik=i*lnpi+c*lnpc+d*lnpd  
  sum(llik)
}

# MLE using optim()
stval=runif(17)*0                
strat.mle=optim(stval,llik,hessian=TRUE,method="BFGS",control=list(fnscale=-1,trace=1),X=indvar,Y=depvar)



# Print results
#u11:znumber_of_provisions,zdelegationratio,zefficient_score
#u13:const
#u14:zmsdis_sum,zmss_sum,zcouconf_sum,typodir,zinterest,numissues
#u23:0
#u24:const,zcomdis_sum,zcomss_sum,zcouconf_sum,typodir,zinterest,dummytime

Variable=c(
  "u11:Number of major provisions","u11:Delegation ratio","u11:Bureaucratic efficiency",
  "u13:Constant",
  "u14:Member state's disagreement","u14:Member State's saliency","u14:Diversity of member states' interests","u14:Type of directive (1 - CM & EP)","u14:Interest group pattern","u14:Number of issues",
  "u24:Constant","u24:Commission's agreement","u24:Commission's saliency","u24:Diversity of member states' interests","u24:Type of directive (1 - CM & EP)","u24:Interest group pattern","u24:Issue of timeliness")
Parameter=as.matrix(strat.mle$par)
StdError=as.matrix(sqrt(diag(solve(-1 * strat.mle$hessian))))
Tstat=as.matrix(Parameter/StdError)
pvalue<-2*pt(abs(Tstat),length(i)-length(Parameter),lower.tail=FALSE)
strat.mle.table=cbind(format(Variable,justify="left",width=20),
                      prettyNum(Parameter,justify="right",digits=4),
                      prettyNum(StdError,justify="right",digits=4),
                      #prettyNum(Tstat,justify="right",digits=4),
                      prettyNum(pvalue,justify="right",digits=4))
cat("\n")
#Results=c("Variable","Parameter","StdError","t","pvalue");
Results=c("Variable","Parameter","StdError","pvalue");
print(format(rbind(Results,strat.mle.table),width=12,justify="right"))
cat("\nLog-Likelihood = ",strat.mle$value,"\n")
cat("N = ",length(i),"\n\n")


# Predicted Probabilities
n=299
const=rep(1,n)
x11=cbind(znumber_of_provisions,zdelegationratio,zefficient_score)
x13=const
x14=cbind(zmsdis_sum,zmss_sum,zcouconf_sum,typodir,zinterest,numissues)
x24=cbind(const,zcomgain_sum,zcomss_sum,zcouconf_sum,typodir,zinterest,dummytime)
b11=matrix(strat.mle$par[1:3],3,1)  
b13=matrix(strat.mle$par[4],1,1)    
b14=matrix(strat.mle$par[5:10],6,1)      
b24=matrix(strat.mle$par[11:17],7,1)     
u11=x11 %*% b11
u13=x13 %*% b13                     
u14=x14 %*% b14
u24=x24 %*% b24
#Probabilities
p4=pnorm(u24/sqrt(2))              
p3=1-p4                            
p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1)) 
p1=1-p2                            
pimpl=p1
pcap=p2*p3
pdis=p2*p4
#Create new variable depend
pdepend <- rep(NA,299)
pdepend[pcap>pimpl & pcap>pdis] <- 3
pdepend[pdis>pimpl & pdis>pcap] <- 4
pdepend[pimpl>pcap & pimpl>pdis] <- 1
#Comparison between observed and predicted
crosstable <- table(pdepend, depend)
crosstable
#Corretly predicted 
n=299
cpredicted=(crosstable[1,1]+crosstable[2,2]+crosstable[3,3])/n
print(cpredicted)



#DIFFERENT OUTCOME PROBABILITIES
n=299
const=rep(1,n)
#Summary of independent variables
summary(zmsdis_sum)
summary(zmss_sum)
summary(znumber_of_provisions)
summary(zdelegationratio)
summary(zefficient_score)
summary(zcouconf_sum)
summary(typodir)
summary(zinterest)
summary(numissues)
summary(zcomgain_sum)
summary(zcomss_sum)
summary(dummytime)

#Low values of indvar (1rd Qu.)
low_zmsdis_sum=rep(-0.6616,n)
low_zmss_sum=rep(-0.6626,n)
low_znumber_of_provisions=rep(-0.6192,n)
low_zdelegationratio=rep(-0.9686,n)
low_zefficient_score=rep(-0.344,n)
low_zcouconf_sum=rep(-0.7757,n)
low_typodir=rep(1,n)
low_zinterest=rep(-0.789,n)
low_numissues=rep(1,n)
low_zcomgain_sum=rep(-0.7757,n)
low_zcomss_sum=rep(-0.7642,n)
low_dummytime=rep(0,n)
#Mean values of indvar
mean_zmsdis_sum=rep(0,n)
mean_zmss_sum=rep(0,n)
mean_znumber_of_provisions=rep(0,n)
mean_zdelegationratio=rep(0,n)
mean_zefficient_score=rep(0,n)
mean_zcouconf_sum=rep(0,n)
mean_typodir=rep(1,n)
mean_zinterest=rep(0,n)
mean_numissues=rep(1.906,n)
mean_zcomgain_sum=rep(0,n)
mean_zcomss_sum=rep(0,n)
mean_dummytime=rep(0,n)
#Moderate values of indvar (3rd Qu.)
mode_zmsdis_sum=rep(0.3854,n)
mode_zmss_sum=rep(0.6348,n)
mode_znumber_of_provisions=rep(0.1397,n)
mode_zdelegationratio=rep(0.6642,n)
mode_zefficient_score=rep(0.7768,n)
mode_zcouconf_sum=rep(0.4183,n)
mode_typodir=rep(1,n)
mode_zinterest=rep(0.8686,n)
mode_numissues=rep(2,n)
mode_zcomgain_sum=rep(0.4586,n)
mode_zcomss_sum=rep(0.7777,n)
mode_dummytime=rep(0,n)



#X-VARIABLE: DIVERSITY OF STATES' INTERESTS
#zcouconf_sum
n=299
const=rep(1,n)
x=seq(-2,2.7,length=n)

# Pr(impl/dis/cap|zcouconf_sum) holding other indvar at min/mean/max values
#u11:znumber_of_provisions,zdelegationratio,zefficient_score
#u13:const
#u14:zmsdis_sum,zmss_sum,zcouconf_sum,typodir,zinterest,numissues
#u23:0
#u24:const,zcomdis_sum,zcomss_sum,zcouconf_sum,typodir,zinterest,dummytime

# Predicted Probabilities for LOW VALUES
n=299
const=rep(1,n)
x11=cbind(low_znumber_of_provisions,low_zdelegationratio,low_zefficient_score)
x13=const
x14=cbind(low_zmsdis_sum,low_zmss_sum,x,low_typodir,low_zinterest,low_numissues)
x24=cbind(const,low_zcomgain_sum,low_zcomss_sum,x,low_typodir,low_zinterest,low_dummytime)
b11=matrix(strat.mle$par[1:3],3,1)  
b13=matrix(strat.mle$par[4],1,1)    
b14=matrix(strat.mle$par[5:10],6,1)      
b24=matrix(strat.mle$par[11:17],7,1)     
u11=x11 %*% b11
u13=x13 %*% b13                     
u14=x14 %*% b14
u24=x24 %*% b24
#Probabilities
p4=pnorm(u24/sqrt(2))              
p3=1-p4                           
p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1))  
p1=1-p2                            

lowpimpl=p1
lowpcap=p2*p3
lowpdis=p2*p4

# Predicted Probabilities for MEAN VALUES
n=299
const=rep(1,n)
x11=cbind(mean_znumber_of_provisions,mean_zdelegationratio,mean_zefficient_score)
x13=const
x14=cbind(mean_zmsdis_sum,mean_zmss_sum,x,mean_typodir,mean_zinterest,mean_numissues)
x24=cbind(const,mean_zcomgain_sum,mean_zcomss_sum,x,mean_typodir,mean_zinterest,mean_dummytime)
b11=matrix(strat.mle$par[1:3],3,1)  
b13=matrix(strat.mle$par[4],1,1)   
b14=matrix(strat.mle$par[5:10],6,1)     
b24=matrix(strat.mle$par[11:17],7,1)     
u11=x11 %*% b11
u13=x13 %*% b13                    
u14=x14 %*% b14
u24=x24 %*% b24
#Probabilities
p4=pnorm(u24/sqrt(2))             
p3=1-p4                           
p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1)) 
p1=1-p2                           

meanpimpl=p1
meanpcap=p2*p3
meanpdis=p2*p4

# Predicted Probabilities for MODERATE VALUES
n=299
const=rep(1,n)
x11=cbind(mode_znumber_of_provisions,mode_zdelegationratio,mode_zefficient_score)
x13=const
x14=cbind(mode_zmsdis_sum,mode_zmss_sum,x,mode_typodir,mode_zinterest,mode_numissues)
x24=cbind(const,mode_zcomgain_sum,mode_zcomss_sum,x,mode_typodir,mode_zinterest,mode_dummytime)
b11=matrix(strat.mle$par[1:3],3,1)  
b13=matrix(strat.mle$par[4],1,1)   
b14=matrix(strat.mle$par[5:10],6,1)      
b24=matrix(strat.mle$par[11:17],7,1)    
u11=x11 %*% b11
u13=x13 %*% b13                     
u14=x14 %*% b14
u24=x24 %*% b24
#Probabilities
p4=pnorm(u24/sqrt(2))              
p3=1-p4                            
p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1))  
p1=1-p2                           

modepimpl=p1
modepcap=p2*p3
modepdis=p2*p4


#Figure 4a
#Plot Pr(enforced compliance|zcouconf_sum) holding other indvar at min/mean/max values
plot(x,meanpdis,ylab="Pr(Enforced Compliance)",xlab="Diversity of states' interests",ylim=c(0,1.0), type="l", lty="solid",lwd="2.0",cex.lab=1.5)
lines(x,lowpdis,lty="F8",lwd="2.0")
lines(x,modepdis,lty="dotted",lwd="2.0")
legend(-2.0, 1.0, c("Mean", "Low", "Moderate"), cex=1, lty=c("solid","F8","dotted"), lwd=1.9)

#Figure 4b
#Plot Pr(compliance deficit|zcouconf_sum) holding other indvar at min/mean/max values
plot(x,meanpcap,ylab="Pr(Deficit)",xlab="Diversity of states' interests",ylim=c(0,1.0), type="l", lty="solid",lwd="2.0",cex.lab=1.5)
lines(x,lowpcap,lty="F8",lwd="2.0")
lines(x,modepcap,lty="dotted",lwd="2.0")
legend(-2.0, 1.0, c("Mean", "Low", "Moderate"), cex=1.5, lty=c("solid","F8","dotted"), lwd=1.0)



#X-VARIABLE: MEMBER STATES' DISAGREEMENT
#zmsdis_sum
n=299
const=rep(1,n)
x=seq(-1.2,4.2,length=n)

# Pr(impl/dis/cap|zcouconf_sum) holding other indvar at min/mean/max values
#u11:znumber_of_provisions,zdelegationratio,zefficient_score
#u13:const
#u14:zmsdis_sum,zmss_sum,zcouconf_sum,typodir,zinterest,numissues
#u23:0
#u24:const,zcomdis_sum,zcomss_sum,zcouconf_sum,typodir,zinterest,dummytime

# Predicted Probabilities for LOW VALUES
n=299
const=rep(1,n)
x11=cbind(low_znumber_of_provisions,low_zdelegationratio,low_zefficient_score)
x13=const
x14=cbind(x,low_zmss_sum,low_zcouconf_sum,low_typodir,low_zinterest,low_numissues)
x24=cbind(const,low_zcomgain_sum,low_zcomss_sum,low_zcouconf_sum,low_typodir,low_zinterest,low_dummytime)
b11=matrix(strat.mle$par[1:3],3,1)  
b13=matrix(strat.mle$par[4],1,1)    
b14=matrix(strat.mle$par[5:10],6,1)    
b24=matrix(strat.mle$par[11:17],7,1)     
u11=x11 %*% b11
u13=x13 %*% b13                     
u14=x14 %*% b14
u24=x24 %*% b24
#Probabilities
p4=pnorm(u24/sqrt(2))              
p3=1-p4                            
p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1))  
p1=1-p2                            

lowpimpl=p1
lowpcap=p2*p3
lowpdis=p2*p4

# Predicted Probabilities for MEAN VALUES
n=299
const=rep(1,n)
x11=cbind(mean_znumber_of_provisions,mean_zdelegationratio,mean_zefficient_score)
x13=const
x14=cbind(x,mean_zmss_sum,mean_zcouconf_sum,mean_typodir,mean_zinterest,mean_numissues)
x24=cbind(const,mean_zcomgain_sum,mean_zcomss_sum,mean_zcouconf_sum,mean_typodir,mean_zinterest,mean_dummytime)
b11=matrix(strat.mle$par[1:3],3,1)  
b13=matrix(strat.mle$par[4],1,1)    
b14=matrix(strat.mle$par[5:10],6,1)      
b24=matrix(strat.mle$par[11:17],7,1)     
u11=x11 %*% b11
u13=x13 %*% b13                    
u14=x14 %*% b14
u24=x24 %*% b24
#Probabilities
p4=pnorm(u24/sqrt(2))              
p3=1-p4                            
p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1))  
p1=1-p2                            

meanpimpl=p1
meanpcap=p2*p3
meanpdis=p2*p4

# Predicted Probabilities for MODERATE VALUES
n=299
const=rep(1,n)
x11=cbind(mode_znumber_of_provisions,mode_zdelegationratio,mode_zefficient_score)
x13=const
x14=cbind(x,mode_zmss_sum,mode_zcouconf_sum,mode_typodir,mode_zinterest,mode_numissues)
x24=cbind(const,mode_zcomgain_sum,mode_zcomss_sum,mode_zcouconf_sum,mode_typodir,mode_zinterest,mode_dummytime)
b11=matrix(strat.mle$par[1:3],3,1)  
b13=matrix(strat.mle$par[4],1,1)    
b14=matrix(strat.mle$par[5:10],6,1)      
b24=matrix(strat.mle$par[11:17],7,1)     
u11=x11 %*% b11
u13=x13 %*% b13                     
u14=x14 %*% b14
u24=x24 %*% b24
#Probabilities
p4=pnorm(u24/sqrt(2))              
p3=1-p4                            
p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1))  
p1=1-p2                           

modepimpl=p1
modepcap=p2*p3
modepdis=p2*p4

#Figure 4c
#Plot Pr(dis|zmsdis_sum) holding other indvar at min/mean/max values
plot(x,meanpdis,ylab="Pr(Enforced Compliance)",xlab="Member state's disagreement",ylim=c(0,1.0), type="l", lty="solid",lwd="2.0",cex.lab=1.5)
lines(x,lowpdis,lty="F8",lwd="2.0")
lines(x,modepdis,lty="dotted",lwd="2.0")
legend(-1.2, 1.0, c("Mean", "Low", "Moderate"), cex=1.2, lty=c("solid","F8","dotted"), lwd=1.0)

#Figure 4d
#Plot Pr(cap|zmsdis_sum) holding other indvar at min/mean/max values
plot(x,meanpcap,ylab="Pr(Deficit)",xlab="Member state's disagreement",ylim=c(0,1.0), type="l", lty="solid",lwd="2.0",cex.lab=1.5)
lines(x,lowpcap,lty="F8",lwd="2.0")
lines(x,modepcap,lty="dotted",lwd="2.0")
legend(-1.2, 1.0, c("Mean", "Low", "Moderate"), cex=1.2, lty=c("solid","F8","dotted"), lwd=1.0)


#Figure 5
#3D Plots outcome probabilies for a change in two explanatory variables
#Outcome function
#x: Diversity of states
range(zcouconf_sum)
x <- seq(-1.347835,2.570057,length=299)
#y: MS disagreement
range(zmsdis_sum)
y <- seq(-0.9607449,4.0349951,length=299)

outcomedis <- function(x,y){
  n=299
  const=rep(1,n)
  x11=cbind(mean_znumber_of_provisions,mean_zdelegationratio,mean_zefficient_score)
  x13=const
  x14=cbind(y,mean_zmss_sum,x,mean_typodir,mean_zinterest,mean_numissues)
  x24=cbind(const,mean_zcomgain_sum,mean_zcomss_sum,x,mean_typodir,mean_zinterest,mean_dummytime)
  b11=matrix(strat.mle$par[1:3],3,1)  
  b13=matrix(strat.mle$par[4],1,1)    
  b14=matrix(strat.mle$par[5:10],6,1)     
  b24=matrix(strat.mle$par[11:17],7,1)     
  u11=x11 %*% b11
  u11 <- as.numeric(u11)
  u13=x13 %*% b13                    
  u13 <- as.numeric(u13)
  u14=x14 %*% b14
  u14 <- as.numeric(u14)
  u24=x24 %*% b24
  u24 <- as.numeric(u24)
  #Probabilities
  p4=pnorm(u24/sqrt(2))
  p3=1-p4
  p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1))  
  p1=1-p2                            
  
  pimpl=p1
  pcap=p2*p3
  pdis=p2*p4
  
  return(pdis)
}

zdis <- outer(x,y,outcomedis)

persp(x, y, zdis, theta = 30, phi = 30, expand = 0.9, col = "gray80", border= NA, shade = 0.8, 
      xlab="Diversity of states' interests",ylab="MS disagreement",zlab="Pr(Enforced Compliance)",
      ticktype="detailed",cex.lab=1,zlim=c(0,1))