﻿##################################################################################################################### ## The test of high-dimensional correlations by (5) in Yata and Aoshima (2016). ## ## INPUT: ## Tecdm(X1,X2); p_1 by n matrix X1 and p_2 by n matrix X2 as X1=(x_{11},...,x_{1n}) and X2=(x_{21},...,x_{2n}), ## where n (\ge 4) is the sample size and p_i (\ge 1) is the dimension of Xi. ## ## ## OUTPUT: ## TestStatistics: The test statistic value for (5) (\hat{T}_n/\hat{\delta}). ## pvalue: The asymptotic p-value for testing (1) by (5) ## (1-\Phi(\hat{T}_n/\hat{\delta}), where \Phi(x) is the c.d.f. of N(0,1)). ## Tn: The estimator of \Delta (\hat{T}_n). ## delta: The estimator of \delta (\hat{delta}). ##################################################################################################################### Tecdm<-function(X1,X2) { n<-dim(X1)[2] n1<-ceiling(n/2) n2<-n-n1 u<-2*n1*n2/((n1-1)*(n2-1)*n*(n-1)) Y1<-rbind(X1,matrix(0,1,n)) Y2<-rbind(X2,matrix(0,1,n)) V1<-function(k,x) { if (floor(k/2)>=n1) { x[,(floor(k/2)-n1+1):floor(k/2)] } else { cbind(x[,1:floor(k/2)],x[,(floor(k/2)+n2+1):n]) } } V2<-function(k,x) { if (floor(k/2)<=n1) { x[,(floor(k/2)+1):(floor(k/2)+n2)] } else { cbind(x[,1:(floor(k/2)-n1)],x[,(floor(k/2)+1):n]) } } H11<-function(k) { apply(V1(k,Y1),1, mean) } H21<-function(k) { apply(V1(k,Y2),1, mean) } H12<-function(k) { apply(V2(k,Y1),1, mean) } H22<-function(k) { apply(V2(k,Y2),1, mean) } S<-c(3:(2*n-1)) M11<-sapply(S,H11) M21<-sapply(S,H21) M12<-sapply(S,H12) M22<-sapply(S,H22) hatt<-function(i,j) { ((Y1[,i]-M11[,i+j-2])%*%(Y1[,j]-M12[,i+j-2])* (Y2[,i]-M21[,i+j-2])%*%(Y2[,j]-M22[,i+j-2])) } hatT<-u*sum(mapply(hatt,sequence(c(1:(n-1))),rep(2:n,1:(n-1)))) w1<-function(i,j) { ((Y1[,i]-M11[,i+j-2])%*%(Y1[,j]-M12[,i+j-2]))^2 } W1<-u*sum(mapply(w1,sequence(c(1:(n-1))),rep(2:n,1:(n-1)))) w2<-function(i,j) { ((Y2[,i]-M21[,i+j-2])%*%(Y2[,j]-M22[,i+j-2]))^2 } W2<-u*sum(mapply(w2,sequence(c(1:(n-1))),rep(2:n,1:(n-1)))) hatdelta<-(2*W1*W2)^{1/2}/n pval<-(1-pnorm((hatT/hatdelta))) return(list(TestStatistics=hatT/hatdelta,pvalue=pval,Tn=hatT,delta=hatdelta)) } ############################################################################################ ## Estimation of tr(\Sigma^2) by the ECDM methodology. ## ## INPUT: ## W(X); p by n matrix X as X=(x_{1},...,x_{n}), ## where n (\ge 4) is the sample size and p (\ge 1) is the dimension. ## ## OUTPUT: ## trSigmaSquare: The estimator of \tr(\Sigma^2) by the ECDM methodology. ############################################################################################ W<-function(X) { n<-dim(X)[2] n1<-ceiling(n/2) n2<-n-n1 u<-2*n1*n2/((n1-1)*(n2-1)*n*(n-1)) Y<-rbind(X,matrix(0,1,n)) V1<-function(k,x) { if (floor(k/2)>=n1) { x[,(floor(k/2)-n1+1):floor(k/2)] } else { cbind(x[,1:floor(k/2)],x[,(floor(k/2)+n2+1):n]) } } V2<-function(k,x) { if (floor(k/2)<=n1) { x[,(floor(k/2)+1):(floor(k/2)+n2)] } else { cbind(x[,1:(floor(k/2)-n1)],x[,(floor(k/2)+1):n]) } } H1<-function(k) { apply(V1(k,Y),1, mean) } H2<-function(k) { apply(V2(k,Y),1, mean) } S<-c(3:(2*n-1)) M1<-sapply(S,H1) M2<-sapply(S,H2) q<-function(i,j) { ((Y[,i]-M1[,i+j-2])%*%(Y[,j]-M2[,i+j-2]))^2 } Q<-u*sum(mapply(q,sequence(c(1:(n-1))),rep(2:n,1:(n-1)))) return(list(trSigmaSquare=Q)) }