Nothing
#' 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)
}
surrogate.data
}
#' Surrogate data testing
#' @details
#' This function tests the null hypothesis (H0) stating that the series
#' is a gaussian 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.
#'
#' To test with a significance level of \eqn{\alpha}{alpha} if the statistic
#' from the original data is smaller than the expected value under the null
#' hypothesis (a one-side test), \eqn{K/\alpha - 1}{K/alpha - 1} surrogates
#' are generated. The null hypothesis is then rejected if the statistic from
#' the data has one of the K smallest values. For a two-sided test,
#' \eqn{2K/\alpha - 1}{2K/alpha - 1} surrogates are generated. The null
#' hypothesis is rejected if the statistic from the data gives one of the K
#' smallest or largest values.
#'
#' 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 one.sided Logical value. If \emph{TRUE}, the routine runs a one-side
#' test. If \emph{FALSE}, a two-side test is applied (default).
#' @param alternative Specifies the concrete type of one-side test that should be
#' performed: If the the user wants to test if the statistic from the original
#' data is smaller (\emph{alternative="smaller"}) or larger
#' (\emph{alternative="larger"}) than the expected value under the
#' null hypothesis.
#' @param K Integer controlling the number of surrogates to be generated (see
#' details).
#' @param FUN The function that computes the discriminating statistic that shall
#' be used for testing.
#' @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 main an overall title for the plot.
#' @param xlab a title for the x axis.
#' @param ylab a title for the y axis.
#' @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
#' @examples
#' \dontrun{
#' lx = lorenz(do.plot=F)$x
#' sdt = surrogateTest(time.series = lx,significance = 0.05,
#' K=5, one.sided = FALSE, FUN=timeAsymmetry)
#' }
surrogateTest = function(time.series, significance = 0.05, one.sided=FALSE,
alternative = c("smaller", "larger"), K = 1, FUN,
verbose = TRUE, do.plot = TRUE,
xlab = "Values of the statistic", ylab = "",
main = "Surrogate data testing",
...) {
ktol = 1e-6
if (K < 1) {
stop("Invalid K (K < 1 not allowed")
}
# make sure that K is an integer
if (abs(K %% 1) > ktol) {
warning("K is not an integer... rounding it!")
K = round(K)
}
if (one.sided) {
alternative = match.arg(alternative)
n.surrogates = ceiling(K / significance - 1 )
} else {
n.surrogates = ceiling((2*K) / significance - 1 )
}
surrogates = FFTsurrogate(time.series, n.samples = n.surrogates)
nltest = list()
if (verbose) {
message("Computing statistics\n")
}
surrogates.statistics = apply(surrogates, MARGIN = 1, FUN = FUN, ...)
FUN = match.fun(FUN)
data.statistic = FUN(time.series,...)
if (verbose) {
if (one.sided) {
sdtMessage1s(data.statistic, surrogates.statistics, alternative, K)
}else{
sdtMessage2s(data.statistic, surrogates.statistics, K)
}
}
nltest$surrogates.statistics = surrogates.statistics
nltest$data.statistic = data.statistic
class(nltest) = "surrogateTest"
if (do.plot) {
plot(nltest, main = main, xlab = xlab, ylab = ylab)
}
nltest
}
#' @export
plot.surrogateTest = function(x, xlab = "Values of the statistic", ylab = "",
main="Surrogate data testing",
type = "h", lwd = 2,
xlim = NULL, ylim = NULL,
add.legend = T, ...) {
surrogates.statistics = x$surrogates.statistics
data.statistic = x$data.statistic
if (is.null(xlim)) {
xlim = range(data.statistic, surrogates.statistics)
}
if (is.null(ylim)) {
ylim = c(0, 3)
}
plot(surrogates.statistics, rep(1, length(surrogates.statistics)),
xlim = xlim, ylim = ylim,
type = type, lwd = lwd, main = main, xlab = xlab, ylab = ylab,
col = 1, lty = 1)
lines(data.statistic,2, type = type, lwd = lwd, col = 2, lty = 2)
if (add.legend) {
legend("top",bty = "n", legend = c("Surrogate data", "Original data"),
col = 1:2, lty = c(1, 2), lwd = lwd)
}
}
# Surrogate Data Testing two-sided Message
sdtMessage2s = function(data.statistic, surrogates.statistics,K){
sorted.values = sort( c(data.statistic,surrogates.statistics) )
# threshold for being significant smaller (being among the K smallest)
min.st = sorted.values[K]
# threshold for being significant larger (being among the K largest)
len = length(sorted.values)
max.st = sorted.values[[len - K + 1]]
message("Null Hypothesis: Data comes from a linear stochastic process")
if (data.statistic <= min.st) {
message(paste("Reject Null hypothesis:\n",
"\tOriginal data's stat is significant smaller",
"than surrogates' stats"))
} else {
if (data.statistic >= max.st) {
message(paste("Reject Null hypothesis:\n",
"\tOriginal data's stat is significant larger",
"than surrogates' stats"))
} else {
message("Accept Null hypothesis")
}
}
}
# Surrogate one sided Data Testing
sdtMessage1s = function(data.statistic, surrogates.statistics, alternative, K) {
sorted.values = sort(c(data.statistic, surrogates.statistics))
# threshold for being significant smaller (being among the K smallest)
min.st = sorted.values[K]
# threshold for being significant larger (being among the K largest)
len = length(sorted.values)
max.st = sorted.values[[len - K + 1]]
message("Null Hypothesis: Original data comes from a linear stochastic process")
if ((alternative == "smaller") && (data.statistic <= min.st)) {
message(paste("Reject Null hypothesis:\n",
"\tOriginal data's stat is significant smaller",
"than surrogates' stats"))
} else {
if ((alternative == "larger") && (data.statistic >= max.st)) {
message(paste("Reject Null hypothesis:\n",
"\tOriginal data's stat is significant larger",
"than surrogates' stats"))
}else{
message("Accept Null hypothesis")
}
}
}
# 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
Re(fft(fourier.transform, inverse = "TRUE")) / length.data
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.