##############################################################################################
# OUTPUTS :
#
# Matrix matstockgam: Matrix of dimension nT x N where T is the number of the different values of (alpha_n,tau_j) considered. This matrix collect the values of the estimations of gamma for each replications (one replication corresponds to one column). The rows (j-1)n+1 to jn corresponds to the estimations of gamma for each X_i, i=1,...,n and for the j-th value of (alpha_n,tau_j)
#
# Matrix matstock: similar to matstockgam but for the extreme quatile estimator
###############################################################################################


#######################################
# SETTINGS                            #
#######################################

# Definition of the Kernel function K #
funcK<-function(t) {
eps<-0.1
((2-eps)+2*(eps-1)*t)*(t <= 1 & t>= 0)
}

# Definition of the kernel function Q #
funcQ<-function(t) {
1/2*(1+2*t-sign(t)*t^2)*(t<1 & t>=-1)+1*(t>=1)
}

# Choice of the distance #
funcdopt<-function(s,t) {
# Distance d_Z #
1/2*abs(sin(4*pi*s)/(4*pi*s)-sin(4*pi*t)/(4*pi*t)) 
# Distance d_X #
#1+1/2*(sin(4*pi*s)/(4*pi*s)+sin(4*pi*t)/(4*pi*t))-sin(2*pi*(s+t))/(2*pi*(s+t))-(sin(2*pi*(s-t))+1*(s==t))/(2*pi*(s-t)+1*(s==t))
}

# Sample size #
n<-500

# Number of replications #
N<-50

# Definition of the function gamma #
funcgam<-function(A) {
2*A-1+0.25
}

# Order of the quantile to estimate #
beta<-5/n

# Values of alpha_n and tau_j # 
rvec<-c(9) # Power of 1/j #
alphamaxvec<-c(5,10,15,20)*log(n)/n
matvar<-matrix(c(rep(rvec,length(alphamaxvec)),sort(rep(alphamaxvec,length(rvec)))),nrow=2,byrow=T)

################################################################################################

##################################
# Simulation of the observations #
##################################

Z<-matrix(runif(N*n,1/4,1),ncol=N)
normeX<-1/2*(1+sin(4*pi*Z)/(4*pi*Z))
gamma<-funcgam(normeX)
rho<-(-2)
U<-matrix(runif(N*n,0,1),ncol=N)
Y<-(U^(rho*gamma)-1)^(-1/rho) # Loi de Burr

###############################################################################################

###################################################
# Selection of h and lambda by cross-validation   #
###################################################

hvec<-seq(0.01,0.1,length=20) 
lambdavec<-c(0.1)
crityao<-function(Y,Z,hvec,lambdavec) {
N<-ncol(Z)
n<-nrow(Z)
mathl<-matrix(c(rep(hvec,length(lambdavec)),sort(rep(lambdavec,length(hvec)))),nrow=2,byrow=T)
ssech<-c(1:n)
crit<-matrix(rep(NA,length(ssech)*N),ncol=N)
critfin<-matrix(rep(NA,ncol(mathl)*N),ncol=N)

for (r in 1:ncol(mathl)) {

h<-mathl[1,r]
lambda<-mathl[2,r]

	for (j in 1:length(ssech)) {
	y<-Y[ssech[j],] # y matrice de dimension 1 x N #
	res<-matrix(rep(NA,n*N),ncol=N)

		for (i in 1:N) {
		Zi<-Z[,i]
		Yi<-Y[,i]
		yi<-y[i]
		funcd<-function(u) { funcdopt(u,Zi) }
		matd<-sapply(Zi,funcd)
		matK<-funcK(matd/h)
		diag(matK)<-rep(0,n)
		matKQ<-matK*matrix(rep(funcQ((Yi-yi)/lambda),n),ncol=n)
		res[,i]<-(1*(Yi>=yi)-apply(matKQ,2,sum)/apply(matK,2,sum))^2
		res[apply(matK,2,sum)==0,i]<-Inf
		}
	crit[j,]<-apply(res,2,sum)
	}
critfin[r,]<-apply(crit,2,sum)
}
critfin
}

res4<-crityao(Y,Z,hvec,lambdavec)

R1<-matrix(rep(apply(res4,2,min),(length(hvec)*length(lambdavec))),ncol=N,byrow=T)
ind<-apply((R1==res4)*matrix(rep(c(1:(length(hvec)*length(lambdavec))),N),ncol=N),2,sum)
mathl<-matrix(c(rep(hvec,length(lambdavec)),sort(rep(lambdavec,length(hvec)))),nrow=2,byrow=T)
hlambdaopt<-matrix(mathl[,ind],ncol=N) 

##############################################################################################

##########################################################
# Computation of the estimators of gamma and q(beta_n|x) #
##########################################################

# Estimation of F #
funcFchap<-function(y,Y,Z,hl) {
funcdiff<-function(s,t) {
s-t
}
N<-ncol(Z)
n<-nrow(Z)
res<-matrix(rep(NA,n*N),ncol=N)

	for (i in 1:N) {

	h<-hl[1,i]
	lambda<-hl[2,i]

	Zi<-Z[,i]
	Yi<-Y[,i]
	yi<-y[,i]
	funcd<-function(u) { funcdopt(u,Zi) }
	matd<-sapply(Zi,funcd)
	matK<-funcK(matd/h)
	funcdiff2<-function(u) { funcdiff(u,yi) }
	matdiff<-sapply(Yi,funcdiff2)
	matQ<-funcQ(matdiff/lambda)
		for (j in 1:n) {
		res[j,i]<-sum(matK[j,]*matQ[j,])/sum(matK[j,])
		}
	}
res
}
#################################################################################

# Estimation of the quantile #
funcqchap<-function(alpha,Y,Z,hl){

tol<-10^(-10) 

A<-2*max(Y)
n<-nrow(Y)
N<-ncol(Y)
Amin<-matrix(rep(0,n*N),ncol=N)
Amax<-matrix(rep(A,n*N),ncol=N)
ecart<-1
niter<-0

	while (ecart>tol) {
	TMP<-funcFchap((Amin+Amax)/2,Y,Z,hl)
	ecart<-max(Amax-Amin)
	matsigne<-(sign(TMP-alpha)+1)/2
	Amaxtmp<-Amax*matsigne+(1-matsigne)*(Amin+Amax)/2
	Amintmp<-matsigne*(Amin+Amax)/2+(1-matsigne)*Amin
	Amax<-Amaxtmp
	Amin<-Amintmp
	}
(Amax+Amin)/2
}

################################################################################

# Estimation of gamma and the extreme quantile #
funcweiss<-function(r,alphamax) {

nbpt<-matrix(rep(NA,n*N),ncol=N)
	for (i in 1:N) {
	Zi<-Z[,i]
	funcd<-function(u) { funcdopt(u,Zi) }
	matd<-sapply(Zi,funcd)
	nbpt[,i]<-apply(matd<=hlambdaopt[1,i],2,sum)
	}
J<-10
alphamin<-1/min(nbpt) 
matqchap<-matrix(rep(NA,n*N*J),nrow=(n*N)) 

# Definition of tau_j #

tauj<-(seq(alphamax,alphamin,length=J)/alphamax)^r

for (j in 1:J) {
matqchap[,j]<-matrix(funcqchap(matrix(rep(alphamax*tauj[j],(n*N)),ncol=N),Y,Z,hlambdaopt),ncol=1)
}

# Estimation of Gamma #

p<-1 
matchapun<-matrix(rep(log(matqchap)[,1],J),ncol=J)
matgamtmp<-(apply(((log(matqchap)-matchapun)^p),1,sum)/sum((log(1/tauj))^p))^(1/p)
matgam<-matrix(matgamtmp,ncol=N)


# Estimation of the extreme quantile #

qweiss<-matrix(matqchap[,1],ncol=N)*exp(log(alphamax/beta)*matgam)

# Computation of the L2-errors delta #

R3<-(beta^(rho*gamma)-1)^(-1/rho) 
classMSE2<-sort(apply((R3-qweiss)^2,2,sum),index.return=T)
c10<-classMSE2$ix[max(1,floor(10*N/100))] # 10 % #
c50<-classMSE2$ix[floor(50*N/100)] # Mediane #
c90<-classMSE2$ix[floor(90*N/100)] # 90 % #
Delta10<-classMSE2$x[c10]
Delta50<-classMSE2$x[c50]
Delta90<-classMSE2$x[c90]
C<-c(c10,c50,c90)
Delta<-c(Delta10,Delta50,Delta90)

return(list(qweiss=qweiss,C=C,matgam=matgam,Delta=Delta))

}

###########################################################################

matstock<-matrix(rep(NA,n*N*length(alphamaxvec)*length(rvec)),ncol=N)
matstockgam<-matrix(rep(NA,n*N*length(alphamaxvec)*length(rvec)),ncol=N)
matC<-matrix(rep(NA,length(alphamaxvec)*length(rvec)*3),ncol=3)
matDelta<-matrix(rep(NA,length(alphamaxvec)*length(rvec)*3),ncol=3)

for (i in 1:(length(alphamaxvec)*length(rvec))) {
cat("compteur",i,"\n")
r<-matvar[1,i]
alphamax<-matvar[2,i]
res<-funcweiss(r,alphamax) 
matstock[((i-1)*n+1):(i*n),]<-res$qweiss
matstockgam[((i-1)*n+1):(i*n),]<-res$matgam
matC[i,]<-res$C
matDelta[i,]<-res$Delta
}