gistfile1.r

##########################
##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)