rm(list=ls())
library(foreach)
library(doMC)
library(np)
library(quantreg)

# Data generating process

gendata = function(n) {
#parameters
mux1<-1
mux2<-5
d0<--1
d1<-5
d2<--5
d3<--0.05
g1<--5
g2<-1
beta<-0 #parameter controlling TE
sige0<-5 #parameter controlling TE
sige1<-5 #parameter controlling TE


# characteristics
x1<-runif(n,mux1-(sqrt(12)/2),mux1+(sqrt(12)/2))
x2<-runif(n,mux2-(sqrt(12)/2),mux2+(sqrt(12)/2))
x1s<-x1^2
x2s<-x2^2
x1x2<-x1*x2

X<-cbind(x1,x2,x1s,x2s,x1x2)
eta<-rlogis(n, location = 0, scale = 10)
  
T<-ifelse(((d0+d1*x1+d2*x2+d3*x1s+eta)>0),1,0)

eps0<-rnorm(n,0,sige0)
eps1<-rnorm(n,beta,sige1)

Y0<-g1*x1+g2*x2+eps0
Y1<-Y0+eps1-eps0
Y<-Y1*T+Y0*(1-T)
  
  simdata = data.frame(Y, X, T)
  
  colnames(simdata) <- c("y","x1","x2","x1s","x2s","x1x2","t")
  
  simdata
} 

MCfun = function(dataset) {
    attach(dataset)
    n<-length(y)
    lvl = seq(0.1, 0.9, 0.02)

    pscore<-glm(t~x1+x2+x1s+x2s+x1x2,family=binomial(link = "logit"))$fitted
    
    w1<-t/(n*pscore)
    w0<-(1-t)/(n*(1-pscore))

    dr.q1 = foreach(i = 1:length(lvl), .combine = c) %do% {
        coef(rq(y ~ 1, tau=lvl[i], weights = w1))
    }
    
    dr.q0 = foreach(i = 1:length(lvl), .combine = c) %do% {
        coef(rq(y ~ 1, tau=lvl[i], weights = w0))
    }
    detach(dataset)
    
    c(dr.q1-dr.q0)
} 


ssfun <- function(dataset, block, smplsz) {
  resamples <- lapply(1:(smplsz-block +1), function(i) dataset[(i:(smplsz-block-1+i)),])
  lapply(resamples, MCfun)
}


#registerDoMC(12)

smplsz = 500

rep = 2000

block<- 350

NumSS<-smplsz-block+1

Numtau<-length(seq(0.1, 0.9, 0.02))

# The following calculate a lot of bootstraps
dup.beta = foreach(j=1:rep)%dopar%{
  dataset = gendata(smplsz)
  beta = MCfun(dataset)
  btbeta = ssfun(dataset, block, smplsz)
  list(beta=beta, btbeta=btbeta)
}


beta.aux<-array(0,dim=c(rep,Numtau))
btbeta.aux<-array(0,dim=c(rep,NumSS,Numtau))
for (j in 1:rep){ 
for (i in 1:NumSS){
for (k in 1:Numtau){
  beta.aux[j,k]<-dup.beta[[j]]$beta[[k]]
  btbeta.aux[j,i,k]<-dup.beta[[j]]$btbeta[[i]][k]
}
}
}

Wb<-array(0,dim=c(rep,NumSS,Numtau))
for (i in 1:rep){
for (j in 1:NumSS){
Wb[i,j,]<-btbeta.aux[i,j,]#-beta.aux[i,]

}
}

Vb<-array(0,dim=c(rep,NumSS))
for (i in 1:rep){
for (j in 1:NumSS){
#Vb[i,j]<-max(abs(Wb[i,j,])) #KS - No SQRT
Vb[i,j]<-sum(abs(Wb[i,j,])) #CVM - No SQRT
}
}

Vf<-array(0,dim=c(rep,1))
cv.01<-array(0,dim=c(rep,1))
cv.05<-array(0,dim=c(rep,1))
cv.10<-array(0,dim=c(rep,1))

for (i in 1:rep){
#Vf[i]<-max(abs(beta.aux[i,])) #KS - No SQRT
Vf[i]<-sum(abs(beta.aux[i,])) #CVM - No SQRT
cv.01[i]<-quantile(Vb[i,],0.99)
cv.05[i]<-quantile(Vb[i,],0.95)
cv.10[i]<-quantile(Vb[i,],0.90)
}


result.aux<-cbind(Vf,cv.01,cv.05,cv.10)


# Size of the tests
rr.01<-array(0,dim=c(rep,1))
rr.05<-array(0,dim=c(rep,1))
rr.10<-array(0,dim=c(rep,1))

for (i in 1:rep){
rr.01[i]<-ifelse(Vf[i]>cv.01[i],1,0)
rr.05[i]<-ifelse(Vf[i]>cv.05[i],1,0)
rr.10[i]<-ifelse(Vf[i]>cv.10[i],1,0)
}

size.01<-sum(rr.01)/rep
size.05<-sum(rr.05)/rep
size.10<-sum(rr.10)/rep

result<-cbind(size.01,size.05,size.10)
result

