simu_RDMonoTest<- function(modelpara,n,numsm=3,bwsel="CCT",Q=5,nbs=200,eta=10^(-6),para=0,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)
  u = rnorm(n)
  y.l =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 =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[Z<0],y.r[Z>0])+u*0.1
  if  (type=="Fuzzy"){
  t.l =(cfs.m[1] + cfs.m[2]*X + cfs.m[3]*Z + cfs.m[4]*Z*X + cfs.m[5]*Z^2 + cfs.m[6]*X^2+u>0)
  t.r =(cfs.p[1] + cfs.p[2]*X + cfs.p[3]*Z + cfs.p[4]*Z*X + cfs.p[5]*Z^2 + cfs.p[6]*X^2+u>0)  
  tdat = c(t.l[Z<0],t.r[Z>0])
  }
  zdat = c(Z[Z<0],Z[Z>0])
  xdat = c(X[Z<0],X[Z>0])

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)
rejinc_int=matrix(NA,nrow=1,ncol=numsm)
rejdec_int=matrix(NA,nrow=1,ncol=numsm)
mrdbw.out = rdbwselect(ydat,zdat, bwselect = bwsel)

for (smind in 1:numsm){
	if (bwsel=="CCT"){
			smooth=4+smind*0.25
	} else {
		smooth=4+(smind-1)*0.25
		print(smooth)
	}

	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]
	if  (type=="Fuzzy"){
		t=tdat[abs(zdat)<=bw]
  	}  else
  	{   t=NULL
  	} 
	z=zdat[abs(zdat)<=bw]
	# test for monotonicity
	if (para==1){
	## parametric test with interaction term
	wght=triawt(z,bw)
	pos=(z>0)
	xpos=x*pos
	zpos=z*pos
	zx=x*z
	zxpos=z*x*pos
	results<-summary(lm(y~x+z+pos+xpos+zpos,weights=wght))$coefficients
#	results<-summary(lm(y~x+z+pos+xpos+zpos+zx+zxpos,weights=wght))$coefficients  // to generate the simulation results reported in footnote 2.
	rejdec_int[smind]<-(results[5,1]/results[5,2]< -1.64)
	rejinc_int[smind]<-(results[5,1]/results[5,2]> 1.64)
	} else
	{	rej_int=NULL
		}
	## proposed test
	results<-RDMonoTest(y,z,x,bw,Q=Q,nbs=nbs,t=t)
	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_int=rejinc_int,rejdec_int=rejdec_int,rejinc_ks=rejinc_ks,rejdec_ks=rejdec_ks,rejinc_ks_gms=rejinc_ks_gms,rejdec_ks_gms=rejdec_ks_gms)
}




