R: Red-Noise Spectra of Time-Series

redfit {dplR}

R Documentation

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)

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.

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.

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`	A `list` with the following items: np number of points in the data. nseg number of points in each segment. nfreq number of frequencies in the results. avgdt average sampling interval. df difference between consecutive frequencies. fnyq average Nyquist frequency. n50 value of the `n50` argument. ofac value of the `ofac` argument. hifac value of the `hifac` argument. segskip the ideal, non-rounded difference between starting indices of consecutive segments. iwin value of the `iwin` argument. If a `character` value was used, this contains the corresponding number used internally. nsim value of the `nsim` argument. mctest value of the `mctest` argument. rhopre value of the `rhopre` argument.
`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`.

Author(s)

Mikko Korpela. Examples by Andy Bunn.

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

# 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()

## Not run: 
# 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)

## End(Not run)
.Random.seed <- rs

[Package dplR version 1.7.7 Index]