rm(list=ls())
library(foreign)
library(locfit)
library(plotrix)

#setwd("/Users/innoshuadmin/Dropbox/RDHeteroTest/RDMonoTest")
setwd("/home/users/innoshu.AD3/RDMonoTest")
## load in robust bandwidth selection functions from CCT(2014)## 
source("./functions/rdbwselect.R")
source("./functions/functions.R")
## load in RD test functions
source("./functions/RDMono_Functions_Generic.R")

## load in data
canadian <- read.dta(file="./APP2/canadiancollege.dta")
canadian<- canadian[c("dist_from_cut","probation_year1","left_school","hsgrade_pct","female","suspended_ever","english")]
#canadian<-subset(canadian,probation_year1==0| dist_from_cut<=0) ## 3 students with GPA above cutoff are on probation 
#canadian<-subset(canadian,probation_year1==1| dist_from_cut>=0) ## 48 students with GPA below cutoff are not on 

nbs=1000
Zeval <- seq(-0.6, 0.6, length=31)  

pdf("./APP2/canadianMRDplot.pdf",width=10,height=4)
par(mfcol=c(1,3))
## first stage
data<-canadian
y=data$probation_year1
z=data$dist_from_cut
xlab="First Year GPA From Cutoff"
ylab1="First Stage"  
ylab2="# of Individuals"
main="Full Sample"
labels<-data.frame(xlab,ylab1,ylab2,main)
results<-MRD_sharp(y,z,bw=0.6,Zeval)
MRDplot(Zeval,y,z,results,labels)

for (HSindex in 1:2){
	## load in data
	if (HSindex==1){
		main="HS grades < median"
		data<-subset(canadian,hsgrade_pct<50)
		} 
	else {
		main="HS grades > median"
		data<-subset(canadian,hsgrade_pct>=50)
		}
## treatment effect
## leave school after first evaluation 
y=data$left_school
z=data$dist_from_cut
xlab="First Year GPA From Cutoff"
ylab1="Leaving School After First Year"  
ylab2="# of Individuals"
labels<-data.frame(xlab,ylab1,ylab2,main)
results<-MRD_sharp(y,z,bw=0.6,Zeval)
MRDplot(Zeval,y,z,results,labels)
}
dev.off()

