Learn R Programming

WaveletComp (version 1.0)

analyze.wavelet: Computation of the wavelet power spectrum of a single time series

Description

The time series is selected from an input data frame by specifying either its name or its column number. Optionally, the time series is detrended, using loess with parameter loess.span. Internally, the series will be further standardized before it undergoes wavelet transformation.

The wavelet power spectrum is computed by applying the Morlet wavelet. P-values to test the null hypothesis that a period (within lowerPeriod and upperPeriod) is irrelevant at a certain time are calculated if desired; this is accomplished with the help of a simulation algorithm. There is a selection of models from which to choose the alternative hypothesis. The selected model will be fitted to the data and simulated according to estimated parameters in order to provide surrogate time series.

Wavelet transformation, as well as p-value computations, are carried out by calling subroutine wt.

The name and parts of the layout of subroutine wt were inspired by a similar function developed by Huidong Tian and Bernard Cazelles (archived R package WaveletCo). The basic concept of the simulation algorithm, and of ridge determination build on ideas developed by these authors. The major part of the code for the computation of the cone of influence, and the code for Fourier-randomized surrogate time series has been adopted from Huidong Tian.

Wavelet computation, the simulation algorithm, and ridge determination build heavily on the use of matrices in order to minimize computation time in R.

This function provides a broad variety of final as well as intermediate results which can be further analyzed in detail.

Usage

analyze.wavelet(my.data, my.series = 1, loess.span = 0.75, dt = 1, dj = 1/20, 
                lowerPeriod = 2*dt, upperPeriod = floor(nrow(my.data)/3)*dt, 
                make.pval = T, method = "white.noise", params = NULL,
                n.sim = 100, verbose = T)

Arguments

my.data
data frame of time series (including header, and dates as rownames or as separate column named "date" if available)
my.series
name or column index indicating the series to be analyzed, e.g. 1, 2, "dji", "ftse". Default: 1
loess.span
parameter alpha in loess controlling the degree of time series smoothing, if the time series is to be be detrended; no detrending if loess.span=0. Default: 0.75
dt
time resolution, i.e. sampling resolution on time domain, 1/dt = number of intervals per time unit. Default: 1.
dj
frequency resolution, i.e. sampling resolution on frequency domain, 1/dj = number of suboctaves (voices per octave). Default: 1/20.
lowerPeriod
lower Fourier period (in time units) for wavelet decomposition. Default: 2*dt.
upperPeriod
upper Fourier period (in time units) for wavelet decomposition. Default: (floor of one third of time series length)*dt
make.pval
Compute p-values? Logical. Default: TRUE
method
the method of generating surrogate time series, select from: rlll{ "white.noise" : white noise "shuffle" : shuffling the given time series "Fourier.rand" : time series with a similar
params
a list of assignments between methods (AR, and ARIMA) and lists of parameter values to apply to surrogates. Default: NULL. Default includes: AR: AR = list(p=1), where: <
n.sim
number of simulations. Default: 100
verbose
Print verbose output on the screen? Logical. Default: TRUE

Value

  • A list of class analyze.wavelet with the following elements:
  • seriesa data frame with the following columns rllll{ date : the calendar date (as given in my.data) : the series which has been analyzed : (detrended, if loess.span != 0; original names retained) .trend : the trend series (if loess.span != 0) } Row names are taken over from my.data, and so are dates if given as row names.
  • loess.spanparameter alpha in loess which controlled the degree of time series smoothing, if the time series was detrended; detrending was omitted if loess.span=0.
  • dttime resolution, i.e. sampling resolution on time domain, 1/dt = number of intervals per time step.
  • djfrequency resolution, i.e. sampling resolution on frequency domain, 1/dj = number of suboctaves (voices per octave).
  • Wavecomplex wavelet transform of the series
  • Phasephases
  • Amplamplitudes
  • Powerwavelet power in the time/frequency domain
  • Power.avgaverage wavelet power in the frequency domain (averages over time)
  • Power.pvalp-values of wavelet power
  • Power.avg.pvalp-values of average wavelet power
  • Ridgepower ridge, in the form of a 0-1 matrix: columns correspond to dt steps, rows correspond to dj steps whose numerical values are given in Period
  • Periodthe Fourier periods (in time units)
  • Scalethe scales
  • ncnumber of columns/time steps
  • nrnumber of rows/scales/Fourier periods
  • coi.1, coi.2borders of the region where the wavelet transforms are not influenced by edge effects (cone of influence)
  • axis.1tick levels corresponding to time steps
  • axis.2tick levels corresponding to Fourier periods (= log2(Period))

Details

Wavelet transformation, as well as p-value computations, are carried out by calling the internal function wt.

References

Aguiar-Conraria L., and Soares M.J., 2011. The Continuous Wavelet Transform: A Primer. NIPE Working Paper Series 16/2011.

Carmona R., Hwang W.-L., and Torresani B., 1998. Practical Time Frequency Analysis. Gabor and Wavelet Transforms with an Implementation in S. Academic Press, San Diego.

Cazelles B., Chavez M., Berteaux, D., Menard F., Vik J.O., Jenouvrier S., and Stenseth N.C., 2008. Wavelet analysis of ecological time series. Oecologia 156, 287--304.

Liu Y., Liang X.S., and Weisberg R.H., 2007. Rectification of the Bias in the Wavelet Power Spectrum. Journal of Atmospheric and Oceanic Technology 24, 2093--2102.

Tian, H., and Cazelles, B., 2012. WaveletCo. Available at http://cran.r-project.org/src/contrib/Archive/WaveletCo/, archived April 2013; accessed July 26, 2013.

Torrence C., and Compo G.P., 1998. A practical guide to wavelet analysis. Bulletin of the American Meteorological Society 79 (1), 61--78.

See Also

wt.avg, wt.image, wt.sel.phases, wt.phase.image, reconstruct

Examples

Run this code
## The following example is adopted from Liu et al, 2007:

series.length = 6*128*24
x1 = periodic.series(start.period = 1*24, length = series.length)
x2 = periodic.series(start.period = 8*24, length = series.length)
x3 = periodic.series(start.period = 32*24, length = series.length)
x4 = periodic.series(start.period = 128*24, length = series.length)
x = x1 + x2 + x3 + x4

plot(ts(x, start=0, frequency=24), type="l", 
  xlab="time (days)", 
  ylab="hourly data", main="a series of hourly data with periods of 1, 8, 32, and 128 days")
     
my.data = data.frame(x=x)

my.wt = analyze.wavelet(my.data, "x", 
                        loess.span=0, 
                        dt=1/24, dj=1/20, 
                        lowerPeriod=1/4, 
                        make.pval=T, n.sim=10)

## Plot of wavelet power spectrum (with equidistant color breakpoints):  
wt.image(my.wt, color.key="interval", legend.params=list(lab="wavelet power levels"))
## Plot of average wavelet power:
wt.avg(my.wt, siglvl=0.05, sigcol="red")

Run the code above in your browser using DataLab