##########################

##code for 1 unit ahead model evaluation of BTSCS data

##James E. Yonamine

##1/3/2013

##data from http://www.apsanet.org/content_29436.cfm

##run time: 12 seconds

##########################

rm(list=ls())

library(foreign)

library(ROCR)

library(lattice)

library(Zelig)

library(Amelia)

library(stats)

data <- read.dta("collier and hoeffler 2004.dta")

data<-data[complete.cases(data$warsa),] #this drops all rows with NA for warsa, i.e. drops all continuations. C&H use this approach, which is reasonable

a.out <- amelia(data, m = 1, ts = "year", cs = "country")

data<-a.out$imputations[[1]]

data<-data[order(data$year),]

data<-cbind(data, data.frame(count=cumsum(c(TRUE,data$year[-1]!=data$year[-length(data$year)]))))

num<-data[nrow(data),ncol(data)] ##figure out how many unique time series you have

start = 4 ##this determine what portion of the data to use as the initial training set

step=1 ##how many steps ahead do you want to forecast

history <- matrix(data=0, nrow=(num-start), ncol=3)

colnames(history) <- c("year", "griev_1.auc", "opp_1.auc")

for (i in start:(num-1)) {

set.seed(i)

train.data<- data[which(data[,ncol(data)] <= i),] #set initial test set

test.data <- data[which(data[,ncol(data)] ==(i+step)),] #set initial trianing set to be one period ahead

#########Grievance################

grievance.1 <- glm(warsa ~ elfo + rf + pol16 + etdo4590 + dem + peace + mount + geogia + lnpop, family=binomial(link="logit"), data = train.data)

grievance.1.predict <- predict(grievance.1, newdata=test.data, type="response")

grievance.1.y.hat<-as.matrix(grievance.1.predict)

## Establish prediction objects from ROCR package

y<-as.matrix(test.data$warsa)

with.predict <-prediction(grievance.1.y.hat,y)

## Calculate and store the AUC

with.auc <- performance(with.predict, measure = "auc")

history[(i-start)+1,2] <- as.numeric(unlist(slot(with.auc,"y.values")))

#########Opportunity###############

opportunity.1 <- glm(warsa ~sxp + sxp2 + coldwar + secm + gy1 + peace + prevwara + mount + geogia + frac + lnpop, family=binomial(link="logit"), data = train.data)

opportunity.1.predict <- predict(opportunity.1, newdata=test.data, type="response")

opportunity.1.y.hat<-as.matrix(opportunity.1.predict)

## Establish prediction objects from ROCR package

y<-as.matrix(test.data$warsa)

without.predict <-prediction(opportunity.1.y.hat,y)

## Calculate and store the AUC

without.auc <- performance(without.predict, measure = "auc")

history[(i-start)+1,3] <- as.numeric(unlist(slot(without.auc,"y.values")))

### store the temporal component

history[(i-start)+1,1] <-test.data[1,4]

}

#####plot it####

history<-as.data.frame(history)

ts(history$year)

ymax=1

ymin=min(history)

yrange=range(ymin, ymax)

plot(history[,2], type="o", col="blue", xaxt="n", ylim=yrange, xlab="Year", ylab="AUC", main="Example Plot")

lines(history[,3], type="o", col="red")

axis(1, 1:nrow(history), lab=history[1:4,1])

legend(1, c("Grievance_1","Opportunity_1"),

col=c("blue","red"), pch=21, lty=1)