##Treatment Group Matched (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)
}

dataT<-read.table("C:/Documents and Settings/Kevin Arceneaux/Desktop/trtdata.txt",header=T)
dataT$ward<-as.factor(dataT$ward)
pscore.glm<-glm(phntrt.one ~ persons + vote96.1 + newvoter + majorpty + age +
                ward + agesq, 
                family=binomial(logit), data=dataT)
D<-dataT$phntrt.one
Y<-dataT$voted98
grp<-dataT$persons
psboot(Y, D, fitted(pscore.glm), sims=500, n.match=1, grp=grp) 
psboot(Y, D, fitted(pscore.glm), sims=500, n.match=2, grp=grp) 
psboot(Y, D, fitted(pscore.glm), sims=500, n.match=3, grp=grp) 
psboot(Y, D, fitted(pscore.glm), sims=500, n.match=4, grp=grp) 
psboot(Y, D, fitted(pscore.glm), sims=500, n.match=5, grp=grp) 


#NOT bootstrap (without replacement)

psboot <- function(y, treat, pscores, sims=1, n.match=3, grp=NULL){ 

  uG<-unique(grp)
  nG<-length(uG)
  ate <- matrix(0,nrow=sims,ncol=nG+1)
  data<-cbind(y,treat,pscores)
  labels <- seq(1:nrow(data))
   for(i in 1:sims){
    bootrows <- labels
    bootsample <- data[bootrows,]
    y.sample <- bootsample[,1]
    treat.sample <- bootsample[,2]
    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[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)
}


dataT<-read.table("C:/Documents and Settings/Kevin Arceneaux/Desktop/trtdata.txt",header=T)
dataT$ward<-as.factor(dataT$ward)
pscore.glm<-glm(phntrt.one ~ persons + vote96.1 + newvoter + majorpty + age +
                ward + agesq, 
                family=binomial(logit), data=dataT)
D<-dataT$phntrt.one
Y<-dataT$voted98
grp<-dataT$persons
psboot(Y, D, fitted(pscore.glm), sims=100, n.match=1, grp=grp) 
psboot(Y, D, fitted(pscore.glm), sims=100, n.match=2, grp=grp) 
psboot(Y, D, fitted(pscore.glm), sims=100, n.match=3, grp=grp) 
psboot(Y, D, fitted(pscore.glm), sims=100, n.match=4, grp=grp) 
psboot(Y, D, fitted(pscore.glm), sims=100, n.match=5, grp=grp)