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)))
}