R/counterfactual.r

In Counterfactual: Estimation and Inference Methods for Counterfactual Analysis

counterfactual <- function(formula,data,weights,na.action=na.exclude,group, treatment =FALSE, decomposition = FALSE, counterfactual_var,transformation=FALSE,quantiles = c(1:9)/10,method="qr",trimming=0.005,nreg=100,scale_variable,counterfactual_scale_variable,censoring=0,right=FALSE,nsteps=3,firstc=0.1,secondc=0.05,noboot=FALSE,weightedboot=FALSE,seed=8,robust=FALSE,reps=100,alpha=0.05,first=0.1,last=0.9,cons_test=0,printdeco=TRUE,sepcore = FALSE,ncore=1){
	nreg = round(nreg)
	if(nreg<1){
		stop( "The option nreg must be a strictly positive integer.")}	
	if(missing(group) & missing(counterfactual_var)){
		stop("One of the options group and counterfactual are required")}		
	taus	  = quantiles
		
	if(method!="qr"& method!="loc"& method!="locsca" & method!="cqr" & method!="cox" & method!="logit" & method!="probit" & method!="lpm"){
		stop("The selected method has not been implemented")}			
	if(method=="locsca" & !missing(scale_variable) & missing(group) & missing(counterfactual_scale_variable)){
		stop("If the location scale estimator is used with the option scale, then either the group option or the counterfactual_scale_variable option must be provided.")}
	if(method =="cqr"){
		if(nsteps<3) 		stop( "The options nsteps must be at least 3")
		if(missing(censoring)) 	stop("The option censoring must be provided to use the censored quantile regression estimator.")
	}
	
	if(trimming>0.1){
		stop("Trimming must be less than 10% of the observations.")
	}

	mf 		= match.call()		
	if(!missing(group)){	 
		m 				        = match(c("formula","data","weights", "group"),names(mf),0) 
		mf  				    = mf[c(1,m)]
		mf$drop.unused.levels   = TRUE 
		mf[[1]] 				= as.name("model.frame")
		mf						= eval.parent(mf)		
		mt						= attr(mf,"terms")		
		dep						= model.response(mf,"numeric")	
		reg						= model.matrix(mt,mf)	 
		reg						= as.matrix(reg[,-1])  		
		weight					= model.weights(mf)
		
		if (!is.null(weight) & !is.numeric(weight)) stop("'weights' must be a numeric vector")	
		if(!is.null(weight)) weight = weight else weight = rep(1,length(dep))
		weight = weight/sum(weight)
		
		groupN	           		= mf[,ncol(mf)]	
		group0             		= groupN==0		
		group1	           		= groupN==1

		dep0	           		= subset(dep,group0)			
		reg0               		= subset(reg,group0)
		dep1               		= subset(dep,group1)
		counterfactual_var 		= subset(reg,group1)		
		weight0            		= subset(weight,group0)
		weight1	           		= subset(weight,group1)	
		obs0	           		= length(dep0)
		obs1	           		= length(dep)-obs0		
		
	}else{
		m 						= match(c("formula","data","weights", "counterfactual_var"),names(mf),0)
		mf  					= mf[c(1,m)]
		mf$drop.unused.levels 	= TRUE
		mf[[1]] 				= as.name("model.frame")
		mf						= eval.parent(mf)		
		mt						= attr(mf,"terms")
		dep						= model.response(mf,"numeric")	
		reg						= model.matrix(mt,mf)
		reg		        		= as.matrix(reg[,-1])  
		weight					= model.weights(mf)
		
		if (!is.null(weight) & !is.numeric(weight)) stop("'weights' must be a numeric vector")	
		if(!is.null(weight)) weight = weight else weight = rep(1,length(dep))
		weight = weight/sum(weight)
		
		counterfactual_var  	= mf[,ncol(mf)]

		dep0			    	= dep		
		reg0			    	= reg	
		dep1                	= dep0			
		weight0			    	= weight		
		weight1			    	= weight
	    obs0 			    	= length(dep)
		obs1 			    	= obs0		
	}    
 
	depeval = sort(unique(dep0))
 	if(!missing(nreg) & (nreg!=-1) & (nreg < length(depeval)) ){
		wq         = (1:(nreg-2)-0.5)/(nreg-2)
		seqU       = floor(wq*length(depeval))
		thredIndex = c(1, seqU, length(depeval)) 
		depeval    = depeval[thredIndex]
	} 

	if(treatment){
		depeval_1 = sort(unique(dep1)) 
		if(!missing(nreg) & (nreg!=-1) & (nreg < length(depeval_1)) ){
			wq           = (1:(nreg-2)-0.5)/(nreg-2)
			seqU_1       = floor(wq*length(depeval_1))
			thredIndex_1 = c(1, seqU_1, length(depeval_1)) 
			depeval_1    = depeval_1[thredIndex_1]
		} 
	}else{
		depeval_1 = depeval
	}

	#################################################### Estimation Begins ########################################## 	   
    if(missing(scale_variable))                scale_variable                = reg0
	if(missing(counterfactual_scale_variable)) counterfactual_scale_variable = counterfactual_var	   
	QE_Results = QteDistEst(dep0,depeval,reg0,weight0,counterfactual_var, dep1, depeval_1, treatment, decomposition, weight1,taus,method,trimming,nreg,scale_variable,counterfactual_scale_variable,censoring,right,nsteps,firstc,secondc)
	#################################################################################################################	 		
	quantile_effect 	   	= QE_Results$quantile_effect
	marginal_obs			= QE_Results$marginal_obs
	marginal_fitted			= QE_Results$marginal_fitted
	marginal_counterfactual = QE_Results$marginal_counterfactual
	if(treatment & decomposition){
		composition_effect  = QE_Results$composition_effect
		total_effect        = QE_Results$total_effect
	}	
	
	if(treatment & decomposition){
    	qte_cov_obs = c(quantile_effect, marginal_obs, marginal_fitted, marginal_counterfactual, composition_effect, total_effect, QE_Results$marginal_obs_1, QE_Results$marginal_fitted_1, QE_Results$F_ms)	
	}else{
		qte_cov_obs = c(quantile_effect, marginal_obs, marginal_fitted, marginal_counterfactual)	
	} 
		
	# Inference & Testing 
	if(!noboot){				
		data0 = list(dep=dep0,reg=reg0,weight=weight0,scale_variable=scale_variable)
		data1 = list(counterfactual_var=counterfactual_var,weight=weight1,scale_variable=counterfactual_scale_variable)
		
		qte_cov_booti     = BootstrapProcedure(transformation,reps,data0,depeval,data1,weightedboot,obs0,obs1, dep1, depeval_1, treatment, decomposition, taus,method,trimming,nreg,censoring,right,nsteps,firstc,secondc,seed,sepcore, ncore)
		Bootstrap_Results = InferenceTestingEval(method,qte_cov_booti,reps,taus,robust,qte_cov_obs,last,first,alpha,cons_test, treatment, decomposition)
	}else{		 
		qte 			= quantile_effect
		distributions 	= cbind(marginal_obs,marginal_fitted,marginal_counterfactual)
		if(treatment & decomposition){
			composition_effect   = composition_effect
			total_effect         = total_effect
		}
	}	

	# Display Confidence Interval
	if(printdeco){
		cat("\n")
		cat(format("Conditional Model: ",width=40))
		if(method=="qr") 		cat("linear quantile regression\n")
		if(method=="loc") 		cat("location model\n")
		if(method=="locsca")	cat("location scale model\n")
		if(method=="cqr") 		cat("linear censored quantile regression\n")
		
		if(method=="logit") 	cat("logit\n")
		if(method=="probit")	cat("probit\n")
		if(method=="lpm") 		cat("linear probability model\n")		
		if(method=="cox") 		cat("cox duration model\n")
		if(method!="cox") 		cat(format("Number of regressions estimated:",width=40),QE_Results$nreg,"\n\n")	
		
		if(!noboot){			
			
			cat("The variance has been estimated by bootstraping the results",reps, "times.\n\n")
			cat(format("No. of obs. in the reference group:",width=40),obs0,"\n")
			cat(format("No. of obs. in the counterfactual group:",width=40),obs1,"\n\n\n")			
			
			if(!treatment & !decomposition){ ## CE
				bootstrap_CI_result_print_ce = Bootstrap_Results$resCE	
				
				cat(format("Quantile Effects -- Composition", width=70,justify = "centre"),"\n")
				cat(paste(rep("-",70),collapse=""),"\n")
				cat(format("",width=20),format("Pointwise",width=8,justify="centre"), format("Pointwise",width=20,justify="centre"),format("Functional",width=18,justify="centre"),"\n")
				if(robust){
					cat(format("Quantile",width=10,justify="centre"),format("Est.",width=10,justify = "centre"), format("Robust S.E.",width=8,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=20,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=18,justify="centre"),"\n")					
				}else{		
					cat(format("Quantile",width=10,justify="centre"),format("Est.",width=10,justify = "centre"), format("Std.Err",width=8,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=20,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=18,justify="centre"),"\n")			
				}
				for(i in 1:length(taus)){
 					cat(format(taus[i],width=8,justify="left"),formatC(bootstrap_CI_result_print_ce[i,],digits=3, width=10),"\n",sep="")
				}
				if(printdeco){		
					Counter_Test_CE = Bootstrap_Results$testCE	
					row_of_test     = nrow(Counter_Test_CE)	
				    constest1       = sort(unique(cons_test))	
				    constestNo0     = setdiff(constest1,0) 
				    nc              = length(constestNo0)		

				    cat("\n\n")
					cat(format("Bootstrap inference on the counterfactual quantile process",width=70,justify="centre"),"\n")
					cat(paste(rep("-",70),collapse=""),"\n")
					cat("\t\t\t\t\t\t    P-values\t\n")
					cat("\t\t\t\t\t       ",paste(rep("=",22),collapse=""),"\n")
					cat(format("NULL-hypthoesis",width=38),"\t", format("KS-statistic",width=12),format("CMS-statistic",width=12),"\n")
					cat(paste(rep("=",70),collapse=""),"\n")
					cat(format("Correct specification of the parametric model",width=40),format(Counter_Test_CE[1,1],digits=2,nsmall =2, width=10,justify="centre"),format(Counter_Test_CE[1,2],width=10,digits=2, nsmall =2),"\n")	
					cat(format("No effect: QE(tau)=0 for all taus            ",width=40),format(Counter_Test_CE[2,1],digits=2,nsmall =2, width=10,justify="centre"),format(Counter_Test_CE[2,2],width=10,digits=2, nsmall =2),"\n")
					if(nc>0){			
						for(i in 1:nc){				
							cat(format(paste("Constant effect: QE(tau)=",constestNo0[i],"for all taus"),width=40),format(Counter_Test_CE[i+2,1],digits=2,nsmall =2,width=10),format(Counter_Test_CE[i+2,2],width=10,digits=2, nsmall =2),"\n")
						}
					}
					cat(format("Constant effect: QE(tau)=QE(0.5) for all taus",width=40),format(Counter_Test_CE[row_of_test-2,1],digits=2, nsmall =2, width=10),format(Counter_Test_CE[row_of_test-2,2],width=10,digits=2, nsmall =2),"\n")
					cat(format("Stochastic dominance: QE(tau)>0 for all taus ",width=40),format(Counter_Test_CE[row_of_test-1,1],digits=2, nsmall =2, width=10),format(Counter_Test_CE[row_of_test-1,2],width=10,digits=2, nsmall =2),"\n")
					cat(format("Stochastic dominance: QE(tau)<0 for all taus ",width=40),format(Counter_Test_CE[row_of_test,1],digits=2, nsmall =2, width=10),format(Counter_Test_CE[row_of_test,2],width=10,digits=2, nsmall =2),"\n")		
				}
			
			}else if(treatment & !decomposition){  
				bootstrap_CI_result_print_se = Bootstrap_Results$resSE	
				
				cat(format("Quantile Effects -- Structure", width=70,justify = "centre"),"\n")
				cat(paste(rep("-",70),collapse=""),"\n")
				cat(format("",width=20),format("Pointwise",width=8,justify="centre"), format("Pointwise",width=20,justify="centre"),format("Functional",width=18,justify="centre"),"\n")
				if(robust){
					cat(format("Quantile",width=10,justify="centre"),format("Est.",width=10,justify = "centre"), format("Robust S.E.",width=8,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=20,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=18,justify="centre"),"\n")					
				}else{		
					cat(format("Quantile",width=10,justify="centre"),format("Est.",width=10,justify = "centre"), format("Std.Err",width=8,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=20,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=18,justify="centre"),"\n")			
				}
				for(i in 1:length(taus)){
					cat(format(taus[i],width=8,justify="left"),formatC(bootstrap_CI_result_print_se[i,],digits=3, width=10),"\n",sep="")
				}
				if(printdeco){		
					Counter_Test_SE = Bootstrap_Results$testSE		
					row_of_test     = nrow(Counter_Test_SE)	
				    constest1       = sort(unique(cons_test))	
				    constestNo0     = setdiff(constest1,0) 
				    nc              = length(constestNo0)			
						
					cat("\n\n")
					cat(format("Bootstrap inference on the counterfactual quantile process",width=70,justify="centre"),"\n")
					cat(paste(rep("-",70),collapse=""),"\n")
					cat("\t\t\t\t\t\t    P-values\t\n")
					cat("\t\t\t\t\t       ",paste(rep("=",22),collapse=""),"\n")
					cat(format("NULL-hypthoesis",width=38),"\t", format("KS-statistic",width=12),format("CMS-statistic",width=12),"\n")
					cat(paste(rep("=",70),collapse=""),"\n")
					cat(format("Correct specification of the parametric model",width=40),format(Counter_Test_SE[1,1],digits=2, nsmall =2, width=10,justify="centre"),format(Counter_Test_SE[1,2],width=10,digits=2, nsmall =2),"\n")	
					cat(format("No effect: QE(tau)=0 for all taus            ",width=40),format(Counter_Test_SE[2,1],digits=2, nsmall =2, width=10,justify="centre"),format(Counter_Test_SE[2,2],width=10,digits=2, nsmall =2),"\n")
					if(nc>0){			
						for(i in 1:nc){				
							cat(format(paste("Constant effect: QE(tau)=",constestNo0[i],"for all taus"),width=40),format(Counter_Test_SE[i+2,1],digits=2, nsmall =2, width=10),format(Counter_Test_SE[i+2,2],width=10,digits=2, nsmall =2),"\n")
						}
					}
					cat(format("Constant effect: QE(tau)=QE(0.5) for all taus",width=40),format(Counter_Test_SE[row_of_test-2,1],digits=2,nsmall =2, width=10),format(Counter_Test_SE[row_of_test-2,2],width=10,digits=2, nsmall =2),"\n")
					cat(format("Stochastic dominance: QE(tau)>0 for all taus ",width=40),format(Counter_Test_SE[row_of_test-1,1],digits=2,nsmall =2, width=10),format(Counter_Test_SE[row_of_test-1,2],width=10,digits=2, nsmall =2),"\n")
					cat(format("Stochastic dominance: QE(tau)<0 for all taus ",width=40),format(Counter_Test_SE[row_of_test,1],digits=2, nsmall =2, width=10),format(Counter_Test_SE[row_of_test,2],width=10,digits=2, nsmall =2),"\n")		
				}

			}else if(treatment & decomposition){
				
				## SE
				bootstrap_CI_result_print_se = Bootstrap_Results$resSE	
			 
				cat(format("Quantile Effects -- Structure", width=70,justify = "centre"),"\n")
				cat(paste(rep("-",70),collapse=""),"\n")
				cat(format("",width=20),format("Pointwise",width=8,justify="centre"), format("Pointwise",width=20,justify="centre"),format("Functional",width=18,justify="centre"),"\n")
				if(robust){
					cat(format("Quantile",width=10,justify="centre"),format("Est.",width=10,justify = "centre"), format("Robust S.E.",width=8,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=20,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=18,justify="centre"),"\n")					
				}else{		
					cat(format("Quantile",width=10,justify="centre"),format("Est.",width=10,justify = "centre"), format("Std.Err",width=8,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=20,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=18,justify="centre"),"\n")			
				}
				for(i in 1:length(taus)){
					cat(format(taus[i],width=8,justify="left"),formatC(bootstrap_CI_result_print_se[i,],digits=3, width=10),"\n",sep="")
				}
				if(printdeco){		
					Counter_Test_SE = Bootstrap_Results$testSE		
					row_of_test     = nrow(Counter_Test_SE)	
				    constest1       = sort(unique(cons_test))	
				    constestNo0     = setdiff(constest1,0) 
				    nc              = length(constestNo0)			
						
					cat("\n\n")
					cat(format("Bootstrap inference on the counterfactual quantile process",width=70,justify="centre"),"\n")
					cat(paste(rep("-",70),collapse=""),"\n")
					cat("\t\t\t\t\t\t    P-values\t\n")
					cat("\t\t\t\t\t       ",paste(rep("=",22),collapse=""),"\n")
					cat(format("NULL-hypthoesis",width=38),"\t", format("KS-statistic",width=12),format("CMS-statistic",width=12),"\n")
					cat(paste(rep("=",70),collapse=""),"\n")
					cat(format("Correct specification of the parametric model",width=40),format(Counter_Test_SE[1,1],digits=2,width=10,justify="centre"),format(Counter_Test_SE[1,2],width=10,digits=2),"\n")	
					cat(format("No effect: QE(tau)=0 for all taus            ",width=40),format(Counter_Test_SE[2,1],digits=2,width=10,justify="centre"),format(Counter_Test_SE[2,2],width=10,digits=2),"\n")
					if(nc>0){			
						for(i in 1:nc){				
							cat(format(paste("Constant effect: QE(tau)=",constestNo0[i],"for all taus"),width=40),format(Counter_Test_SE[i+2,1],digits=2,width=10),format(Counter_Test_SE[i+2,2],width=10,digits=2),"\n")
						}
					}
					cat(format("Constant effect: QE(tau)=QE(0.5) for all taus",width=40),format(Counter_Test_SE[row_of_test-2,1],digits=2,width=10),format(Counter_Test_SE[row_of_test-2,2],width=10,digits=2),"\n")
					cat(format("Stochastic dominance: QE(tau)>0 for all taus ",width=40),format(Counter_Test_SE[row_of_test-1,1],digits=2,width=10),format(Counter_Test_SE[row_of_test-1,2],width=10,digits=2),"\n")
					cat(format("Stochastic dominance: QE(tau)<0 for all taus ",width=40),format(Counter_Test_SE[row_of_test,1],digits=2,width=10),format(Counter_Test_SE[row_of_test,2],width=10,digits=2),"\n")		
				}

				## CE 
				bootstrap_CI_result_print_ce = Bootstrap_Results$resCE	
				cat("\n\n")
				cat(format("Quantile Effects -- Composition", width=70,justify = "centre"),"\n")
				cat(paste(rep("-",70),collapse=""),"\n")
				cat(format("",width=20),format("Pointwise",width=8,justify="centre"), format("Pointwise",width=20,justify="centre"),format("Functional",width=18,justify="centre"),"\n")
				if(robust){
					cat(format("Quantile",width=10,justify="centre"),format("Est.",width=10,justify = "centre"), format("Robust S.E.",width=8,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=20,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=18,justify="centre"),"\n")					
				}else{		
					cat(format("Quantile",width=10,justify="centre"),format("Est.",width=10,justify = "centre"), format("Std.Err",width=8,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=20,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=18,justify="centre"),"\n")			
				}
				for(i in 1:length(taus)){
					cat(format(taus[i],width=8,justify="left"),formatC(bootstrap_CI_result_print_ce[i,],digits=3, width=10),"\n",sep="")
				}
				if(printdeco){		
					Counter_Test_CE = Bootstrap_Results$testCE			
				    cat("\n\n")
					cat(format("Bootstrap inference on the counterfactual quantile process",width=70,justify="centre"),"\n")
					cat(paste(rep("-",70),collapse=""),"\n")
					cat("\t\t\t\t\t\t    P-values\t\n")
					cat("\t\t\t\t\t       ",paste(rep("=",22),collapse=""),"\n")
					cat(format("NULL-hypthoesis",width=38),"\t", format("KS-statistic",width=12),format("CMS-statistic",width=12),"\n")
					cat(paste(rep("=",70),collapse=""),"\n")
					cat(format("Correct specification of the parametric model",width=40),format(Counter_Test_CE[1,1],digits=2,width=10,justify="centre"),format(Counter_Test_CE[1,2],width=10,digits=2),"\n")	
					cat(format("No effect: QE(tau)=0 for all taus            ",width=40),format(Counter_Test_CE[2,1],digits=2,width=10,justify="centre"),format(Counter_Test_CE[2,2],width=10,digits=2),"\n")
					if(nc>0){			
						for(i in 1:nc){				
							cat(format(paste("Constant effect: QE(tau)=",constestNo0[i],"for all taus"),width=40),format(Counter_Test_CE[i+2,1],digits=2,width=10),format(Counter_Test_CE[i+2,2],width=10,digits=2),"\n")
						}
					}
					cat(format("Constant effect: QE(tau)=QE(0.5) for all taus",width=40),format(Counter_Test_CE[row_of_test-2,1],digits=2,width=10),format(Counter_Test_CE[row_of_test-2,2],width=10,digits=2),"\n")
					cat(format("Stochastic dominance: QE(tau)>0 for all taus ",width=40),format(Counter_Test_CE[row_of_test-1,1],digits=2,width=10),format(Counter_Test_CE[row_of_test-1,2],width=10,digits=2),"\n")
					cat(format("Stochastic dominance: QE(tau)<0 for all taus ",width=40),format(Counter_Test_CE[row_of_test,1],digits=2,width=10),format(Counter_Test_CE[row_of_test,2],width=10,digits=2),"\n")		
				}

				## TE
				bootstrap_CI_result_print_te = Bootstrap_Results$resTE	
				cat("\n\n")
				cat(format("Quantile Effects -- Total", width=70,justify = "centre"),"\n")
				cat(paste(rep("-",70),collapse=""),"\n")
				cat(format("",width=20),format("Pointwise",width=8,justify="centre"), format("Pointwise",width=20,justify="centre"),format("Functional",width=18,justify="centre"),"\n")
				if(robust){
					cat(format("Quantile",width=10,justify="centre"),format("Est.",width=10,justify = "centre"), format("Robust S.E.",width=8,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=20,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=18,justify="centre"),"\n")					
				}else{		
					cat(format("Quantile",width=10,justify="centre"),format("Est.",width=10,justify = "centre"), format("Std.Err",width=8,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=20,justify="centre"),format(paste((1-alpha)*100,"% Conf.Interval",sep=""),width=18,justify="centre"),"\n")			
				}
				for(i in 1:length(taus)){
					cat(format(taus[i],width=8,justify="left"),formatC(bootstrap_CI_result_print_te[i,],digits=3, width=10),"\n",sep="")
				}
				if(printdeco){		
					Counter_Test_TE = Bootstrap_Results$testTE		
					cat("\n\n")
					cat(format("Bootstrap inference on the counterfactual quantile process",width=70,justify="centre"),"\n")
					cat(paste(rep("-",70),collapse=""),"\n")
					cat("\t\t\t\t\t\t    P-values\t\n")
					cat("\t\t\t\t\t       ",paste(rep("=",22),collapse=""),"\n")
					cat(format("NULL-hypthoesis",width=38),"\t", format("KS-statistic",width=12),format("CMS-statistic",width=12),"\n")
					cat(paste(rep("=",70),collapse=""),"\n")
					cat(format("Correct specification of the parametric model",width=40),format(Counter_Test_TE[1,1],digits=2,width=10,justify="centre"),format(Counter_Test_TE[1,2],width=10,digits=2),"\n")	
					cat(format("No effect: QE(tau)=0 for all taus            ",width=40),format(Counter_Test_TE[2,1],digits=2,width=10,justify="centre"),format(Counter_Test_TE[2,2],width=10,digits=2),"\n")
					if(nc>0){			
						for(i in 1:nc){				
							cat(format(paste("Constant effect: QE(tau)=",constestNo0[i],"for all taus"),width=40),format(Counter_Test_TE[i+2,1],digits=2,width=10),format(Counter_Test_TE[i+2,2],width=10,digits=2),"\n")
						}
					}
					cat(format("Constant effect: QE(tau)=QE(0.5) for all taus",width=40),format(Counter_Test_TE[row_of_test-2,1],digits=2,width=10),format(Counter_Test_TE[row_of_test-2,2],width=10,digits=2),"\n")
					cat(format("Stochastic dominance: QE(tau)>0 for all taus ",width=40),format(Counter_Test_TE[row_of_test-1,1],digits=2,width=10),format(Counter_Test_TE[row_of_test-1,2],width=10,digits=2),"\n")
					cat(format("Stochastic dominance: QE(tau)<0 for all taus ",width=40),format(Counter_Test_TE[row_of_test,1],digits=2,width=10),format(Counter_Test_TE[row_of_test,2],width=10,digits=2),"\n")		
				}
			}	

		}else{			
			cat("The variance has not been computed.\n")
			cat("Do not turn the option noboot on if you want to compute it.\n\n")
			cat(format("No. of obs. in the reference group:",width=40),obs0,"\n")
			cat(format("No. of obs. in the counterfactual group:",width=40),obs1,"\n\n\n")	
			
			if(treatment){
				cat(format("Quantile Effects -- Structure", width=50,justify = "centre"),"\n")
			}else{
				cat(format("Quantile Effects -- Composition", width=50,justify = "centre"),"\n")
			}
			cat(paste(rep("-",50),collapse=""),"\n")
			cat(format("Quantile",width=18,justify="right"),format("Est.",width=20,justify = "centre"),"\n")	
			for(i in 1:length(quantile_effect)){
				cat(format(taus[i],width=15,justify="centre"),"\t",format(quantile_effect[i],justify="centre",width=15,digits=2),"\n",sep="")
			}	
			
			if(treatment & decomposition){

				cat("\n\n")
				cat(format("Quantile Effects -- Composition", width=50,justify = "centre"),"\n")
				cat(paste(rep("-",50),collapse=""),"\n")
				cat(format("Quantile",width=18,justify="right"),format("Est.",width=20,justify = "centre"),"\n")	
				for(i in 1:length(composition_effect)){
					cat(format(taus[i],width=15,justify="centre"),"\t",format(composition_effect[i],justify="centre",width=15,digits=2),"\n",sep="")
				}	

				cat("\n\n")
				cat(format("Quantile Effects -- Total", width=50,justify = "centre"),"\n")
				cat(paste(rep("-",50),collapse=""),"\n")
				cat(format("Quantile",width=18,justify="right"),format("Est.",width=20,justify = "centre"),"\n")	
				for(i in 1:length(total_effect)){
					cat(format(taus[i],width=15,justify="centre"),"\t",format(total_effect[i],justify="centre",width=15,digits=2),"\n",sep="")
				}	
			}
		}
	}
	
	## returns
	results = NULL 
	
 	if(treatment & decomposition){
		results = list(quantiles = taus, structral_effect = quantile_effect, composition_effect = composition_effect, total_effect = total_effect, nreg=nreg)
	}else if(treatment & !decomposition){
		results = list(quantiles = taus, structral_effect = quantile_effect, nreg=nreg)
	}else{
		results = list(quantiles = taus, composition_effect = quantile_effect, nreg=nreg)
	}

	if(!noboot){
		results = c(results, Bootstrap_Results)
		class(results) <- c("counterfactual", class(results)) 	
	}
	return(results)
}

QteDistEst <- function(dep0,depeval,reg0,weight0,counterfactual_var, dep1, depeval_1, treatment, decomposition, weight1,taus,method,trimming,nreg,scale_variable,counterfactual_scale_variable,censoring,right,nsteps,firstc,secondc){	
	reg1 = counterfactual_var
	if(method == "qr"){		
		wq	 					= (1:nreg-0.5)/nreg	
        wqNew                   = wq[(wq>=trimming)&wq<=(1-trimming)]	

		wqW                     = rep(1,length(wqNew))
        wei0                    = as.matrix(weight0)%*%wqW		
		wei1                    = as.matrix(weight1)%*%wqW

		conditional_coef	    = coef(rq(dep0~reg0,tau=wqNew,weights=weight0))	
		marginal_fitted  		= wtd.quantile(as.vector(trimming+cbind(1,reg0)%*%conditional_coef), as.vector(wei0), probs=taus, na.rm=TRUE)
		marginal_counterfactual = wtd.quantile(as.vector(trimming+cbind(1,counterfactual_var)%*%conditional_coef), as.vector(wei1), probs=taus, na.rm=TRUE)	
		
		if(treatment){
			conditional_coef_1	= coef(rq(dep1~reg1,tau=wqNew,weights=weight1))	
        	marginal_fitted_1  	= wtd.quantile(as.vector(trimming+cbind(1,reg1)%*%conditional_coef_1), as.vector(wei1), probs=taus,na.rm=TRUE)
		}	
		if(treatment&decomposition){

			y_thred         = c(depeval, depeval_1)
			y_te            = c(dep0, dep1)
			 
			conditionTaus_0 =  colMeans(apply(as.matrix(y_thred), 1, function(xxx) {rowMeans(ifelse((trimming+cbind(1,reg0)%*%conditional_coef) <= xxx, 1, 0))}))
			conditionTaus_1 =  colMeans(apply(as.matrix(y_thred), 1, function(xxx) {rowMeans(ifelse((trimming+cbind(1,reg1)%*%conditional_coef_1) <= xxx, 1, 0))}))
			
			F_0hat          = length(dep0)/(length(dep0) + length(dep1))
			F_est           = conditionTaus_0*F_0hat + conditionTaus_1*(1- F_0hat)
			 
			F_obs           = colMeans(apply(as.matrix(y_thred), 1, function(xxx) {ifelse( y_te <= xxx, 1, 0) }) ) 
			F_ms            = F_obs - F_est
		}
	}

	if(method == "loc"){		
		wq                  	= (1:nreg-0.5)/nreg

		wqW                     = rep(1,length(wq))
        wei0                    = as.matrix(weight0)%*%wqW		
		wei1                    = as.matrix(weight1)%*%wqW
		 
		loc                 	= lm(dep0~reg0,weights=weight0)		
		conditional_coef    	= kronecker(as.matrix(coef(loc)),matrix(1,1,nreg)) 			
		conditional_coef[1,]	= conditional_coef[1,]+wtd.quantile(as.matrix(resid(loc)),weight0,wq,na.rm=TRUE)	
		
        marginal_fitted  		= wtd.quantile(as.vector(cbind(1,reg0)%*%conditional_coef),as.vector(wei0),taus,na.rm=TRUE)
		marginal_counterfactual = wtd.quantile(as.vector(cbind(1,counterfactual_var)%*%conditional_coef),as.vector(wei1),taus,na.rm=TRUE)	
		
		if(treatment){
			loc_1                 	= lm(dep1~reg1, weights=weight1)		
			conditional_coef_1    	= kronecker(as.matrix(coef(loc_1)), matrix(1,1,nreg)) 			
			conditional_coef_1[1,]	= conditional_coef_1[1,] + wtd.quantile(as.matrix(resid(loc_1)), weight1, wq, na.rm=TRUE)	
	
			marginal_fitted_1  	    = wtd.quantile(as.vector(cbind(1,reg1)%*%conditional_coef_1),as.vector(wei1),taus,na.rm=TRUE)
		}			
		if(treatment&decomposition){

			y_thred         = c(depeval, depeval_1)
			y_te            = c(dep0, dep1)
			 
			conditionTaus_0 =  colMeans(apply(as.matrix(y_thred), 1, function(xxx) {rowMeans(ifelse((cbind(1,reg0)%*%conditional_coef) <= xxx, 1, 0))}))
			conditionTaus_1 =  colMeans(apply(as.matrix(y_thred), 1, function(xxx) {rowMeans(ifelse((cbind(1,reg1)%*%conditional_coef_1) <= xxx, 1, 0))}))
			
			F_0hat          = length(dep0)/(length(dep0) + length(dep1))
			F_est           = conditionTaus_0*F_0hat + conditionTaus_1*(1- F_0hat)
			 
			F_obs           = colMeans(apply(as.matrix(y_thred), 1, function(xxx) {ifelse( y_te <= xxx, 1, 0) }) ) 
			F_ms            = F_obs - F_est
		}
	}

	if(method == "locsca"){				
		wq	 	        		= (1:nreg-0.5)/nreg	

		wqW                     = rep(1,length(wq))
        wei0                    = as.matrix(weight0)%*%wqW		
		wei1                    = as.matrix(weight1)%*%wqW
		 
		loc		        		= lm(dep0~reg0,weights=weight0)	
		conditional_resid       = resid(loc)
		lmR                     = lm(log(conditional_resid^2)~scale_variable, weights=weight0)		
		predsca         		= sqrt(exp(predict(lmR)))			
		residN          		= wtd.quantile(conditional_resid/predsca, weight0, wq, na.rm=TRUE)	
		
		predall	        		= kronecker(predict(loc),matrix(1,1,nreg))+kronecker(predsca,t(residN)) 
		marginal_fitted			= wtd.quantile(as.vector(predall),wei0,probs=taus,na.rm=TRUE)	
				
		predCounter  			= cbind(1,counterfactual_var)%*%coef(loc)
		predscaCounter  		= sqrt(exp(cbind(1,counterfactual_scale_variable)%*%coef(lmR)))
		predallCounter			= kronecker(predCounter,matrix(1,1,nreg))+kronecker(predscaCounter,t(residN)) 		
		marginal_counterfactual	= wtd.quantile(as.vector(predallCounter),wei1,probs=taus,na.rm=TRUE)
		
		if(treatment){
			
			loc_1		        	= lm(dep1~reg1, weights=weight1)	
			conditional_resid_1     = resid(loc_1)
			lmR_1                   = lm(log(conditional_resid_1^2)~counterfactual_scale_variable, weights=weight1)		
			predsca_1         		= sqrt(exp(predict(lmR_1)))			
			residN_1          		= wtd.quantile(conditional_resid_1/predsca_1, weight1, wq, na.rm=TRUE)	
			
			predall_1	        	= kronecker(predict(loc_1),matrix(1,1,nreg))+kronecker(predsca_1,t(residN_1)) 
			marginal_fitted_1		= wtd.quantile(as.vector(predall_1), wei1, probs=taus, na.rm=TRUE)	
		}	

		if(treatment&decomposition){

			y_thred         = c(depeval, depeval_1)
			y_te            = c(dep0, dep1)
			 
			conditionTaus_0 =  colMeans(apply(as.matrix(y_thred), 1, function(xxx) {rowMeans(ifelse(predall <= xxx, 1, 0))}))
			conditionTaus_1 =  colMeans(apply(as.matrix(y_thred), 1, function(xxx) {rowMeans(ifelse(predall_1 <= xxx, 1, 0))}))
			
			F_0hat          = length(dep0)/(length(dep0) + length(dep1))
			F_est           = conditionTaus_0*F_0hat + conditionTaus_1*(1- F_0hat)
			 
			F_obs           = colMeans(apply(as.matrix(y_thred), 1, function(xxx) {ifelse( y_te <= xxx, 1, 0) }) ) 
			F_ms            = F_obs - F_est
		}
    }	

	if(method == "cqr"){						
		wq              = (1:nreg-0.5)/nreg	
				
		est_cqrl <- function(dep,reg,weight,wq,censoring,nsteps,c1,c2){					
			coef_logit	= coef(glm((dep>censoring)~reg,weights=weight,family=binomial(link=logit)))
			pred		= plogis(cbind(1,reg)%*%coef_logit)
			
			coef_wq		= NULL		
			for(i in wq){		
				delta1   = quantile(pred[pred>1-i],c1)
				select1  = (pred>=delta1) 
				temp     = coef(rq(dep[select1]~as.matrix(reg)[select1,],i,weights=weight[select1]))
								
				steps	= 3
				while(steps <= nsteps){	
					pred1   = cbind(1,reg)%*%temp
					delta2  = quantile((pred1-censoring)[pred1>censoring],c2)
					select2 = (pred1>=delta2) 
					temp    = coef(rq(dep[select2]~as.matrix(reg)[select2,],i,weights=weight[select2]))						
					steps	= steps + 1						
					select1 = select2
				}
				coef_wq = cbind(coef_wq,temp)		
			}	           		
			return(coef_wq)
		}			
		if(!right){
			conditional_coef    = est_cqrl(dep0, reg0, weight0, wq, censoring, nsteps, firstc, secondc)
		}else{	
			conditional_coef	= -est_cqrl(-dep0, reg0, weight0, 1-wq, -censoring, nsteps, firstc, secondc)			
		}			  
       	
		marginal_fitted         = quantile(as.vector(cbind(1,reg0)%*%conditional_coef), probs=taus, na.rm=TRUE)
		marginal_counterfactual = quantile(as.vector(cbind(1,counterfactual_var)%*%conditional_coef), probs=taus, na.rm=TRUE)	
		
		if(treatment){
			if(!right){
				conditional_coef_1  = est_cqrl(dep1, reg1, weight1, wq, censoring, nsteps, firstc, secondc)
			}else{	
				depR_1              = -dep1
				censoringR_1        = -censoring 
				conditional_coef_1	= -est_cqrl(depR_1, reg1, weight1, 1-wq, censoringR_1, nsteps, firstc, secondc)			
			}			         		
			marginal_fitted_1  	    = quantile(as.vector(cbind(1,reg1)%*%conditional_coef_1), probs=taus, na.rm=TRUE)
		}

		if(treatment&decomposition){

			y_thred         = c(depeval, depeval_1)
			y_te            = c(dep0, dep1)
			 
			conditionTaus_0 =  colMeans(apply(as.matrix(y_thred), 1, function(xxx) {rowMeans(ifelse((cbind(1,reg0)%*%conditional_coef) <= xxx, 1, 0))}))
			conditionTaus_1 =  colMeans(apply(as.matrix(y_thred), 1, function(xxx) {rowMeans(ifelse((cbind(1,reg1)%*%conditional_coef_1) <= xxx, 1, 0))}))
			
			F_0hat          = length(dep0)/(length(dep0) + length(dep1))
			F_est           = conditionTaus_0*F_0hat + conditionTaus_1*(1- F_0hat)

			F_obs           = colMeans(apply(as.matrix(y_thred), 1, function(xxx) {ifelse( y_te <= xxx, 1, 0) }) ) 
			F_ms            = F_obs - F_est
		}

	}

	if(method == "cox"|method=="logit"|method=="probit"|method=="lpm"){
		
		getquantile <- function(depevalg,pred,tausg){	
			Q		= NULL
			for(i in 1:length(tausg)){
				Q	= rbind(Q,depevalg[min(sum(pred<=tausg[i])+1,length(depevalg))]) # left inverse of depevalg, from pred <= tausg[i]					
			}
			return(Q)
		}	

	    if(method=="cox"){
	    	
	    	if(min(dep0)<0) stop("Only positive dependent variable allowed")
	    	
			fit0	                = coxph(Surv(dep0)~reg0, weights=weight0, ties = "breslow")
			depeval                 = basehaz(fit0)$time
			coef0                   = coef(fit0)
			S0                      = (survfit(fit0)$surv)^exp(-colMeans(reg0)%*%coef0)   
	        NS0                     = length(S0)

			reg0                    = as.matrix(reg0)
			counterfactual_var      = as.matrix(counterfactual_var)
			index0                  = exp(reg0%*%coef(fit0))
			index1                  = exp(counterfactual_var%*%coef0)
			
			NReg0                   = nrow(reg0)
			Sc0                     = kronecker(matrix(S0,1,NS0), matrix(1,NReg0,1))^kronecker(matrix(index0,NReg0,1), matrix(1,1,NS0))
			F0                      = colMeans(1-Sc0)
			marginal_fitted         = getquantile(depeval, F0, taus) 
			
			NReg1                   = nrow(counterfactual_var)
			Sc1                     = kronecker(matrix(S0,1,NS0), matrix(1,NReg1,1))^kronecker(matrix(index1,NReg1,1), matrix(1,1,NS0))
			F1                      = colMeans(1-Sc1)
			marginal_counterfactual = getquantile(depeval, F1, taus) 
			
			if(treatment){
				fit_1	            = coxph(Surv(dep1)~reg1, weights=weight1, ties = "breslow")
				depeval_1           = basehaz(fit_1)$time
				coef_1              = coef(fit_1)
				reg1                = as.matrix(reg1)
				S_1                 = (survfit(fit_1)$surv)^exp(-colMeans(reg1)%*%coef_1)   
		        NS_1                = length(S_1)

				index_1             = exp(reg1%*%coef(fit_1))
				
				NReg_1              = nrow(reg1)
				Sc0_1               = kronecker(matrix(S_1,1,NS_1), matrix(1,NReg_1,1))^kronecker(matrix(index_1, NReg_1, 1), matrix(1,1,NS_1))
				F0_1                = colMeans(1-Sc0_1)
				marginal_fitted_1  	= getquantile(depeval_1, F0_1, taus) 
			}	

			if(treatment&decomposition){

				y_thred        = c(depeval, depeval_1)
				thredIndex     = c(1, floor((1:(nreg-2)-0.5)/(nreg-2)*length(y_thred)), length(y_thred)) 
				y_thred_approx = y_thred[thredIndex]
				y_te           = c(dep0, dep1)
				F_obs          = colMeans(apply(as.matrix(y_thred_approx), 1, function(xxx) {ifelse( y_te <= xxx, 1, 0) }) ) 
 
			 	F_approx_fun0  = approxfun(depeval,  F0, yleft=0, yright=1)
				F_approx_fun1  = approxfun(depeval_1, F0_1, yleft=0, yright=1)

			 	F_0hat         = length(dep0)/(length(dep0) + length(dep1))
				F_est          = F_approx_fun0(y_thred_approx)*F_0hat + F_approx_fun1(y_thred_approx)*(1-F_0hat)

				F_ms           = F_obs - F_est
			}
		}	
		
		if(method=="logit"|method=="probit"|method=="lpm"){	
			
			dep2_0          = depeval 
			nreg			= length(dep2_0)	
		 		
		 	if(treatment){
				 dep2_1 = depeval_1
			}

			if(method == "logit"){		
				
				conditional_coef 		= apply(as.matrix(dep2_0), 1, function(xxx) coef(glm((dep0<=xxx)~reg0, weights=weight0, family=binomial(link="logit"))))
				conditional_coef        = conditional_coef[, colSums(!is.nan(conditional_coef))>0]
				
				pred	     			= plogis(cbind(1,reg0)%*%conditional_coef) 
				conditional_fitted 		= sort(apply(pred, 2, function(a) wtd.mean(a, weights=weight0)))
				if(max(conditional_fitted) == 1 ){
					conditional_fitted_s = conditional_fitted
				}else{
					conditional_fitted_s = c(conditional_fitted,1)
				}
				marginal_fitted	        = getquantile(dep2_0, conditional_fitted_s, taus) 

				predCounter 			  = plogis(cbind(1,counterfactual_var)%*%conditional_coef) 
				conditional_fittedCounter = sort(apply(predCounter, 2, function(a) wtd.mean(a, weights=weight1)))
				if(max(conditional_fittedCounter) == 1){
					predCounter_s       = conditional_fittedCounter
				}else{
					predCounter_s       = c(conditional_fittedCounter, 1)
				}
				marginal_counterfactual = getquantile(dep2_0, predCounter_s, taus) 
				
				if(treatment){
				
					conditional_coef_1 	= apply(as.matrix(dep2_1), 1, function(xxx) coef(glm((dep1<=xxx)~reg1, weights=weight1, family=binomial(link="logit"))))
					conditional_coef_1  = conditional_coef_1[, colSums(!is.nan(conditional_coef_1))>0]
					
					pred_1	     		       = plogis(cbind(1,reg1)%*%conditional_coef_1) 
					conditional_fitted_1       = sort(apply(pred_1, 2, function(a) wtd.mean(a, weights=weight1)))
					if(max(conditional_fitted_1) ==1 ){
						conditional_fitted_1_s = conditional_fitted_1
					}else{
						conditional_fitted_1_s = c(conditional_fitted_1,1)
					}
					marginal_fitted_1  	= getquantile(dep2_1, conditional_fitted_1_s, taus) 
				}	
				if(treatment&decomposition){

					y_thred        = c(depeval, depeval_1)
					thredIndex     = c(1, floor((1:(nreg-2)-0.5)/(nreg-2)*length(y_thred)), length(y_thred)) 
					y_thred_approx = y_thred[thredIndex]
					y_te           = c(dep0, dep1)
					F_obs          = colMeans(apply(as.matrix(y_thred_approx), 1, function(xxx) {ifelse( y_te <= xxx, 1, 0) }) ) 
 
			 		F_approx_fun0  = approxfun(dep2_0, conditional_fitted, yleft=0, yright=1)
					F_approx_fun1  = approxfun(dep2_1, conditional_fitted_1, yleft=0, yright=1)

			 		F_0hat         = length(dep0)/(length(dep0) + length(dep1))
					F_est          = F_approx_fun0(y_thred_approx)*F_0hat + F_approx_fun1(y_thred_approx)*(1-F_0hat)

					F_ms           = F_obs - F_est
				}
		 	}	

			if(method == "probit"){	

				conditional_coef 		= apply(as.matrix(dep2_0), 1, function(xxx) coef(glm((dep0<=xxx)~reg0, weights=weight0, family=binomial(link="probit"))))
				conditional_coef        = conditional_coef[, colSums(!is.nan(conditional_coef))>0]
				pred 					= pnorm(cbind(1,reg0)%*%conditional_coef)
				predCounter 			= pnorm(cbind(1,counterfactual_var)%*%conditional_coef)

				conditional_fitted      = sort(apply(pred, 2, function(a) wtd.mean(a, weights=weight0)))
				if(max(conditional_fitted) ==1 ){
					conditional_fitted_s = conditional_fitted
				}else{
					conditional_fitted_s = c(conditional_fitted,1)
				}

				conditional_fittedCounter = sort(apply(predCounter, 2, function(a) wtd.mean(a, weights=weight1)))
				if(max(conditional_fittedCounter) == 1){
					predCounter_s       = conditional_fittedCounter
				}else{
					predCounter_s       = c(conditional_fittedCounter, 1)
				}

				marginal_fitted	        = getquantile(dep2_0, conditional_fitted_s, taus) 
				marginal_counterfactual = getquantile(dep2_0, predCounter_s, taus) 
				
				if(treatment){
					conditional_coef_1 	= apply(as.matrix(dep2_1), 1, function(xxx) coef(glm((dep1<=xxx)~reg1, weights=weight1, family=binomial(link="probit"))))
					conditional_coef_1  = conditional_coef_1[, colSums(!is.nan(conditional_coef_1))>0]
					pred_1	     		= pnorm(cbind(1,reg1)%*%conditional_coef_1) 
					
					conditional_fitted_1= sort(apply(pred_1, 2, function(a) wtd.mean(a, weights=weight1)))
					if(max(conditional_fitted_1) ==1 ){
						conditional_fitted_1_s = conditional_fitted_1
					}else{
						conditional_fitted_1_s = c(conditional_fitted_1,1)
					}
					marginal_fitted_1  	= getquantile(dep2_1, conditional_fitted_1_s, taus) 
				}	

				if(treatment&decomposition){

					y_thred        = c(depeval, depeval_1)
					thredIndex     = c(1, floor((1:(nreg-2)-0.5)/(nreg-2)*length(y_thred)), length(y_thred)) 
					y_thred_approx = y_thred[thredIndex]
					y_te           = c(dep0, dep1)
					F_obs          = colMeans(apply(as.matrix(y_thred_approx), 1, function(xxx) {ifelse( y_te <= xxx, 1, 0) }) ) 
 
			 		F_approx_fun0  = approxfun(dep2_0, conditional_fitted, yleft=0, yright=1)
					F_approx_fun1  = approxfun(dep2_1, conditional_fitted_1, yleft=0, yright=1)

			 		F_0hat         = length(dep0)/(length(dep0) + length(dep1))
					F_est          = F_approx_fun0(y_thred_approx)*F_0hat + F_approx_fun1(y_thred_approx)*(1-F_0hat)

					F_ms           = F_obs - F_est
				}

			}			
			if(method == "lpm"){	

				conditional_coef 		= apply(as.matrix(dep2_0), 1, function(xxx) coef(glm((dep0<=xxx)~reg0, weights=weight0)))
				conditional_coef        = conditional_coef[, colSums(!is.nan(conditional_coef))>0]
				pred 					= cbind(1, reg0)%*%conditional_coef
				predCounter 			= cbind(1, counterfactual_var)%*%conditional_coef

				conditional_fitted      = sort(apply(pred, 2, function(a) wtd.mean(a, weights=weight0)))
				if(max(conditional_fitted) ==1 ){
					conditional_fitted_s = conditional_fitted
				}else{
					conditional_fitted_s = c(conditional_fitted,1)
				}

				conditional_fittedCounter = sort(apply(predCounter, 2, function(a) wtd.mean(a, weights=weight1)))
				if(max(conditional_fittedCounter) == 1){
					predCounter_s       = conditional_fittedCounter
				}else{
					predCounter_s       = c(conditional_fittedCounter, 1)
				}

				marginal_fitted	        = getquantile(dep2_0, conditional_fitted_s, taus) 
				marginal_counterfactual = getquantile(dep2_0, predCounter_s, taus) 

				if(treatment){
					conditional_coef_1 	= apply(as.matrix(dep2_1), 1, function(xxx) coef(glm((dep1<=xxx)~reg1, weights=weight1)))
				    conditional_coef_1  = conditional_coef_1[, colSums(!is.nan(conditional_coef_1))>0]
					pred_1	     		= cbind(1,reg1)%*%conditional_coef_1
					
					conditional_fitted_1= sort(apply(pred_1, 2, function(a) wtd.mean(a, weights=weight1)))
					if(max(conditional_fitted_1) ==1 ){
						conditional_fitted_1_s = conditional_fitted_1
					}else{
						conditional_fitted_1_s = c(conditional_fitted_1,1)
					}
					marginal_fitted_1  	= getquantile(dep2_1, conditional_fitted_1_s, taus) 
				}
				if(treatment&decomposition){

					y_thred        = c(depeval, depeval_1)
					thredIndex     = c(1, floor((1:(nreg-2)-0.5)/(nreg-2)*length(y_thred)), length(y_thred)) 
					y_thred_approx = y_thred[thredIndex]
					y_te           = c(dep0, dep1)
					F_obs          = colMeans(apply(as.matrix(y_thred_approx), 1, function(xxx) {ifelse( y_te <= xxx, 1, 0) }) ) 
 
			 		F_approx_fun0  = approxfun(dep2_0, conditional_fitted, yleft=0, yright=1)
					F_approx_fun1  = approxfun(dep2_1, conditional_fitted_1, yleft=0, yright=1)

			 		F_0hat         = length(dep0)/(length(dep0) + length(dep1))
					F_est          = F_approx_fun0(y_thred_approx)*F_0hat + F_approx_fun1(y_thred_approx)*(1-F_0hat)

					F_ms           = F_obs - F_est
		 		}
	
			}
		}		
	}	
	
	if(treatment & decomposition){
		structral_effect            = marginal_fitted_1 - marginal_counterfactual
		quantile_effect 			= structral_effect 
		composition_effect          = marginal_counterfactual - marginal_fitted
		total_effect                = marginal_fitted_1 - marginal_fitted
	}else if(treatment & !decomposition){
		structral_effect            = marginal_fitted_1 - marginal_counterfactual
		quantile_effect             = structral_effect
	}else{
		composition_effect          = marginal_counterfactual - marginal_fitted
		quantile_effect 			= composition_effect         
	}

	marginal_obs    				= wtd.quantile(dep0, weights = weight0*length(weight0), probs=taus, na.rm=TRUE)
	names(marginal_fitted)  		= NULL
	names(marginal_counterfactual) 	= NULL
	names(marginal_obs) 			= NULL
	if(treatment){ marginal_obs_1 = wtd.quantile(dep1, weights = weight1*length(weight1), probs=taus, na.rm=TRUE) }
	
	if(treatment & decomposition){
    	
    	return(list(taus = taus, quantile_effect=quantile_effect, composition_effect = composition_effect, total_effect = total_effect, marginal_obs= marginal_obs,marginal_fitted=marginal_fitted,marginal_counterfactual=marginal_counterfactual,nreg=nreg, F_ms = F_ms, marginal_obs_1 = marginal_obs_1, marginal_fitted_1 = marginal_fitted_1))
	
	} else if(treatment & !decomposition){

    	return(list(taus = taus, quantile_effect=quantile_effect, marginal_obs= marginal_obs,marginal_fitted=marginal_fitted,marginal_counterfactual=marginal_counterfactual,nreg=nreg, F_ms = F_ms, marginal_obs_1 = marginal_obs_1, marginal_fitted_1 = marginal_fitted_1))
	
	} else {
		return(list(taus = taus, quantile_effect=quantile_effect, marginal_obs= marginal_obs,marginal_fitted=marginal_fitted,marginal_counterfactual=marginal_counterfactual,nreg=nreg))
	}
}

BootstrapProcedure <- function(transformation,reps,data0,depeval,data1,weightedboot,obs0,obs1, dep1, depeval_1, treatment, decomposition, taus,method,trimming,nreg,censoring,right,nsteps,firstc,secondc,seed,sepcore,ncore){
	subsamplek <- function(transformation,data0,data1,weightedboot,obs0,obs1,dep1, depeval_1, treatment, decomposition, taus,method,trimming,nreg,censoring,right,nsteps,firstc,secondc){				
		bootstrap_weight                  = NULL
		bootstrap_counterfactual_weight   = NULL
		if(weightedboot){
			bootstrap_weight                  = rexp(obs0,rate=1)
			bootstrap_counterfactual_weight   = rexp(obs1,rate=1)
		}			
		
		subsample_Index0           = sample(obs0,obs0,replace=TRUE,prob=bootstrap_weight)	
		subsampled_dep0            = (data0$dep)[subsample_Index0]
		subsampled_reg0            = as.matrix(data0$reg)[subsample_Index0,]
		subsampled_weight0         = (data0$weight)[subsample_Index0]
		subsampled_scale_variable  = as.matrix(data0$scale_variable)[subsample_Index0,]			
		if(transformation){					
			subsample_Index1 = subsample_Index0			
		}else{
			subsample_Index1 = sample(obs1,obs1,replace=TRUE,prob=bootstrap_counterfactual_weight)			
		}
		subsampled_counterfactual_var       = as.matrix(data1$counterfactual)[subsample_Index1,]
		subsampled_weight1                  = (data1$weight)[subsample_Index1]
		subsampled_counterfactual_scale_var = as.matrix(data1$scale_variable)[subsample_Index1,]	
		
		if(treatment){ subsampled_dep1 = dep1[subsample_Index1] }
		sub_PE_res = QteDistEst(subsampled_dep0,depeval,subsampled_reg0,subsampled_weight0,subsampled_counterfactual_var, subsampled_dep1, depeval_1, treatment, decomposition, subsampled_weight1,taus,method,trimming,nreg,subsampled_scale_variable,subsampled_counterfactual_scale_var,censoring,right,nsteps,firstc,secondc)	
		
		if(treatment & decomposition){
	    	qte_cov_boot = c(sub_PE_res$quantile_effect,sub_PE_res$marginal_obs,sub_PE_res$marginal_fitted,sub_PE_res$marginal_counterfactual, sub_PE_res$composition_effect, sub_PE_res$total_effect, sub_PE_res$marginal_obs_1, sub_PE_res$marginal_fitted_1, sub_PE_res$F_ms)
		} else if(treatment & !decomposition){
			qte_cov_boot = c(sub_PE_res$quantile_effect,sub_PE_res$marginal_obs,sub_PE_res$marginal_fitted,sub_PE_res$marginal_counterfactual, sub_PE_res$marginal_obs_1, sub_PE_res$marginal_fitted_1, sub_PE_res$F_ms)
	
		} else{
			qte_cov_boot = c(sub_PE_res$quantile_effect,sub_PE_res$marginal_obs,sub_PE_res$marginal_fitted,sub_PE_res$marginal_counterfactual)
		}
		return(qte_cov_boot)
	}	
	
	if(sepcore){
		if(ncore>=1){	
		    cores <- detectCores()
			if(ncore>=cores){
			    stop("The number of cores specified is bigger than or equal to the computer has.")
			}else{				
			    if(ncore>1) {cat(format(paste("cores used=",ncore,"\n"),width=50))}
			    cl <- makeCluster(ncore) 
				registerDoParallel(cl)			 
				qte_cov_booti <- foreach(i=1:reps, .combine = cbind, .packages =c("quantreg", "Hmisc", "survival"), .export = c("QteDistEst"), .options.RNG = seed) %dorng%{
				    tryCatch({subsamplek(transformation,data0,data1,weightedboot,obs0,obs1, dep1, depeval_1, treatment, decomposition, taus,method,trimming,nreg,censoring,right,nsteps,firstc,secondc)},
					  error = function(e) NULL 
					)
				}
				stopCluster(cl)
	        }		
		}else{	
			stop("The number of cores specified needs to be an integer.")
		}
	}else{
		cores <- detectCores()
		cl    <- makeCluster(cores-1) 
		registerDoParallel(cl)			 
		qte_cov_booti <- foreach(i=1:reps, .combine = cbind, .packages =c("quantreg", "Hmisc", "survival"), .export = c("QteDistEst"), .options.RNG = seed) %dorng%{
			tryCatch({subsamplek(transformation,data0,data1,weightedboot,obs0,obs1, dep1, depeval_1, treatment, decomposition, taus,method,trimming,nreg,censoring,right,nsteps,firstc,secondc)},
			error = function(e) NULL 
			)
		}
		stopCluster(cl)
	}			
		
	qte_cov_booti = t(qte_cov_booti)
	return(qte_cov_booti=qte_cov_booti)
}
	
InferenceTestingEval<- function(method,qte_cov_booti,reps,taus,robust,qte_cov_obs,last,first,alpha,cons_test, treatment, decomposition){
		
	nqs               = length(taus)
 	reps              = nrow(qte_cov_booti)
	sel0              = ((sum(taus<first))+1):(sum(taus<=last)) 
	
	## deal with median
	if(is.element(0.5,taus)){
		median_sel    = which(taus==0.5)	    
	}else{
		print("Q_0.5 is not evaluated. This is the median of the specified rather than the Q_0.5")
		median_sel    = ceiling(nqs/2) 			
	}
	sel_median                        = c((sum(taus<first)+1):(median_sel-1),(median_sel+1):sum(taus<=last)) 
	selall                            = 1:length(sel_median)
	
	quantile_effect_bootstrap	      = qte_cov_booti[,1:nqs]	
	marginal_obs_bootstrap	          = qte_cov_booti[,(nqs+1):(2*nqs)]	
	marginal_fitted_bootstrap         = qte_cov_booti[,(2*nqs+1):(3*nqs)]
	marginal_counterfactual_bootstrap = qte_cov_booti[,(3*nqs+1):(4*nqs)]		

    quantile_effect 	   	          = qte_cov_obs[1:nqs]
	marginal_obs				      = qte_cov_obs[(nqs+1):(2*nqs)]
	marginal_fitted				      = qte_cov_obs[(2*nqs+1):(3*nqs)]
	marginal_counterfactual           = qte_cov_obs[(3*nqs+1):(4*nqs)]	

	if(treatment & decomposition){
		composition_effect_bootstrap  = qte_cov_booti[,(4*nqs+1):(5*nqs)]
		total_effect_bootstrap	      = qte_cov_booti[,(5*nqs+1):(6*nqs)]
		marginal_obs_1_bootstrap      = qte_cov_booti[,(6*nqs+1):(7*nqs)]
		marginal_fitted_1_bootstrap   = qte_cov_booti[,(7*nqs+1):(8*nqs)]
		F_ms_bootstrap                = qte_cov_booti[,-(1:8*nqs)]
		
		composition_effect 	   	      = qte_cov_obs[(4*nqs+1):(5*nqs)]
		total_effect 	   	          = qte_cov_obs[(5*nqs+1):(6*nqs)]
		marginal_obs_1                = qte_cov_obs[(6*nqs+1):(7*nqs)]
		marginal_fitted_1             = qte_cov_obs[(7*nqs+1):(8*nqs)]
		F_ms_obs                 	  = qte_cov_obs[-(1:8*nqs)]
		sel_ms                        = 1:length(F_ms_obs)
	}else if(treatment & !decomposition){
		marginal_obs_1_bootstrap      = qte_cov_booti[,(4*nqs+1):(5*nqs)]
		marginal_fitted_1_bootstrap   = qte_cov_booti[,(5*nqs+1):(6*nqs)]
		F_ms_bootstrap                = qte_cov_booti[,-(1:6*nqs)]
		
		marginal_obs_1                = qte_cov_obs[(4*nqs+1):(5*nqs)]
		marginal_fitted_1             = qte_cov_obs[(5*nqs+1):(6*nqs)]
		F_ms_obs                 	  = qte_cov_obs[-(1:6*nqs)]
		sel_ms                        = 1:length(F_ms_obs)
	}
   
	#***************************************************** Variance ****************************************		
	## variance for quantile_effect
	pe_boot         = VarianceEval(quantile_effect_bootstrap,quantile_effect,reps,nqs,alpha,robust,sel0)		
	pe_lb_point 	= pe_boot$qte_cov - pe_boot$se*qnorm(1-alpha/2)
	pe_ub_point 	= pe_boot$qte_cov + pe_boot$se*qnorm(1-alpha/2)
	
	if(!treatment&!decomposition){
		
		resCE           = cbind(pe_boot$qte_cov, pe_boot$se, pe_lb_point, pe_ub_point, pe_boot$lb, pe_boot$ub)
		colnames(resCE) = c("composition.effect","se.ce","lb.ce.point","ub.ce.point","lb.ce","ub.ce")
	
	}else if(treatment&!decomposition){
		
		resSE           = cbind(pe_boot$qte_cov, pe_boot$se, pe_lb_point, pe_ub_point, pe_boot$lb, pe_boot$ub)
		colnames(resSE) = c("structure.effect","se.se","lb.se.point","ub.se.point","lb.se","ub.se")
	
	}else if(treatment & decomposition){
		
		## variance for structure_effect
		resSE           = cbind(pe_boot$qte_cov, pe_boot$se, pe_lb_point, pe_ub_point, pe_boot$lb, pe_boot$ub)
		colnames(resSE) = c("structure.effect","se.se","lb.se.point","ub.se.point","lb.se","ub.se")

		## variance for composition_effect
		ce_boot         = VarianceEval(composition_effect_bootstrap, composition_effect, reps, nqs, alpha, robust, sel0)		
		ce_lb_point 	= ce_boot$qte_cov - ce_boot$se*qnorm(1-alpha/2)
		ce_ub_point 	= ce_boot$qte_cov + ce_boot$se*qnorm(1-alpha/2)
		resCE			= cbind(ce_boot$qte_cov, ce_boot$se, ce_lb_point, ce_ub_point, ce_boot$lb, ce_boot$ub)		
		colnames(resCE) = c("composition.effect","se.ce","lb.ce.point","ub.ce.point","lb.ce","ub.ce")

		## variance for total_effect
		te_boot         = VarianceEval(total_effect_bootstrap, total_effect, reps, nqs, alpha, robust, sel0)		
		te_lb_point 	= te_boot$qte_cov - te_boot$se*qnorm(1-alpha/2)
		te_ub_point 	= te_boot$qte_cov + te_boot$se*qnorm(1-alpha/2)
		resTE			= cbind(te_boot$qte_cov, te_boot$se, te_lb_point, te_ub_point, te_boot$lb, te_boot$ub)		
		colnames(resTE) = c("total.effect", "se.te", "lb.te.point", "ub.te.point", "lb.te", "ub.te")
	}
	
	## Distribution	
	sample_quantile_ref0             = VarianceEval(marginal_obs_bootstrap, marginal_obs,reps,nqs,alpha,robust,sel0) 
	colnames(sample_quantile_ref0)   = c("ME.ref0","se.ref0","lb.ref0","ub.ref0")

	model_quantile_ref0              = VarianceEval(marginal_fitted_bootstrap, marginal_fitted,reps,nqs,alpha,robust,sel0)	
	colnames(model_quantile_ref0)    = c("ME.fitted0","se.fitted0","lb.fitted0","ub.fitted0")

	model_quantile_counter           = VarianceEval(marginal_counterfactual_bootstrap, marginal_counterfactual,reps,nqs,alpha,robust,sel0)	
	colnames(model_quantile_counter) = c("ME.counter","se.counter","lb.counter","ub.counter")

	if(treatment){
		sample_quantile_ref1   = VarianceEval(marginal_obs_1_bootstrap, marginal_obs_1,reps,nqs,alpha,robust,sel0) 
	    model_quantile_ref1    = VarianceEval(marginal_fitted_1_bootstrap,marginal_fitted_1,reps,nqs,alpha,robust,sel0)	

	    colnames(sample_quantile_ref1)   = c("ME.ref1","se.ref1","lb.ref1","ub.ref1")
	    colnames(model_quantile_ref1)    = c("ME.fitted1","se.fitted1","lb.fitted1","ub.fitted1")
	}
		
	#************************************************** Testing  *******************************************
	## test of no misspecification, obs-fitted		
	if(method=="logit"| method=="lpm"){	
		test_MS = c(NA,NA)
 	}else{
		qte_cov_boot_ms = marginal_obs_bootstrap - marginal_fitted_bootstrap	
	    qte_cov_ms      = marginal_obs - marginal_fitted
		test_MS         = TestingEval((qte_cov_boot_ms-kronecker(matrix(qte_cov_ms,1,nqs),matrix(1,reps,1)))[,sel0],qte_cov_ms[sel0],qte_cov_boot_ms,sel0,reps,robust)
	}
	# test of no effect, test_const == 0	
	test_0 = TestingEval((quantile_effect_bootstrap-kronecker(matrix(quantile_effect,1,nqs),matrix(1,reps,1)))[,sel0],quantile_effect[sel0],quantile_effect_bootstrap,sel0,reps,robust)	
	# test of const effect	
	constest1      = sort(unique(cons_test))	
	constestNo0    = setdiff(constest1,0) 
	nc             = length(constestNo0)	
	if(nc>0){
	    test_const     = matrix(0,nc,2)
		for(i in 1:nc){						
			test_const[i,] = TestingEval((quantile_effect_bootstrap-kronecker(matrix(quantile_effect,1,nqs),matrix(1,reps,1)))[,sel0],quantile_effect[sel0]-constestNo0[i],quantile_effect_bootstrap,sel0,reps,robust)	
		}	
	}else test_const = NULL
	# test of median		
	qte_cov_boot_median = quantile_effect_bootstrap[,sel_median]-quantile_effect_bootstrap[,median_sel] 
	qte_cov_def_median  = quantile_effect[sel_median]-quantile_effect[median_sel]	
	test_median         = TestingEval(qte_cov_boot_median-kronecker(matrix(qte_cov_def_median,1,ncol(qte_cov_boot_median)),matrix(1,reps,1)),qte_cov_def_median,qte_cov_boot_median,selall,reps,robust)	
	# test of stochastic dominance	
	qte_cov_boot_SD_ex = (quantile_effect_bootstrap-kronecker(matrix(quantile_effect,1,nqs),matrix(1,reps,1)))[,sel0]   	
	test_SD            = TestingEval((qte_cov_boot_SD_ex*(qte_cov_boot_SD_ex<=0)),(quantile_effect*(quantile_effect<=0))[sel0],quantile_effect_bootstrap,sel0,reps,robust)
	# test of being stochastically dominated		
	test_SDD  = TestingEval((qte_cov_boot_SD_ex*(qte_cov_boot_SD_ex>=0)),(quantile_effect*(quantile_effect>=0))[sel0],quantile_effect_bootstrap,sel0,reps,robust)
	

	## Tests
	if(!treatment&!decomposition){ 
		testCE  = rbind(test_MS, test_0, test_const, test_median, test_SD, test_SDD)
	}else if(treatment&!decomposition){ 
		testSE  = rbind(test_MS, test_0, test_const, test_median, test_SD, test_SDD)
	}else if(treatment&decomposition){
		
		testSE  = rbind(test_MS, test_0, test_const, test_median, test_SD, test_SDD)

		## **** CE *****
		## test of no misspecification, obs-fitted		
		if(method=="logit"|method=="lpm"){
			test_MS_CE = c(NA,NA)
	 	}else{
	 		qte_cov_boot_ms_CE = marginal_obs_1_bootstrap - marginal_fitted_1_bootstrap
		    qte_cov_ms_CE      = marginal_obs_1 - marginal_fitted_1
		    test_MS_CE         = TestingEval((qte_cov_boot_ms_CE-kronecker(matrix(qte_cov_ms_CE,1,length(qte_cov_ms_CE)),matrix(1,reps,1)))[,sel0], qte_cov_ms_CE[sel0], qte_cov_boot_ms_CE, sel0, reps, robust)
		}
		# # test of no effect, test_const == 0
		test_0_CE = TestingEval((composition_effect_bootstrap-kronecker(matrix(composition_effect,1,nqs),matrix(1,reps,1)))[,sel0],composition_effect[sel0],composition_effect_bootstrap,sel0,reps,robust)	
		# test of const effect	
		if(nc>0){
		    test_const_CE     = matrix(0,nc,2)
			for(i in 1:nc){						
				test_const_CE[i,] = TestingEval((composition_effect_bootstrap-kronecker(matrix(composition_effect,1,nqs),matrix(1,reps,1)))[,sel0],composition_effect[sel0]-constestNo0[i],composition_effect_bootstrap,sel0,reps,robust)	
			}	
		}else test_const_CE = NULL
		# test of median		
		ce_cov_boot_median = composition_effect_bootstrap[,sel_median]-composition_effect_bootstrap[,median_sel] 
		ce_cov_def_median  = composition_effect[sel_median]-composition_effect[median_sel]	
		test_median_CE     = TestingEval(ce_cov_boot_median-kronecker(matrix(ce_cov_def_median,1,ncol(ce_cov_boot_median)),matrix(1,reps,1)),ce_cov_def_median,ce_cov_boot_median,selall,reps,robust)	
		# test of stochastic dominance	
		ce_cov_boot_SD_ex = (composition_effect_bootstrap-kronecker(matrix(composition_effect,1,nqs),matrix(1,reps,1)))[,sel0]   	
		test_SD_CE        = TestingEval((ce_cov_boot_SD_ex*(ce_cov_boot_SD_ex<=0)),(composition_effect*(composition_effect<=0))[sel0],composition_effect_bootstrap,sel0,reps,robust)
		# test of being stochastically dominated		
		test_SDD_CE  = TestingEval((ce_cov_boot_SD_ex*(ce_cov_boot_SD_ex>=0)),(composition_effect*(composition_effect>=0))[sel0],composition_effect_bootstrap,sel0,reps,robust)
		# testCE
		testCE  = rbind(test_MS_CE, test_0_CE,test_const_CE,test_median_CE,test_SD_CE,test_SDD_CE)

		## **** TE *****
		## test of no misspecification, obs-fitted		
		if(method=="logit"|method=="lpm"){
			test_MS_TE = c(NA,NA)
	 	}else{
			qte_cov_boot_ms_TE = F_ms_bootstrap 
		    qte_cov_ms_TE      = F_ms_obs
		    test_MS_TE         = TestingEval((qte_cov_boot_ms_TE-kronecker(matrix(qte_cov_ms_TE,1,length(qte_cov_ms_TE)),matrix(1,reps,1))), qte_cov_ms_TE, qte_cov_boot_ms_TE, sel_ms,reps, robust)
		}
		# test of no effect, test_const == 0
		test_0_TE = TestingEval((total_effect_bootstrap-kronecker(matrix(total_effect,1,nqs),matrix(1,reps,1)))[,sel0],total_effect[sel0],total_effect_bootstrap,sel0,reps,robust)	
		# test of const effect				
		if(nc>0){
		    test_const_TE     = matrix(0,nc,2)
			for(i in 1:nc){						
				test_const_TE[i,] = TestingEval((total_effect_bootstrap-kronecker(matrix(total_effect,1,nqs),matrix(1,reps,1)))[,sel0],total_effect[sel0]-constestNo0[i],total_effect_bootstrap,sel0,reps,robust)	
			}	
		}else test_const_TE = NULL
		# test of median		
		te_cov_boot_median = total_effect_bootstrap[,sel_median]-total_effect_bootstrap[,median_sel] 
		te_cov_def_median  = total_effect[sel_median]-total_effect[median_sel]	
		test_median_TE     = TestingEval(te_cov_boot_median-kronecker(matrix(te_cov_def_median,1,ncol(te_cov_boot_median)),matrix(1,reps,1)),te_cov_def_median,te_cov_boot_median,selall,reps,robust)	
		# test of stochastic dominance	
		te_cov_boot_SD_ex  = (total_effect_bootstrap-kronecker(matrix(total_effect,1,nqs),matrix(1,reps,1)))[,sel0]   	
		test_SD_TE         = TestingEval((te_cov_boot_SD_ex*(te_cov_boot_SD_ex<=0)),(total_effect*(total_effect<=0))[sel0],total_effect_bootstrap,sel0,reps,robust)
		# test of being stochastically dominated		
		test_SDD_TE  = TestingEval((te_cov_boot_SD_ex*(te_cov_boot_SD_ex>=0)),(total_effect*(total_effect>=0))[sel0],total_effect_bootstrap,sel0,reps,robust)
		#testTE
		testTE  = rbind(test_MS_TE,test_0_TE,test_const_TE,test_median_TE,test_SD_TE,test_SDD_TE)
	}	

	res_boot 	                    = NULL	 
	res_boot$sample_quantile_ref0   = sample_quantile_ref0 	
	res_boot$model_quantile_ref0    = model_quantile_ref0  
	res_boot$model_quantile_counter = model_quantile_counter 

	if(!treatment&!decomposition){
		res_boot$resCE 	  = resCE
    	res_boot$testCE   = testCE 	

	}else if(treatment&!decomposition){
		res_boot$sample_quantile_ref1  = sample_quantile_ref1
	    res_boot$model_quantile_ref1   = model_quantile_ref1

		res_boot$resSE 	               = resSE	
		res_boot$testSE                = testSE

	}else if(treatment & decomposition){
		res_boot$sample_quantile_ref1  = sample_quantile_ref1
	    res_boot$model_quantile_ref1   = model_quantile_ref1

		res_boot$resSE 	               = resSE	
		res_boot$testSE                = testSE

    	res_boot$resCE 	               = resCE
    	res_boot$testCE                = testCE	

    	res_boot$resTE                 = resTE 
		res_boot$testTE                = testTE
	}
	
	return(res_boot)
}	

VarianceEval <- function(qte_cov_boot,qte_cov,reps,nqs,alpha,robust,sel){
	if(robust){
		seuqf = apply(qte_cov_boot,2,function(x) {(quantile(x,.75,na.rm=TRUE)-quantile(x,.25,na.rm=TRUE))/1.34} )
	}else{
		seuqf = sqrt(diag(var(qte_cov_boot))) #nqs, hence this is the variance for each tau
	}
	Vuqf    = seuqf^2

	##
	Kuqf    = sqrt((qte_cov_boot-kronecker(matrix(qte_cov,1,nqs),matrix(1,reps,1)))^2/kronecker(matrix(Vuqf,1,nqs),matrix(1,reps,1)))[,sel]
	Kuqfsel	= Kuqf[,apply(Kuqf, 2, function(x) all(is.finite(x)))]
	Kmaxuqf = apply(Kuqfsel,1,max)  
	##

	Kalpha	    = quantile(Kmaxuqf,1-alpha,na.rm=TRUE,names=FALSE)			
	lb		    = qte_cov - seuqf*Kalpha
	ub		    = qte_cov + seuqf*Kalpha	
	simbootres  = data.frame(qte_cov=qte_cov,se=seuqf,lb=lb,ub=ub)	
	
	return(simbootres)
}	

TestingEval<- function(boot_test_numerator,obs_test_numerator,variable_for_variance,sel,reps,robust){			
	if(robust){
		seboot	= apply(variable_for_variance,2,function(x) {(quantile(x,.75,na.rm=TRUE)-quantile(x,.25,na.rm=TRUE))/1.34} )		
		Vuqf    = (seboot)^2
	}else{		
		Vuqf	= diag(var(variable_for_variance))
	}		
	Vuqf        = Vuqf[sel]+1e-9
	nqs         = ncol(boot_test_numerator)
	Kuqf        = sqrt(boot_test_numerator^2/kronecker(matrix(Vuqf,1,nqs),matrix(1,reps,1)))
	
	Kmaxuqf     = apply(Kuqf,1,max)
	KSstat      = max(sqrt(obs_test_numerator^2/Vuqf))
 
	Kuqf2       = Kuqf^2
	Kmeanuqf    = rowMeans(Kuqf2)		
	CMSstat     = mean(obs_test_numerator^2/Vuqf)

	testboot    = c(mean(Kmaxuqf>KSstat),mean(Kmeanuqf>CMSstat))
	return(testboot)
}