## Version: July 10, 2019
triawt<-function(z,bw,ctff=0){
	kern=(1-abs((z-ctff)/bw))*(abs(z-ctff)<bw)
	S0<-sum(kern)
	S1<-sum(kern*(z-ctff))
	S2<-sum(kern*(z-ctff)^2)
	wt<-kern*(S2-S1*(z-ctff))/(S0*S2-S1^2)
	wt
}

gausswt<-function(z,bw,ctff=0){
	kern=dnorm((z-ctff)/bw)
	S0<-sum(kern)
	S1<-sum(kern*(z-ctff))
	S2<-sum(kern*(z-ctff)^2)
	wt<-kern*(S2-S1*(z-ctff))/(S0*S2-S1^2)
	wt
}

estim_func<-function(wtLT,wtRT,cubeyLT,cubeyRT,z,ctff=0){
	locfitLT=sum(wtLT*cubeyLT)
	locfitRT=sum(wtRT*cubeyRT)
	hatrho=locfitRT-locfitLT
	resrho=rep(NA,length(z))
	resrho[z<ctff]=-(cubeyLT-locfitLT)*wtLT
	resrho[z>=ctff]=(cubeyRT-locfitRT)*wtRT
	list(hatrho=hatrho,resrho=resrho)
}

estim_SRDFSfunc<-function(wt,cubeind){
	hatp=sum(wt*cubeind)
	resp=(cubeind-hatp)*wt
	list(hatp=hatp,resp=resp)
} 


RDMono<-function(y,z,x,bw,t=NULL,s_d=NULL,s_c=NULL,kernel="tria",ctff=0,Q=5){
	Q=Q-1 ## there is a discrepency in the Q definition in the codes below and the that in the paper
	nobs=length(y)
	if (kernel=="tria"){
		wt=triawt(z,bw)
		wtLT=triawt(z[z<ctff],bw)
		wtRT=triawt(z[z>=ctff],bw)
	} else 
	if (kernel=="gaussian"){
		wt=gausswt(z,bw)
		wtLT=gausswt(z[z<ctff],bw)
		wtRT=gausswt(z[z>=ctff],bw)
	} else {
		print("only triangular or gaussian kernels are allowed")
	}

	## split data based on RD running variable and cutoff
	zLT=z[z<ctff]
	zRT=z[z>=ctff]
	yLT=y[z<ctff]
	yRT=y[z>=ctff]
	if (is.null(t)){
		type="Sharp"
	} else
	{	type="Fuzzy"
		tLT=t[z<ctff]
		tRT=t[z>=ctff]
	}

## check s_d and s_c which are optional
	if (is.null(s_d)){
		s_d=rep(1,nobs)
	}
	sdLT=s_d[z<ctff]
  	sdRT=s_d[z>=ctff]
  	sdunqval=unique(s_d, fromLast = TRUE)	
  	if (length(sdunqval)==1){
		nsdcube=1
	} else
	{	nsdcube=length(sdunqval)+1
	}
	if (is.null(s_c)) {
		nell=(Q*(Q+1)*(2*Q+1)/6+Q*(Q+1)/2)/2*nsdcube
	} else
	{	nell=(Q*(Q+1)*(2*Q+1)/6+Q^2*(Q+1)^2/4)/2*nsdcube
	} 

	infl=matrix(NA,nrow=nell,ncol=nobs)
	hatnu=matrix(NA,ncol=nell,nrow=1)
	temp=0
	## loop over different values of s_d
	for (sdindex in 1:nsdcube){
			if (sdindex==nsdcube){
				sdLTind=rep(1,length(sdLT))
				sdRTind=rep(1,length(sdRT))
				sdind=rep(1,nobs)
			} else
			{
			sdLTind=(sdLT==sdunqval[sdindex])
			sdRTind=(sdRT==sdunqval[sdindex])	
			sdind=(s_d==sdunqval[sdindex])
			}
		for (q in 1:Q){
			# divide x to cubes
			xgridtemp=seq(min(x),max(x),length.out=q+2)
			catx=cut(x,xgridtemp,labels=FALSE)
			catx[is.na(catx)==TRUE]<- 1
			catxLT=catx[z<ctff]
			catxRT=catx[z>=ctff]
			# divide continuous s_c to cubes
			if (is.null(s_c)) {	
				catsc=rep(1,nobs)
				qsc=1
			} else
			{	qsc=q
				scgridtemp=seq(min(s_c),max(s_c),length.out=qsc+1) 
				catsc=cut(s_c,scgridtemp,labels=FALSE)
				catsc[is.na(catsc)==TRUE]<- 1
			}
			catscLT=catsc[z<ctff]
			catscRT=catsc[z>=ctff]
			# loop over different cubes of s_c
			for (scindex in 1:qsc){
				hatrho1=rep(NA,q)
				hatp1=rep(NA,q)
				resrho1=matrix(NA,nrow=q,ncol=nobs)
				resp1=matrix(NA,nrow=q,ncol=nobs)
				scLTind=(catscLT==scindex)
				scRTind=(catscRT==scindex)	
				scind=(catsc==scindex)
				## estimate elements of the influence functions
				# loop over different cube choices of x1
				for (x1 in 2:(q+1)){
					cubeyLT=yLT*(catxLT==x1)*scLTind*sdLTind
					cubeyRT=yRT*(catxRT==x1)*scRTind*sdRTind
					results=estim_func(wtLT,wtRT,cubeyLT,cubeyRT,z,ctff)
					hatrho1[x1-1]=results$hatrho
					resrho1[x1-1,]=results$resrho
					if (is.null(t)){
						cubeind=(catx==x1)*scind*sdind
						results=estim_SRDFSfunc(wt,cubeind)
						hatp1[x1-1]=results$hatp
						resp1[x1-1,]=results$resp
					} else
					{	cubetLT=tLT*(catxLT==x1)*scLTind*sdLTind
						cubetRT=tRT*(catxRT==x1)*scRTind*sdRTind
						results=estim_func(wtLT,wtRT,cubetLT,cubetRT,z,ctff)
						hatp1[x1-1]=results$hatrho
						resp1[x1-1,]=results$resrho
					}
				}
				# loop over different cube choices of x2
				hatrho2=rep(NA,q)
				hatp2=rep(NA,q)
				resrho2=matrix(NA,nrow=q,ncol=nobs)
				resp2=matrix(NA,nrow=q,ncol=nobs)
				for (x2 in 1:q){
					cubeyLT=yLT*(catxLT==x2)*scLTind*sdLTind
					cubeyRT=yRT*(catxRT==x2)*scRTind*sdRTind
					results=estim_func(wtLT,wtRT,cubeyLT,cubeyRT,z,ctff)
					hatrho2[x2]=results$hatrho
					resrho2[x2,]=results$resrho
					if (is.null(t)){
						cubeind=(catx==x2)*scind*sdind
						results=estim_SRDFSfunc(wt,cubeind)
						hatp2[x2]=results$hatp
						resp2[x2,]=results$resp
					} else
					{	cubetLT=tLT*(catxLT==x2)*scLTind*sdLTind
						cubetRT=tRT*(catxRT==x2)*scRTind*sdRTind
						results=estim_func(wtLT,wtRT,cubetLT,cubetRT,z,ctff)
						hatp2[x2]=results$hatrho
						resp2[x2,]=results$resrho
					}
				}
				## obtain the influence functions
				for (x1 in 2:(q+1)){
					for (x2 in 1:(x1-1)){	
						temp=temp+1
						infl[temp,]=hatp1[x1-1]*resrho2[x2,]+hatrho2[x2]*resp1[x1-1,]-hatp2[x2]*resrho1[x1-1,]-hatrho1[x1-1]*resp2[x2,]
						hatnu[temp]=hatrho2[x2]*hatp1[x1-1]-hatrho1[x1-1]*hatp2[x2]
					}
				}
			}
		} 
	}
list(hatnu=hatnu,infl=infl)
}


RDMono_multiplier<-function(infl,nbs=1000){
	nell=nrow(infl)
	nobs=ncol(infl)
	hatnu_bs=matrix(NA,ncol=nell,nrow=nbs)
	for (bs in 1:nbs){
		inflbs=infl*((matrix(1,nrow=nell,ncol=1)%*%matrix(rnorm(nobs),nrow=1,ncol=nobs)))
		hatnu_bs[bs,]=rowSums(inflbs)
	}
	hatnu_bs
}

RDMonoTest<-function(y,z,x,bw,t=NULL,s_d=NULL,s_c=NULL,kernel="tria",ctff=0,Q=5,nbs=1000,eps=0.005){
	## calculate test statistics
	results<-RDMono(y,z,x,bw,t=t,s_d=s_d,s_c=s_c,kernel=kernel,ctff=ctff,Q=Q)
	infl<-results$infl
	hatnu<-results$hatnu
	nell=nrow(results$infl)
	nobs=ncol(results$infl)
	sigma_nu<-pmax(apply(infl^2, 1, sum),rep(eps,nell)*sum(infl[1,]^2))^0.5
	hatnu_stu<-hatnu/sigma_nu
	statinc_ks<-max(hatnu_stu)
	statdec_ks<-max(-hatnu_stu)

	### define nuisance parameters
	an=(0.3*log(nobs*bw))^{1/2}
	Bn=(0.4*log(nobs*bw)/log(log(nobs*bw)))^{1/2}
	psi_nu_inc=-Bn*(hatnu_stu < -an)
	psi_nu_dec=-Bn*(-hatnu_stu < -an)

	# conduct multiplier bootstrap
	hatnu_bs<-RDMono_multiplier(infl=infl,nbs=nbs)
	statinc_ks_bs<-apply(hatnu_bs/(matrix(1,nrow=nbs,ncol=1)%*%sigma_nu), 1, max) # null: monotonically increasing (or non-deacreasing)
	statdec_ks_bs<-apply(-hatnu_bs/(matrix(1,nrow=nbs,ncol=1)%*%sigma_nu), 1, max) # null: monotonically decreasing (or non-inacreasing)
	statinc_ks_bs_gms<-apply(hatnu_bs/(matrix(1,nrow=nbs,ncol=1)%*%sigma_nu) +matrix(1,nrow=nbs,ncol=1)%*%psi_nu_inc, 1, max) 
	statdec_ks_bs_gms<-apply(-hatnu_bs/(matrix(1,nrow=nbs,ncol=1)%*%sigma_nu) +matrix(1,nrow=nbs,ncol=1)%*%psi_nu_dec, 1, max) 
list(hatnu=hatnu,hatnu_stu=hatnu_stu,statinc_ks=statinc_ks,statinc_ks_bs=statinc_ks_bs,statdec_ks=statdec_ks,statdec_ks_bs=statdec_ks_bs,statinc_ks_bs_gms=statinc_ks_bs_gms,statdec_ks_bs_gms=statdec_ks_bs_gms)
}

MRD_sharp<-function(y,z,bw,Zeval=0,kernel="tria",ctff=0,weight=1){
if (length(Zeval)==1) {
	Zeval.C=Zeval
	Zeval.T=Zeval
} else {
Zeval.C<-Zeval[Zeval<=ctff]
Zeval.T<-Zeval[Zeval>=ctff]
}
if (length(weight)==1) {
	sss=z*0+1
} else {
	sss<-weight
}
n<-length(y)
yC=y[z<ctff]
yT=y[z>=ctff]
zC=z[z<ctff]
zT=z[z>=ctff]
sssC=sss[z<ctff]
sssT=sss[z>=ctff]
yCtrim=y[z<ctff&z>=-bw]
yTtrim=y[z>=ctff&z<=bw]
zCtrim=z[z<ctff&z>=-bw]
zTtrim=z[z>=ctff&z<=bw]
sssCtrim=sss[z<ctff&z>=-bw]
sssTtrim=sss[z>=ctff&z<=bw]
locfitC<-predict(locfit(yC~zC,deg=1,ev=Zeval.C,alpha=c(0,bw),kern=kernel,weight=sssC),se=TRUE,band="local")
evalC=locfitC$fit
evalCse=locfitC$se.fit
residualC<-residuals(locfit(yCtrim~zCtrim,deg=1,alpha=c(0,bw),kern=kernel,weight=sssCtrim))
locfitT<-predict(locfit(yT~zT,deg=1,ev=Zeval.T,alpha=c(0,bw),kern=kernel,weight=sssT),se=TRUE,band="local")
evalT=locfitT$fit
evalTse=locfitT$se.fit
residualT<-residuals(locfit(yTtrim~zTtrim,deg=1,alpha=c(0,bw),kern=kernel,weight=sssTtrim))
estimate=evalT[1]-evalC[length(evalC)]
tstat=estimate/sqrt(evalTse[1]^2+evalCse[length(evalC)]^2)
list(estimate=estimate,tstat=tstat,evalC=evalC,evalT=evalT,evalCse=evalCse,evalTse=evalTse,bw=bw)
}

MRDplot<-function(Zeval,daty,datz,results,labels){
daty<-daty[datz!=0]
datz<-datz[datz!=0]
ones<-daty*0+1
aa= cut(datz,breaks=Zeval, right = FALSE)
aana=is.na(aa)
agg=aggregate(daty[aana==0], by=list(aa[aana==0]),  FUN=mean, na.rm=TRUE)
results$avgy<-agg$x
agg=aggregate(ones[aana==0], by=list(aa[aana==0]),  FUN=sum, na.rm=TRUE)
results$numy<-agg$x
bw=results$bw
plot(Zeval[Zeval<=0],results$evalC,type="l",xlim=c(min(Zeval),max(Zeval)),ylim=c(2*min(results$avgy)-max(results$avgy),1.05*max(results$avgy)-0.05*min(results$avgy)),xlab=labels$xlab, ylab=labels$ylab1,main=labels$main, sub=paste("bandwidth =",as.character(format(bw,digits=3)),sep=""),cex.lab=0.75, cex.sub=0.75,cex.main=0.75)
lines(Zeval[Zeval>=0],results$evalT)
points(Zeval[1:(length(Zeval)-1)],results$avgy)
par(new=TRUE)
barplot(results$numy,xaxt="n",yaxt="n",xlab="",ylab="",ylim=c(0,max(results$numy)*2.5))
axis(4)
mtext(labels$ylab2,side=4,line=0,cex=0.75)
}

