rm(list=ls())
library(foreign)
library(AER)
#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 probation, some are suspended in the first year 
canadian<-subset(canadian,abs(dist_from_cut)<=min(abs(min(dist_from_cut)),max(dist_from_cut))) ## avoid any chance of having unbalanced estimation window
n<-length(canadian$dist_from_cut)
bwsel="CCT"
set.seed(123)
mrdbw.out = rdbwselect(canadian$left_school,canadian$dist_from_cut, bwselect = bwsel)
	
## following are the main codes
#################################################################
nbs=1000
Qgrid=c(30,40)
ngrid=length(Qgrid)
nsm=3
pvalue1=matrix(NA,nrow=4*nsm,ncol=ngrid)
pvalue2=matrix(NA,nrow=4*nsm,ncol=ngrid)
pvalue3=matrix(NA,nrow=4*nsm,ncol=ngrid)
for (gridind in 1:ngrid){
	Q=Qgrid[gridind]
	print(Q)
	for (smind in 1:nsm){
        smooth=4+0.25*smind
# leave school after first year
bw=mrdbw.out$bws[1,1]*n^(1/5-1/smooth)
data<-subset(canadian,abs(canadian$dist_from_cut)<=bw)
## test
set.seed(123)
y=data$left_school
z=data$dist_from_cut
x=data$hsgrade_pct
print(c(bw,c(length(y[z<0]),length(y[z>=0]))/Q))
results<-RDMonoTest(y,z,x,bw,Q=Q,nbs=nbs)
bscdf_inc=ecdf(results$statinc_ks_bs)
pvalue1[smind,gridind]=1-bscdf_inc(results$statinc_ks)
bscdf_inc_gms=ecdf(results$statinc_ks_bs_gms)
pvalue1[nsm+smind,gridind]=1-bscdf_inc_gms(results$statinc_ks)
bscdf_dec=ecdf(results$statdec_ks_bs)
pvalue1[2*nsm+smind,gridind]=1-bscdf_dec(results$statdec_ks)
bscdf_dec_gms=ecdf(results$statdec_ks_bs_gms)
pvalue1[3*nsm+smind,gridind]=1-bscdf_dec_gms(results$statdec_ks)
set.seed(123)
s_d=data$english
print(c(bw,c(length(y[z<0&s_d==0]),length(y[z>=0&s_d==0]),length(y[z<0&s_d==1]),length(y[z>=0&s_d==1]))/Q))
results<-RDMonoTest(y,z,x,bw,s_d=s_d,Q=Q,nbs=nbs)
bscdf_inc=ecdf(results$statinc_ks_bs)
pvalue2[smind,gridind]=1-bscdf_inc(results$statinc_ks)
bscdf_inc_gms=ecdf(results$statinc_ks_bs_gms)
pvalue2[nsm+smind,gridind]=1-bscdf_inc_gms(results$statinc_ks)
bscdf_dec=ecdf(results$statdec_ks_bs)
pvalue2[2*nsm+smind,gridind]=1-bscdf_dec(results$statdec_ks)
bscdf_dec_gms=ecdf(results$statdec_ks_bs_gms)
pvalue2[3*nsm+smind,gridind]=1-bscdf_dec_gms(results$statdec_ks)
s_d=data$female
print(c(bw,c(length(y[z<0&s_d==0]),length(y[z>=0&s_d==0]),length(y[z<0&s_d==1]),length(y[z>=0&s_d==1]))/Q))
results<-RDMonoTest(y,z,x,bw,s_d=s_d,Q=Q,nbs=nbs)
bscdf_inc=ecdf(results$statinc_ks_bs)
pvalue3[smind,gridind]=1-bscdf_inc(results$statinc_ks)
bscdf_inc_gms=ecdf(results$statinc_ks_bs_gms)
pvalue3[nsm+smind,gridind]=1-bscdf_inc_gms(results$statinc_ks)
bscdf_dec=ecdf(results$statdec_ks_bs)
pvalue3[2*nsm+smind,gridind]=1-bscdf_dec(results$statdec_ks)
bscdf_dec_gms=ecdf(results$statdec_ks_bs_gms)
pvalue3[3*nsm+smind,gridind]=1-bscdf_dec_gms(results$statdec_ks)
}
print(pvalue1[,gridind])
print(pvalue2[,gridind])
print(pvalue3[,gridind])
}

save.image("./APP2/canadian.RData")

require(stargazer) 
pval<-rbind(t(pvalue1),t(pvalue2),t(pvalue3))
stargazer(pval,digits=2,out = paste("APP2table.tex",sep="")) 	
