#These are the codes that were used to generate the tables 1 to 5 (Monte Carlo simulations)

library(MASS)

#parametric estimator
obs<-100
resultsp1<-matrix(,10000,21)
for(i in 1:10000){
set.seed(i)
z0<-mvrnorm(n=obs,mu=c(0,0,0),Sigma=matrix(c(1,0.5,0.5,0.5,1,0.5,0.5,0.5,1),3,3))
z0[,1]<-pnorm(z0[,1])
z0[,2]<-(z0[,2]>0)
z0<-cbind(rep(1,obs)+z0%*%c(1,1,1)+rcauchy(obs)*(1+z0[,1]),z0)
z1<-mvrnorm(n=obs,mu=c(0,0,0),Sigma=matrix(c(1,0.5,0.5,0.5,1,0.5,0.5,0.5,1),3,3))
z1[,1]<-pnorm(z1[,1])
z1[,2]<-(z1[,2]>0)
z1<-cbind(rep(1,obs)+z1%*%c(1,1,1)+rnorm(obs)*(1+z1[,1]),z1)
if(i>4000){
temp<-matrix(,0,3)
while(dim(temp)[1]<100){
boot<-sample(1:obs,obs,replace=TRUE)
z0boot<-z0[boot,]
z1boot<-z1[boot,]
temp1<-rqp(z0boot[,1],z0boot[,2:4])
temp2<-evalu(z1boot[,2:4],temp1$beta,temp1$wq,prob=c(0.05,0.25,0.5))
temp1<-rqp(z1boot[,1],z1boot[,2:4])
temp3<-evalu(z1boot[,2:4],temp1$beta,temp1$wq,prob=c(0.05,0.25,0.5))
temp<-rbind(temp,temp3-temp2)}
} else temp<-matrix(0,2,3)
temp1<-rqdeco(z0[,1],z1[,1],z0[,2:4],z1[,2:4],prob=c(0.05,0.25,0.5),est.cov=TRUE)
resultsp1[i,]<-c(quantile(z0[,1],c(0.05,0.25,0.5)),quantile(z1[,1],c(0.05,0.25,0.5)),temp1$res[,1],temp1$res[,2],temp1$deco[,3],temp1$deco[,4],sd(temp))
}

obs<-400
resultsp2<-matrix(,5000,22)
for(i in 1:5000){
set.seed(i)
z0<-mvrnorm(n=obs,mu=c(0,0,0),Sigma=matrix(c(1,0.5,0.5,0.5,1,0.5,0.5,0.5,1),3,3))
z0[,1]<-pnorm(z0[,1])
z0[,2]<-(z0[,2]>0)
z0<-cbind(rep(1,obs)+z0%*%c(1,1,1)+rcauchy(obs)*(1+z0[,1]),z0)
z1<-mvrnorm(n=obs,mu=c(0,0,0),Sigma=matrix(c(1,0.5,0.5,0.5,1,0.5,0.5,0.5,1),3,3))
z1[,1]<-pnorm(z1[,1])
z1[,2]<-(z1[,2]>0)
z1<-cbind(rep(1,obs)+z1%*%c(1,1,1)+rnorm(obs)*(1+z1[,1]),z1)
if(i<=2000){
temp<-matrix(,0,3)
while(dim(temp)[1]<100){
boot<-sample(1:obs,obs,replace=TRUE)
z0boot<-z0[boot,]
z1boot<-z1[boot,]
temp1<-rqp(z0boot[,1],z0boot[,2:4])
temp2<-evalu(z1boot[,2:4],temp1$beta,temp1$wq,prob=c(0.05,0.25,0.5))
temp1<-rqp(z1boot[,1],z1boot[,2:4])
temp3<-evalu(z1boot[,2:4],temp1$beta,temp1$wq,prob=c(0.05,0.25,0.5))
temp<-rbind(temp,temp3-temp2)}
} else temp<-matrix(0,2,3)
temp1<-rqdeco(z0[,1],z1[,1],z0[,2:4],z1[,2:4],prob=c(0.05,0.25,0.5),est.cov=TRUE)
resultsp2[i,]<-c(quantile(z0[,1],c(0.05,0.25,0.5)),quantile(z1[,1],c(0.05,0.25,0.5)),temp1$res[,1],temp1$res[,2],temp1$deco[,3],temp1$deco[,4],sd(temp),temp1$res[3,3])
}

obs<-1600
resultsp3<-matrix(,2500,21)
for(i in 1:2500){
set.seed(i)
z0<-mvrnorm(n=obs,mu=c(0,0,0),Sigma=matrix(c(1,0.5,0.5,0.5,1,0.5,0.5,0.5,1),3,3))
z0[,1]<-pnorm(z0[,1])
z0[,2]<-(z0[,2]>0)
z0<-cbind(rep(1,obs)+z0%*%c(1,1,1)+rcauchy(obs)*(1+z0[,1]),z0)
z1<-mvrnorm(n=obs,mu=c(0,0,0),Sigma=matrix(c(1,0.5,0.5,0.5,1,0.5,0.5,0.5,1),3,3))
z1[,1]<-pnorm(z1[,1])
z1[,2]<-(z1[,2]>0)
z1<-cbind(rep(1,obs)+z1%*%c(1,1,1)+rnorm(obs)*(1+z1[,1]),z1)
if(i<=1000){
temp<-matrix(,0,3)
while(dim(temp)[1]<100){
boot<-sample(1:obs,obs,replace=TRUE)
z0boot<-z0[boot,]
z1boot<-z1[boot,]
temp1<-rqpl(z0boot[,1],z0boot[,2:4],nrq=200)
temp2<-evalu(z1boot[,2:4],temp1$beta,temp1$wq,prob=c(0.05,0.25,0.5))
temp1<-rqpl(z1boot[,1],z1boot[,2:4],nrq=200)
temp3<-evalu(z1boot[,2:4],temp1$beta,temp1$wq,prob=c(0.05,0.25,0.5))
temp<-rbind(temp,temp3-temp2)}
} else temp<-matrix(0,2,3)
temp1<-rqdeco(z0[,1],z1[,1],z0[,2:4],z1[,2:4],prob=c(0.05,0.25,0.5),est.cov=TRUE,nrq=200,method2=2)
resultsp3[i,]<-c(quantile(z0[,1],c(0.05,0.25,0.5)),quantile(z1[,1],c(0.05,0.25,0.5)),temp1$res[,1],temp1$res[,2],temp1$deco[,3],temp1$deco[,4],sd(temp))
}

#machado and mata
obs<-400
resultsm<-matrix(,5000,3*505)
m<-c(seq(100,50000,100),60000,70000,80000,90000,100000)
for(i in 1:5000){
set.seed(i)
z0<-mvrnorm(n=obs,mu=c(0,0,0),Sigma=matrix(c(1,0.5,0.5,0.5,1,0.5,0.5,0.5,1),3,3))
z0[,1]<-pnorm(z0[,1])
z0[,2]<-(z0[,2]>0)
z0<-cbind(rep(1,obs)+z0%*%c(1,1,1)+rcauchy(obs)*(1+z0[,1]),z0)
z1<-mvrnorm(n=obs,mu=c(0,0,0),Sigma=matrix(c(1,0.5,0.5,0.5,1,0.5,0.5,0.5,1),3,3))
z1[,1]<-pnorm(z1[,1])
z1[,2]<-(z1[,2]>0)
z1<-cbind(rep(1,obs)+z1%*%c(1,1,1)+rnorm(obs)*(1+z1[,1]),z1)
temp<-rqp(z0[,1],z0[,2:4])
boot<-sample(1:obs,100000,replace=TRUE)
a<-sample(1:length(temp$wq),100000,replace=TRUE,prob=temp$wq)
ypred<-temp$beta[1,a]+temp$beta[2,a]*z1[boot,2]+temp$beta[3,a]*z1[boot,3]+temp$beta[4,a]*z1[boot,4]
for(mm in 1:505) resultsm[i,((mm-1)*3+1):(mm*3)]<-quantile(ypred[1:m[mm]],c(0.05,0.25,0.5))
}

#nonparametric estimator
obs1<-100
obs0<-100
resultsn1<-matrix(,8000,18)
for(i in 1:8000){
set.seed(i)
z0<-sort(rnorm(obs0,sd=2))
z0<-cbind(5+z0+4*cos(z0)+2*sin(z0*3)+rcauchy(obs0)*(0.5+abs(z0)),z0)
z1<-sort(rnorm(obs1))
z1<-cbind(6+z1+4*cos(z1)+2*sin(z1*3)+rnorm(obs1)*(0.5+abs(z1)),z1)
temp<-qtet(z0[,1],z0[,2],z1[,1],z1[,2],prob=c(0.05,0.25,0.5),hmed0=obs0^-0.25,hmed1=obs1^-0.25,method0=1,method1=1,q0=c(-2,2),q1=c(-2,2),step0=0.1,step1=0.1)
resultsn1[i,]<-c(quantile(z0[,1],c(0.05,0.25,0.5)),quantile(z1[,1],c(0.05,0.25,0.5)),quantile(npirq(z0[,2],z0[,1],z0[,2],100,obs0^-0.25),c(0.05,0.25,0.5)),temp$res[,2],temp$qtet[,1],temp$qtet[,2])
}

obs1<-400
obs0<-400
resultsn2<-matrix(,4000,18)
for(i in 1:4000){
set.seed(i)
z0<-sort(rnorm(obs0,sd=2))
z0<-cbind(5+z0+4*cos(z0)+2*sin(z0*3)+rcauchy(obs0)*(0.5+abs(z0)),z0)
z1<-sort(rnorm(obs1))
z1<-cbind(6+z1+4*cos(z1)+2*sin(z1*3)+rnorm(obs1)*(0.5+abs(z1)),z1)
temp<-qtet(z0[,1],z0[,2],z1[,1],z1[,2],prob=c(0.05,0.25,0.5),hmed0=obs0^-0.25,hmed1=obs1^-0.25,method0=2,method1=2,q0=c(-2,2),q1=c(-2,2),step0=0.1,step1=0.1)
resultsn2[i,]<-c(quantile(z0[,1],c(0.05,0.25,0.5)),quantile(z1[,1],c(0.05,0.25,0.5)),quantile(npirq2(z0[,2],z0[,1],z0[,2],100,obs0^-0.25,q=c(-4,4),step=0.1),c(0.05,0.25,0.5)),temp$res[,2],temp$qtet[,1],temp$qtet[,2])
}

obs1<-1600
obs0<-1600
resultsn3<-matrix(,2000,18)
for(i in 1:2000){
set.seed(i)
z0<-sort(rnorm(obs0,sd=2))
z0<-cbind(5+z0+4*cos(z0)+2*sin(z0*3)+rcauchy(obs0)*(0.5+abs(z0)),z0)
z1<-sort(rnorm(obs1))
z1<-cbind(6+z1+4*cos(z1)+2*sin(z1*3)+rnorm(obs1)*(0.5+abs(z1)),z1)
temp<-qtet(z0[,1],z0[,2],z1[,1],z1[,2],prob=c(0.05,0.25,0.5),hmed0=obs0^-0.25,hmed1=obs1^-0.25,method0=2,method1=2,q0=c(-4,4),q1=c(-2,2),step0=0.1,step1=0.1)
resultsn3[i,]<-c(quantile(z0[,1],c(0.05,0.25,0.5)),quantile(z1[,1],c(0.05,0.25,0.5)),quantile(npirq2(z0[,2],z0[,1],z0[,2],100,obs0^-0.25,q=c(-5,5),step=0.05),c(0.05,0.25,0.5)),temp$res[,2],temp$qtet[,1],temp$qtet[,2])
}

obs1<-6400
obs0<-6400
resultsn4<-matrix(,1000,18)
for(i in 1:1000){
set.seed(i)
z0<-sort(rnorm(obs0,sd=2))
z0<-cbind(5+z0+4*cos(z0)+2*sin(z0*3)+rcauchy(obs0)*(0.5+abs(z0)),z0)
z1<-sort(rnorm(obs1))
z1<-cbind(6+z1+4*cos(z1)+2*sin(z1*3)+rnorm(obs1)*(0.5+abs(z1)),z1)
temp<-qtet(z0[,1],z0[,2],z1[,1],z1[,2],prob=c(0.05,0.25,0.5),hmed0=obs0^-0.25,hmed1=obs1^-0.25,method0=2,method1=2,q0=c(-4.5,4.5),q1=c(-2.25,2.25),step0=0.05,step1=0.05)
resultsn4[i,]<-c(quantile(z0[,1],c(0.05,0.25,0.5)),quantile(z1[,1],c(0.05,0.25,0.5)),quantile(npirq2(z0[,2],z0[,1],z0[,2],100,obs0^-0.25,q=c(-6,6),step=0.05),c(0.05,0.25,0.5)),temp$res[,2],temp$qtet[,1],temp$qtet[,2])
}

#true values (obtained by the quantiles of a sample of size 10^9)
obs<-1000000; tau<-c(0.05,0.25,0.5); res1<-rep(list(c()),length(tau)); res0<-rep(list(c()),length(tau)); resc<-rep(list(c()),length(tau)); a0<-b0<-a1<-b1<-ac<-bc<-matrix(,length(tau),1000)
set.seed(1)
z0<-mvrnorm(n=obs,mu=c(0,0,0),Sigma=matrix(c(1,0.5,0.5,0.5,1,0.5,0.5,0.5,1),3,3))
z0[,1]<-pnorm(z0[,1])
z0[,2]<-(z0[,2]>0)
z0<-rep(1,obs)+z0%*%c(1,1,1)+rcauchy(obs)*(1+z0[,1])
z1<-mvrnorm(n=obs,mu=c(0,0,0),Sigma=matrix(c(1,0.5,0.5,0.5,1,0.5,0.5,0.5,1),3,3))
z1[,1]<-pnorm(z1[,1])
z1[,2]<-(z1[,2]>0)
zc<-rep(1,obs)+z1%*%c(1,1,1)+rcauchy(obs)*(1+z1[,1])
z1<-rep(1,obs)+z1%*%c(1,1,1)+rnorm(obs)*(1+z1[,1])
q0<-q1<-qc<-matrix(,length(tau),2)
for(j in 1:length(tau)){
q0[j,]<-c(quantile(z0,tau[j]-0.001),quantile(z0,tau[j]+0.001))
q1[j,]<-c(quantile(z1,tau[j]-0.001),quantile(z1,tau[j]+0.001))
qc[j,]<-c(quantile(zc,tau[j]-0.001),quantile(zc,tau[j]+0.001))}
for(i in 1:1000){
z0<-mvrnorm(n=obs,mu=c(0,0,0),Sigma=matrix(c(1,0.5,0.5,0.5,1,0.5,0.5,0.5,1),3,3))
z0[,1]<-pnorm(z0[,1])
z0[,2]<-(z0[,2]>0)
z0<-rep(1,obs)+z0%*%c(1,1,1)+rcauchy(obs)*(1+z0[,1])
z1<-mvrnorm(n=obs,mu=c(0,0,0),Sigma=matrix(c(1,0.5,0.5,0.5,1,0.5,0.5,0.5,1),3,3))
z1[,1]<-pnorm(z1[,1])
z1[,2]<-(z1[,2]>0)
zc<-rep(1,obs)+z1%*%c(1,1,1)+rcauchy(obs)*(1+z1[,1])
z1<-rep(1,obs)+z1%*%c(1,1,1)+rnorm(obs)*(1+z1[,1])
for(j in 1:length(tau)){
res1[[j]]<-c(res1[[j]],z1[z1>=q1[j,1] & z1<=q1[j,2]])
res0[[j]]<-c(res0[[j]],z0[z0>=q0[j,1] & z0<=q0[j,2]])
resc[[j]]<-c(resc[[j]],zc[zc>=qc[j,1] & zc<=qc[j,2]])
b1[j,i]<-sum(z1<q1[j,1])
b0[j,i]<-sum(z0<q0[j,1])
bc[j,i]<-sum(zc<qc[j,1])
a1[j,i]<-sum(z1>q1[j,2])
a0[j,i]<-sum(z0>q0[j,2])
ac[j,i]<-sum(zc>qc[j,2])
}
}
n<-1000*obs
truep<-c()
for(j in 1:length(tau)) if(sum(b1[j,1:1000])<n*tau[j] & sum(a1[j,1:1000])<n*(1-tau[j])) truep<-c(truep,quantile(res1[[j]],(tau[j]*n-sum(b1[j,1:1000]))/(n-sum(a1[j,1:1000]+b1[j,1:1000])))) else break
for(j in 1:length(tau)) if(sum(b0[j,1:1000])<n*tau[j] & sum(a0[j,1:1000])<n*(1-tau[j])) truep<-c(truep,quantile(res0[[j]],(tau[j]*n-sum(b0[j,1:1000]))/(n-sum(a0[j,1:1000]+b0[j,1:1000])))) else break
for(j in 1:length(tau)) if(sum(bc[j,1:1000])<n*tau[j] & sum(ac[j,1:1000])<n*(1-tau[j])) truep<-c(truep,quantile(resc[[j]],(tau[j]*n-sum(bc[j,1:1000]))/(n-sum(ac[j,1:1000]+bc[j,1:1000])))) else break

obs<-1000000; tau<-c(0.05,0.25,0.5); res1<-rep(list(c()),length(tau)); res0<-rep(list(c()),length(tau)); resc<-rep(list(c()),length(tau)); a0<-b0<-a1<-b1<-ac<-bc<-matrix(,length(tau),1000)
set.seed(1)
z0<-rnorm(obs,sd=2)
z0<-5+z0+4*cos(z0)+2*sin(z0*3)+rcauchy(obs)*(0.5+abs(z0))
z1<-rnorm(obs)
zc<-5+z1+4*cos(z1)+2*sin(z1*3)+rcauchy(obs)*(0.5+abs(z1))
z1<-6+z1+4*cos(z1)+2*sin(z1*3)+rnorm(obs)*(0.5+abs(z1))
q0<-q1<-qc<-matrix(,length(tau),2)
for(j in 1:length(tau)){
q0[j,]<-c(quantile(z0,tau[j]-0.001),quantile(z0,tau[j]+0.001))
q1[j,]<-c(quantile(z1,tau[j]-0.001),quantile(z1,tau[j]+0.001))
qc[j,]<-c(quantile(zc,tau[j]-0.001),quantile(zc,tau[j]+0.001))}
for(i in 1:1000){
z0<-rnorm(obs,sd=2)
z0<-5+z0+4*cos(z0)+2*sin(z0*3)+rcauchy(obs)*(0.5+abs(z0))
z1<-rnorm(obs)
zc<-5+z1+4*cos(z1)+2*sin(z1*3)+rcauchy(obs)*(0.5+abs(z1))
z1<-6+z1+4*cos(z1)+2*sin(z1*3)+rnorm(obs)*(0.5+abs(z1))
for(j in 1:length(tau)){
res1[[j]]<-c(res1[[j]],z1[z1>=q1[j,1] & z1<=q1[j,2]])
res0[[j]]<-c(res0[[j]],z0[z0>=q0[j,1] & z0<=q0[j,2]])
resc[[j]]<-c(resc[[j]],zc[zc>=qc[j,1] & zc<=qc[j,2]])
b1[j,i]<-sum(z1<q1[j,1])
b0[j,i]<-sum(z0<q0[j,1])
bc[j,i]<-sum(zc<qc[j,1])
a1[j,i]<-sum(z1>q1[j,2])
a0[j,i]<-sum(z0>q0[j,2])
ac[j,i]<-sum(zc>qc[j,2])
}
}
n<-1000*obs
truen<-c()
for(j in 1:length(tau)) if(sum(b1[j,1:1000])<n*tau[j] & sum(a1[j,1:1000])<n*(1-tau[j])) truen<-c(truen,quantile(res1[[j]],(tau[j]*n-sum(b1[j,1:1000]))/(n-sum(a1[j,1:1000]+b1[j,1:1000])))) else break
for(j in 1:length(tau)) if(sum(b0[j,1:1000])<n*tau[j] & sum(a0[j,1:1000])<n*(1-tau[j])) truen<-c(truen,quantile(res0[[j]],(tau[j]*n-sum(b0[j,1:1000]))/(n-sum(a0[j,1:1000]+b0[j,1:1000])))) else break
for(j in 1:length(tau)) if(sum(bc[j,1:1000])<n*tau[j] & sum(ac[j,1:1000])<n*(1-tau[j])) truen<-c(truen,quantile(resc[[j]],(tau[j]*n-sum(bc[j,1:1000]))/(n-sum(ac[j,1:1000]+bc[j,1:1000])))) else break

#tables
skweness<-function(x) sum((x-mean(x))^3)/(length(x)-1)/sd(x)^3
kurtosis<-function(x) sum((x-mean(x))^4)/(length(x)-1)/sd(x)^4
table1<-matrix(,0,6)
for(i in 1:3){
temp<-c()
temp<-c(temp,round(mean(resultsp1[,6+i])-truep[3+i],4))
temp<-c(temp,round(sd(resultsp1[,6+i]),4))
temp<-c(temp,round(mean((truep[3+i]-resultsp1[,6+i])^2),4))
temp<-c(temp,round(skweness(resultsp1[,6+i]),4))
temp<-c(temp,round(kurtosis(resultsp1[,6+i]),4))
temp<-c(temp,round(mean((truep[3+i]-resultsp1[,i])^2)/mean((truep[3+i]-resultsp1[,6+i])^2),4))
table1<-rbind(table1,temp)
temp<-c()
temp<-c(temp,round(mean(resultsp2[,6+i])-truep[3+i],4))
temp<-c(temp,round(sd(resultsp2[,6+i]),4))
temp<-c(temp,round(mean((truep[3+i]-resultsp2[,6+i])^2),4))
temp<-c(temp,round(skweness(resultsp2[,6+i]),4))
temp<-c(temp,round(kurtosis(resultsp2[,6+i]),4))
temp<-c(temp,round(mean((truep[3+i]-resultsp2[,i])^2)/mean((truep[3+i]-resultsp2[,6+i])^2),4))
table1<-rbind(table1,temp)
temp<-c()
temp<-c(temp,round(mean(resultsp3[,6+i])-truep[3+i],4))
temp<-c(temp,round(sd(resultsp3[,6+i]),4))
temp<-c(temp,round(mean((truep[3+i]-resultsp3[,6+i])^2),4))
temp<-c(temp,round(skweness(resultsp3[,6+i]),4))
temp<-c(temp,round(kurtosis(resultsp3[,6+i]),4))
temp<-c(temp,round(mean((truep[3+i]-resultsp3[,i])^2)/mean((truep[3+i]-resultsp3[,6+i])^2),4))
table1<-rbind(table1,temp)
}
for(i in 1:3){
temp<-c()
temp<-c(temp,round(mean(resultsp1[,9+i])-truep[i],4))
temp<-c(temp,round(sd(resultsp1[,9+i]),4))
temp<-c(temp,round(mean((truep[i]-resultsp1[,9+i])^2),4))
temp<-c(temp,round(skweness(resultsp1[,9+i]),4))
temp<-c(temp,round(kurtosis(resultsp1[,9+i]),4))
temp<-c(temp,round(mean((truep[i]-resultsp1[,3+i])^2)/mean((truep[i]-resultsp1[,9+i])^2),4))
table1<-rbind(table1,temp)
temp<-c()
temp<-c(temp,round(mean(resultsp2[,9+i])-truep[i],4))
temp<-c(temp,round(sd(resultsp2[,9+i]),4))
temp<-c(temp,round(mean((truep[i]-resultsp2[,9+i])^2),4))
temp<-c(temp,round(skweness(resultsp2[,9+i]),4))
temp<-c(temp,round(kurtosis(resultsp2[,9+i]),4))
temp<-c(temp,round(mean((truep[i]-resultsp2[,3+i])^2)/mean((truep[i]-resultsp2[,9+i])^2),4))
table1<-rbind(table1,temp)
temp<-c()
temp<-c(temp,round(mean(resultsp3[,9+i])-truep[i],4))
temp<-c(temp,round(sd(resultsp3[,9+i]),4))
temp<-c(temp,round(mean((truep[i]-resultsp3[,9+i])^2),4))
temp<-c(temp,round(skweness(resultsp3[,9+i]),4))
temp<-c(temp,round(kurtosis(resultsp3[,9+i]),4))
temp<-c(temp,round(mean((truep[i]-resultsp3[,3+i])^2)/mean((truep[i]-resultsp3[,9+i])^2),4))
table1<-rbind(table1,temp)
}
for(i in 1:3){
temp<-c()
temp<-c(temp,round(mean(resultsp1[,9+i]-resultsp1[,12+i])-truep[i]+truep[6+i],4))
temp<-c(temp,round(sd(resultsp1[,9+i]-resultsp1[,12+i]),4))
temp<-c(temp,round(mean((truep[i]-truep[6+i]-resultsp1[,9+i]+resultsp1[,12+i])^2),4))
temp<-c(temp,round(skweness(resultsp1[,9+i]-resultsp1[,12+i]),4))
temp<-c(temp,round(kurtosis(resultsp1[,9+i]-resultsp1[,12+i]),4),NA)
table1<-rbind(table1,temp)
temp<-c()
temp<-c(temp,round(mean(resultsp2[,9+i]-resultsp2[,12+i])-truep[i]+truep[6+i],4))
temp<-c(temp,round(sd(resultsp2[,9+i]-resultsp2[,12+i]),4))
temp<-c(temp,round(mean((truep[i]-truep[6+i]-resultsp2[,9+i]+resultsp2[,12+i])^2),4))
temp<-c(temp,round(skweness(resultsp2[,9+i]-resultsp2[,12+i]),4))
temp<-c(temp,round(kurtosis(resultsp2[,9+i]-resultsp2[,12+i]),4),NA)
table1<-rbind(table1,temp)
temp<-c()
temp<-c(temp,round(mean(resultsp3[,9+i]-resultsp3[,12+i])-truep[i]+truep[6+i],4))
temp<-c(temp,round(sd(resultsp3[,9+i]-resultsp3[,12+i]),4))
temp<-c(temp,round(mean((truep[i]-truep[6+i]-resultsp3[,9+i]+resultsp3[,12+i])^2),4))
temp<-c(temp,round(skweness(resultsp3[,9+i]-resultsp3[,12+i]),4))
temp<-c(temp,round(kurtosis(resultsp3[,9+i]-resultsp3[,12+i]),4),NA)
table1<-rbind(table1,temp)
}
table2<-matrix(,0,5)
for(i in 1:3){
temp<-c()
for(alpha in c(0.01,0.05,0.1)) temp<-c(temp,sum((truep[i]-truep[6+i])<(resultsp1[,12+i]+qnorm(alpha/2)*resultsp1[,15+i]))+sum((truep[i]-truep[6+i])>(resultsp1[,12+i]+qnorm(1-alpha/2)*resultsp1[,15+i])))
temp<-c(temp/10000,median(resultsp1[,15+i])-sd(resultsp1[,12+i]))
temp<-c(temp,median(abs(resultsp1[,15+i]-sd(resultsp1[,12+i]))))
table2<-rbind(table2,temp)
temp<-c()
for(alpha in c(0.01,0.05,0.1)) temp<-c(temp,sum((truep[i]-truep[6+i])<(resultsp1[1:4000,12+i]+qnorm(alpha/2)*resultsp1[1:4000,18+i]))+sum((truep[i]-truep[6+i])>(resultsp1[1:4000,12+i]+qnorm(1-alpha/2)*resultsp1[1:4000,18+i])))
temp<-c(temp/4000,median(resultsp1[1:4000,18+i])-sd(resultsp1[,12+i]))
temp<-c(temp,median(abs(resultsp1[1:4000,18+i]-sd(resultsp1[,12+i]))))
table2<-rbind(table2,temp)
temp<-c()
for(alpha in c(0.01,0.05,0.1)) temp<-c(temp,sum((truep[i]-truep[6+i])<(resultsp2[,12+i]+qnorm(alpha/2)*resultsp2[,15+i]))+sum((truep[i]-truep[6+i])>(resultsp2[,12+i]+qnorm(1-alpha/2)*resultsp2[,15+i])))
temp<-c(temp/5000,median(resultsp2[,15+i])-sd(resultsp1[,12+i]))
temp<-c(temp,median(abs(resultsp2[,15+i]-sd(resultsp1[,12+i]))))
table2<-rbind(table2,temp)
temp<-c()
for(alpha in c(0.01,0.05,0.1)) temp<-c(temp,sum((truep[i]-truep[6+i])<(resultsp2[3001:5000,12+i]+qnorm(alpha/2)*resultsp2[3001:5000,18+i]))+sum((truep[i]-truep[6+i])>(resultsp2[3001:5000,12+i]+qnorm(1-alpha/2)*resultsp2[3001:5000,18+i])))
temp<-c(temp/2000,median(resultsp2[3001:5000,18+i])-sd(resultsp2[,12+i]))
temp<-c(temp,median(abs(resultsp2[3001:5000,18+i]-sd(resultsp2[,12+i]))))
table2<-rbind(table2,temp)
temp<-c()
for(alpha in c(0.01,0.05,0.1)) temp<-c(temp,sum((truep[i]-truep[6+i])<(resultsp3[,12+i]+qnorm(alpha/2)*resultsp3[,15+i]))+sum((truep[i]-truep[6+i])>(resultsp3[,12+i]+qnorm(1-alpha/2)*resultsp3[,15+i])))
temp<-c(temp/2500,median(resultsp3[,15+i])-sd(resultsp3[,12+i]))
temp<-c(temp,median(abs(resultsp3[,15+i]-sd(resultsp3[,12+i]))))
table2<-rbind(table2,temp)
temp<-c()
for(alpha in c(0.01,0.05,0.1)) temp<-c(temp,sum((truep[i]-truep[6+i])<(resultsp3[1:1000,12+i]+qnorm(alpha/2)*resultsp3[1:1000,18+i]))+sum((truep[i]-truep[6+i])>(resultsp3[1:1000,12+i]+qnorm(1-alpha/2)*resultsp3[1:1000,18+i])))
temp<-c(temp/1000,median(resultsp3[1:1000,18+i])-sd(resultsp3[,12+i]))
temp<-c(temp,median(abs(resultsp3[1:1000,18+i]-sd(resultsp3[,12+i]))))
table2<-rbind(table2,temp)
}
table3<-matrix(,0,7)
for(i in c(3,12,30,300,1515)){
temp<-c()
temp<-c(temp,round(mean(resultsm[,i])-truep[9],4))
temp<-c(temp,round(sd(resultsm[,i]),4))
temp<-c(temp,round(mean((truep[9]-resultsm[,i])^2),4))
temp<-c(temp,round(skweness(resultsm[,i]),4))
temp<-c(temp,round(kurtosis(resultsm[,i]),4))
temp<-c(temp,round(mean((truep[9]-resultsm[,i])^2)/mean((truep[9]-resultsp2[,22])^2),4))
temp<-c(temp,round(cor(resultsm[,i],resultsp2[,22]),4))
table3<-rbind(table3,temp)
}
temp<-c()
temp<-c(temp,round(mean(resultsp2[,22])-truep[9],4))
temp<-c(temp,round(sd(resultsp2[,22]),4))
temp<-c(temp,round(mean((truep[9]-resultsp2[,22])^2),4))
temp<-c(temp,round(skweness(resultsp2[,22]),4))
temp<-c(temp,round(kurtosis(resultsp2[,22]),4),1,1)
table3<-rbind(table3,temp)
table4<-matrix(,0,6)
for(i in 1:3){
temp<-c()
temp<-c(temp,round(mean(resultsn1[,6+i])-truen[i+3],4))
temp<-c(temp,round(sd(resultsn1[,6+i]),4))
temp<-c(temp,round(mean((truen[i+3]-resultsn1[,6+i])^2),4))
temp<-c(temp,round(skweness(resultsn1[,6+i]),4))
temp<-c(temp,round(kurtosis(resultsn1[,6+i]),4))
temp<-c(temp,round(mean((truen[i+3]-resultsn1[,i])^2)/mean((truen[i+3]-resultsn1[,6+i])^2),4))
table4<-rbind(table4,temp)
temp<-c()
temp<-c(temp,round(mean(resultsn2[,6+i])-truen[i+3],4))
temp<-c(temp,round(sd(resultsn2[,6+i]),4))
temp<-c(temp,round(mean((truen[i+3]-resultsn2[,6+i])^2),4))
temp<-c(temp,round(skweness(resultsn2[,6+i]),4))
temp<-c(temp,round(kurtosis(resultsn2[,6+i]),4))
temp<-c(temp,round(mean((truen[i+3]-resultsn2[,i])^2)/mean((truen[i+3]-resultsn2[,6+i])^2),4))
table4<-rbind(table4,temp)
temp<-c()
temp<-c(temp,round(mean(resultsn3[,6+i])-truen[i+3],4))
temp<-c(temp,round(sd(resultsn3[,6+i]),4))
temp<-c(temp,round(mean((truen[i+3]-resultsn3[,6+i])^2),4))
temp<-c(temp,round(skweness(resultsn3[,6+i]),4))
temp<-c(temp,round(kurtosis(resultsn3[,6+i]),4))
temp<-c(temp,round(mean((truen[i+3]-resultsn3[,i])^2)/mean((truen[i+3]-resultsn3[,6+i])^2),4))
table4<-rbind(table4,temp)
temp<-c()
temp<-c(temp,round(mean(resultsn4[,6+i])-truen[i+3],4))
temp<-c(temp,round(sd(resultsn4[,6+i]),4))
temp<-c(temp,round(mean((truen[i+3]-resultsn4[,6+i])^2),4))
temp<-c(temp,round(skweness(resultsn4[,6+i]),4))
temp<-c(temp,round(kurtosis(resultsn4[,6+i]),4))
temp<-c(temp,round(mean((truen[i+3]-resultsn4[,i])^2)/mean((truen[i+3]-resultsn4[,6+i])^2),4))
table4<-rbind(table4,temp)
}
for(i in 1:3){
temp<-c()
temp<-c(temp,round(mean(resultsn1[,9+i])-truen[i],4))
temp<-c(temp,round(sd(resultsn1[,9+i]),4))
temp<-c(temp,round(mean((truen[i]-resultsn1[,9+i])^2),4))
temp<-c(temp,round(skweness(resultsn1[,9+i]),4))
temp<-c(temp,round(kurtosis(resultsn1[,9+i]),4))
temp<-c(temp,round(mean((truen[i]-resultsn1[,3+i])^2)/mean((truen[i]-resultsn1[,9+i])^2),4))
table4<-rbind(table4,temp)
temp<-c()
temp<-c(temp,round(mean(resultsn2[,9+i])-truen[i],4))
temp<-c(temp,round(sd(resultsn2[,9+i]),4))
temp<-c(temp,round(mean((truen[i]-resultsn2[,9+i])^2),4))
temp<-c(temp,round(skweness(resultsn2[,9+i]),4))
temp<-c(temp,round(kurtosis(resultsn2[,9+i]),4))
temp<-c(temp,round(mean((truen[i]-resultsn2[,3+i])^2)/mean((truen[i]-resultsn2[,9+i])^2),4))
table4<-rbind(table4,temp)
temp<-c()
temp<-c(temp,round(mean(resultsn3[,9+i])-truen[i],4))
temp<-c(temp,round(sd(resultsn3[,9+i]),4))
temp<-c(temp,round(mean((truen[i]-resultsn3[,9+i])^2),4))
temp<-c(temp,round(skweness(resultsn3[,9+i]),4))
temp<-c(temp,round(kurtosis(resultsn3[,9+i]),4))
temp<-c(temp,round(mean((truen[i]-resultsn3[,3+i])^2)/mean((truen[i]-resultsn3[,9+i])^2),4))
table4<-rbind(table4,temp)
temp<-c()
temp<-c(temp,round(mean(resultsn4[,9+i])-truen[i],4))
temp<-c(temp,round(sd(resultsn4[,9+i]),4))
temp<-c(temp,round(mean((truen[i]-resultsn4[,9+i])^2),4))
temp<-c(temp,round(skweness(resultsn4[,9+i]),4))
temp<-c(temp,round(kurtosis(resultsn4[,9+i]),4))
temp<-c(temp,round(mean((truen[i]-resultsn4[,3+i])^2)/mean((truen[i]-resultsn4[,9+i])^2),4))
table4<-rbind(table4,temp)
}
for(i in 1:3){
temp<-c()
temp<-c(temp,round(mean(resultsn1[,12+i])-truen[i]+truen[6+i],4))
temp<-c(temp,round(sd(resultsn1[,12+i]),4))
temp<-c(temp,round(mean((truen[i]-truen[6+i]-resultsn1[,12+i])^2),4))
temp<-c(temp,round(skweness(resultsn1[,12+i]),4))
temp<-c(temp,round(kurtosis(resultsn1[,12+i]),4),NA)
table4<-rbind(table4,temp)
temp<-c()
temp<-c(temp,round(mean(resultsn2[,12+i])-truen[i]+truen[6+i],4))
temp<-c(temp,round(sd(resultsn2[,12+i]),4))
temp<-c(temp,round(mean((truen[i]-truen[6+i]-resultsn2[,12+i])^2),4))
temp<-c(temp,round(skweness(resultsn2[,12+i]),4))
temp<-c(temp,round(kurtosis(resultsn2[,12+i]),4),NA)
table4<-rbind(table4,temp)
temp<-c()
temp<-c(temp,round(mean(resultsn3[,12+i])-truen[i]+truen[6+i],4))
temp<-c(temp,round(sd(resultsn3[,12+i]),4))
temp<-c(temp,round(mean((truen[i]-truen[6+i]-resultsn3[,12+i])^2),4))
temp<-c(temp,round(skweness(resultsn3[,12+i]),4))
temp<-c(temp,round(kurtosis(resultsn3[,12+i]),4),NA)
table4<-rbind(table4,temp)
temp<-c()
temp<-c(temp,round(mean(resultsn4[,12+i])-truen[i]+truen[6+i],4))
temp<-c(temp,round(sd(resultsn4[,12+i]),4))
temp<-c(temp,round(mean((truen[i]-truen[6+i]-resultsn4[,12+i])^2),4))
temp<-c(temp,round(skweness(resultsn4[,12+i]),4))
temp<-c(temp,round(kurtosis(resultsn4[,12+i]),4),NA)
table4<-rbind(table4,temp)
}
table5<-matrix(,0,6)
for(i in 1:3){
temp<-c()
for(alpha in c(0.01,0.05,0.1)) temp<-c(temp,sum((truen[i]-truen[6+i])<(resultsn1[,12+i]+qnorm(alpha/2)*resultsn1[,15+i]))+sum((truen[i]-truen[6+i])>(resultsn1[,12+i]+qnorm(1-alpha/2)*resultsn1[,15+i])))
temp<-temp/8000
temp<-c(temp,sd(resultsn1[,12+i]))
temp<-c(temp,median(resultsn1[,15+i])-sd(resultsn1[,12+i]))
temp<-c(temp,median(abs(resultsn1[,15+i]-sd(resultsn1[,12+i]))))
table5<-rbind(table5,temp)
temp<-c()
for(alpha in c(0.01,0.05,0.1)) temp<-c(temp,sum((truen[i]-truen[6+i])<(resultsn2[,12+i]+qnorm(alpha/2)*resultsn2[,15+i]))+sum((truen[i]-truen[6+i])>(resultsn2[,12+i]+qnorm(1-alpha/2)*resultsn2[,15+i])))
temp<-temp/4000
temp<-c(temp,sd(resultsn2[,12+i]))
temp<-c(temp,median(resultsn2[,15+i])-sd(resultsn2[,12+i]))
temp<-c(temp,median(abs(resultsn2[,15+i]-sd(resultsn2[,12+i]))))
table5<-rbind(table5,temp)
temp<-c()
for(alpha in c(0.01,0.05,0.1)) temp<-c(temp,sum((truen[i]-truen[6+i])<(resultsn3[,12+i]+qnorm(alpha/2)*resultsn3[,15+i]))+sum((truen[i]-truen[6+i])>(resultsn3[,12+i]+qnorm(1-alpha/2)*resultsn3[,15+i])))
temp<-temp/2000
temp<-c(temp,sd(resultsn3[,12+i]))
temp<-c(temp,median(resultsn3[,15+i])-sd(resultsn3[,12+i]))
temp<-c(temp,median(abs(resultsn3[,15+i]-sd(resultsn3[,12+i]))))
table5<-rbind(table5,temp)
temp<-c()
for(alpha in c(0.01,0.05,0.1)) temp<-c(temp,sum((truen[i]-truen[6+i])<(resultsn4[,12+i]+qnorm(alpha/2)*resultsn4[,15+i]))+sum((truen[i]-truen[6+i])>(resultsn4[,12+i]+qnorm(1-alpha/2)*resultsn4[,15+i])))
temp<-temp/1000
temp<-c(temp,sd(resultsn4[,12+i]))
temp<-c(temp,median(resultsn4[,15+i])-sd(resultsn4[,12+i]))
temp<-c(temp,median(abs(resultsn4[,15+i]-sd(resultsn4[,12+i]))))
table5<-rbind(table5,temp)
}