setwd("P:/Arbeitsverzeichnis/Promotion/R/Projekt_1/Multivariat")
library(CDVine)
library(mvtnorm)
# function to simulate copula data from a Gaussian copula
gausscop_sim <- function(n, d, rho) {
sigma <- matrix(rho, ncol=d, nrow=d)
diag(sigma) <- 1
x <- rmvnorm(n, mean = rep(0,d), sigma = sigma)
u <- pnorm(x)
#out <- list("u"=u, "x"=x)
return(u)
}
# test the function
dat <- gausscop_sim(n=1000, d=6, rho=0.7)
dim(dat)
head(dat)
### INPUT
# copula family
fam <- 1 # 1=Gaussian, 2=t, 3=Clayton, 4=Gumbel, 5=Frank
fam_name <- "gauss"
# elliptical or not
elliptical <- 1 # 1=elliptical, 0=not elliptical
# second parameter for copula (if necessary)
p2 <- 0
# 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 <- gausscop_sim(100,d,rho[1]) # 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 <- gausscop_sim(n[j], d=d, rho=rho[k])
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))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.