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


datacps<-read.table("data_cps.txt",header=TRUE)

attach(datacps)
n<-length(real_hrwage)
NumBoot = 500
lvl = seq(0.01, 0.99, 0.008)
Numtau<-length(lvl)
age2<-age^2
wgt<-wgt/1e+03
for(level in unique(datacps$year)){
  datacps[paste("T")] <- ifelse(datacps$year == 1995, 0, 1)
}
detach(datacps)


#####
##MTE
#####
#MTEfun returns: gamma.factual, gamma.counterfactual, mu_t, mu_t-1

MTEfun = function(dataset) {
attach(dataset)
y<-real_hrwage
n<-length(real_hrwage)

pscore<-glm(T~dummy(state)+female+dummy(age)+dummy(race)+metropolitan+rural+dummy(educ)+married+public+dummy(codocup)+dummy(industry),family=binomial(link = "logit"),weights=wgt)$fitted

    p<-sum(wgt*T)/sum(wgt)
    wT11<-wgt*T
    wT10<-wgt*T*(1-pscore)/pscore
    wT00<-wgt*(1-T)
    
    mte11 = coef(lm(y ~ 1, weights = wT11))
    mte10 = coef(lm(y ~ 1, weights = wT10))
    mte00 = coef(lm(y ~ 1, weights = wT00))

detach(dataset)
mte<-(mte11-mte00)/mte00
mtes<-(mte10-mte00)/mte00
cbind(mte, mtes, mte11, mte00, mte10)
}

#####
##QTE
#####
#QTEfun returns: gic.factual, gic.counterfactual, q_t, q_t-1

QTEfun = function(dataset) {
attach(dataset)
y<-real_hrwage
n<-length(real_hrwage)
lvl = seq(0.01, 0.99, 0.008)
Numtau<-length(lvl)

pscore<-glm(T~dummy(state)+female+dummy(age)+dummy(race)+metropolitan+rural+dummy(educ)+married+public+dummy(codocup)+dummy(industry),family=binomial(link = "logit"),weights=wgt)$fitted

    p<-sum(wgt*T)/sum(wgt)
    wT11<-wgt*T
    wT10<-wgt*T*(1-pscore)/pscore
    wT00<-wgt*(1-T)

    qtt11 = foreach(i = 1:length(lvl), .combine = c) %do% {
        coef(rq(y ~ 1, tau=lvl[i], weights = wT11))
    }

   qtt10 = foreach(i = 1:length(lvl), .combine = c) %do% {
        coef(rq(y ~ 1, tau=lvl[i], weights = wT10))
    }
    
    qtt00 = foreach(i = 1:length(lvl), .combine = c) %do% {
        coef(rq(y ~ 1, tau=lvl[i], weights = wT00))
    }

detach(dataset)  
  gic<-((qtt11-qtt00)/qtt00)
  gics<-((qtt10-qtt00)/qtt00)
cbind(gic, gics, qtt11, qtt00, qtt10)
}


# Compute the beta (GICF, GICC) and alpha (mu_t, mu_t-1, gammaF, gammaC) using the entire sample

beta = QTEfun(datacps)
alpha = MTEfun(datacps)

gic.f<-beta[,1]
gic.c<-beta[,2]

gic.avg.f<-mean(beta[,1])
gic.avg.c<-mean(beta[,2])

#####
#PLOT GIC
#####
plot(lvl, beta[,1], type="l", xlab = "quantiles", ylab = "GIC", ylim=c(0,0.4),lwd=2, lty =1, cex = 0.5)
title(paste("GIC for US 1995-2007"), cex = 0.75)
lines(lvl,beta[,2],type="l", lwd=2, lty=4, cex = 0.5)
lines(lvl,rep(gic.avg.f,Numtau),type="l", lwd=2, lty =6, cex = 0.5)
lines(lvl,rep(gic.avg.c,Numtau),type="l", lwd=2, lty=3, cex = 0.5)
legend(0, 0.4, c("GIC", "Counterfactual GIC", "Average", "Counterfactual Average"), lwd = c(2,2,2,2), lty=c(1,4,6,3), cex=0.85)

plot(ecdf(beta[,3]), xlab = "quantiles", ylab = "CDF Income", main= "CDF Income for US 1995-2007", lty =1, lwd=1, do.points=FALSE, verticals=TRUE)
plot(ecdf(beta[,4]), lty =2, lwd=2, add=TRUE, do.points=FALSE, verticals=TRUE)
plot(ecdf(beta[,5]), lty=3, lwd=3, add=TRUE, do.points=FALSE, verticals=TRUE)
legend(50, 0.2, c("Income 2007", "Income 1995", "Counterfactual Income"), lwd = c(1, 2, 3), lty=c(1,2,3), cex=0.85)

# Lorenz Curves
#Factual
loc.t.aux.f<-beta[,3]
o.t.f <- order(loc.t.aux.f)
loc.t.aux.f <- loc.t.aux.f[o.t.f]
loc.t.f<-cumsum(loc.t.aux.f)/sum(loc.t.aux.f)

loc.lt.aux.f<-beta[,4]
o.lt.f <- order(loc.lt.aux.f)
loc.lt.aux.f <- loc.lt.aux.f[o.lt.f]
loc.lt.f<-cumsum(loc.lt.aux.f)/sum(loc.lt.aux.f)

delta.loc.f<-loc.t.f-loc.lt.f

#CounterFactual

loc.t.aux.c<-beta[,5]
o.t.c <- order(loc.t.aux.c)
loc.t.aux.c <- loc.t.aux.c[o.t.c]
loc.t.c<-cumsum(loc.t.aux.c)/sum(loc.t.aux.c)

loc.lt.aux.c<-beta[,4]
o.lt.c <- order(loc.lt.aux.c)
loc.lt.aux.c <- loc.lt.aux.c[o.lt.c]
loc.lt.c<-cumsum(loc.lt.aux.c)/sum(loc.lt.aux.c)

delta.loc.c<-loc.t.c-loc.lt.c

#####
#PLOT Lorenz
#####
plot(lvl, loc.t.f, type="l", xlab = "quantiles", ylab = "Cumulative share", ylim=c(0,1),lwd=2, cex = 0.5)
title(paste("Lorenz Curves for US 1995-2007"), cex = 0.75)
lines(lvl,loc.lt.f,type="l", lwd=2,cex = 0.5, lty=2)
lines(lvl, loc.t.c, type="l", lwd=2, lty=3, cex = 0.5)
legend(0, 1, c("Factual 2007", "Factual 1995", "Counterfactual 2007"), lwd = c(2, 2,2), lty=c(1,2,3), cex=0.85)

#####
#PLOT Delta Lorenz
#####
plot(lvl, delta.loc.f, type="l", xlab = "quantiles", ylab = "Delta", ylim=c(-0.017,0.001),lwd=2, cex = 0.5)
title(paste("Delta Lorenz Curves for US 1995-2007"), cex = 0.75)
lines(lvl, delta.loc.c, type="l", lwd=2, cex = 0.5, lty=2)
legend(0, -0.015, c("Delta Lorenz Factual", "Delta Lorenz Counterfactual"), lwd = c(2, 2), lty=c(1,2), cex=0.85)


########
#Testing
########


bootQTE <- function(data, B, n) {
  x<-1:n	
  resamples <- lapply(1:B, function(i) data[sample(x, replace = TRUE),])
  lapply(resamples, QTEfun)
}

bootMTE <- function(data, B, n) {
  x<-1:n	
  resamples <- lapply(1:B, function(i) data[sample(x, replace = TRUE),])
  lapply(resamples, MTEfun)
}



# Compute the bootstrap for alpha and beta

bootbeta = bootQTE(datacps, NumBoot, n)

#Organize the alpha and beta vectors
# f for factural and c for couterfactual
# bootbeta.aux is bootbeta organized

boot.mte.f.aux<-array(0,dim=c(NumBoot,1))
boot.mte.c.aux<-array(0,dim=c(NumBoot,1))
boot.mu.t.aux<-array(0,dim=c(NumBoot,1))
boot.mu.lt.aux<-array(0,dim=c(NumBoot,1))
boot.mu.ts.aux<-array(0,dim=c(NumBoot,1))

gic.f.aux<-array(0,dim=c(Numtau,1))
gic.c.aux<-array(0,dim=c(Numtau,1))
q.t.aux<-array(0,dim=c(Numtau,1))
q.lt.aux<-array(0,dim=c(Numtau,1))
q.ts.aux<-array(0,dim=c(Numtau,1))

boot.gic.f.aux<-array(0,dim=c(NumBoot,Numtau))
boot.gic.c.aux<-array(0,dim=c(NumBoot,Numtau))
boot.q.t.aux<-array(0,dim=c(NumBoot,Numtau))
boot.q.lt.aux<-array(0,dim=c(NumBoot,Numtau))
boot.q.ts.aux<-array(0,dim=c(NumBoot,Numtau))

for (i in 1:NumBoot){
for (k in 1:Numtau){
  gic.f.aux[k]<-beta[k,1]
  gic.c.aux[k]<-beta[k,2]
  q.t.aux[k]<-beta[k,3]
  q.lt.aux[k]<-beta[k,4]
  q.ts.aux[k]<-beta[k,5]

  boot.gic.f.aux[i,k]<-bootbeta[[i]][k,1]
  boot.gic.c.aux[i,k]<-bootbeta[[i]][k,2]
  boot.q.t.aux[i,k]<-bootbeta[[i]][k,3]
  boot.q.lt.aux[i,k]<-bootbeta[[i]][k,4]
  boot.q.ts.aux[i,k]<-bootbeta[[i]][k,5]

  boot.mte.f.aux[i]<-mean(bootbeta[[i]][,1])
  boot.mte.c.aux[i]<-mean(bootbeta[[i]][,2])
  boot.mu.t.aux[i]<-mean(bootbeta[[i]][,3])
  boot.mu.lt.aux[i]<-mean(bootbeta[[i]][,4])
  boot.mu.ts.aux[i]<-mean(bootbeta[[i]][,5])
}
}

boot.delta.loc.f<-array(0,dim=c(NumBoot, Numtau))
boot.delta.loc.c<-array(0,dim=c(NumBoot, Numtau))

for (i in 1:NumBoot){
	loc.t.aux.1<-boot.q.t.aux[i,]
	o.t.1 <- order(loc.t.aux.1)
	loc.t.aux.1 <- loc.t.aux.1[o.t.1]
	loc.t.1<-cumsum(loc.t.aux.1)/sum(loc.t.aux.1)
	
	loc.lt.aux.1<-boot.q.lt.aux[i,]
	o.lt.1 <- order(loc.lt.aux.1)
	loc.lt.aux.1 <- loc.lt.aux.1[o.lt.1]
	loc.lt.1<-cumsum(loc.lt.aux.1)/sum(loc.lt.aux.1)

	loc.ts.aux.2<-boot.q.ts.aux[i,]
	o.ts.2 <- order(loc.ts.aux.2)
	loc.ts.aux.2 <- loc.ts.aux.2[o.ts.2]
	loc.ts.2<-cumsum(loc.ts.aux.2)/sum(loc.ts.aux.2)

	boot.delta.loc.f[i,]<-loc.t.1-loc.lt.1
	boot.delta.loc.c[i,]<-loc.ts.2-loc.lt.1
}

# Compute MTE Factual and Counterfactual and their CI
avg5.f<-array(0,dim=c(1,3))
avg10.f<-array(0,dim=c(1,3))
avg5.f[1,1]<-gic.avg.f
avg10.f[1,1]<-gic.avg.f
avg5.f[1,2]<-quantile(boot.mte.f.aux,0.025)
avg5.f[1,3]<-quantile(boot.mte.f.aux,0.975)
avg10.f[1,2]<-quantile(boot.mte.f.aux,0.05)
avg10.f[1,3]<-quantile(boot.mte.f.aux,0.95)

avg5.c<-array(0,dim=c(1,3))
avg10.c<-array(0,dim=c(1,3))
avg5.c[1,1]<-gic.avg.c
avg10.c[1,1]<-gic.avg.c
avg5.c[1,2]<-quantile(boot.mte.c.aux,0.025)
avg5.c[1,3]<-quantile(boot.mte.c.aux,0.975)
avg10.c[1,2]<-quantile(boot.mte.c.aux,0.05)
avg10.c[1,3]<-quantile(boot.mte.c.aux,0.95)


# Compute GIC Factual and Counterfactual and their CI
gic5.f<-array(0,dim=c(Numtau,3))
gic10.f<-array(0,dim=c(Numtau,3))
gic5.f[,1]<-beta[,1]
gic10.f[,1]<-beta[,1]

gic5.c<-array(0,dim=c(Numtau,3))
gic10.c<-array(0,dim=c(Numtau,3))
gic5.c[,1]<-beta[,2]
gic10.c[,1]<-beta[,2]


for (j in 1:Numtau){
gic5.f[j,2]<-quantile(boot.gic.f.aux[,j],0.025)
gic5.f[j,3]<-quantile(boot.gic.f.aux[,j],0.975)
gic10.f[j,2]<-quantile(boot.gic.f.aux[,j],0.05)
gic10.f[j,3]<-quantile(boot.gic.f.aux[,j],0.95)

gic5.c[j,2]<-quantile(boot.gic.c.aux[,j],0.025)
gic5.c[j,3]<-quantile(boot.gic.c.aux[,j],0.975)
gic10.c[j,2]<-quantile(boot.gic.c.aux[,j],0.05)
gic10.c[j,3]<-quantile(boot.gic.c.aux[,j],0.95)
}


# Compute Delta Factual and Counterfactual and their CI
delta5.f<-array(0,dim=c(Numtau,3))
delta10.f<-array(0,dim=c(Numtau,3))
delta5.f[,1]<-delta.loc.f
delta10.f[,1]<-delta.loc.f

delta5.c<-array(0,dim=c(Numtau,3))
delta10.c<-array(0,dim=c(Numtau,3))
delta5.c[,1]<-delta.loc.c
delta10.c[,1]<-delta.loc.c


for (j in 1:Numtau){
delta5.f[j,2]<-quantile(boot.delta.loc.f[,j],0.025)
delta5.f[j,3]<-quantile(boot.delta.loc.f[,j],0.975)
delta10.f[j,2]<-quantile(boot.delta.loc.f[,j],0.05)
delta10.f[j,3]<-quantile(boot.delta.loc.f[,j],0.95)

delta5.c[j,2]<-quantile(boot.delta.loc.c[,j],0.025)
delta5.c[j,3]<-quantile(boot.delta.loc.c[,j],0.975)
delta10.c[j,2]<-quantile(boot.delta.loc.c[,j],0.05)
delta10.c[j,3]<-quantile(boot.delta.loc.c[,j],0.95)
}



# Compute the variance over the bootstraps

var.mte.f<-var(boot.mte.f.aux)
var.mte.c<-var(boot.mte.c.aux)
var.mu.t<-var(boot.mu.t.aux)
var.mu.lt<-var(boot.mte.f.aux)

var.gic.f<-array(0,dim=c(Numtau))
var.gic.c<-array(0,dim=c(Numtau))
var.delta.f<-array(0,dim=c(Numtau))
var.delta.c<-array(0,dim=c(Numtau))
var.q.t<-array(0,dim=c(Numtau))
var.q.lt<-array(0,dim=c(Numtau))
for (j in 1:Numtau){
var.gic.f[j]<-var(boot.gic.f.aux[,j])
var.gic.c[j]<-var(boot.gic.c.aux[,j])
var.delta.f[j]<-var(boot.delta.loc.f[,j])
var.delta.c[j]<-var(boot.delta.loc.c[,j])
var.q.t[j]<-var(boot.q.t.aux[,j])
var.q.lt[j]<-var(boot.q.lt.aux[,j])
}

# Compute the Wb() over bootstraps for the statistic of test
# r(\tau)=0

Wb.gic.f<-array(0,dim=c(NumBoot,Numtau))
Wb.gic.c<-array(0,dim=c(NumBoot,Numtau))
Wb.delta.f<-array(0,dim=c(NumBoot,Numtau))
Wb.delta.c<-array(0,dim=c(NumBoot,Numtau))
for (j in 1:NumBoot){
Wb.gic.f[j,]<-(boot.gic.f.aux[j,]-gic.f.aux)
Wb.gic.c[j,]<-(boot.gic.c.aux[j,]-gic.c.aux)
Wb.delta.f[j,]<-(boot.delta.loc.f[j,]-delta.loc.f)
Wb.delta.c[j,]<-(boot.delta.loc.c[j,]-delta.loc.c)
}

# r(\tau)=rhat

Wb.gic.f.rh<-array(0,dim=c(NumBoot,Numtau))
Wb.gic.c.rh<-array(0,dim=c(NumBoot,Numtau))
Wb.gic.gic.rh<-array(0,dim=c(NumBoot,Numtau))
for (j in 1:NumBoot){
Wb.gic.f.rh[j,]<-(boot.gic.f.aux[j,]-boot.mte.f.aux[j,])-(gic.f.aux-gic.avg.f)
Wb.gic.c.rh[j,]<-(boot.gic.c.aux[j,]-boot.mte.c.aux[j,])-(gic.c.aux-gic.avg.c)
Wb.gic.gic.rh[j,]<-((boot.gic.f.aux[j,]-boot.mte.f.aux[j,])-(gic.f.aux-gic.avg.f))-((boot.gic.c.aux[j,]-boot.mte.c.aux[j,])-(gic.c.aux-gic.avg.c))
}


# Compute the statistic of tests over the bootstraps
# r(\tau)=0

VbKS.gic.f<-array(0,dim=c(NumBoot))
VbCVM.gic.f<-array(0,dim=c(NumBoot))
VbKS.gic.c<-array(0,dim=c(NumBoot))
VbCVM.gic.c<-array(0,dim=c(NumBoot))
VbKS.delta.f<-array(0,dim=c(NumBoot))
VbCVM.delta.f<-array(0,dim=c(NumBoot))
VbKS.delta.c<-array(0,dim=c(NumBoot))
VbCVM.delta.c<-array(0,dim=c(NumBoot))
for (j in 1:NumBoot){
VbKS.gic.f[j]<-max(abs(Wb.gic.f[j,])) #KS - No SQRT
VbCVM.gic.f[j]<-sum(abs(Wb.gic.f[j,])) #CVM - No SQRT
VbKS.gic.c[j]<-max(abs(Wb.gic.c[j,])) #KS - No SQRT
VbCVM.gic.c[j]<-sum(abs(Wb.gic.c[j,])) #CVM - No SQRT
VbKS.delta.f[j]<-max(abs(Wb.delta.f[j,])) #KS - No SQRT
VbCVM.delta.f[j]<-sum(abs(Wb.delta.f[j,])) #CVM - No SQRT
VbKS.delta.c[j]<-max(abs(Wb.delta.c[j,])) #KS - No SQRT
VbCVM.delta.c[j]<-sum(abs(Wb.delta.c[j,])) #CVM - No SQRT
}

# r(\tau)=rhat
VbKS.gic.f.rh<-array(0,dim=c(NumBoot))
VbCVM.gic.f.rh<-array(0,dim=c(NumBoot))
VbKS.gic.c.rh<-array(0,dim=c(NumBoot))
VbCVM.gic.c.rh<-array(0,dim=c(NumBoot))
VbKS.gic.gic.rh<-array(0,dim=c(NumBoot))
VbCVM.gic.gic.rh<-array(0,dim=c(NumBoot))
for (j in 1:NumBoot){
VbKS.gic.f.rh[j]<-max(abs(Wb.gic.f.rh[j,])) #KS - No SQRT
VbCVM.gic.f.rh[j]<-sum(abs(Wb.gic.f.rh[j,])) #CVM - No SQRT
VbKS.gic.c.rh[j]<-max(abs(Wb.gic.c.rh[j,])) #KS - No SQRT
VbCVM.gic.c.rh[j]<-sum(abs(Wb.gic.c.rh[j,])) #CVM - No SQRT
VbKS.gic.gic.rh[j]<-max(abs(Wb.gic.c.rh[j,])) #KS - No SQRT
VbCVM.gic.gic.rh[j]<-sum(abs(Wb.gic.c.rh[j,])) #CVM - No SQRT
}

# Compute critical values by taking empirical quantile
# r(\tau)=0
#Factual
cvKS.gic.f.01<-quantile(VbKS.gic.f,0.99)
cvKS.gic.f.05<-quantile(VbKS.gic.f,0.95)
cvKS.gic.f.10<-quantile(VbKS.gic.f,0.90)
cvCVM.gic.f.01<-quantile(VbCVM.gic.f,0.99)
cvCVM.gic.f.05<-quantile(VbCVM.gic.f,0.95)
cvCVM.gic.f.10<-quantile(VbCVM.gic.f,0.90)

cvKS.delta.f.01<-quantile(VbKS.delta.f,0.99)
cvKS.delta.f.05<-quantile(VbKS.delta.f,0.95)
cvKS.delta.f.10<-quantile(VbKS.delta.f,0.90)
cvCVM.delta.f.01<-quantile(VbCVM.delta.f,0.99)
cvCVM.delta.f.05<-quantile(VbCVM.delta.f,0.95)
cvCVM.delta.f.10<-quantile(VbCVM.delta.f,0.90)

#Counter
cvKS.gic.c.01<-quantile(VbKS.gic.c,0.99)
cvKS.gic.c.05<-quantile(VbKS.gic.c,0.95)
cvKS.gic.c.10<-quantile(VbKS.gic.c,0.90)
cvCVM.gic.c.01<-quantile(VbCVM.gic.c,0.99)
cvCVM.gic.c.05<-quantile(VbCVM.gic.c,0.95)
cvCVM.gic.c.10<-quantile(VbCVM.gic.c,0.90)

cvKS.delta.c.01<-quantile(VbKS.delta.c,0.99)
cvKS.delta.c.05<-quantile(VbKS.delta.c,0.95)
cvKS.delta.c.10<-quantile(VbKS.delta.c,0.90)
cvCVM.delta.c.01<-quantile(VbCVM.delta.c,0.99)
cvCVM.delta.c.05<-quantile(VbCVM.delta.c,0.95)
cvCVM.delta.c.10<-quantile(VbCVM.delta.c,0.90)

# r(\tau)=rhat
#Factual
cvKS.gic.f.rh.01<-quantile(VbKS.gic.f.rh,0.99)
cvKS.gic.f.rh.05<-quantile(VbKS.gic.f.rh,0.95)
cvKS.gic.f.rh.10<-quantile(VbKS.gic.f.rh,0.90)
cvCVM.gic.f.rh.01<-quantile(VbCVM.gic.f.rh,0.99)
cvCVM.gic.f.rh.05<-quantile(VbCVM.gic.f.rh,0.95)
cvCVM.gic.f.rh.10<-quantile(VbCVM.gic.f.rh,0.90)

#Counter
cvKS.gic.c.rh.01<-quantile(VbKS.gic.c.rh,0.99)
cvKS.gic.c.rh.05<-quantile(VbKS.gic.c.rh,0.95)
cvKS.gic.c.rh.10<-quantile(VbKS.gic.c.rh,0.90)
cvCVM.gic.c.rh.01<-quantile(VbCVM.gic.c.rh,0.99)
cvCVM.gic.c.rh.05<-quantile(VbCVM.gic.c.rh,0.95)
cvCVM.gic.c.rh.10<-quantile(VbCVM.gic.c.rh,0.90)

cvKS.gic.gic.rh.01<-quantile(VbKS.gic.gic.rh,0.99)
cvKS.gic.gic.rh.05<-quantile(VbKS.gic.gic.rh,0.95)
cvKS.gic.gic.rh.10<-quantile(VbKS.gic.gic.rh,0.90)
cvCVM.gic.gic.rh.01<-quantile(VbCVM.gic.gic.rh,0.99)
cvCVM.gic.gic.rh.05<-quantile(VbCVM.gic.gic.rh,0.95)
cvCVM.gic.gic.rh.10<-quantile(VbCVM.gic.gic.rh,0.90)


#Compute the statistic of test over the entire sample
# r(\tau)=0

VfKS.gic.f<-max(abs(gic.f.aux-0)) #KS - No SQRT
VfCVM.gic.f<-sum(abs(gic.f.aux-0)) #CVM - No SQRT
VfKS.gic.c<-max(abs(gic.c.aux-0)) #KS - No SQRT
VfCVM.gic.c<-sum(abs(gic.c.aux-0)) #CVM - No SQRT

VfKS.delta.f<-max(abs(delta.loc.f-0)) #KS - No SQRT
VfCVM.delta.f<-sum(abs(delta.loc.f-0)) #CVM - No SQRT
VfKS.delta.c<-max(abs(delta.loc.c-0)) #KS - No SQRT
VfCVM.delta.c<-sum(abs(delta.loc.c-0)) #CVM - No SQRT

# r(\tau)=rhat

VfKS.f.rh<-max(abs(gic.f.aux-gic.avg.f)) #KS - No SQRT
VfCVM.f.rh<-sum(abs(gic.f.aux-gic.avg.f)) #CVM - No SQRT
VfKS.c.rh<-max(abs(gic.c.aux-gic.avg.c)) #KS - No SQRT
VfCVM.c.rh<-sum(abs(gic.c.aux-gic.avg.c)) #CVM - No SQRT

VfKS.gic.gic.rh<-max(abs(gic.f.aux-gic.avg.f-(gic.c.aux-gic.avg.c))) #KS - No SQRT
VfCVM.gic.gic.rh<-sum(abs(gic.f.aux-gic.avg.f-(gic.c.aux-gic.avg.c))) #CVM - No SQRT


#Final results
# Coeffs
resultcoeffs<-cbind(lvl,gic.f,var.gic.f,gic.c,var.gic.c,delta.loc.f,var.delta.f, delta.loc.c,var.delta.c)

# r(\tau)=0
resultKS.gic.f<-cbind(VfKS.gic.f,cvKS.gic.f.01,cvKS.gic.f.05,cvKS.gic.f.10)
resultKS.gic.c<-cbind(VfKS.gic.c,cvKS.gic.c.01,cvKS.gic.c.05,cvKS.gic.c.10)
resultCVM.gic.f<-cbind(VfCVM.gic.f,cvCVM.gic.f.01,cvCVM.gic.f.05,cvCVM.gic.f.10)
resultCVM.gic.c<-cbind(VfCVM.gic.c,cvCVM.gic.c.01,cvCVM.gic.c.05,cvCVM.gic.c.10)

resultKS.delta.f<-cbind(VfKS.delta.f,cvKS.delta.f.01,cvKS.delta.f.05,cvKS.delta.f.10)
resultKS.delta.c<-cbind(VfKS.delta.c,cvKS.delta.c.01,cvKS.delta.c.05,cvKS.delta.c.10)
resultCVM.delta.f<-cbind(VfCVM.delta.f,cvCVM.delta.f.01,cvCVM.delta.f.05,cvCVM.delta.f.10)
resultCVM.delta.c<-cbind(VfCVM.delta.c,cvCVM.delta.c.01,cvCVM.delta.c.05,cvCVM.delta.c.10)

# r(\tau)=rhat

resultKS.gic.f.rh<-cbind(VfKS.f.rh,cvKS.gic.f.rh.01,cvKS.gic.f.rh.05,cvKS.gic.f.rh.10)
resultKS.gic.c.rh<-cbind(VfKS.c.rh,cvKS.gic.c.rh.01,cvKS.gic.c.rh.05,cvKS.gic.c.rh.10)
resultCVM.gic.f.rh<-cbind(VfCVM.f.rh,cvCVM.gic.f.rh.01,cvCVM.gic.f.rh.05,cvCVM.gic.f.rh.10)
resultCVM.gic.c.rh<-cbind(VfCVM.c.rh,cvCVM.gic.c.rh.01,cvCVM.gic.c.rh.05,cvCVM.gic.c.rh.10)

resultKS.gic.gic.rh<-cbind(VfKS.gic.gic.rh,cvKS.gic.gic.rh.01,cvKS.gic.gic.rh.05,cvKS.gic.gic.rh.10)
resultCVM.gic.gic.rh<-cbind(VfCVM.gic.gic.rh,cvCVM.gic.gic.rh.01,cvCVM.gic.gic.rh.05,cvCVM.gic.gic.rh.10)


