R/GoF_MC_simulation_Mult_mvtnorm_t.R

In MiriamJaser/elliptest:

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


library(CDVine)
library(mvtnorm)


# function to simulate copula data from a Gaussian copula

tcop_sim <- function(n, d, rho, nu) {
  
  sigma       <- matrix(rho, ncol=d, nrow=d)
  diag(sigma) <- 1
  
  x <- rmvt(n, delta=rep(0,d), sigma=sigma, df=nu)
  #cov(x) ### cov = rho*nu/(nu-2)
  
  u <- pnorm(x)
  
  return(u)
}

# test the function
dat <- tcop_sim(d=6, n=1000, rho=0.7, nu=5)
dim(dat)
head(dat)







### INPUT

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


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


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


# number of replications
nr <- 1000


# dimension
d <- 6
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
z <- tcop_sim(100,d,rho[1],p2) # Gaussian copula data
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))


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)))
  
  #for(j in 1:2){
  for(j in 1:length(n)){
    
    print(n[j])
    
    for (i in 1:nr){
      z <- tcop_sim(n[j], d=d, rho=rho[k], nu=p2)
      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,"_mvtnorm.RData\")")
#eval(parse(text=eq))