R/surrogate.data.R

In nonlinearAnalysis: R package for nonlinear time series analysis

#' Generate surrogate data using the Fourier transform
#' @description
#' Generates surrogate samples from the original time series.
#' @details
#' This function uses the phase randomization procedure for 
#' generating the surrogated data. This algorithm generates surrogate data with the same
#' mean and autocorrelation function (and thus, the same power spectrum because of the 
#' Wiener-Khinchin theorem) as the original time series.
#' 
#' The phase randomization algorithm is often used when the null hypothesis being tested
#' consist on the assumption that the time.series data comes from a stationary linear
#' stochastic process with Gaussian inputs. The phase randomization preserves the Gaussian
#' distribution.
#' @param time.series The original time.series from which the surrogate data is generated.
#' @param n.samples The number of surrogate data sets to generate,
#' @return A matrix containing the generated surrogate data (one time series per row).
#' @references H. Kantz  and T. Schreiber: Nonlinear Time series Analysis (Cambridge university press)
#' @author Constantino A. Garcia
#' @examples
#' \dontrun{
#' # generate 20 surrogate sets using as original time series
#' # an arma(1,1) simulation
#' time.series = arima.sim(list(order = c(1,0,1), ar = 0.6, ma = 0.5), n = 200)
#' surrogate = FFTsurrogate(time.series, 20)
#' }
#' @export FFTsurrogate
FFTsurrogate = function(time.series,n.samples=1){
  fourier.transform = fft(time.series)
  length.data = length(time.series)
  surrogate.data=matrix(0, nrow = n.samples, ncol=length.data)
  for (counter in 1:n.samples){
    surrogate.data[counter,] = generateSingleFFTsurr(fourier.transform, length.data)
  }
  return(surrogate.data)
  
}

#' Surrogate data testing 
#' @details
#' This function tests the null hypothesis (H0) stating that the series describes a linear process. The
#' test is performed by generating several surrogate data according to H0 and comparing the values
#' of a discriminating statistic between both original data and the surrogate data. If the value of 
#' the statistic is significantly different for the original series than for the surrogate set,
#'  the null hypothesis is rejected and nonlinearity assumed.  The  surrogate data is generated by using a phase randomization procedure.
#' @param time.series The original time.series from which the surrogate data is generated.
#' @param significance Significance of the test
#' @param verbose Logical value. If TRUE, a brief summary of the test is shown.
#' @param do.plot Logical value. If TRUE, a graphical representation of the statistic value for
#' both surrogates and original data is shown.
#' @param FUN The function that computes the discriminating statistic that shall be used for testing.
#' @param ... Additional arguments for the FUN function.
#' @return  A list containing the values of the statistic for the surrogates (\emph{surrogates.statistics}
#' field) and the value for the original time series (\emph{data.statistic})
#' @references SCHREIBER, Thomas; SCHMITZ, Andreas. Surrogate time series. Physica D: 
#' Nonlinear Phenomena, 2000, vol. 142, no 3, p. 346-382.
#' @author Constantino A. Garcia
#' @export surrogateTest
surrogateTest <- function(time.series, significance = 0.05,
                          verbose = TRUE, do.plot = TRUE,FUN, ...){
  
  surrogates = FFTsurrogate(time.series, n.samples = ceiling(2/significance - 1 ))
  nltest = list()
  if (verbose){
    cat("      Computing statistics\n")  
  }
  surrogates.statistics = apply(surrogates,MARGIN=1,FUN=FUN,...)  
  FUN = match.fun(FUN)
  data.statistic = FUN(time.series,...)
  
  if (do.plot){
    plot(surrogates.statistics,rep(1,length(surrogates.statistics)),
         xlim=range(data.statistic,surrogates.statistics),ylim=c(0,3),
         type="h",col="blue",lty=1,xlab="Statistic",ylab="",main="Surrogate data testing")
    lines(data.statistic,1,type="h",col="red",lty=2)
    legend("center",legend=c("Surrogate data","Original data"),col=c("blue","red"),
           lty=c(1,2))
  }
  
  if (verbose){
    
    cat("      Null Hypothesis: Original data comes from a linear stochastic process\n")
    if (data.statistic < min(surrogates.statistics)){
      cat("      Reject Null hypothesis: Original data's statistic is the smallest one\n")
    }else{
      if (data.statistic > max(surrogates.statistics)){
        cat("      Reject Null hypothesis: Original data's statistic is the greatest one\n")
      }else{
        cat("      Accept Null hypothesis\n")
      }  
    }
  }
  nltest$surrogates.tatistics = surrogates.statistics
  nltest$data.statistic = data.statistic
  
  return(nltest)
}

# private function
# generates one surrogate data set using the fourier.transform Fourier transform of length length.data
generateSingleFFTsurr = function(fourier.transform, length.data){
  # we must respect the hermitian simmetry to generate real time series
  # phase[2] = - phase[N]; phase[3] = -phase[N-1]... phase[k] = -phase[N-k+2]
  
  # Differentiate the odd and the even case
  if (length.data%%2 == 0){ 
    #even case
    for (k in 1:(length.data/2-1) ){
      phase = runif(1,min=0,2*pi)
      fourier.transform[[k+1]]=abs(fourier.transform[[k+1]])*exp(1i*phase)
      fourier.transform[[length.data-k+1]]=abs(fourier.transform[[length.data-k+1]])*exp(-1i*phase)
    }
  }else{
    # odd case
    for (k in 1:( (length.data-1)/2 ) ){
      phase = runif(1,min=0,2*pi)
      fourier.transform[[k+1]]=abs(fourier.transform[[k+1]])*exp(1i*phase)
      fourier.transform[[length.data-k+1]]=abs(fourier.transform[[length.data-k+1]])*exp(-1i*phase)
    }
  }
  # give warning if some imaginary part is too big 
  return(Re(fft(fourier.transform,inverse="TRUE"))/length.data)
  
}