library(expm); library(MASS); library(vars); library(Matrix); ######################################################### # Function to vectorize a matrix without k diagonals# ######################################################### vec.k <- function(X,k){ n = dim(X)[1]; p = dim(X)[2]; Y = matrix(0,n,p*(p+1)/2); for(i in 1:{n-k}){ Y[i,] = (X[i,]%*%t(X[(i+k),]))[upper.tri(diag(p),diag = TRUE)] } return(Y); } ######################################## # Function performing the mean testing # # Y: raw data, nxp m:trimming parameter# ######################################## SN <- function(Y,m = 10){ n = dim(Y)[1]; p = dim(Y)[2]; tilde.n = n - m; theta = rep(0,tilde.n); theta[1] = sum(Y[(1+m),]*Y[1,])/p; for(i in 2:tilde.n){ theta[i] = theta[i-1]+ sum(Y[1:i,]%*%Y[i+m,])/p; } theta = theta/((1:tilde.n)*(2:(tilde.n+1))/2); W2 = sum((((1:tilde.n)*(2:(tilde.n+1))/2)*(theta - theta[tilde.n]))^2)/(tilde.n^2*(tilde.n+1)/2); Q =(tilde.n*(tilde.n+1)/2)*(theta[tilde.n])^2/(W2); if(Q > 54.70){return(1)} #For different significant level, change 54.70 to other critical values else{return(0)} } ############################################# # Function performing White Noise Testing # # M: different DGP # ############################################## wn <- function(n,p,M){ if(M == 1){ S <- matrix(0,p,p); S <- 0.995**abs(row(S)-col(S)); A <- sqrtm(S); A <- (A + t(A))/2; z <- mvrnorm(n,rep(0,p),diag(p)); x <- z %*% t(A); } if(M == 2){ A <- matrix(runif(p^2,-1,1),p,p); z <- mvrnorm(n,rep(0,p),diag(p)); x <- z %*% t(A); } if(M == 3){ A <- matrix(0,p,p); A[subset(row(A) <= 12 & subset(col(A) <= 12))] <- runif(144,-0.25,0.25) x <- matrix(0,n,p); z <- matrix(rt(n*p,8),n,p); x[1,] <- z[1,]; for(i in 2:n){ x[i,] <- A%*%x[i-1,]+z[i,] } } if(M == 4){ A <- matrix(runif(p^2,-1,1)*(runif(p^2,0,1)<1/3),p,p); diag(A) <- 0.8 Sigma <- matrix(0,n,n); Sigma <- 0.5*(abs(row(Sigma) - col(Sigma))^(-0.6)*(abs(row(Sigma) - col(Sigma)) <= 7 & abs(row(Sigma) - col(Sigma)) >= 1)) diag(Sigma) <- 1 z <- matrix(rnorm(n*p),n,p) z[,1:12] <- sqrtm(Sigma) %*% z[,1:12] x <- z %*% t(A) } if(M == 5){ A <- matrix(0,p,p); A <- 0.9^(abs(row(A)-col(A))) A <- 0.9*A/norm(A,type = "2") x <- matrix(0,n,p); z <- matrix(rt(n*p,8),n,p); x[1,] <- z[1,]; for(i in 2:n){ x[i,] <- A%*%x[i-1,]+z[i,] } } lag <- seq(2,10,2); result.SN <-NULL; for(k in lag){ Y <- vec.k(x,k); result.SN <- c(result.SN,SN(Y)) } #result.SN <- ifelse(sum(result.SN) >=1 , 1,0); result.WN <- wntest(t(segmentTS(x,k0 = 10)\$X),1000,k_max = 10, kk = lag)\$res; #result.WN <- ifelse(sum(result.WN) >= 1,1,0) return(c(result.SN,result.WN)) } ######################################## # Function for bandness testing # # k: different bandedness level # # k0: bandness level under H0 # ######################################## band <- function(n,p,k0,k,delta){ X <- matrix(0,n,p); gamma = 0; if(k == 1){gamma = 1}; if(k == 2){gamma = c(0.5,0.25)} if(k == 5){gamma = c(0.4,0.4,0.4,0.4,0.4)} X[1,] <- arima.sim(model=list(ma=gamma),p); for(i in 2:n){ X[i,] <- arima.sim(model=list(ma=gamma),p)+delta*X[i-1,] } Y <- vec.band(X,k0); return(c(adstat(X,k0) > 3.289, SN(Y))) }