evaluate.pm: Evaluate the periodogram maximum of a time series

Description

Calculate the Lomb-Scargle periodogram of a time series and estimate the statistical significance of the periodogram maximum based on the null hypothesis of an Ornstein-Uhlenbeck state space (OUSS) process.

Usage

evaluate.pm(times, signal, minPeakFreq = 0, minFitFreq = 0, 
            iterations = 100, accuracy = 0.005, startRadius = 5, 
            verbose = FALSE)

Value

A list with the following entries:

error: Will be TRUE if an error occured, FALSE otherwise.
errorMessage: A short error message if error==TRUE.

If error==FALSE, the returned list also includes:

frequencies: Available periodogram frequencies as a numerical vector.
periodogram: Periodogram powers corresponding to the returned frequencies, as a numerical vector.
fittedPS: Maximum-likelihood fitted OUSS power spectrum corresponding to frequencies, as a numerical vector.
power_o: Estimated power spectrum at zero-frequency generated by the underlying OU process.
lambda: Estimated resilience of the OU process (in inverse time units).
power_e: Estimated power spectrum at large frequencies due to the random measurement errors.
sigma: Estimated stationary standard deviation of the underlying OU process.
epsilon: Estimated standard deviation of measurement errors.
time_step: The average time step of the time series, as used to fit the OUSS power spectrum.
peakMode: An integer indicating the position of the periodogram maximum in the vector frequencies.
minPeakMode: The minimum periodogram mode considered for determining the periodogram maximum. This will be 1 if minPeakFreq==0.
minFitMode: The minimum periodogram mode considered for estimating the white noise power. This will be 1 if minFitFreq==0.
MLL: Log-likelihood value at calculated maximum.
P: Statistical significance of the periodogram peak against the null hypothesis of the estimated OUSS process. This is the probability that the estimated OUSS process would generate a periodogram with global maximum at least as ``extreme'' as the observed peak (among the considered frequencies). See details above.
Plocal: Statistical significance of the relative power of the periodogram peak. Mainly used for comparison to the statistical significance of secondary peaks. See significanceOfLocalPeak.

Arguments

times: Numerical vector. Time points of the time series.
signal: Numerical vector of same size as times. Values of the time series.
minPeakFreq: Non-negative number. Minimum considered frequency when determining periodogram peak. Use this if you want to ignore low-frequency components in the spectrum when searching for the periodogram maximum.
minFitFreq: Non-negative number. Minimum considered frequency when fitting the OUSS model to the periodogram. Use this to ignore low-frequency components in the spectrum when estimating the OUSS parameters, e.g. if you suspect your time series to be trended. minFitFreq might differ from minPeakFreq.
iterations: Number of iterations for the maximum-likelihood fitter. Increasing this will result in greater estimation accuracy but also reduced performance.
startRadius: Single integer. The number of initial guesses for each OUSS parameter during optimization of the likelihood function, per direction (up and down). Increasing this will improve estimation accuracy. However, execution time scales with startRadius^2.
accuracy: An upper bound for the standard deviation of the P-value estimator using repeated Bernoulli trials. The number of trials scales as 1/accuracy^2. If you are using 0.05 as a nominal significance threshold, then accuracy=0.005 is probably good enough.
verbose: Single logical. If TRUE, the function prints out occasional progress reports.

Author

Stilianos Louca

Details

The OUSS model describes the measurement of an Ornstein-Uhlenbeck (OU) stochastic process at discrete times with additional uncorrelated Gaussian measurement errors of fixed variance. The OU process itself is a continuous-time random walk with linear stabilizing forces, described by the stochastic differential equation $$ dX = \lambda(\mu-X) dt + s dW, $$ where $W$ is the standard Wiener process. The OUSS process is a parsimonious model for describing stochastically fluctuating variables subject to linear stabilizing forces and uncorrelated measurement errors.

Due to temporal correlations, the OUSS power spectrum increases gradually towards lower frequencies, as opposed to the white noise spectrum, which is flat. Using white noise as a null hypothesis for the evaluation of cyclicity in time series, particularly for systems with long correlation times, may result in increased false cycle detection rates because of the increased low-frequency power. The OUSS model is an attempt to account for these correlations.

The OUSS model parameters are estimated using maximum-likelihood fitting to the periodogram. The likelihood function, and therefore the OUSS parameter estimates, are only approximations that become exact for long regular time series. The statistical significance of the periodogram peak (power $S$ at frequency $F$) under the null-hypothesis of an OUSS process is defined as the probability that the same OUSS process would generate a periodogram whose maximum (power $s$ at frequency $f$) satisfies $$ s^2/e\ge S^2/E, $$ where $e$ and $E$ are the expected periodogram powers at the frequencies $f$ and $F$, respectively. The P-value is estimated via repeated Bernoulli trials in which random OUSS periodograms are emulated by exponentially distributed numbers.

If you want to evaluate secondary peaks (i.e. non-global periodogram maxima), you will need to either (a) adjust the parameters minPeakFreq and minFitFreq to omit low-frequency modes or (b) use significanceOfLocalPeak after using evaluate.pm.

References

Louca, S., Doebeli, M. (2015) Detecting cyclicity in ecological time series, Ecology 96: 1724--1732

Examples

Run this code

# In this example we generate a cyclic time series
# and analyse its periodogram using evaluate.pm

# Parameters
lambda         = 1;    # inverse correlation time of OU process
cyclePeriod    = 1;    # Cycle period
cycleAmplitude = 0.6;
times          = seq(0,20,0.25);

# generate cyclic time series by adding a periodic signal to an OUSS process
signal = cycleAmplitude * cos(2*pi*times/cyclePeriod) +
         generate_ouss(times, mu=0, sigma=1, lambda=lambda, epsilon=0.5);

# Find periodogram peak and estimate statistical significance
# Ignore frequencies lower than a pre-defined threshold
# to avoid masking by low-frequency maximum
report = evaluate.pm(times=times, signal=signal, 
                     minPeakFreq=lambda/3, 
                     minFitFreq=lambda/3,
                     startRadius=2);

# plot overview of periodogram peak analysis
plotReport(sprintf("Cyclic at frequency %.2g",1/cyclePeriod), 
           times=times, signal=signal, report=report);

Run the code above in your browser using DataLab