### Function to compute the test statistic for a given sample ###
#' Elliptical test for multivariate data
#'
#' @param z n times d matrix
#' @return The values of the test statistic
ellipticalTestMult <- function(z) {
n <- dim(z)[1]
d <- dim(z)[2]
FIRST <- c()
SECOND <- c()
for(i in 1:(d-1))
{
for(j in (i+1):d)
{
FIRST <- cbind(FIRST, z[, i])
SECOND <- cbind(SECOND, z[, j])
}
}
dd <- dim(FIRST)[2]
A <- matrix(nrow=n,ncol=2*dd)
TT <- numeric(dd)
for(i in 1:dd)
{
X <- FIRST[,i]
Y <- SECOND[,i]
# Estimator for Blomquvist's beta
#hatBeta <- mean( sign(X-median(X))*sign(Y-median(Y)) )
hatBeta <- mean( sign(X-0.5)*sign(Y-0.5) ) # for copula data the median of the margins is 0.5
# Estimator for Kendall's tau [Def. 3.5]
hatTau <- cor(x=X, y=Y, method="kendall")
# The test statistic
T <- hatBeta - hatTau
TT[i] <- T
# Vector of empirical cdf of (X,Y) evaluated at (X_i, Y_i)
hatF_xy <- numeric(n)
for(j in 1:n) {
hatF_xy[j] <- sum((X<=X[j]) & (Y<=Y[j])) / n
}
# Hoeffding projection according to [Lemma 3.10]
hatPhi <- 1 - 2 * X - 2 * Y + 4 * hatF_xy
# The following computes the estimated covariance matrix V
A[,i] <- sign(X-0.5)*sign(Y-0.5)-hatBeta
A[,dd+i] <- 2*hatPhi-mean(2*hatPhi)
}
V <- t(A) %*% A / n # Asymptoric Variance from the CLT
GRAD <- cbind(diag(dd), -diag(dd))
VV <- GRAD %*% V %*% t(GRAD)
out <- n*t(TT) %*% solve(VV) %*% (TT)
return(out)
}
# Test the function
library(CDVine)
d <- 4
dd <- d*(d-1)/2
z <- CDVineSim(100, family=rep(5,dd), par=rep(10,dd), type=2)
ellipticalTestMult(z)
# save the function to load in other files
#save("ellipticalTestMult", file="ellipticalTestMult.Rdata")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.