elliptest:

setwd("P:/Arbeitsverzeichnis/Promotion/R/Projekt_1/Multivariat")


library(CDVine)
library(copula)



### INPUT

# copula family
fam <- 5 # 1=Gaussian, 2=t, 3=Clayton, 4=Gumbel, 5=Frank
fam_name <- "frank"


# elliptical or not
elliptical <- 0 # 1=elliptical, 0=not elliptical


# second parameter for copula (if necessary)
p2 <- 0


# number of replications
nr <- 1000


# dimension
d <- 3
dd <- d*(d-1)/2


# sample size
n <- c(100, 250, 500, 1000, 5000)
#n <- c(50, 100, 250, 500, 1000, 5000, 10000)
#n <- c(250, 500)
#n <- 1000

# different values for Kendall's tau
ttau <- c(0.1, 0.25, 0.5, 0.75, 0.9)
#ttau <- c(0.25,0.5,0.75)
#ttau <- 0.25


# corresponding values for the correlation
if(elliptical==1){
  rho  <- sin(pi*ttau/2) # holds for elliptical copulas
}else{
  rho <- rep(NA,length(ttau))
  for(i in 1:length(ttau))
    rho[i]  <- BiCopTau2Par(family=fam, tau=ttau[i]) 
}

# for BB1 copula
if(fam==7){
  rho <- 2 / (p2*(1-ttau)) - 2 
}
# test if correct parameters for given ttau values
BiCopPar2Tau(fam=7, par=rho, par2=p2)


# function to compute values of the test statistic for a given sample
load("ellipticalTestMult.RData")


# test the function
cop_test <- archmCopula(fam_name, param=rho[5] , dim=6)
z <- rCopula(100, cop_test)
dim(z)
head(z)
ellipticalTestMult(z)




### Compute the values of the test statistic 
# for different values of Kendall's tau and the sample size


# List to store the results (values of the test statistic)
# List elements for different values of Kendall's tau/rho

test <- vector("list", length(ttau))

# lhs  <- paste0("test_", fam_name)
# rhs  <- paste0("vector(\"list\",length(ttau))")
# eq   <- paste(paste(lhs, rhs, sep="<-"), collapse=";")
# eval(parse(text=eq))


for(k in 1:length(ttau)){
  
  # array to store the results 
  # for nr replications (rows) and different sample sizes n (cols)
  test[[k]] <- array(NA, c(nr, length(n)))
  
  print(Sys.time())
  print(paste0("k: ",k, " of ", length(ttau)))
  
  # set the copula object
  cop <- archmCopula(fam_name, param=rho[k] , dim=d)
  
  #for(j in 1:2){
  for(j in 1:length(n)){
    
    print(n[j])
    
    for (i in 1:nr){
      
      z <-  rCopula(n[j], cop)
      test[[k]][i,j] <-  ellipticalTestMult(z)
      
    }
    
    emp_power <- sum(test[[k]][,j] >= qchisq(0.95,df=dd)) / nr
    print(paste0("empirical level/power: ", emp_power))
    qqplot(qchisq(ppoints(nr), df = dd),test[[k]][,j], sub=paste0("tau=", ttau[k], ", n=", n[j]))
    abline(0, 1, col=2, lwd=2)
    
  }
  
  print(Sys.time())
  
}




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


### determine the level 
# for different values of Kendall's tau (rows) and different sample sizes n (cols)

if(elliptical==1){
  
  level_alpha05 <- array(NA, c(length(ttau), length(n)))
  
  for(i in 1:length(ttau)){
    for(j in 1:length(n)){
      
      level_alpha05[i,j] <- sum(test[[i]][,j] >= qchisq(0.95,df=dd)) / nr
      
    }
  }
  
  level_alpha05
  
}else{
  
  power_alpha05 <- array(NA, c(length(ttau), length(n)))
  
  for(i in 1:length(ttau)){
    for(j in 1:length(n)){
      
      power_alpha05[i,j] <- sum(test[[i]][,j] >= qchisq(0.95,df=dd)) / nr
      
    }
  }
  
  power_alpha05
  
}




### check the qq-plots

for(i in 1:length(ttau)){
  for(j in 1:length(n)){
    
    if(sum(is.na(test[[i]][,j]))==0){
      qqplot(qchisq(ppoints(nr), df = dd),test[[i]][,j], sub=paste0("tau=", ttau[i], ", n=", n[j]))
      #qqnorm(test[[i]][,j], sub=paste0("tau=", ttau[i], ", n=", n[j]))
      abline(0, 1, col=2, lwd=2)
    }
  }
}




### save results

#eq <- paste0("save.image(\"erg_",fam_name,"_mult_",d,"_copula.RData\")")
#eval(parse(text=eq))