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
romanian <- read.dta(file="./APP1/romanian_twoschooltown.dta")
romanian<- romanian[c("bct_Dmean","dzag","meanzag_jk","attend")]
romanian<-subset(romanian,abs(dzag)<=min(abs(min(dzag)),max(dzag))) ## avoid the possibility of having unbalanced estimation window
n<-length(romanian$dzag)
bwsel="CCT"
mrdbw1.out = rdbwselect(romanian$attend,romanian$dzag, bwselect = bwsel)
mrdbw2.out = rdbwselect(romanian$bct_Dmean,romanian$dzag, bwselect = bwsel)

## following are the main codes
#################################################################
nbs=1000
ngrid=2
Qgrid=c(50,75)
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)
para1=matrix(NA,nrow=nsm,ncol=3)
para2=matrix(NA,nrow=nsm,ncol=3)
para3=matrix(NA,nrow=nsm,ncol=3)
for (gridind in 1:ngrid){
	Q=Qgrid[gridind]
	print(Q)
	for (smind in 1:nsm){
        smooth=4+0.25*smind

## first stage
bw=mrdbw1.out$bws[1,1]*n^(1/5-1/smooth)
data<-subset(romanian,abs(romanian$dzag)<=bw)
## test
set.seed(123)
y=data$attend
z=data$dzag
x=data$meanzag_jk
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)
wght=triawt(z,bw)
pos=(z>0)
xpos=x*pos
zpos=z*pos
results<-summary(lm(y~x+z+pos+xpos+zpos,weights=wght))$coefficients
para1[smind,]<-results[5,1:3]

## reduced-form
bw=mrdbw2.out$bws[1,1]*n^(1/5-1/smooth)
data<-subset(romanian,abs(romanian$dzag)<=bw)
# test
set.seed(123)
y=data$bct_Dmean
z=data$dzag
x=data$meanzag_jk
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)
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)
wght=triawt(z,bw)
pos=(z>0)
xpos=x*pos
zpos=z*pos
results<-summary(lm(y~x+z+pos+xpos+zpos,weights=wght))$coefficients
para2[smind,]<-results[5,1:3]

## treatment effect
set.seed(123)
t=data$attend
results<-RDMonoTest(y,z,x,bw,Q=Q,nbs=nbs,t=t)
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)
wght=triawt(z,bw)
pos=(z>0)
xpos=x*pos
zpos=z*pos
xt=x*t
results<-summary(ivreg(y~x+z+t+xt+zpos|x+z+pos+xpos+zpos,weights=wght))$coefficients
para3[smind,]<-results[5,1:3]
}
print(pvalue1[,gridind])
print(pvalue2[,gridind])
print(pvalue3[,gridind])
print(para1)
print(para2)
print(para3)
}

save.image("./APP1/romanian.RData")

require(stargazer) 
pval<-rbind(t(pvalue1[1:6,]),t(pvalue1[7:12,]), t(pvalue2[1:6,]),t(pvalue2[7:12,]),t(pvalue3[1:6,]),t(pvalue3[7:12,]))
stargazer(pval,digits=3,out = paste("APP1table_np.tex",sep="")) 	

para<-rbind(para1,para2,para3)
stargazer(para,digits=3,out = paste("APP1table_p.tex",sep="")) 	
