##Phone:  Noncompliance (Data constructed from the "apsr_corrected_latest.dat"
##file located on Imai's website and associated with his replication file)
##We use code written by Imai to conduct matching.  Computer errors in the program were corrected by us (see
##corrected_code.R).

#Bootstrapped (without replacement)

psboot <- function(y, treat, pscores, sims=500, n.match=3, grp=NULL){ 
  uG<-unique(grp)
  nG<-length(uG)
  ate <- matrix(0,nrow=sims,ncol=nG+1)
  data<-cbind(y,treat,pscores,grp)
  labels <- seq(1:nrow(data))
  for(i in 1:sims){
    bootrows <- sample(labels,nrow(data),replace=TRUE)
    bootsample <- data[bootrows,]
    y.sample <- bootsample[,1]
    treat.sample <- bootsample[,2]
    grp.sample <- bootsample[,4]
    match <- psmatch(treat.sample, bootsample[,3],
                     n.match=n.match)
    diff <- y.sample[treat.sample==1]-y.sample[c(match$match.treat)]
    ate[i,1] <- mean(diff)
    for(j in 1:nG)
      ate[i,j+1] <- mean(diff[grp.sample[treat.sample==1]==uG[j]])
    if(i==sims){
      cat(i,"iter.(ate=",mean(ate[1:i,1]),"sd=",sd(ate[1:i,1]),")...\n")
      for(j in 1:nG)
        cat("**",uG[j],"grp.(ate=",mean(ate[1:i,j+1]),"sd=",sd(ate[1:i,j+1]),")...\n")
    }
  }
}

psmatch <- function(treat, pscores, n.match=3){
  n <- length(treat)
  n1 <- length(treat[treat==1])
  labels <- seq(1,n)
  clabels <- labels[treat==0]
  pscore1 <- pscores[treat==1]
  pscore0 <- pscores[treat==0]
  matchedT <- matrix(0, nrow=length(pscore1), ncol=n.match)
  p0<-pscore0
  for(i in 1:n1){
    needed<-n.match
    matched<-NULL
    for(j in 1:n.match)
      if(needed>0){
        deviation <- abs(p0-pscore1[i])
        mindev <- clabels[min(deviation)==deviation]
        if(length(mindev)>needed)
          temp <- sample(mindev,needed,replace=F)
        else 
          temp <- mindev
        for(k in 1:n.match)
          p0[clabels==temp[k]] <- 1000
        matched<-c(matched, temp)
        needed<-n.match-length(matched)
      }
    matchedT[i,] <- matched
  }
  list(match.treat = matchedT, match.pscore=pscores)
}

data2<-read.table("C:/Documents and Settings/Kevin Arceneaux/Desktop/noncomp.txt",header=T)
data2$ward<-as.factor(data2$ward)
data2$new<-(data2$vote96.0==0 & data2$vote96.1==0)*1
data2$phone<-((data2$phonegrp+data2$mailcall)>0)*1
data2$age2<-data2$age^2/100
D<-(data2$appeal>1 & data2$noncomp==1 & data2$mailcall==0 &
    data2$persngrp==0 & data2$mailings==0)*1 
C<-(data2$phonegrp==0 & data2$mailcall==0 & data2$persngrp==0 &
    data2$mailings==0)*1
dataP<-subset(data2, D+C>0)
pscore.glm<-glm(noncomp ~ persons + vote96.1 + new + majorpty + age +
                ward + persons:vote96.1 + persons:new + persons:majorpty + persons:age + 
		persons:ward + vote96.1:majorpty + vote96.1:age + vote96.1:ward + new:majorpty +
		new:age + new:ward + majorpty:age + majorpty:ward + age:ward, 
                family=binomial(logit), data=dataP)
D<-dataP$noncomp
Y<-dataP$voted98
grp<-dataP$persons
psboot(Y, D, fitted(pscore.glm), sims=500, n.match=5, grp=grp)