#######################################
## Replication R-code:               ##
## Parties in the Council?           ##
## Hagemann & Hoyland (2008)         ##
## Journal of European Public Policy ##
## Vol 15, Issue 8, 1205 -- 1221     ##
####################################### 
###########
rm(list=ls())
data <- "http://folk.uio.no/bjornkho/JEPPData.txt"
Votes<-read.table(url(data))
party <- "http://folk.uio.no/bjornkho/parties.txt"
Party <-data.frame(t(read.table(url(party))))
library(MCMCpack)
allresults<- MCMCirtKdRob(t(Votes), dimensions=1,
burnin=1e4, mcmc=1e5, thin=10,
B0=.25, store.item=FALSE, store.ability=TRUE,
verbose=1e3)

t2 <- allresults[,3:43]
tra<- data.frame(NA)
for( i in 1:dim(allresults)[1]){
tr<-as.data.frame(rank(t2[i,]))
tra <- cbind(tra,tr)
}
tra <- tra[,-1]
summary(t(tra))
ttra <-t(tra)
rankings <- data.frame(rep(0,41),rep(0,41),rep(0,41),rep(0,41),
rep(0,41))
row.names(rankings) <-colnames(Votes)
for (i in 1:dim(ttra)[2]){
rankings[i,] <- quantile(ttra[,i],probs=c(.5,.1,.9,.25,.75))
}

colnames(rankings) <-c("med","lower","higher","low25","high75")
rankings$party <- Party
New<-data.frame(rankings[order(rankings$med),])
names(New) <- c("med","lower","higher","low25","high75","party")
PES <- New[New$party=="pes ",]
EPP <- New[New$party=="epp ",]
index <- seq(dim(New)[1]:1)

pdf("Rankings.pdf")
par(mfrow=c(1,1))
plot(New$med,index,type="p",xlim=c(0,60),bty="n",ylab="",
main="1999 - Oct 2007",col=1,pch=16,yaxt="n",
xlab="Rankings in the Council")
for (j in 1:dim(New)[1]){
segments(x0=New$low25[j],y0=index[j],x1=New$high75[j],y1=index[j],lwd=1)
#segments(x0=New$lower[j],y0=index[j],x1=New$higher[j],y1=index[j],lwd=1,col="grey")
}
abline(v=median(EPP[,1]),col="gray",lwd=2,lty=2)
abline(v=median(PES[,1]),col="gray",lwd=2,lty=3)
text(46,y=index,labels=row.names(New),adj=0,cex=.75)
text(43,10,labels=("Median PES"),adj=1,cex=.75)
text(43,24,labels=("Median EPP"),adj=1,cex=.75)
#legend(18,4,legend=c("95% credibility of rank order"),cex=.75,
#lwd=1,bty="n",adj=0,col="grey")
legend(18,3,legend=c("Central Tendency of rank order"),cex=.75,
lwd=1,bty="n",adj=0)
legend(20.5,2,legend=c("Median ranking"),cex=.75,
pch=16,bty="n",adj=0)
dev.off()
#######################
notvoting <- NA
for (i in 1:dim(Votes)[2]){
  nonvoting <- summary(Votes[,i])[7]
  notvoting <- c(notvoting,nonvoting)
}
notvoting <- notvoting[-1]
notvoting[is.na(notvoting)] <- 0
nrvotes <- dim(Votes)[1] - notvoting
diff <- New$high75 - New$low25
cor(nrvotes,diff)
#################
# The code below just calculates the difference
# between selected governments
difference <- allresults[,3]- allresults[,4]
count <- 0
for(i in 1:length(difference)){
   if (difference[i]<= 0){ 
        count <<- count + 1    }
}
count/length(difference)