R/statiotest_function.R

In D-Roberts/statiotest: A Second-Order Cumulant Spectrum Test for Strict Stationarity

#' Statiotest - Test for Strict Stationarity of a Time Series.
#'
#' The function tests the null hypothesis of strict stationarity of
#' the numerical input vector using a second order cumulant spectrum
#' test, statiotest. The alternatives can be violation via unit root and
#' violations of time invariant second moment of the time series.
#'
#' Reference: Hinich M., Patterson D., and Roberts D., Statiotest. A Second
#' Order Cumulant Spectrum Test for Strict Stationarity in the Frequency Domain,
#' to be published.
#'
#' @param x A vector.
#' @param L A number.
#' @param idemean A boolean value indicating if the mean is removed.
#' @param idetrend A boolean value indicating if the trend is removed.
#' @param itrim An integer value for percentage each tail trimming.
#' @return A list of statistics, including a p-value.
#' statistic: a vector of values that are tested for normality using the
#' Kolmogorov-Smirnoff test. If normal, then fail to reject null. The
#' elements are the real and imaginary part of the complex valued test
#' statistic calculated using the cumulant spectrum of the data.
#'
#' p.value:The p-value returned by the K-S test.
#'
#' alternative: "non-stationary"
#'
#' method: "statiotest"
#'
#' data.name: Name of the input vector.
#'
#' data: Returns a copy of the input vector after preprocessing by demeaning,
#' detrending or trimming.
#'
#' orig_data: Returns a copy of the input vector.
#'
#' bifreqs: Returns the tuple of frequencies in the principal domain
#' where the cumulant is
#' calculated based on the sample data and then used to calculate the
#' test statistic.
#'
#'desc_orig: Returns the numerical summaries of the original data.

#' @examples
#' statiotest(rnorm(1000))
#'
#' @export
statiotest <- function(x = rnorm(1000),
                      L = round(sqrt(length(x))),
                      idemean = FALSE,
                      idetrend = FALSE,
                      itrim = 0) {


  # depending on preprocessing, N
  # might be different from the initial length

  DNAME <- deparse(substitute(x))

  x2 <- as.vector(x, mode = "double")

  if (NCOL(x2) > 1)
    stop("x is not a vector or univariate time series")


  x1 <- preprocess(x = x2,
                  idemean = idemean, idetrend = idetrend, itrim = itrim)
  N <- length(x1)

  if (N < 20)
    stop("Sample size should be at least 20.")

  if (L < 4)
    stop("Frame length should be at least 4.")

  ##################################End data pre-process
  il <- 1
  iu <- floor(L / 2)
  P <- floor(N / L)

  if (P < 4)
    stop("Number of frames should be at least 4.")

  # the (k1,k2) and (k1, -k2) pairs
  tup2 <- bifreqs(L, N)

  #number of bifrequencies in the cum2 principal domain
  kc <- dim(tup2)[1]

  #half kc to compute cm for k2 <0 and k2 > 0 sepparately
  hkc <- kc / 2

  #Other counters initialization

  ks <- il - 1
  jf <- iu - il + 1
  i <- 1:L

  k <- il:iu
  jk <- k - il + 1

  kh <- 1:kc

  x <- rep(0, N)

  #Initialize arrays
  sp <- rep(0, jf)
  cm <- matrix(rep(0,2*kc), nrow = kc, ncol = 2)

  #Start loop over frames

  for (kb in 1 : P){

    # fast fourier transform of each frame
    f <- fft( x1[ i + {kb - 1} * L ] )[ 2 : L ]

    #######################


    #Periodogram(sp) and Cum2 il < k1=k2 <= iu

    #vectorized sp

    sp[ jk ] <- sp[ jk ] +  Re( (f[ k ] * f[ L - k ]) )

    ##########compute Cum2

    #cm contains both the imaginary and the real parts

    # for k2 > 0
    cm[1:hkc, 1 ] <- cm[1:hkc, 1] +
      Re(f[tup2[1:hkc, 1]] * f[tup2[1:hkc, 2]])
    cm[1:hkc, 2] <- cm[1:hkc, 2] +
      Im(f[tup2[1:hkc, 1]] * f[tup2[1:hkc, 2 ]])

    # for k2 < 0
    cm[(hkc + 1):kc, 1] <- cm[(hkc + 1):kc, 1] +
      Re(f[tup2[(hkc + 1):kc, 1]] * f[L - tup2[(hkc + 1):kc, 2]])
    cm[(hkc + 1): kc, 2] <- cm[ (hkc + 1): kc, 2] +
      Im(f[tup2[ (hkc + 1):kc, 1]] * f[L - tup2[(hkc + 1):kc, 2]])

  } #end loop over P


  #normalization factor for cum2 in vectorized format

  h <- rep(0, kc)

  h <- sqrt(sp[tup2[ kh , 1 ] - ks] * sp[abs(tup2[kh, 2]) - ks])

  # Getting the Y-HIPARO pivotal quantity

  cm <- cm * sqrt(P) * sqrt(2)/h


  ##########################next testing for normality with ks

  y <- c(cm[ , 1 ], cm[ , 2 ])

  #ks.test is part of base R
  # test for normality of the imaginary and real parts (stacked)
  # of the Y test statistic
  # other tests for normality may be used

  kstest <- stats::ks.test(y, "pnorm", exact = F)

  desc1 <- descriptives(x1) # desc of preprocessed data
  desc2 <- descriptives(x2) # desc of preprocessed data

  structure(list(statistic = y,
                 p.value = kstest$p.value,
                 alternative = "NON-STATIONARY",
                 method = "statiotest",
                 data.name = DNAME,
                 data = x1,
                 orig_data = x2,
                 bifreqs = tup2,
                 desc_orig = desc1,
                 desc_processed = desc2,
                 framel = L),
            class = "htest")

} #end statiotest function