TimeseriesSurrogates.jl

TimeseriesSurrogates is a Julia package for generating surrogate timeseries. It is part of JuliaDynamics, a GitHub organization dedicated to creating high quality scientific software.

If you are new to this method of surrogate time series, feel free to read the What is a timeseries surrogate? page. Please note that timeseries surrogates should not be confused with surrogate models, such as those provided by Surrogates.jl.

Installation

TimeseriesSurrogates is a registered Julia package. To install the latest version, run the following code:

import Pkg; Pkg.add("TimeseriesSurrogates")

API

TimeseriesSurrogates.jl exports two main functions. Both of them dispatch on the chosen method, a subtype of Surrogate. It is recommended to standardize the signal before using these functions, i.e. subtract mean and divide by standard deviation.

TimeseriesSurrogates.surrogateFunction
surrogate(x, method::Surrogate [, rng]) → s

Create a single surrogate timeseries s from x based on the given method. If you want to generate multiple surrogates from x, you should use surrogenerator for better performance.

source
TimeseriesSurrogates.surrogeneratorFunction
surrogenerator(x, method::Surrogate [, rng]) → sg::SurrogateGenerator

Initialize a generator that creates surrogates of x on demand, based on the given method. This is efficient, because for most methods some things can be initialized and reused for every surrogate. Optionally you can provide an rng::AbstractRNG object that will control the random number generation and hence establish reproducibility of the generated surrogates. By default Random.default_rng() is used.

Notice that the generated surrogates overwrite, in-place, a common vector container. Use copy if you need to actually store multiple surrogates.

To generate a surrogate, call sg as a function with no arguments, e.g.:

sg = surrogenerator(x, method)
for i in 1:1000
    s = sg()
    # do stuff with s and or x
    result[i] = stuff
end
source

Surrogate methods

TimeseriesSurrogates.RandomShuffleType
RandomShuffle() <: Surrogate

A random constrained surrogate, generated by shifting values around.

Random shuffle surrogates preserve the mean, variance and amplitude distribution of the original signal. Properties not preserved are any temporal information, such as the power spectrum and hence linear correlations.

The null hypothesis this method can test for is whether the data are uncorrelated noise, possibly measured via a nonlinear function. Specifically, random shuffle surrogate can test the null hypothesis that the original signal is produced by independent and identically distributed random variables[^Theiler1991, ^Lancaster2018].

Beware: random shuffle surrogates do not cover the case of correlated noise[Lancaster2018].

source
TimeseriesSurrogates.BlockShuffleType
BlockShuffle(n::Int; shift = false)

A block shuffle surrogate constructed by dividing the time series into n blocks of roughly equal width at random indices (end blocks are wrapped around to the start of the time series).

If shift is true, then the input signal is circularly shifted by a random number of steps prior to picking blocks.

Block shuffle surrogates roughly preserve short-range temporal properties in the time series (e.g. correlations at lags less than the block length), but break any long-term dynamical information (e.g. correlations beyond the block length).

Hence, these surrogates can be used to test any null hypothesis aimed at comparing short-range dynamical properties versus long-range dynamical properties of the signal.

source
TimeseriesSurrogates.CycleShuffleType
CycleShuffle(n::Int = 7, σ = 0.5)

Cycle shuffled surrogates[Theiler1994] that identify successive local peaks in the data and shuffle the cycles in-between the peaks. Similar to BlockShuffle, but here the "blocks" are defined as follows:

  1. The timeseries is smoothened via convolution with a Gaussian (DSP.gaussian(n, σ)).
  2. Local maxima of the smoothened signal define the peaks, and thus the blocks in between them.
  3. The first and last index of timeseries can never be peaks and thus signals that should have peaks very close to start or end of the timeseries may not perform well. In addition, points before the first or after the last peak are never shuffled.
  4. The defined blocks are randomly shuffled as in BlockShuffle.

CSS are used to test the null hypothesis that the signal is generated by a periodic oscillator with no dynamical correlation between cycles, i.e. the evolution of cycles is not deterministic.

See also PseudoPeriodic.

source
TimeseriesSurrogates.CircShiftType
CircShift(n)

Surrogates that are circularly shifted versions of the original timeseries.

n can be an integer (the surrogate is the original time series shifted by n indices), or any vector of integers, which which means that each surrogate is shifted by an integer selected randomly among the entries in n.

source
TimeseriesSurrogates.RandomFourierType
RandomFourier(phases = true)

A surrogate that randomizes the Fourier components of the signal in some manner. If phases==true, the phases are randomized, otherwise the amplitudes are randomized. FT is an alias for RandomFourier.

Random Fourier phase surrogates[Theiler1991] preserve the autocorrelation function, or power spectrum, of the original signal. Random Fourier amplitude surrogates preserve the mean and autocorrelation function but do not preserve the variance of the original. Random amplitude surrogates are not common in the literature, but are provided for convenience.

Random phase surrogates can be used to test the null hypothesis that the original signal was produced by a linear Gaussian process [Theiler1991].

source
TimeseriesSurrogates.TFTDRandomFourierType
TFTD(phases::Bool = true, fϵ = 0.05)

The TFTDRandomFourier (or just TFTD for short) surrogate was proposed by Lucio et al. (2012)[Lucio2012] as a combination of truncated Fourier surrogates[Nakamura2006] (TFTS) and detrend-retrend surrogates.

The TFTD part of the name comes from the fact that it uses a combination of truncated Fourier transforms (TFT) and de-trending and re-trending (D) the time series before and after surrogate generation. Hence, it can be used to generate surrogates also from (strongly) nonstationary time series.

Implementation details

Here, a best-fit linear trend is removed/added from the signal prior to and after generating the random Fourier signal. In principle, any trend can be removed, but so far, we only provide the linear option.

See also: TFTDAAFT, TFTDIAAFT.

source
TimeseriesSurrogates.PartialRandomizationType
PartialRandomization(α = 0.5)

PartialRandomization surrogates[Ortega1998] are similar to RandomFourier phase surrogates, but during the phase randomization step, instead of drawing phases from [0, 2π], phases are drawn from [0, 2π]*α, where α ∈ [0, 1]. The authors refers to α as the "degree" of phase randomization, where α = 0 means 0 % randomization and α = 1 means 100 % randomization.

source
TimeseriesSurrogates.PartialRandomizationAAFTType
PartialRandomizationAAFT(α = 0.5)

PartialRandomizationAAFF surrogates are similar to PartialRandomization surrogates[Ortega1998], but adds a rescaling step, so that the surrogate has the same values as the original time series (analogous to the rescaling done for AAFT surrogates). Partial randomization surrogates have, to the package authors' knowledge, not been published in scientific literature.

source
TimeseriesSurrogates.TFTSType
TFTS(fϵ::Real)

A truncated Fourier transform surrogate[Nakamura2006] (TFTS).

TFTS surrogates are generated by leaving some frequencies untouched when performing the phase shuffling step (as opposed to randomizing all frequencies, like for RandomFourier surrogates).

These surrogates were designed to deal with data with irregular fluctuations superimposed over long term trends (by preserving low frequencies)[Nakamura2006]. Hence, TFTS surrogates can be used to test the null hypothesis that the signal is a stationary linear system generated the irregular fluctuations part of the signal[Nakamura2006].

Controlling the truncation of the spectrum

The truncation parameter fϵ ∈ [-1, 0) ∪ (0, 1] controls which parts of the spectrum are preserved.

  • If fϵ > 0, then indicates the ratio of high frequency domain to the entire frequency domain. For example, fϵ = 0.5 preserves 50% of the frequency domain (randomizing the higher frequencies, leaving low frequencies intact).
  • If fϵ < 0, then indicates ratio of low frequency domain to the entire frequency domain. For example, fϵ = -0.2 preserves 20% of the frequency domain (leaving higher frequencies intact, randomizing the lower frequencies).
  • If fϵ ± 1, then all frequencies are randomized. The method is then equivalent to RandomFourier.

The appropriate value of strongly depends on the data and time series length, and must be manually determined[Nakamura2006], for example by comparing periodograms for the time series and the surrogates.

source
TimeseriesSurrogates.TFTDType
TFTD(phases::Bool = true, fϵ = 0.05)

The TFTDRandomFourier (or just TFTD for short) surrogate was proposed by Lucio et al. (2012)[Lucio2012] as a combination of truncated Fourier surrogates[Nakamura2006] (TFTS) and detrend-retrend surrogates.

The TFTD part of the name comes from the fact that it uses a combination of truncated Fourier transforms (TFT) and de-trending and re-trending (D) the time series before and after surrogate generation. Hence, it can be used to generate surrogates also from (strongly) nonstationary time series.

Implementation details

Here, a best-fit linear trend is removed/added from the signal prior to and after generating the random Fourier signal. In principle, any trend can be removed, but so far, we only provide the linear option.

See also: TFTDAAFT, TFTDIAAFT.

source
TimeseriesSurrogates.TFTDIAAFTType
TFTDIAAFT(fϵ = 0.05; M::Int = 100, tol::Real = 1e-6, W::Int = 75)

TFTDIAAFT[Lucio2012] are similar to TFTDAAFT, but adds an iterative procedure to better match the periodograms of the surrogate and the original time series, analogously to how IAAFT improves upon AAFT.

fϵ ∈ (0, 1] is the fraction of the powerspectrum corresponding to the lowermost frequencies to be preserved. M is the maximum number of iterations. tol is the desired maximum relative tolerance between power spectra. W is the number of bins into which the periodograms are binned when comparing across iterations.

See also: TFTD, TFTDAAFT.

source
TimeseriesSurrogates.AAFTType
AAFT()

An amplitude-adjusted-fourier-transform (AAFT) surrogate[Theiler1991].

AAFT surrogates have the same linear correlation, or periodogram, and also preserves the amplitude distribution of the original data.

AAFT surrogates can be used to test the null hypothesis that the data come from a monotonic nonlinear transformation of a linear Gaussian process (also called integrated white noise)[Theiler1991].

source
TimeseriesSurrogates.TAAFTType
TAAFT(fϵ)

An truncated version of the amplitude-adjusted-fourier-transform surrogate[Theiler1991][Nakamura2006].

The truncation parameter and phase randomization procedure is identical to TFTS, but here an additional step of rescaling back to the original data is performed. This preserves the amplitude distribution of the original data.

source
TimeseriesSurrogates.IAAFTType
IAAFT(M = 100, tol = 1e-6, W = 75)

An iteratively adjusted amplitude-adjusted-fourier-transform surrogate[SchreiberSchmitz1996].

IAAFT surrogates have the same linear correlation, or periodogram, and also preserves the amplitude distribution of the original data, but are improved relative to AAFT through iterative adjustment (which runs for a maximum of M steps). During the iterative adjustment, the periodograms of the original signal and the surrogate are coarse-grained and the powers are averaged over W equal-width frequency bins. The iteration procedure ends when the relative deviation between the periodograms is less than tol (or when M is reached).

IAAFT, just as AAFT, can be used to test the null hypothesis that the data come from a monotonic nonlinear transformation of a linear Gaussian process.

source
TimeseriesSurrogates.AutoRegressiveType
AutoRegressive(n, method = LPCLevinson())

Autoregressive surrogates of order-n. The autoregressive coefficients φ are estimated using DSP.lpc(x, n, method), and thus see the documentation of DSP.jl for possible methods.

While these surrogates are obviously suited to test the null hypothesis whether the data are coming from a autoregressive process, the Fourier Transform-based surrogates are probably a better option. The current method is more like an explicit way to produce surrogates for the same hypothesis by fitting a model. It can be used as convient way to estimate autoregressive coefficients and automatically generate surrogates based on them.

The coefficients φ of the autoregressive fit can be found by doing

sg = surrogenerator(x, AutoRegressive(n))
φ = sg.init.φ
source
TimeseriesSurrogates.PseudoPeriodicType
PseudoPeriodic(d, τ, ρ, shift = true)

Create surrogates suitable for pseudo-periodic signals. They retain the periodic structure of the signal, while inter-cycle dynamics that are either deterministic or correlated noise are destroyed (for appropriate ρ choice). Therefore these surrogates are suitable to test the null hypothesis that the signal is a periodic orbit with uncorrelated noise[Small2001].

Arguments d, τ, ρ are as in the paper, the embedding dimension, delay time and noise radius. The method works by performing a delay coordinates ambedding via the library DynamicalSystems.jl. See its documentation for choosing appropriate values for d, τ. For ρ, we have implemented the method proposed in the paper in the function noiseradius.

The argument shift is not discussed in the paper. If shift=false we adjust the algorithm so that there is little phase shift between the periodic component of the original and surrogate data.

See also CycleShuffle.

source
TimeseriesSurrogates.PseudoPeriodicTwinType
PseudoPeriodicTwin(d::Int, τ::Int, δ = 0.2, ρ = 0.1, metric = Euclidean())
PseudoPeriodicTwin(δ = 0.2, ρ = 0.1, metric = Euclidean())

A pseudoperiodic twin surrogate[Miralles2015], which is a fusion of the twin surrogate[Thiel2006] and the pseudo-periodic surrogate[Small2001].

Input parameters

A delay reconstruction of the input timeseries is constructed using embedding dimension d and embedding delay τ. The threshold δ ∈ (0, 1] determines which points are "close" (neighbors) or not, and is expressed as a fraction of the attractor diameter, as determined by the input data. The authors of the original twin surrogate paper recommend 0.05 ≤ δ ≤ 0.2[Thiel2006].

If you have pre-embedded your timeseries, and timeseries is already a ::Dataset, use the three-argument constructor (so that no delay reconstruction is performed). If you want a surrogate for a scalar-valued timeseries, use the five-argument constructor to also provide the embedding delay τ and embedding dimension d.

Null hypothesis

Pseudo-periodic twin surrogates generate signals similar to the original data if the original signal is (quasi-)periodic. If the original signal is not (quasi-)periodic, then these surrogates will have different recurrence plots than the original signal, but preserve the overall shape of the attractor. Thus, PseudoPeriodicTwin surrogates can be used to test null hypothesis that the observed timeseries (or orbit) is consistent with a quasi-periodic orbit[Miralles2015].

Returns

A d-dimensional surrogate orbit (a Dataset) is returned. Sample the first column of this dataset if a scalar-valued surrogate is desired.

source
TimeseriesSurrogates.WLSType
WLS(shufflemethod::Surrogate = IAAFT();
    f::Union{Nothing, Function} = Statistics.cor,
    rescale::Bool = true,
    wt::Wavelets.WT.OrthoWaveletClass = Wavelets.WT.Daubechies{16}())

A wavelet surrogate generated by the following procedure:

  1. Compute the wavelet transform of the signal. This results in a set of detail coefficients over a set of dyadic scales. As in Keylock (2006), we here use the maximal overlap discrete wavelet transform, or MODWT, so that the number of coefficients at each scale are the same.
  2. Shuffle the detail coefficients at each dyadic scale using the provided shufflemethod. See "Shuffling methods" below for alternatives.
  3. Apply the inverse wavelet transform to the shuffled detail coefficients to obtain a surrogate time series.

Shuffling methods

You may choose to use any surrogate from this package to perform the randomization of the detail coefficients at each dyadic scale.

The following methods have been discussed in the literature (more may exist):

  • Random permutations of wavelet coefficients within each scale (Breakspear et al., 2003). To get this behaviour, use WLS(x, RandomShuffle(), rescale = false, f = nothing).
  • Cyclic rotation of wavelet coefficients within each scale (Breakspear et al., 2003). To get this behaviour, use WLS(x, Circshift(1:length(x)), rescale = false, f = nothing).
  • Block resampling of wavelet coefficients within each scale (Breakspear et al., 2003). To get this behaviour, use WLS(x, BlockShuffle(nblocks, randomize = true), rescale = false, f = nothing).
  • IAAFT resampling of wavelet coefficients within each scale (Keylock, 2006). To get this behaviour, use WLS(x, IAAFT(), rescale = true, f = Statistics.cor). This method preserves the local mean and variance structure of the signal, but randomises nonlinear properties of the signal (i.e. Hurst exponents)[Keylock2006]. These surrogates can therefore be used to test for changes in nonlinear properties of the original signal. In contrast to IAAFT surrogates, the IAAFT-wavelet surrogates also preserves nonstationarity. Using other shufflemethods does not necessarily preserve nonstationarity. To deal with nonstationary signals, Keylock (2006) recommends using a wavelet with a high number of vanishing moments. Thus, our default is to use a Daubechies wavelet with 16 vanishing moments. Note: The iterative procedure after the rank ordering step (step [v] in [Keylock2006]) is not performed in this implementation.

The default method and parameters replicate the behaviour of Keylock (2006)'s IAAFT wavelet surrogates.

Error minimization

For the IAAFT approach introduced in Keylock (2006), detail coefficients at each level are circularly rotated to minimize an error function. The methods introduced in Breakspear et al. (2003) do not apply this error minimization.

In our implementation, you can turn this option on/off using the f parameter of the WLS constructor. If f = nothing turns off error minization. If f is set to a two-argument function that computes some statistic, for example f = Statistics.cor, then detail coefficients at each scale are circularly rotated until that function is maximized (and hence the "error" minimized). If you want to minimize some error function, then instead provide an appropriate transform of your function. For example, if using the root mean squared deviation, define rmsd_inv(x, y) = 1 - StatsBase.rmsd(x, y) and set f = rmsd_inv.

Rescaling

If rescale == true, then surrogate values are mapped onto the values of the original time series, as in the AAFT algorithm. If rescale == false, surrogate values are not constrained to the original time series values. If AAFT or IAAFT shuffling is used, rescale should be set to true. For other methods, it does not necessarily need to be.

source
TimeseriesSurrogates.RandomCascadeType
RandomCascade(paddingmode::String = "zeros")

A random cascade multifractal wavelet surrogate (Paluš, 2008)[Paluš2008].

If the input signal length is not a power of 2, the signal must be padded before the surrogate is constructed. paddingmode determines how the signal is padded. Currently supported padding modes: "zeros". The final surrogate (constructed from the padded signal) is subset to match the length of the original signal.

Random cascade surrogate preserve multifractal properties of the input time series, that is, interactions among dyadic scales and nonlinear dependencies[Paluš2008].

source
TimeseriesSurrogates.ShuffleDimensionsType
ShuffleDimensions()

Multidimensional surrogates of input datasets (DelayEmbeddings.Dataset, which are also multidimensional) that have shuffled dimensions in each point.

These surrogates destroy the state space structure of the dataset and are thus suited to distinguish deterministic datasets from high dimensional noise.

source
TimeseriesSurrogates.IrregularLombScargleType
IrregularLombScargle(t; tol = 1, n_total = 100000, n_acc = 50000, q = 1)

IrregularLombScargle surrogates for unevenly sampled time series with supporting time steps t, generated using the simulated annealing algorithm described in [SchreiberSchmitz1999].

IrregularLombScargle surrogates (given enough iterations and a low enough tolerance) preserve the periodogram and the amplitude distribution of the original signal. For time series with equidistant time steps, surrogates generated by this method result in surrogates similar to those produced by the IAAFT method.

This algorithm starts with a random permutation of the original data. Then it iteratively approaches the power spectrum of the original data by swapping two randomly selected values in the surrogate data if the Minkowski distance of order q between the power spectrum of the surrogate data and the original data is less than before. The iteration procedure ends when the relative deviation between the periodograms is less than tol or when n_total number of tries or n_acc number of actual swaps is reached.

source

Utilities

TimeseriesSurrogates.noiseradiusFunction
noiseradius(x::AbstractVector, d::Int, τ, ρs, n = 1) → ρ

Use the proposed* algorithm of[Small2001] to estimate optimal ρ value for PseudoPeriodic surrogates, where ρs is a vector of possible ρ values. *The paper is ambiguous about exactly what to calculate. Here we count how many times we have pairs of length-2 that are identical in x and its surrogate, but are not also part of pairs of length-3.

This function directly returns the arg-maximum of the evaluated distribution of these counts versus ρ, use TimeseriesSurrogates._noiseradius with same arguments to get the actual distribution. n means to repeat τhe evaluation n times, which increases accuracy.

source

Visualization

TimeseriesSurrogates.jl has defined a simple function surroplot(x, s). This comes into scope when using Makie (you also need a plotting backend).

To load the function, do:

using TimeseriesSurrogates
using CairoMakie, Makie
using TimeseriesSurrogates, CairoMakie, Makie
ts = AR1() # create a realization of a random AR(1) process
s = surrogate(ts, AAFT())
fig = surroplot(ts, s)

Citing

Please use the following BiBTeX entry, or DOI, to cite TimeseriesSurrogates.jl:

DOI: https://doi.org/10.21105/joss.04414

BiBTeX:

@article{TimeseriesSurrogates.jl,
    doi = {10.21105/joss.04414},
    url = {https://doi.org/10.21105/joss.04414},
    year = {2022},
    publisher = {The Open Journal},
    volume = {7},
    number = {77},
    pages = {4414},
    author = {Kristian Agasøster Haaga and George Datseris},
    title = {TimeseriesSurrogates.jl: a Julia package for generating surrogate data},
    journal = {Journal of Open Source Software}
}
  • Theiler1991J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, J. Farmer, Testing for nonlinearity in time series: The method of surrogate data, Physica D 58 (1–4) (1992) 77–94.
  • Theiler1994J. Theiler, On the evidence for low-dimensional chaos in an epileptic electroencephalogram, Phys. Lett. A 196
  • Theiler1991J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, J. Farmer, Testing for nonlinearity in time series: The method of surrogate data, Physica D 58 (1–4) (1992) 77–94.
  • Nakamura2006Nakamura, Tomomichi, Michael Small, and Yoshito Hirata. "Testing for nonlinearity in irregular fluctuations with long-term trends." Physical Review E 74.2 (2006): 026205.
  • Lucio2012Lucio, J. H., Valdés, R., & Rodríguez, L. R. (2012). Improvements to surrogate data methods for nonstationary time series. Physical Review E, 85(5), 056202.
  • Ortega1998Ortega, Guillermo J.; Louis, Enrique (1998). Smoothness Implies Determinism in Time Series: A Measure Based Approach. Physical Review Letters, 81(20), 4345–4348. doi:10.1103/PhysRevLett.81.4345
  • Ortega1998Ortega, Guillermo J.; Louis, Enrique (1998). Smoothness Implies Determinism in Time Series: A Measure Based Approach. Physical Review Letters, 81(20), 4345–4348. doi:10.1103/PhysRevLett.81.4345
  • Nakamura2006Nakamura, Tomomichi, Michael Small, and Yoshito Hirata. "Testing for nonlinearity in irregular fluctuations with long-term trends." Physical Review E 74.2 (2006): 026205.
  • Nakamura2006Nakamura, Tomomichi, Michael Small, and Yoshito Hirata. "Testing for nonlinearity in irregular fluctuations with long-term trends." Physical Review E 74.2 (2006): 026205.
  • Lucio2012Lucio, J. H., Valdés, R., & Rodríguez, L. R. (2012). Improvements to surrogate data methods for nonstationary time series. Physical Review E, 85(5), 056202.
  • Lucio2012Lucio, J. H., Valdés, R., & Rodríguez, L. R. (2012). Improvements to surrogate data methods for nonstationary time series. Physical Review E, 85(5), 056202.
  • Lucio2012Lucio, J. H., Valdés, R., & Rodríguez, L. R. (2012). Improvements to surrogate data methods for nonstationary time series. Physical Review E, 85(5), 056202.
  • Theiler1991J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, J. Farmer, Testing for nonlinearity in time series: The method of surrogate data, Physica D 58 (1–4) (1992) 77–94.
  • Theiler1991J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, J. Farmer, Testing for nonlinearity in time series: The method of surrogate data, Physica D 58 (1–4) (1992) 77–94.
  • Nakamura2006Nakamura, Tomomichi, Michael Small, and Yoshito Hirata. "Testing for nonlinearity in irregular fluctuations with long-term trends." Physical Review E 74.2 (2006): 026205.
  • SchreiberSchmitz1996T. Schreiber; A. Schmitz (1996). "Improved Surrogate Data for Nonlinearity Tests". Phys. Rev. Lett. 77 (4)
  • Small2001Small et al., Surrogate test for pseudoperiodic time series data, Physical Review Letters, 87(18)
  • Small2001Small et al., Surrogate test for pseudoperiodic timeseries data, Physical Review Letters, 87(18)
  • Thiel2006Thiel, Marco, et al. "Twin surrogates to test for complex synchronisation." EPL (Europhysics Letters) 75.4 (2006): 535.
  • Miralles2015Miralles, R., et al. "Characterization of the complexity in short oscillating timeseries: An application to seismic airgun detonations." The Journal of the Acoustical Society of America 138.3 (2015): 1595-1603.
  • Breakspear2003Breakspear, M., Brammer, M., & Robinson, P. A. (2003). Construction of multivariate surrogate sets from nonlinear data using the wavelet transform. Physica D: Nonlinear Phenomena, 182(1-2), 1-22.
  • Keylock2006C.J. Keylock (2006). "Constrained surrogate time series with preservation of the mean and variance structure". Phys. Rev. E. 73: 036707. doi:10.1103/PhysRevE.73.036707.
  • Paluš2008Paluš, Milan (2008). Bootstrapping Multifractals: Surrogate Data from Random Cascades on Wavelet Dyadic Trees. Physical Review Letters, 101(13), 134101–. doi:10.1103/PhysRevLett.101.134101
  • SchmitzSchreiber1999A.Schmitz T.Schreiber (1999). "Testing for nonlinearity in unevenly sampled time series" Phys. Rev E
  • Small2001Small et al., Surrogate test for pseudoperiodic time series data, Physical Review Letters, 87(18)