#The dCov.R file must be sourced before doing the following examples. #Examples of Section 10.1 #Student copula require(mvtnorm) #for multivariate Gaussian and Student random number generator n <- 100 d <- 2 rw <- 0.5 #correlation within rb <- 0.3 #correlation between D <- (1-rw)*diag(d)+rw*matrix(1,ncol=d,nrow=d) E <- matrix(rb,ncol=d,nrow=d) R <- rbind(cbind(D,E,E),cbind(E,D,E),cbind(E,E,D)) #correlation matrix Z <- pt(rmvt(n, sigma = R, df = 2),df=2) #sample from the Student copula system.time(out.dCov <- dCov.test(Z,rep(d,3))) system.time(out.kh <- multIndepTest(Z,rep(d,3))) #Kojadinovic and Holmes system.time(out.csorgo <- csorgo.test(Z,rep(d,3))) #Csorgo, eqn (16) out.dCov out.kh out.csorgo #Frank copula require(copula) n <- 100 frank.cop <- frankCopula(param=2, dim = 3) Z <- rCopula(n, frank.cop) par(mfrow=c(1,2)) dCov.test(Z,c(1,1,1)) dCov.test(Z,c(1,1,1),method="hsic",cte=rep(.001,3)) par(mfrow=c(1,1)) #Example of Section 10.2 n <- 100 d1 <- 2 d2 <- 3 R1 <- (1-.5)*diag(d1)+.5*matrix(1,ncol=d1,nrow=d1) R2 <- (1-.3)*diag(d2)+.3*matrix(1,ncol=d2,nrow=d2) theta <- 0.4 X1 <- rmvnorm(n, sigma=R1) Z2 <- rmvnorm(n, sigma=R1) Z3 <- rmvnorm(n, sigma=R1) Z1 <-(1-theta)*abs(X1)+ theta*cbind(abs(X1[,1])*sign(Z2[,1]*Z3[,1]),abs(X1[,2])*sign(Z2[,1]*Z3[,1])) Z4 <- rmvnorm(n, sigma=R2) X5 <- rmvnorm(n, sigma=R2) Z5 <- theta*Z4+X5 Z <- cbind(Z1,Z2,Z3,Z4,Z5) par(mfrow=c(1,2)) dCov.test(Z,c(2,2,2,3,3)) dCov.test(Z,c(2,2,2,3,3),method="hsic", cte=rep(.0001,5)) par(mfrow=c(1,1)) #Example of Section 10.3 require(mAr) w <- c(0,0) C <- rbind(c(1,0.5),c(0.5,1)) theta <- 0.4 A <- rbind(c(0,theta),c(theta,0)) m <- 100 Y <- mAr.sim(w,A,C,N=m) par(mfrow=c(1,2)) dCovSerial.test(Y,3,dependogram=TRUE) dCovSerial.test(Y,3,method="hsic") par(mfrow=c(1,1)) csorgoSerial.test(Y,3) #Application of Section 11.1 require(MVN) #for the multinormality test of Henze-Zirkler Z <- read.table(file="/home/bilodeau/etudiants/Guetsop/JMRL/code/TEMPERAT.txt", header = TRUE, sep = "", skip = 4) hzTest(Z, cov = TRUE, qqplot = TRUE) #Henze-Zirkler test of multivariate normality Z.scaled <- scale(Z) par(mfrow=c(1,2)) dCov.test(Z.scaled, c(3,3,3,1,1), index=1, ord=3) dCov.test(Z.scaled, c(3,3,3,1,1), index=1, ord=3, method="hsic") par(mfrow=c(1,1)) #Application of Section 11.2 #Verification of 5% significance levels of various tests for white noise trivariate sequences of length m=200 with p=5 #This simulation will take several minutes because levels are evaluated with 1000 tests power.dcov.fisher <- 0 power.dcov.tippett <- 0 power.hsic.fisher <- 0 power.hsic.tippett <- 0 m <- 200 p <- 5 r <- 1000 for(i in 1:r){ Z <- matrix(rnorm(3*m), ncol=3) #trivariate white noise sequence of length m=200 out.dcov <- dCovSerial.test(Z, p, ord=2,dependogram=FALSE) power.dcov.fisher <- power.dcov.fisher +(out.dcov$Fisher.pvalue < .05) power.dcov.tippett <- power.dcov.tippett +(out.dcov$Tippett.pvalue < .05) out.hsic <- dCovSerial.test(Z, p, ord=2,method="hsic",dependogram=FALSE) power.hsic.fisher <- power.hsic.fisher +(out.hsic$Fisher.pvalue < .05) power.hsic.tippett <- power.hsic.tippett +(out.hsic$Tippett.pvalue < .05) if(i%%10==0) print(i) } power.dcov.fisher/r power.dcov.tippett/r power.hsic.fisher/r power.hsic.tippett/r # Application of Section 11.2 TSX <- read.table(file="/home/bilodeau/etudiants/Guetsop/JMRL/code/TSXComp.csv", header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE) nbrow.TSX <- nrow(TSX) rate.TSX <- (-1+TSX[2:nbrow.TSX,5]/TSX[1:(nbrow.TSX-1),5])*100 m <- length(rate.TSX) par(mfrow=c(2,1)) plot.ts(TSX[,5],main=expression("TSX"),ylab="Index") plot.ts(rate.TSX,main=expression("TSX"),ylab="DPIR") p <- 10 dCovSerial.test(rate.TSX,p,ord=2) n <- m-4 TSX1 <- rate.TSX[1:n] TSX2 <- rate.TSX[(1+4):m] cor.test(TSX1,TSX2,method="pearson") cor.test(TSX1,TSX2,method="kendall") require(energy) dcov.test(TSX1,TSX2,R=1000) # dcov test of package energy out <- dCovSerial.test(rate.TSX,5,ord=2) out$statistics[4]*n #same statistic as dcov of package energy x1 <- TSX1*(TSX1<=0) x2 <- TSX1*(TSX1>0) summary(lm(TSX2~x1+x2))