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)

Arguments

a numeric vector representing a possibly unevenly spaced time-series.

a numeric vector of times or ages corresponding to 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",

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.

a numeric or bigq 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"

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

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

an integral value giving the length of the sequence in the
    number of runs test.

`Value`

Function runcrit returns a list containing
  rcritlo, rcrithi and rcritexact
  (see below). Function redfit returns a list with the
  following elements:
varxvariance of x estimated from its spectrum.
    A numeric value.
rhoaverage autocorrelation coefficient, either estimated
    from the data or prescribed (rhopre). A numeric
    value.
tauthe time scale of an AR1 process corresponding to
    rho. A numeric value.
rcnta 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.
rcritloa numeric vector of critical low values for
    the number of runs, i.e. the lowest value for accepting the null
    hyphothesis at each level of significance p.  When
    returned from redfit, NULL when rcnt is
    NULL.
rcrithia numeric vector of critical high values for
    the number of runs, i.e. the highest value for accepting the null
    hyphothesis at each level of significance p.  When
    returned from redfit, NULL when rcnt is
    NULL.
rcritexacta 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.
freqthe frequencies used.  A numeric vector.  The
    other numeric vectors have the same length, i.e. one value for each
    frequency studied.
gxxestimated spectrum of the data (t, x).  A
    numeric vector.
gxxcred noise corrected spectrum of the data.  A
    numeric vector.
grravgaverage AR1 spectrum over nsim
    simulations.  A numeric vector.
gredththeoretical AR1 spectrum.  A numeric vector.
corra 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.
ci80a numeric vector containing the bias-corrected
    80th percentile (by frequency) red noise spectrum.  Only if
    mctest is TRUE.
ci90a numeric vector containing the 90th percentile
    red noise spectrum.
ci9595th percentile red noise spectrum.
ci9999th percentile red noise spectrum.
callthe call to the function.  See match.call.
paramsA list with the following items:
    [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object]
versversion of dplR.  See packageVersion.
seedvalue of the seed argument.
tif duplicated values of t are given, the
    non-duplicated numeric time or age vector (see
    tType) is returned here.  Otherwise NULL.  See
    txOrdered.
xif duplicated values of t are given, the
    averaged numeric data vector is returned here.  Otherwise
    NULL.

`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
    betweenredfitand 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 ofruncritbelow) differ betweenredfitand 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.
  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, the normal distribution is used for
  approximation.

`References`

Function redfit is based on the Fortran program
  http://www.geo.uni-bremen.de/geomod/staff/mschulz/{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.

`See Also`

print.redfit

`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)
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)

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 (dB)", 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()

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

data(cana157)
yrs <- as.numeric(rownames(cana157))
x <- cana157[, 1]
redf.x <- redfit(x, nsim = 1000)

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

plot(redf.x[["freq"]], redf.x[["gxxc"]],
     ylim = range(redf.x[["ci99"]], redf.x[["gxxc"]]),
     type = "n", ylab = "Spectrum (dB)", 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)
Run the code above in your browser using DataLab