rm(list=ls())
#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")
source("./simulationR2/simu_functions.R") 

load("./simulationR2/DGP_MonoTest.RData")
numn=3
nummodel=4
model=1
numsm=3
numsim=1000
nbs=1000
Q=10
bwsel="CCT"

source("./simulationR2/SRD_DGPs.R")
modelpara=data.frame(crd.p,crd.m,cfs.p,cfs.m)

avgbw=matrix(NA,nrow=numn*nummodel,ncol=numsm)
rejincprop_ks=matrix(NA,nrow=numn*nummodel,ncol=numsm)
rejdecprop_ks=matrix(NA,nrow=numn*nummodel,ncol=numsm)
rejincprop_ks_gms=matrix(NA,nrow=numn*nummodel,ncol=numsm)
rejdecprop_ks_gms=matrix(NA,nrow=numn*nummodel,ncol=numsm)
rejincprop_int=matrix(NA,nrow=numn*nummodel,ncol=numsm)
rejdecprop_int=matrix(NA,nrow=numn*nummodel,ncol=numsm)
ptm <- proc.time()
for (nind in 1:numn){
	n=2000*2^(nind-1)
	set.seed(123)
	bw.sim=matrix(NA,nrow=numsim,ncol=numsm)
	rejinc_ks.sim=matrix(NA,nrow=numsim,ncol=numsm)
	rejdec_ks.sim=matrix(NA,nrow=numsim,ncol=numsm)
	rejinc_ks_gms.sim=matrix(NA,nrow=numsim,ncol=numsm)
	rejdec_ks_gms.sim=matrix(NA,nrow=numsim,ncol=numsm)
	rejinc_int.sim=matrix(NA,nrow=numsim,ncol=numsm)
	rejdec_int.sim=matrix(NA,nrow=numsim,ncol=numsm)
	for (s in 1:numsim){
		results=simu_RDMonoTest(modelpara,n,numsm,bwsel=bwsel,Q=Q,nbs=nbs,type=type,para=1)
		bw.sim[s,]=results$bw
		rejinc_ks.sim[s,]=results$rejinc_ks
		rejdec_ks.sim[s,]=results$rejdec_ks
		rejinc_ks_gms.sim[s,]=results$rejinc_ks_gms
		rejdec_ks_gms.sim[s,]=results$rejdec_ks_gms
		rejinc_int.sim[s,]=results$rejinc_int
		rejdec_int.sim[s,]=results$rejdec_int
		}
	avgbw[(nind-1)*nummodel+model,]=colMeans(bw.sim)
	rejincprop_ks[(nind-1)*nummodel+model,]=colMeans(rejinc_ks.sim)
	rejdecprop_ks[(nind-1)*nummodel+model,]=colMeans(rejdec_ks.sim)
	rejincprop_ks_gms[(nind-1)*nummodel+model,]=colMeans(rejinc_ks_gms.sim)
	rejdecprop_ks_gms[(nind-1)*nummodel+model,]=colMeans(rejdec_ks_gms.sim)
	rejincprop_int[(nind-1)*nummodel+model,]=colMeans(rejinc_int.sim)
	rejdecprop_int[(nind-1)*nummodel+model,]=colMeans(rejdec_int.sim)	
	print(c(n,model,proc.time() - ptm))
	print(avgbw[1:(nind*nummodel),])
	print(rejincprop_ks[1:(nind*nummodel),])
	print(rejdecprop_ks[1:(nind*nummodel),])
	print(rejincprop_ks_gms[1:(nind*nummodel),])
	print(rejdecprop_ks_gms[1:(nind*nummodel),])
	print(rejincprop_int[1:(nind*nummodel),])
	print(rejdecprop_int[1:(nind*nummodel),])
	save.image(file=paste("./simulationR2/mktable/SRDpanelA_model",model,".RData",sep=""))
}

proc.time() - ptm