R/T_stationary.R

In ftsa: Functional Time Series Analysis

T_stationary <- function(sample, L = 49, J = 500, MC_rep=1000, cumulative_var = .90, Ker1 = FALSE, Ker2 = TRUE, h = ncol(sample)^.5, pivotal=FALSE, use_table = FALSE, significance = c("10%", "5%", "1%")){
  xrefine = N = ncol(sample)
  trefine = nrow(sample)
  if(Ker1)
  {
    K=function(x)
    {
      output = min(1, max(1.1-abs(x),0))
      return(output)
    }
  }
  if(Ker2)
  {
    K=function(x)
    {
      output = min(1, max(2-2*abs(x),0))
      return(output)
    }
  }
  
  
  basis = create.fourier.basis(c(0,1),L)
  ld = length(cumulative_var)
  h1 = h
  
  if(use_table & (ld == 1)){
    possible_levels = c("10%", "5%", "1%")
    match.arg(significance, possible_levels)
    significance = match(significance, possible_levels)
    text_level = possible_levels[significance]
    
    critical_values = generate_critical_values()
  } else if(ld > 1){
    stop("Table procedure is only available for one cumulative variance value")
  }
  
  
  
  X1_bar = rowMeans(sample)
  mean_subtracted_X1 = sample - X1_bar
  gamma_hat = lapply(0:(N-1), gamma_matrix, process = mean_subtracted_X1)
  cov_sample1 = gamma_hat[[1]]
  for(index in 1:(N-1))
  {
    cov_sample1 = cov_sample1 + K(index/h1) *(gamma_hat[[index+1]] + t(gamma_hat[[index+1]]))
  }
  Z_matrix = cov_sample1
  
  e1 = list()
  for(index in 1:L)
  {
    e1[index] = list(as.matrix(eval.basis(evalarg =(1:trefine)/trefine, basisobj = basis,Lfdobj=0)[,index]))
  }
  eigenvalues1 = (eigen(Z_matrix)$values)/trefine
  D = matrix(0,L,L)
  for(k in 1:L)
  {
    for(ell in 1:L)
    {
      Integrand_matrix = Z_matrix * (e1[[k]] %*% t(e1[[ell]]))
      D[k,ell] = 1/(trefine^2)*sum(Integrand_matrix)
    }
  }
  eigenpairs = eigen(D)
  eigenvectors = eigenpairs$vec
  eigenvectors2 = eigen(Z_matrix)$vectors
  eigenvalues = eigenpairs$val
  evals = eigenvalues
  if(pivotal)
  {
    d = c(1:ld)
    switch = 0
    stoper = 1
    spot = 1
    while(switch==0)
    {
      while((sum(eigenvalues[c(1:spot)])/sum(eigenvalues)) < cumulative_var[stoper])
      {
        spot = spot+1
      }
      d[stoper] = spot
      stoper = stoper+1
      if(stoper == (length(d)+1))
      {
        switch = 1
      }
    }
    T_N0=1:ld
    for(r in 1:ld)
    {
      ds=d[r]
      inp.matrix=matrix(0,ds,N)
      eig.v.norm=((trefine)^.5)*eigenvectors2
      for(j in (1:ds))
      {
        for(k in (1:N))
        {
          inp.matrix[j,k]=t(sample[,k])%*%(eig.v.norm[,j])/trefine
        }
      }
      T_Nsum=rep(0,ds)
      for(j in (1:ds))
      {
        s.0=sum(inp.matrix[j,(1:xrefine)])
        for(x in (1:xrefine))
        {
          T_Nsum[j]=T_Nsum[j]+(1/xrefine)*((1/N^.5)*(sum(inp.matrix[j,(1:x)])-(x/xrefine)*s.0))^2
        }
      }
      T_N0[r]=sum(T_Nsum/eigenvalues[1:ds])
    }
    
    if(use_table){
      
      if(d>30){
        warning("The table only allows for 30 FPC, the conclusion might be wrong")
        d = 30
      }
      T_N0 = round(T_N0, 6)
      is_rejected = 1*(T_N0 > critical_values[significance, d])
      if(is_rejected){
        conclusion = 'Reject the null hypothesis\n'
      } else {
        conclusion = "Do not reject the null hypothesis\n"
      }
      
      cat("\n")
      cat("Pivotal test of stationarity for a functional time series\n")
      cat("\n")
      cat("null hypothesis: the series is stationary\n")
      cat("\n")
      
      cat(paste("Critical Value at ", text_level, " = ", critical_values[significance, d], sep=""),"\n")
      cat(paste('Test Statistic:', T_N0, '\n'))
      cat(paste("N (number of functions) = ", N, sep=""),"\n")
      cat(paste("Conclusion:", conclusion))
      
      cat("\n")
      
    } else {
      
      T = vector(, MC_rep)
      T_array = matrix(0, length(d), MC_rep)
      lambda = eigenvalues
      for(dd in 1:length(d))
      {
        for(k in 1:MC_rep)
        {
          z=rnorm(d[dd]*J)
          tot=0
          for(n in c(1:d[dd]))
          {
            sum1 = 0
            sum1 =sum((z[c(((n-1)*d[dd]+1):((n-1)*d[dd]+J))]/(pi*c(1:J)))^2)
            tot = tot+sum1
          }
          T_array[dd,k] = T[k] = tot
        }
      }
      p_values=(1:ld)
      for(dd in 1:length(d))
      {
        p_values[dd] = round(1 - ecdf(T_array[dd,])(T_N0[dd]),4)
      }
      cat("\n")
      cat("Pivotal test of stationarity for a functional time series\n")
      cat("\n")
      cat("null hypothesis: the series is stationary\n")
      cat("\n")
      
      cat(paste("p-value = ", p_values, sep=""),"\n")
      cat(paste("N (number of functions) = ", N, sep=""),"\n")
      cat(paste("number of MC replications = ", MC_rep, sep=""))
      cat("\n")
    }
  } else {
    T = vector(, MC_rep)
    T_array = (1: MC_rep)
    lambda = eigenvalues
    d = min(c(length(which(lambda>0)),15))
    for(k in 1:MC_rep)
    {
      z=rnorm(d*J)
      tot=0
      for(n in c(1:d))
      {
        sum1 = 0
        sum1 =sum((z[c(((n-1)*d+1):((n-1)*d+J))]/(pi*c(1:J)))^2)
        tot = tot+lambda[n]*sum1
      }
      T_array[k] = T[k] = tot
    }
    int = sum(((1/sqrt(N)) *((sample[,1])-(1/xrefine)*rowSums(sample)))^2/(xrefine * trefine))
    for(x in 2:xrefine)
    {
      int = int + sum(((1/sqrt(N)) * (rowSums(sample[,1:x]) -(x/xrefine) * rowSums(sample)))^2/(xrefine * trefine))
    }
    statT_N = int
    p_values=(1:ld)
    for(dd in 1:length(d))
    {
      p_values = round(1 - ecdf(T_array)(statT_N),4)
    }
    cat("\n")
    cat("Monte Carlo test of stationarity of a functional time series\n")
    cat("\n")
    cat("null hypothesis: the series is stationary\n")
    cat("\n")
    
    cat(paste("p-value = ", p_values, sep=""),"\n")
    cat(paste("N (number of functions) = ", N, sep=""),"\n")
    cat(paste("number of MC replications = ", MC_rep, sep=""))
    cat("\n")
  }
}