analyze.wavelet: Computation of the wavelet power spectrum of a single time...
In WaveletComp: Computational Wavelet Analysis

Description Usage Arguments Details Value Author(s) References See Also Examples

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.

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 = TRUE, method = "white.noise", params = NULL,
                n.sim = 100, 
                date.format = NULL, date.tz = NULL, 
                verbose = TRUE)

my.data

data frame of time series (including header, and dates as row names 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 detrended; no detrending if loess.span = 0.

Default: 0.75.

dt

time resolution, i.e. sampling resolution in the time domain, 1/dt = number of observations per time unit. For example: a natural choice of dt in case of hourly data is dt = 1/24, resulting in one time unit equaling one day. This is also the time unit in which periods are measured. If dt = 1, the time interval between two consecutive observations will equal one time unit.

Default: 1.

dj

frequency resolution, i.e. sampling resolution in the frequency domain, 1/dj = number of suboctaves (voices per octave).

Default: 1/20.

lowerPeriod

lower Fourier period (measured in time units determined by dt, see the explanations concerning dt) for wavelet decomposition.
If dt = 1, the minimum admissible value is 2.

Default: 2*dt.

upperPeriod

upper Fourier period (measured in time units determined by dt, see the explanations concerning dt) 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:

`"white.noise"`	:	white noise
`"shuffle"`	:	shuffling the given time series
`"Fourier.rand"`	:	time series with a similar spectrum
`"AR"`	:	AR(p)
`"ARIMA"`	:	ARIMA(p,0,q)

Default: "white.noise".

params

a list of assignments between methods (AR, and ARIMA) and lists of parameter values applying to surrogates. Default: NULL.

Default includes two lists named AR and ARIMA:

AR = list(...), a list containing one single element:

p : AR order.

Default: 1.

ARIMA = list(...), a list of six elements:

`p`	:	AR order.
		Default: `1`.
`q`	:	MA order.
		Default: `1`.
`include.mean`	:	Include a mean/intercept term?
		Default: `TRUE`.
`sd.fac`	:	magnification factor to boost the
		residual standard deviation.
		Default: `1`.
`trim`	:	Simulate trimmed data?
		Default: `FALSE`.
`trim.prop`	:	high/low trimming proportion.
		Default: `0.01`.

n.sim

number of simulations.

Default: 100.

date.format

optional, and for later reference: the format of calendar date (if available in the input data frame) given as a character string, e.g. "%Y-%m-%d", or equivalently "%F"; see strptime for a list of implemented date conversion specifications. Explicit information given here will be overwritten by any later specification given in e.g. wt.image. If unspecified, date formatting will be attempted according to as.Date.

Default: NULL.

date.tz

optional, and for later reference: a character string specifying the time zone of calendar date (if available in the input data frame); see strptime. Explicit information given here will be overwritten by any specification given in e.g. wt.image. If unspecified, "" (the local time zone) will be used.

Default: NULL.

verbose

Print verbose output on the screen? Logical.

Default: TRUE

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

A list of class "analyze.wavelet" with elements of different dimensions. The elements of matrix type (namely, Wave, Phase, Ampl, Power, Power.pval, Ridge) have the following structure:
columns correspond to observations (observation epochs; "epoch" meaning point in time), rows correspond to scales (Fourier periods) whose values are given in Scale (Period). Here is a detailed list of all elements:

series

a data frame with the following columns:

date	:	the calendar date
		(if available as column in `my.data`)
<x>	:	the series which has been analyzed
		(detrended, if `loess.span != 0`;
		original name retained)
<x>.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.span

parameter alpha in loess controlling the degree of time series smoothing if the time series was detrended; no detrending if loess.span = 0

`dt`	time resolution, i.e. sampling resolution in the time domain, `1/dt` = number of observations per time unit
`dj`	frequency resolution, i.e. sampling resolution in the frequency domain, `1/dj` = number of suboctaves (voices per octave)

`Wave`	complex wavelet transform of the series
`Phase`	phases
`Ampl`	amplitudes

`Power`	wavelet power in the time/frequency domain
`Power.avg`	average wavelet power in the frequency domain (averages over time)
`Power.pval`	p-values of wavelet power
`Power.avg.pval`	p-values of average wavelet power

Ridge

wavelet power ridge, in the form of a matrix of 0s and 1s

`Period`	the Fourier periods (measured in time units determined by `dt`, see the explanations concerning `dt`)
`Scale`	the scales (the Fourier periods divided by the Fourier factor)

`nc`	number of columns = number of observations = number of observation epochs; "epoch" meaning point in time
`nr`	number of rows = number of scales (Fourier periods)

coi.1, coi.2

borders of the region where the wavelet transforms are not influenced by edge effects (cone of influence). The coordinates of the borders are expressed in terms of internal axes axis.1 and axis.2.

`axis.1`	tick levels corresponding to the time steps used for (cross-)wavelet transformation: `1, 1+dt, 1+2dt, ...`. The default time axis in plot functions provided by `WaveletComp` is determined by observation epochs, however; "epoch" meaning point in time.
`axis.2`	tick levels corresponding to the log of Fourier periods: `log2(Period)`. This determines the period axis in plot functions provided by `WaveletComp`.

`date.format`	the format of calendar date (if available)
`date.tz`	the time zone of calendar date (if available)

Angi Roesch and Harald Schmidbauer; credits are also due to Huidong Tian, and Bernard Cazelles.

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 https://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.

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

## Not run: 
## 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(x, type = "l", xlab = "index", ylab = "", xaxs = "i",
     main = "hourly series with periods of 1, 8, 32, 128 days")
   
## The following dates refer to the local time zone
## (possibly allowing for daylight saving time):      
my.date <- seq(as.POSIXct("2014-10-14 00:00:00", format = "%F %T"), 
               by = "hour", 
               length.out = series.length)     
my.data <- data.frame(date = my.date, x = x)

## Computation of wavelet power:
## a natural choice of 'dt' in the case of hourly data is 'dt = 1/24',
## resulting in one time unit equaling one day. 
## This is also the time unit in which periods are measured.
## There is an option to store the date format and time zone as additional 
## parameters within object 'my.wt' for later reference.    

my.wt <- analyze.wavelet(my.data, "x", 
                         loess.span = 0, 
                         dt = 1/24, dj = 1/20, 
                         lowerPeriod = 1/4, 
                         make.pval = TRUE, n.sim = 10,
                         date.format = "%F %T", date.tz = "")

## Plot of wavelet power spectrum (with equidistant color breakpoints):  
wt.image(my.wt, color.key = "interval", main = "wavelet power spectrum",
   legend.params = list(lab = "wavelet power levels"),
   periodlab = "period (days)")

## Plot of average wavelet power:
wt.avg(my.wt, siglvl = 0.05, sigcol = "red", 
   periodlab = "period (days)")

## Please see our guide booklet for further examples:
## URL http://www.hs-stat.com/projects/WaveletComp/WaveletComp_guided_tour.pdf.


## End(Not run)