simu_sd<- function(modelpara,n,numsm=3,bwsel="CCT",Q=5,nbs=200,eta=10^(-6),type="Sharp"){
## DGP
crd.m=modelpara$crd.m
crd.p=modelpara$crd.p
cfs.m=modelpara$cfs.m
cfs.p=modelpara$cfs.p

  Z = 2*rbeta(n,2,2)-1 
  X = runif(n)
  S = round(runif(n))
  u = rnorm(n)
  y.l.s0 =crd.m[1] + crd.m[3]*Z + crd.m[5]*Z^2 
  y.r.s0 =crd.p[1] + crd.p[3]*Z + crd.p[5]*Z^2 
  y.l.s1 =crd.m[1] + crd.m[2]*X+ crd.m[3]*Z + crd.m[4]*Z*X + crd.m[5]*Z^2 + crd.m[6]*X^2
  y.r.s1 =crd.p[1] + crd.p[2]*X + crd.p[3]*Z + crd.p[4]*Z*X + crd.p[5]*Z^2 + crd.p[6]*X^2 
  ydat = c(y.l.s0[Z<0&S==0],y.l.s1[Z<0&S==1],y.r.s0[Z>0&S==0],y.r.s1[Z>0&S==1])+u*0.1
  zdat = c(Z[Z<0&S==0],Z[Z<0&S==1],Z[Z>0&S==0],Z[Z>0&S==1])
  xdat = c(X[Z<0&S==0],X[Z<0&S==1],X[Z>0&S==0],X[Z>0&S==1])
  sdat = c(S[Z<0&S==0],S[Z<0&S==1],S[Z>0&S==0],S[Z>0&S==1])
bwmat=matrix(NA,nrow=1,ncol=numsm)
rejinc_ks=matrix(NA,nrow=1,ncol=numsm)
rejdec_ks=matrix(NA,nrow=1,ncol=numsm)
rejinc_ks_gms=matrix(NA,nrow=1,ncol=numsm)
rejdec_ks_gms=matrix(NA,nrow=1,ncol=numsm)
mrdbw.out = rdbwselect(ydat,zdat, bwselect = bwsel)

for (smind in 1:numsm){
	smooth=4+smind*0.25
	bw=mrdbw.out$bws[1,1]*n^(1/5-1/smooth)
	## trim sample to speed up simulation
	y=ydat[abs(zdat)<=bw]
	x=xdat[abs(zdat)<=bw]
	s_d=sdat[abs(zdat)<=bw]
	if  (type=="Fuzzy"){
		t=tdat[abs(zdat)<=bw]
  	}  else
  	{   t=NULL
  	} 
	z=zdat[abs(zdat)<=bw]
	results<-RDMonoTest(y,z,x,bw,s_d=s_d,Q=Q,nbs=nbs,t=t)
	#print(results$hatnu_st)
	# test for monotonicity
	bwmat[smind]<-bw
	rejinc_ks[smind]<-(results$statinc_ks>quantile(results$statinc_ks_bs,0.95))
	rejdec_ks[smind]<-(results$statdec_ks>quantile(results$statdec_ks_bs,0.95))
	rejinc_ks_gms[smind]<-(results$statinc_ks>(quantile(results$statinc_ks_bs_gms,0.95+eta)+eta))
	rejdec_ks_gms[smind]<-(results$statdec_ks>(quantile(results$statdec_ks_bs_gms,0.95+eta)+eta))
}
list(bw=bwmat,rejinc_ks=rejinc_ks,rejdec_ks=rejdec_ks,rejinc_ks_gms=rejinc_ks_gms,rejdec_ks_gms=rejdec_ks_gms)
}


