R/ellipticalTestMult_v2.R

In MiriamJaser/elliptest:

### 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")