TimeseriesSurrogates.jl

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.

surrogenerator(x, method::Surrogate [, rng]) → sgen::SurrogateGenerator

Initialize a generator that creates surrogates of x on demand, based on the given method. This is more efficient than surrogate, 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.

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 sgen as a function with no arguments, e.g.:

sgen = surrogenerator(x, method)
s = sgen()

You can use the generator syntax of Julia to map over surrogates generated by sg. For example, let q be a function returning a discriminatory statistic. To test some null hypothesis with TimeseriesSurrogates.jl you'd do

using TimeseriesSurrogates
q, x # inputs
method = RandomFourier() # some example method
sgen = surrogenerator(x, method)
siter = (sgen() for _ in 1:1000)
qx = q(x)
qs = map(q, siter)
# compare `qx` with quantiles
using Statistics: quantile
q01, q99 = quantile(qs, [0.01, 0.99])
q01 ≤ qx ≤ q99 # if false, hypothesis can be rejected!

SurrogateTest(f::Function, x, method::Surrogate; kwargs...) → test

Initialize a surrogate test for input data x, which can be used in pvalue. The tests requires as input a function f that given a timeseries (like x) it outputs a real number, and a method of how to generate surrogates. f is the function that computes the discriminatory statistic.

Once called with pvalue, the test estimates and then stores the real value rval and surrogate values vals of the discriminatory statistic in the fields rval, vals respectively. Alternatively, you can use fill_surrogate_test! directly if you don't care about the p-value.

SurrogateTest automates the process described in the documentation page Performing surrogate hypothesis tests.

SurrogateTest subtypes HypothesisTest and is part of the StatsAPI.jl interface.

Keywords

• rng = Random.default_rng(): a random number generator.
• n::Int = 10_000: how many surrogates to generate and compute f on.
• threaded = true: Whether to parallelize looping over surrogate computations in to the available threads (Threads.nthreads()).

fill_surrogate_test!(test::SurrgateTest) → rval, vals

Perform the computations foreseen by test and return the value of the discriminatory statistic for the real data rval and the distribution of values for the surrogates vals.

This function is called by pvalue.

pvalue(test::SurrogateTest; tail = :left)

Return the p-value corresponding to the given SurrogateTest, optionally specifying what kind of tail test to do (one of :left, :right, :both).

For SurrogateTest, the p-value is simply the proportion of surrogate statistics that exceed (for tail = :right) or subseed (tail = :left) the discriminatory statistic computed from the input data.

The default value of tail assumes that the surrogate data are expected to have higher discriminatory statistic values. This is the case for statistics that quantify entropy. For statistics that quantify autocorrelation, use tail = :right instead.

Supertype of all surrogate methods.

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]}.

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.

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.

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.

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]}.

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.

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.

See RelativePartialRandomization and SpectralPartialRandomization for alternative partial-randomization algorithms

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.

RelativePartialRandomization(α = 0.5)

RelativePartialRandomization surrogates are similar to PartialRandomization phase surrogates, but instead of drawing phases uniformly from [0, 2π], phases are drawn from ϕ + [0, 2π]*α, where α ∈ [0, 1] and ϕ is the original Fourier phase.

See the documentation for a detailed comparison between partial randomization algorithms.

RelativePartialRandomizationAAFT(α = 0.5)

RelativePartialRandomizationAAFT surrogates are similar to RelativePartialRandomization surrogates, but add a rescaling step, so that the surrogate has the same values as the original time series (analogous to the rescaling done for AAFT surrogates).

SpectralSpectralPartialRandomization(α = 0.5)

SpectralPartialRandomization surrogates are similar to PartialRandomization phase surrogates, but instead of drawing phases uniformly from [0, 2π], phases of the highest frequency components responsible for a proportion α of power are replaced by random phases drawn from [0, 2π]

See the documentation for a detailed comparison between partial randomization algorithms.

SpectralPartialRandomizationAAFT(α = 0.5)

SpectralPartialRandomizationAAFT surrogates are similar to PartialRandomization surrogates, but add a rescaling step, so that the surrogate has the same values as the original time series (analogous to the rescaling done for AAFT surrogates).

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]}.

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.

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.

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 fϵ 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 fϵ 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 fϵ 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.

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.

TFTDAAFT(fϵ = 0.05)

TFTDAAFT^[Lucio2012] are similar to TFTD surrogates, but also re-scales back to the original values of the time series. fϵ ∈ (0, 1] is the fraction of the powerspectrum corresponding to the lowermost frequencies to be preserved.

See also: TFTD, TFTDIAAFT.

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.

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 embedding from DelayEmbeddings.jl (see that docs for choosing appropriate 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.

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 ::StateSpaceSet, 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 StateSpaceSet) is returned. Sample the first column of this dataset if a scalar-valued surrogate is desired.

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.

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].

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 a convenient 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.φ

ShuffleDimensions()

Multidimensional surrogates of input StateSpaceSets from StateSpaceSets.jl. Each point in the set is individually shuffled, but the points themselves are not shuffled.

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

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.

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.