redfit: Red-Noise Spectra of Time-Series

Description

Estimate red-noise spectra from a possibly unevenly spaced time-series.

Usage

redfit(x, t, tType = c("time", "age"), nsim = 1000, mctest = TRUE,
       ofac = 4, hifac = 1, n50 = 3, rhopre = NULL,
       p = c(0.10, 0.05, 0.02), iwin = 2,
       txOrdered = FALSE, verbose = FALSE, seed = NULL,
       maxTime = 10, nLimit = 10000)
runcrit(n, p = c(0.10, 0.05, 0.02), maxTime = 10, nLimit = 10000)

Value

Function runcrit returns a list containing

rcritlo, rcrithi and rcritexact

(see below). Function redfit returns a list with the following elements:

varx: variance of x estimated from its spectrum. A numeric value.
rho: average autocorrelation coefficient, either estimated from the data or prescribed (rhopre). A numeric value.
tau: the time scale of an AR1 process corresponding to rho. A numeric value.
rcnt: a numeric value giving the number of runs to be used for a statistical test studying the difference between a theoretical AR1 spectrum and the bias-corrected spectrum estimated from the data. Null hypothesis: the two spectra agree, i.e. the probability of either being larger than the other is 0.5 at every point. Requires that iwin == 0 ("rectangular"), ofac == 1 and n50 == 1. Otherwise the test is not performed and rcnt is NULL.
rcritlo: a numeric vector of critical low values for the number of runs, i.e. the lowest value for accepting the null hypothesis at each level of significance p. When returned from redfit, NULL when rcnt is NULL.
rcrithi: a numeric vector of critical high values for the number of runs, i.e. the highest value for accepting the null hypothesis at each level of significance p. When returned from redfit, NULL when rcnt is NULL.
rcritexact: a logical vector specifying whether each pair of rcritlo and rcrithi are exact values (TRUE) or approximated from a normal distribution (FALSE). When returned from redfit, NULL when rcnt is NULL.
freq: the frequencies used. A numeric vector. The other numeric vectors have the same length, i.e. one value for each frequency studied.
gxx: estimated spectrum of the data (t, x). A numeric vector.
gxxc: red noise corrected spectrum of the data. A numeric vector.
grravg: average AR1 spectrum over nsim simulations. A numeric vector.
gredth: theoretical AR1 spectrum. A numeric vector.
corr: a numeric vector containing the by-frequency correction: gxxc equals gxx divided by corr (or multiplied by the inverse correction). Also used for computing ci80, ci90, ci95 and ci99.
ci80: a numeric vector containing the bias-corrected 80th percentile (by frequency) red noise spectrum. Only if mctest is TRUE.
ci90: a numeric vector containing the 90th percentile red noise spectrum.
ci95: 95th percentile red noise spectrum.
ci99: 99th percentile red noise spectrum.
call: the call to the function. See match.call.
params
vers: version of dplR. See packageVersion.
seed: value of the seed argument.
t: if duplicated values of t are given, the non-duplicated numeric time or age vector (see tType) is returned here. Otherwise NULL. See txOrdered.
x: if duplicated values of t are given, the averaged numeric data vector is returned here. Otherwise NULL.

Arguments

x: a numeric vector representing a possibly unevenly spaced time-series.
t: a numeric vector of times or ages corresponding to x. If missing, this is set to an evenly spaced increasing integer sequence (1, 2, ...) along x. See txOrdered.
tType: a character string indicating the type of the t vector: either times or ages.
nsim: a numeric value giving the number of simulated AR1 spectra to compute.
mctest: a logical flag. If TRUE, performs a Monte Carlo test for computing red noise false-alarm levels. In that case, the result list will contain non-NULL elements "ci80", "ci90", "ci95" and "ci99".
ofac: oversampling factor for Lomb-Scargle Fourier transform. A numeric value.
hifac: maximum frequency to analyze relative to the Nyquist frequency. A numeric value.
n50: number of segments. The segments overlap by about 50 percent.
rhopre: a numeric value giving the prescribed autocorrelation coefficient. If NULL or negative, the autocorrelation coefficient will be estimated from the data.
p: a numeric or bigq (if package "gmp" is installed) vector of significance levels for a statistical test considering the number of runs in a sequence. See ‘Details’.
iwin: the type of window used for scaling the values of each segment of data. A numeric value or one of "rectangular", "welch i", "hanning", "triangular" and "blackman-harris". Integers 0:4 correspond to the character values, in this order. The default iwin = 2 corresponds to the "hanning" window.
txOrdered: a logical flag. If TRUE, it is assumed that t is in ascending order without duplicates. If FALSE (the default), t will be sorted and x reordered accordingly. Values of x at identical values of t are averaged. If duplicates are found, the non-duplicated t and averaged x will be included in the return value of the function.
verbose: a logical flag. If TRUE, the function prints some information while operating.
seed: a value to be given to set.seed at the start of the function. The value is also recorded in the list returned. If not NULL, this can be used for reproducible results.
maxTime: a numeric value giving the approximate maximum time in seconds to use for computing the exact acceptance regions of the number of runs test. See also nLimit.
nLimit: a numeric value giving the maximum n for which runcrit will try to compute an exact solution to the acceptance regions of the number of runs test. Precomputed solutions may exist for larger n). This limit is in place because a part of the exact solution always needs to be computed in order to roughly estimate the total time and whether it would exceed maxTime. If nLimit is very large, maxTime may be (greatly) exceeded while computing the aforementioned part of the exact solution.
n: an integral value giving the length of the sequence in the number of runs test.

Author

Mikko Korpela. Examples by Andy Bunn.

Details

Function redfit computes the spectrum of a possibly unevenly sampled time-series by using the Lomb-Scargle Fourier transform. The spectrum is bias-corrected using spectra computed from simulated AR1 series and the theoretical AR1 spectrum.

The function duplicates the functionality of program REDFIT by Schulz and Mudelsee. See the manual of that program for more information. The results of this function should be very close to REDFIT. However, some changes have been made:

More precision is used in some constants and computations.
All the data are used: the last segment always contains the last pair of (t, x). There may be small differences between redfit and REDFIT with respect to the number of points per segment and the overlap of consecutive segments.
The critical values of the runs test (see the description of runcrit below) differ between redfit and REDFIT. The approximate equations in REDFIT produce values quite far off from the exact values when the number of frequencies is large.
The user can select the significance levels of the runs test.
Most of the window functions have been adjusted.
6 dB bandwidths have been computed for discrete-time windows.

Function runcrit computes the limits of the acceptance region of a number of runs test: assuming a sequence of n i.i.d. discrete random variables with two possible values a and b of equal probability (0.5), we are examining the distribution of the number of runs. A run is an uninterrupted sequence of only a or only b. The minimum number of runs is 1 (a sequence with only a or only b) while the maximum number is n (alternating a and b). See Bradley, p. 253–254, 259–263. The function is also called from redfit; see rcnt in ‘Value’ for the interpretation. In this case the arguments p, maxTime and nLimit are passed from redfit to runcrit, and n is the number of output frequencies.

The results of runcrit have been essentially precomputed for some values of p and n. If a precomputed result is not found and n is not too large (nLimit, maxTime), the exact results are computed on-demand. Otherwise, or if package "gmp" is not installed, the normal distribution is used for approximation.

References

Function redfit is based on the Fortran program REDFIT (version 3.8e), which is in the public domain.

Bradley, J. V. (1968) Distribution-Free Statistical Tests. Prentice-Hall.

Schulz, M. and Mudelsee, M. (2002) REDFIT: estimating red-noise spectra directly from unevenly spaced paleoclimatic time series. Computers & Geosciences, 28(3), 421–426.

Examples

Run this code

# Create a simulated tree-ring width series that has a red-noise
# background ar1=phi and sd=sigma and an embedded signal with 
# a period of 10 and an amplitude of have the rednoise sd.
library(graphics)
library(stats)
runif(1)
rs <- .Random.seed
set.seed(123)
nyrs <- 500
yrs <- 1:nyrs

# Here is an ar1 time series with a mean of 2mm,
# an ar1 of phi, and sd of sigma
phi <- 0.7
sigma <- 0.3
sigma0 <- sqrt((1 - phi^2) * sigma^2)
x <- arima.sim(list(ar = phi), n = nyrs, sd = sigma0) + 2

# Here is a sine wave at f=0.1 to add in with an amplitude
# equal to half the sd of the red noise background
per <- 10
amp <- sigma0 / 2
wav <- amp * sin(2 * pi / per * yrs)

# Add them together so we have signal and noise
x <- x + wav

# Here is the redfit spec
redf.x <- redfit(x, nsim = 500)

# Acceptance region of number of runs test
# (not useful with default arguments of redfit())
runcrit(length(redf.x[["freq"]]))

op <- par(no.readonly = TRUE) # Save to reset on exit
par(tcl = 0.5, mar = rep(2.2, 4), mgp = c(1.1, 0.1, 0))

plot(redf.x[["freq"]], redf.x[["gxxc"]],
     ylim = range(redf.x[["ci99"]], redf.x[["gxxc"]]),
     type = "n", ylab = "Spectrum", xlab = "Frequency (1/yr)",
     axes = FALSE)
grid()
lines(redf.x[["freq"]], redf.x[["gxxc"]], col = "#1B9E77")
lines(redf.x[["freq"]], redf.x[["ci99"]], col = "#D95F02")
lines(redf.x[["freq"]], redf.x[["ci95"]], col = "#7570B3")
lines(redf.x[["freq"]], redf.x[["ci90"]], col = "#E7298A")
freqs <- pretty(redf.x[["freq"]])
pers <- round(1 / freqs, 2)
axis(1, at = freqs, labels = TRUE)
axis(3, at = freqs, labels = pers)
mtext(text = "Period (yr)", side = 3, line = 1.1)
axis(2); axis(4)
legend("topright", c("x", "CI99", "CI95", "CI90"), lwd = 2,
       col = c("#1B9E77", "#D95F02", "#7570B3", "#E7298A"),
       bg = "white")
box()

if (FALSE) {
# Second example with tree-ring data
# Note the long-term low-freq signal in the data. E.g.,
# crn.plot(cana157)

library(utils)
data(cana157)
yrs <- time(cana157)
x <- cana157[, 1]
redf.x <- redfit(x, nsim = 1000)

plot(redf.x[["freq"]], redf.x[["gxxc"]],
     ylim = range(redf.x[["ci99"]], redf.x[["gxxc"]]),
     type = "n", ylab = "Spectrum", xlab = "Frequency (1/yr)",
     axes = FALSE)
grid()
lines(redf.x[["freq"]], redf.x[["gxxc"]], col = "#1B9E77")
lines(redf.x[["freq"]], redf.x[["ci99"]], col = "#D95F02")
lines(redf.x[["freq"]], redf.x[["ci95"]], col = "#7570B3")
lines(redf.x[["freq"]], redf.x[["ci90"]], col = "#E7298A")
freqs <- pretty(redf.x[["freq"]])
pers <- round(1 / freqs, 2)
axis(1, at = freqs, labels = TRUE)
axis(3, at = freqs, labels = pers)
mtext(text = "Period (yr)", side = 3, line = 1.1)
axis(2); axis(4)
legend("topright", c("x", "CI99", "CI95", "CI90"), lwd = 2,
       col = c("#1B9E77", "#D95F02", "#7570B3", "#E7298A"),
       bg = "white")
box()
par(op)
}
.Random.seed <- rs

Run the code above in your browser using DataLab