Computational Wavelet Analysis

Description

Wavelet analysis and reconstruction of time series, cross-wavelets and phase difference (with filtering options), significance with bootstrap algorithms.

Details

Periodic phenomena of a single time series can be analyzed with function analyze.wavelet. Results of the analysis (a time/period image of the wavelet power spectrum, plots of the average power, and phase plots for selected periods and a time/period image of phases) can be accessed through various plot functions (wt.image, wt.avg, wt.sel.phases, wt.phase.image). Function reconstruct returns the reconstructed time series where reconstruction is according to constraints on significance, period specification, and cone of influence.

The cross-wavelet spectrum and coherency spectrum of two time series can be analyzed with function analyze.coherency. Results (a time/period image of cross-wavelet power or coherency, plots of averages, plots of phases and phase differences for selected periods and the time/period image of phase differences) can be accessed through corresponding functions (wc.image, wc.avg, wc.sel.phases, wc.phasediff.image).

Detrending of the time series at hand is offered as an option. Wavelet transformations are computed using the Morlet wavelet. Smoothing filters are provided in the case of cross-wavelet transformation to compute wavelet coherency.

Significance is assessed with simulation algorithms, a variety of alternative hypotheses to test is available, for which surrogate time series are provided: white noise, shuffling the given time series, time series with a similar spectrum, AR, and ARIMA.

Names and parts of the layout of some routines were inspired by similar functions 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. The implementation of a choice of filtering windows for the computation of the wavelet coherence was inspired by Luis Aguiar-Conraria and Maria Joana Soares (GWPackage).

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

What is new in WaveletComp version 1.2?

Nothing really new, just some formal changes in the URLs, as required. Fully downward compatible. In addition, we have paved the way for the next version.

What is new in WaveletComp version 1.1?

Tools for displaying and analyzing periodic phenomena across time have been extended. The main innovations are:

• All functions of family wt.<> (showing results concerning a single time series) can now also be applied to extract univariate outcomes from cross-wavelet and coherence analysis (objects of class "analyze.coherency").

• It is possible to control the color gradation of time-period spectrum plots, and accentuate the contrast, by raising the wavelet power values to any (positive) exponent before plotting.

• Setting a maximum level for the color bar facilitates the visual comparison of time-period spectrum plots. Maximum and minimum plot levels are options for plots of averages too.

• The time and period axes are now easier to individualize by specifying tick marks and labels. Coordinates on the time axis can be conveniently addressed via an index or a "POSIXct" object.

• Graphical parameters of global coverage (cex.axis, font.axis, cex.lab, font.lab, mgp etc., see par) as well as parameters of local coverage (within axis specification options) help fine-tune plots.

• Two more real-world data sets have been included in WaveletComp, namely:

  • Data set "weather.radiation.Mannheim", containing daily weather and ambient radiation readings from Mannheim (Germany).

  • Data set "USelection2016.Instagram", containing hourly numbers of candidate-related media uploads to Instagram right before the 2016 US presidential election.

Author(s)

Angi Roesch and Harald Schmidbauer; credits are also due to Huidong Tian, Bernard Cazelles, Luis Aguiar-Conraria, and Maria Joana Soares.

References

Aguiar-Conraria L., and Soares M.J., 2011. Business cycle synchronization and the Euro: A wavelet analysis. Journal of Macroeconomics 33 (3), 477–489.

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

Aguiar-Conraria L., and Soares M.J., 2012. GWPackage.
https://sites.google.com/site/aguiarconraria/joanasoares-wavelets (as of April 29, 2025, no longer online, package defunct); accessed September 4, 2013.

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 P.C., 1994. Wavelet spectrum analysis and ocean wind waves. In: Foufoula-Georgiou E., and Kumar P., (eds.), Wavelets in Geophysics, Academic Press, San Diego, 151–166.

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.

Schmidbauer H., Roesch A., Stieler F., 2018. The 2016 US presidential election and media on Instagram: Who was in the lead? Computers in Human Behavior 81, 148–160. doi:10.1016/j.chb.2017.11.021

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.

Veleda D., Montagne R., and Araujo M., 2012. Cross-Wavelet Bias Corrected by Normalizing Scales. Journal of Atmospheric and Oceanic Technology 29, 1401–1408.

Simulation of an AR(p) surrogate for a given time series x

Description

It simulates an AR(p) surrogate for the time series x analyzed by wavelet transformation using either function analyze.wavelet or function analyze.coherency. A set of surrogates is used for significance assessment to test the hypothesis of equal periodic components.

Simulation is subject to coefficient estimates referring to the given AR order.

Usage

AR(x, params = list(AR = list(p = 1)))

Arguments

Details

Basics of model estimation and simulation are provided by application of the R functions pcf and arima.sim.

Value

A surrogate series for x is returned which has the same length and properties according to estimates resulting from the given AR order.

Author(s)

Angi Roesch and Harald Schmidbauer

See Also

analyze.wavelet, analyze.coherency, SurrogateData, ARIMA, FourierRand

Simulation of an ARIMA(p,0,q) surrogate for a given time series x

Description

It simulates an ARIMA(p,0,q) surrogate for the time series x analyzed by wavelet transformation using either function analyze.wavelet or function analyze.coherency. A set of surrogates is used for significance assessment to test the hypothesis of equal periodic components.

Simulation is subject to parameters concerning: AR and MA order, an indicator of mean/intercept term inclusion, a boosting factor which is applied to the residual standard deviation, an indicator to trim the data according to the given trimming proportion (i.e. to restrict them to the interval between the lower and upper quantile corresponding to the given probability).

Usage

ARIMA(x, params = 
      list(ARIMA = list(p = 1, q = 1, 
                        include.mean = TRUE, sd.fac = 1, 
                        trim = FALSE, trim.prop = 0.01)))

Arguments

Details

Basics of model estimation and simulation are provided by application of the R functions arima and arima.sim.

Value

A surrogate series for x is returned which has the same length and properties according to estimates resulting from the given ARIMA setting.

Author(s)

Angi Roesch and Harald Schmidbauer

See Also

analyze.wavelet, analyze.coherency, SurrogateData, AR, FourierRand

Computation of the cone of influence

Description

Wavelet transformation of a time series includes padding with zeros at the edges. This function determines the spectrum borders influenced by edge effects.

This code has been developed by Huidong Tian (archived R package WaveletCo).

Usage

COI(start = start, dt = dt, nc = nc, nr = nr, Period = Period)

Arguments

Value

a list with the following elements:

Author(s)

Code developed by Huidong Tian.

References

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.

See Also

analyze.wavelet, analyze.coherency, wt.image, wc.image

Series of FX trade transactions

Description

USD/euro FX (foreign exchange) trade: number of worldwide transactions recorded per 5-minute intervals in July 2010. The data set contains four full weekly cycles (plus three days at the beginning of July 2010), a weekly cycle lasting from Sunday, 21:00, to Friday, 20:55. The number of transactions between Friday, 21:00, and Sunday, 20:55, is 0 or close to 0. For these intervals, variable "active" is FALSE, otherwise TRUE.

Derived from data delivered by Morning Star.

Usage

data("FXtrade.transactions")

Format

A data frame of two columns:

Source

Morning Star, https://www.morningstar.com (limited access)

Examples

data(FXtrade.transactions)
plot(as.POSIXct(FXtrade.transactions$date, format = "%F %T", tz = "GMT"),
     FXtrade.transactions$transactions, 
     type = "l",
     xlab = "day", ylab = "transactions in 5-minute intervals")

Simulation of a surrogate for a given time series x by Fourier randomization

Description

It simulates a surrogate with similar Fourier spectrum for a given time series x analyzed by wavelet transformation using either function analyze.wavelet or function analyze.coherency. A set of surrogates is used for significance assessment to test the hypothesis of equal periodic components.

This code has been developed by Huidong Tian (archived R package WaveletCo).

Usage

FourierRand(x)

Arguments

Value

A surrogate series for x is returned which has the same length and a similar Fourier spectrum.

Author(s)

Code developed by Huidong Tian.

References

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.

See Also

analyze.wavelet, analyze.coherency, SurrogateData, AR, ARIMA

Simulation of surrogates for a given time series x, subject to the specified method and parameters

Description

It simulates a surrogate for the time series x to be analyzed by wavelet transformation using either function analyze.wavelet or function analyze.coherency. A set of surrogates is used for significance assessment to test the hypothesis of equal periodic components.

Simulation is subject to model/method specification and parameter setting: Currently, one can choose from a variety of 6 methods (white noise, series shuffling, Fourier randomization, AR, and ARIMA) with respective lists of parameters to set.

The name and layout were inspired by a similar function developed by Huidong Tian (archived R package WaveletCo).

Usage

SurrogateData(x, method = "white.noise", params = list( 
         AR    = list(p = 1), 
         ARIMA = list(p = 1, q = 1, include.mean = TRUE, sd.fac = 1, 
                      trim = FALSE, trim.prop = 0.01)))

Arguments

Value

A surrogate series for x is returned which has the same length and properties according to estimates resulting from the model/method specification and parameter setting.

Author(s)

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

References

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.

See Also

analyze.wavelet, analyze.coherency, AR, ARIMA, FourierRand

Hourly time series of the number of candidate-related media posted on Instagram during the week before the 2016 US presidential election

Description

One week (Sunday, 2016-10-30 23:00:00 EDT through Sunday, 2016-11-06 23:00:00 EST; 170 hours) of hourly readings of the number of media posted on Instagram, supposedly (according to hashtags with which they were annotated) positive/neutral or in opposition towards candidates Trump and Clinton. The 2016 US presidential election took place on Tuesday, 2016-11-08.

Usage

data("USelection2016.Instagram")

Format

A data frame of four columns:

Details

Media posted on Instagram are usually annotated with hashtags. A hashtag can be used to determine whether a candidate-related posting is positive/neutral towards a candidate (for example, #makeamericagreatagain for Trump, #hillary2016 for Clinton) or opposing a candidate (for example, #dumptrump for Trump, #neverhillary for Clinton). In this way, four hourly time series are obtained: Trump vs. Clinton, supporters vs. opponents. For further details, see the reference below. The date column also contains a label EDT (Eastern Daylight Time) or EST (Eastern Standard Time); daylight saving time ended 2016-11-06 at 1:00 a.m. when clocks were moved back to 1:00 a.m. EST. The time stamp "2016-11-06 02:00:00" therefore occurs twice, once with EDT and once with EST.

Source

Hourly readings of the number of media posted on Instagram were collected using Instagram's built-in API.

References

Schmidbauer H., Roesch A., Stieler F., 2018. The 2016 US presidential election and media on Instagram: Who was in the lead? Computers in Human Behavior 81, 148–160. doi:10.1016/j.chb.2017.11.021

Examples

data(USelection2016.Instagram)
attach(USelection2016.Instagram)

my.date <- as.POSIXct(date, format = "%F %T", tz = "EST5EDT")

plot(my.date, trump.pos, type = "l", col = 1, lwd = 2,
     ylab = "number of media posted on Instagram", ylim = c(0,6e+6),
     xlab = "the week before the Election Day (Tuesday, 2016-11-08)")
lines(my.date, clinton.pos, col = 2, lwd = 2) 
lines(my.date, trump.neg, col = 3, lwd = 2)
lines(my.date, clinton.neg, col = 4, lwd = 2) 
legend("topleft", legend=names(USelection2016.Instagram[-1]), 
       lty = 1, lwd = 2, col = 1:4) 

detach(USelection2016.Instagram)

Computation of the cross-wavelet power and wavelet coherence of two time series

Description

Given two time series x and y having the same length and timestamp, this function computes the cross-wavelet power and wavelet coherence applying the Morlet wavelet, subject to criteria concerning: the time and frequency resolution, an (optional) lower and/or upper Fourier period, and filtering method for the coherence computation.

The output is further processed by the higher-order function wc and can be retrieved from analyze.coherency.

The name and layout were inspired by a similar function developed by Huidong Tian and Bernard Cazelles (archived R package WaveletCo). The implementation of a choice of filtering windows for the computation of the wavelet coherence was inspired by Luis Aguiar-Conraria and Maria Joana Soares (GWPackage).

Usage

WaveletCoherency(x, y, dt = 1, dj = 1/20, 
                 lowerPeriod = 2*dt, upperPeriod = floor(length(x)*dt/3), 
                 window.type.t = 1, window.type.s = 1, 
                 window.size.t = 5, window.size.s = 1/4)

Arguments

Value

A list with the following elements:

Author(s)

Angi Roesch and Harald Schmidbauer; credits are also due to Huidong Tian, Bernard Cazelles, Luis Aguiar-Conraria, and Maria Joana Soares.

References

Aguiar-Conraria L., and Soares M.J., 2011. Business cycle synchronization and the Euro: A wavelet analysis. Journal of Macroeconomics 33 (3), 477–489.

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

Aguiar-Conraria L., and Soares M.J., 2012. GWPackage.
https://sites.google.com/site/aguiarconraria/joanasoares-wavelets (as of April 29, 2025, no longer online, package defunct); accessed September 4, 2013.

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 P.C., 1994. Wavelet spectrum analysis and ocean wind waves. In: Foufoula-Georgiou E., and Kumar P., (eds.), Wavelets in Geophysics, Academic Press, San Diego, 151–166.

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.

Veleda D., Montagne R., and Araujo M., 2012. Cross-Wavelet Bias Corrected by Normalizing Scales. Journal of Atmospheric and Oceanic Technology 29, 1401–1408.

See Also

WaveletTransform, wc, analyze.coherency

Morlet wavelet transformation of a time series

Description

It computes the Morlet wavelet transformation of a given time series, subject to criteria concerning: the time and frequency resolution, an (optional) lower and/or upper Fourier period.

The output is further processed by higher-order functions wt, WaveletCoherency and wc, and can be retrieved from analyze.wavelet and analyze.coherency.

The name and layout were inspired by a similar function developed by Huidong Tian and Bernard Cazelles (archived R package WaveletCo).

Usage

WaveletTransform(x, dt = 1, dj = 1/20, 
                 lowerPeriod = 2*dt, upperPeriod = floor(length(x)*dt/3))

Arguments

Value

A list of class analyze.wavelet with the following elements:

Author(s)

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

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 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.

See Also

wt, WaveletCoherency, wc, analyze.wavelet, analyze.coherency

Computation of the cross-wavelet power and wavelet coherence spectrum of two time series

Description

The two time series are selected from an input data frame by specifying either their names or their column numbers. Optionally, the time series are detrended, using loess with parameter loess.span. Internally, the series will be standardized before they undergo wavelet transformation.

The cross-wavelet power spectrum is computed 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.

For the computation of wavelet coherence, a variety of filtering methods is provided, with flexible window parameters.

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

The name and parts of the layout of subroutine wc 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. The implementation of a choice of filtering windows for the computation of the wavelet coherence was inspired by Luis Aguiar-Conraria and Maria Joana Soares (GWPackage).

Cross-wavelet and coherence 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.coherency(my.data, my.pair = c(1, 2), loess.span = 0.75, 
                  dt = 1, dj = 1/20, 
                  lowerPeriod = 2*dt, 
                  upperPeriod = floor(nrow(my.data)/3)*dt, 
                  window.type.t = 1, window.type.s = 1, 
                  window.size.t = 5, window.size.s = 1/4, 
                  make.pval = TRUE, method = "white.noise", params = NULL, 
                  n.sim = 100, 
                  date.format = NULL, date.tz = NULL, 
                  verbose = TRUE)

Arguments

Value

A list of class "analyze.coherency" with elements of different dimensions. The elements of matrix type, namely:

• Wave.xy, Angle, sWave.xy, sAngle,

• Power.xy, Power.xy.pval,

• Coherency, Coherence, Coherence.pval,

• Wave.x, Wave.y, Phase.x, Phase.y, Ampl.x, Ampl.y,

• Power.x, Power.y, Power.x.pval, Power.y.pval, sPower.x, sPower.y,

• Ridge.xy, Ridge.co, Ridge.x, Ridge.y,

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:

Author(s)

Angi Roesch and Harald Schmidbauer; credits are also due to Huidong Tian, Bernard Cazelles, Luis Aguiar-Conraria, and Maria Joana Soares.

References

Aguiar-Conraria L., and Soares M.J., 2011. Business cycle synchronization and the Euro: A wavelet analysis. Journal of Macroeconomics 33 (3), 477–489.

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

Aguiar-Conraria L., and Soares M.J., 2012. GWPackage.
https://sites.google.com/site/aguiarconraria/joanasoares-wavelets (as of April 29, 2025, no longer online, package defunct); accessed September 4, 2013.

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 P.C., 1994. Wavelet spectrum analysis and ocean wind waves. In: Foufoula-Georgiou E., and Kumar P., (eds.), Wavelets in Geophysics, Academic Press, San Diego, 151–166.

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.

Veleda D., Montagne R., and Araujo M., 2012. Cross-Wavelet Bias Corrected by Normalizing Scales. Journal of Atmospheric and Oceanic Technology 29, 1401–1408.

See Also

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

Examples

## Not run: 
## The following example is modified from Veleda et al, 2012:

series.length <- 3*128*24
x1 <- periodic.series(start.period = 1*24, length = series.length)
x2 <- periodic.series(start.period = 2*24, length = series.length)
x3 <- periodic.series(start.period = 4*24, length = series.length)
x4 <- periodic.series(start.period = 8*24, length = series.length)
x5 <- periodic.series(start.period = 16*24, length = series.length)
x6 <- periodic.series(start.period = 32*24, length = series.length)
x7 <- periodic.series(start.period = 64*24, length = series.length)
x8 <- periodic.series(start.period = 128*24, length = series.length)

x <- x1 + x2 + x3 + x4 + 3*x5 + x6 + x7 + x8 + rnorm(series.length)
y <- x1 + x2 + x3 + x4 - 3*x5 + x6 + 3*x7 + x8 + rnorm(series.length)

matplot(data.frame(x, y), type = "l", lty = 1, xaxs = "i", col = 1:2, 
 xlab = "index", ylab = "",
 main = "hourly series with periods of 1, 2, 4, 8, 16, 32, 64, 128 days", 
 sub = "(out of phase at period 16, different amplitudes at period 64)")
legend("topright", legend = c("x","y"), col = 1:2, lty = 1)

## 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, y = y)

## Computation of cross-wavelet power and wavelet coherence, x over y:
## 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.wc' for later reference. 

my.wc <- analyze.coherency(my.data, c("x","y"), 
                           loess.span = 0, 
                           dt = 1/24, dj = 1/20, 
                           window.size.t = 1, window.size.s = 1/2, 
                           lowerPeriod = 1/4,
                           make.pval = TRUE, n.sim = 10,
                           date.format = "%F %T", date.tz = "")
## Note:                           
## By default, Bartlett windows are used for smoothing in order to obtain
## the wavelet coherence spectrum; window lengths in this example:
## 1*24 + 1 = 25 observations in time direction,
## (1/2)*20 + 1 = 11 units in scale (period) direction.                             
                         
## Plot of cross-wavelet power 
## (with color breakpoints according to quantiles):
wc.image(my.wc, main = "cross-wavelet power spectrum, x over y",
   legend.params = list(lab = "cross-wavelet power levels"),
   periodlab = "period (days)")
   
## The same plot, now with calendar axis
## (according to date format stored in 'my.wc'):
wc.image(my.wc, main = "cross-wavelet power spectrum, x over y",
   legend.params = list(lab = "cross-wavelet power levels"),
   periodlab = "period (days)", show.date = TRUE)   
         
## Plot of average cross-wavelet power:
wc.avg(my.wc, siglvl = 0.05, sigcol = 'red', 
   periodlab = "period (days)")

## Plot of wavelet coherence 
## (with color breakpoints according to quantiles):
wc.image(my.wc, which.image = "wc",  main = "wavelet coherence, x over y", 
   legend.params = list(lab = "wavelet coherence levels", 
                        lab.line = 3.5, label.digits = 3),
   periodlab = "period (days)")
         
## plot of average coherence:
wc.avg(my.wc, which.avg = "wc", 
   siglvl = 0.05, sigcol = 'red', 
   legend.coords = "topleft", 
   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)

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

Arguments

Details

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

Value

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:

Author(s)

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

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 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.

See Also

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

Examples

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

Plot an arrow

Description

It plots a polygon arrow subject to criteria concerning position, length, width, direction and color.

In plots of the cross-wavelet power (or wavelet coherency) spectrum of two time series as performed by function wc.image, arrows are used as a symbolic indicator of phase differences between the series in the time/frequency domain.

This code has been developed by Huidong Tian (archived R package WaveletCo).

Usage

arrow(x, y, l = 0.1, w = 0.3*l, alpha, col.arrow = "black")

Arguments

Author(s)

Code developed by Huidong Tian.

References

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.

Series of monthly marriages in Turkey

Description

Series of monthly marriages in Turkey (1988-2013), as reported by DIE (Devlet Istatistik Enstitusu) / TUIK (Turkiye Istatistik Kurumu).

Usage

data("marriages.Turkey")

Format

A data frame of three columns:

Source

DIE (Devlet Istatistik Enstitusu) / TUIK (Turkiye Istatistik Kurumu)

Jan 1988 to Dec 2000:
"Evlenme istatistikleri", DIE (Devlet Istatistik Enstitusu, Ankara), ISSN: 1300-1086; several issues.

Jan 2001 to Dec 2013:
https://www.tuik.gov.tr/; accessed Oct 1, 2014.

Examples

data(marriages.Turkey)
plot(ts(marriages.Turkey$marriages, frequency = 12, start = c(1988,1)), 
     type = "l", 
     xlab = "", ylab = "", 
     main = "monthly marriages in Turkey (1988-2013)")     

Computation of a (deterministic) periodic time series of linearly changing period.

Description

It computes and returns a sinusoid of a specified length, which has the given initial phase, and linearly changing periods (if requested) starting from a given period length through the given length at the end. There is an option to plot the time series.

Usage

  periodic.series(start.period = 100, end.period = start.period, 
                  phase = 0, 
                  length = 600, 
                  make.plot = FALSE)

Arguments

Details

This function can be used for illustrating methods and functions.

Producing a sinusoid, periodic.series will work best if start.period (and end.period, if different from start.period) is not too small.

Value

the series as vector

Author(s)

Angi Roesch and Harald Schmidbauer

See Also

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

Examples

# The following time series involves periods from 100 through 50:
x <- periodic.series(start.period = 100, end.period = 50, make.plot = TRUE)
title("time series with period changing linearly from 100 to 50")

# The following three time series involve three different types of period evolution 
# starting from period 100:
x1 <- 0.8*periodic.series(start.period = 100, end.period = 95, phase = 0, length = 1000)
x2 <-     periodic.series(start.period = 100, end.period = 100, phase = 0, length = 1000)
x3 <- 1.2*periodic.series(start.period = 100, end.period = 105, phase = 0, length = 1000)

ts.plot(x2, ylim = c(-2, +2), xlab = "time", ylab = "series with variable period")
lines(x1, col = "blue")
lines(x3, col = "red")
legend("topleft", 
   legend = c("speeding up (end period = 95)", "period = 100", 
              "slowing down (end period = 105)"), 
   lty = 1, col = c("blue", "black", "red"))

Reconstruction of a (detrended) time series from output provided by an object of class "analyze.wavelet" or "analyze.coherency"

Description

This function reconstructs a (detrended) time series analyzed by wavelet transformation using either function analyze.wavelet or function analyze.coherency, subject to optional criteria concerning: minimum wavelet power, significance of wavelet power at a given significance level, specification of (Fourier) periods or period bands, exclusive use of the power ridge and/or the cone of influence. An option is provided to prevent the reconstructed series from final rescaling (applying the original (detrended) series' mean and standard deviation).

(If the object provided as input is of class "analyze.coherency", then the number or name of the time series must be specified.)

Optional: plot of wavelets used for reconstruction, plot of reconstructed series against original (detrended) series. An option is given to individualize the time axis by specifying tick marks and labels.

Output includes the original (detrended) and the reconstructed time series, along with reconstruction waves and parameters.

Usage

reconstruct(WT, my.series = 1, lvl = 0, 
      only.coi = FALSE, 
      only.sig = TRUE, siglvl = 0.05, 
      only.ridge = FALSE, 
      sel.period = NULL, sel.lower = NULL, sel.upper = NULL,  
      rescale = TRUE,
      plot.waves = FALSE, plot.rec = TRUE, 
      lty = 1, lwd = 1, col = 1:2, ylim = NULL,
      show.legend = TRUE, 
      legend.coords = "topleft", legend.horiz = FALSE, legend.text = NULL,
      label.time.axis = TRUE, 
      show.date = FALSE, date.format = NULL, date.tz = NULL,
      timelab = NULL, timetck = 0.02, timetcl = 0.5,
      spec.time.axis = list(at = NULL, labels = TRUE, 
                            las = 1, hadj = NA, padj = NA),
      main.waves = NULL, main.rec = NULL, main = NULL, 
      lwd.axis = 1, 
      verbose = TRUE)

Arguments

Value

A list of class reconstruct with the following elements:

Author(s)

Angi Roesch and Harald Schmidbauer

References

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.

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.

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

analyze.wavelet, wt.image, wt.avg, wt.sel.phases, wt.phase.image, analyze.coherency, wc.image, wc.avg, wc.sel.phases, wc.phasediff.image

Examples

## 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")
  
my.data <- data.frame(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.
my.w <- analyze.wavelet(my.data, "x", 
                       loess.span = 0, 
                       dt = 1/24, dj = 1/20, 
                       lowerPeriod = 1/4, 
                       make.pval = TRUE, n.sim = 10)

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

## Reconstruction of the time series, 
## including significant components only: 
reconstruct(my.w)

## The same reconstruction, but showing wave components first:
reconstruct(my.w, plot.waves = TRUE)

## Reconstruction, including all components whether significant or not: 
reconstruct(my.w, only.sig = FALSE)

## Reconstruction, including significant components, 
## but selected periods only (e.g. ignoring period 8):  
reconstruct(my.w, sel.period = c(1,32,128))

## Reconstruction, including significant components,
## but the ridge only:
reconstruct(my.w, only.ridge = TRUE)

## Alternate styles of the time axis:

## The plot with time elapsed in days, starting from 0 and proceeding 
## in steps of 50 days (50*24 hours),
## instead of the (default) time index:
index.ticks  <- seq(1, series.length, by = 50*24)
index.labels <- (index.ticks-1)/24

## Insert your specification of time axis: 
reconstruct(my.w, only.ridge = TRUE,
   timelab = "time elapsed (days)",
   spec.time.axis = list(at = index.ticks, labels = index.labels))
            
## See the periods involved:
my.rec <- reconstruct(my.w, only.ridge = TRUE)
print(my.rec$Period[my.rec$rnum.used])

## The original and reconstructed time series can be retrieved:
plot(my.rec$series$x, type = "l", xlab = "index", ylab = "")
lines(my.rec$series$x.r, col="red")
legend("topleft", legend = c("original","reconstructed"), 
       lty = 1, col = c("black","red"))

## Please see also the examples in our guide booklet,
## URL http://www.hs-stat.com/projects/WaveletComp/WaveletComp_guided_tour.pdf.


## End(Not run)

Ridge of a wavelet power or coherence spectrum

Description

It finds the ridge of a (cross-)wavelet power or coherence spectrum, subject to criteria concerning: the band of neighboring values to be included in the search area for local maximization, and a scale factor applied to the global maximum defining a minimum level constraint.

The basic concept of ridge determination builds on ideas developed by Huidong Tian and Bernard Cazelles (archived R package WaveletCo). This code is an efficient implementation of this concept, building heavily on the use of matrices in order to minimize computation time in R.

Usage

ridge(wavelet.spectrum, band = 5, scale.factor = 0.1)

Arguments

Value

A matrix of 0s and 1s is returned which determines the course of (cross-)wavelet power or coherence ridge in the time/frequency domain.

Note

Credits are also due to Huidong Tian and Bernard Cazelles whose basic concept of ridge determination has been adopted.

Author(s)

Angi Roesch and Harald Schmidbauer

References

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.

Cross-wavelet transformation, wavelet coherence computation, and a simulation algorithm

Description

This function provides Morlet cross-wavelet transformation results of the given two time series, performed within the lower-order functions WaveletCoherency and WaveletTransform subject to criteria concerning the time and frequency resolution, an (optional) lower and/or upper Fourier period, and a variety of filtering methods for the coherence computation. It performs a simulation algorithm to test against a specified alternative hypothesis, which can be chosen from a variety of options, and provides p-values. The selected model will be fitted to the data and simulated according to estimated parameters in order to provide surrogate time series. This function is called by function analyze.coherency.

The name and parts of the layout were inspired by a similar function developed by Huidong Tian and Bernard Cazelles (archived R package WaveletCo). The major part of the code for the computation of the cone of influence and the code for Fourier-randomized surrogate time series have been adopted from Huidong Tian. The implementation of a choice of filtering windows for the computation of the wavelet coherence was inspired by Luis Aguiar-Conraria and Maria Joana Soares (GWPackage).

Usage

wc(x, y, start = 1, dt = 1, dj = 1/20, 
   lowerPeriod = 2*dt, upperPeriod = floor(length(x)/3)*dt, 
   window.type.t = 1, window.type.s = 1, 
   window.size.t = 5, window.size.s = 1/4,
   make.pval = TRUE, method = "white.noise", params = NULL, 
   n.sim = 100, save.sim = FALSE)

Arguments

Value

A list with the following elements:

Author(s)

Angi Roesch and Harald Schmidbauer; credits are also due to Huidong Tian, Bernard Cazelles, Luis Aguiar-Conraria, and Maria Joana Soares.

References

Aguiar-Conraria L., and Soares M.J., 2011. Business cycle synchronization and the Euro: A wavelet analysis. Journal of Macroeconomics 33 (3), 477–489.

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

Aguiar-Conraria L., and Soares M.J., 2012. GWPackage.
https://sites.google.com/site/aguiarconraria/joanasoares-wavelets (as of April 29, 2025, no longer online, package defunct); accessed September 4, 2013.

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 P.C., 1994. Wavelet spectrum analysis and ocean wind waves. In: Foufoula-Georgiou E., and Kumar P., (eds.), Wavelets in Geophysics, Academic Press, San Diego, 151–166.

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.

Veleda D., Montagne R., and Araujo M., 2012. Cross-Wavelet Bias Corrected by Normalizing Scales. Journal of Atmospheric and Oceanic Technology 29, 1401–1408.

Plot arrows into a cross-wavelet or wavelet coherency spectrum

Description

It adds arrows as symbolic indicators of phase-differences into a cross-wavelet spectrum or wavelet coherency spectrum of two time series.

In particular, the area to be filled with arrows can be determined in several ways: to reflect significance (at a given level) with respect to cross-wavelet power, wavelet coherency, or individual wavelet power, and/or to flag a high-value region ("high" according to a given level) with respect to cross-wavelet power (coherency values, respectively).

There is an option to choose "smoothed" arrows (phase differences) for plotting; these have been computed from smoothing filters as defined in analyze.coherency.

The name and layout were inspired by a similar function developed by Huidong Tian and Bernard Cazelles (archived R package WaveletCo). The code for the arrow design to reflect phase differences has been adopted from Huidong Tian.

Usage

wc.angle(WC = WC, exponent = 1, use.sAngle = TRUE, p = 1, 
         which.lvl = "wp", lvl = 0, 
         which.sig = which.lvl, siglvl = 0.05, 
         col.arrow = "black")

Arguments

Author(s)

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

References

Aguiar-Conraria L., and Soares M.J., 2011. Business cycle synchronization and the Euro: A wavelet analysis. Journal of Macroeconomics 33 (3), 477–489.

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

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 P.C., 1994. Wavelet spectrum analysis and ocean wind waves. In: Foufoula-Georgiou E., and Kumar P., (eds.), Wavelets in Geophysics, Academic Press, San Diego, 151–166.

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.

Veleda D., Montagne R., and Araujo M., 2012. Cross-Wavelet Bias Corrected by Normalizing Scales. Journal of Atmospheric and Oceanic Technology 29, 1401–1408.

See Also

analyze.coherency, wc.image, arrow

Plot cross-wavelet power averages and wavelet coherence averages across time of two time series

Description

This function plots cross-wavelet power averages across time, or alternatively wavelet coherence averages, of two time series, which are provided by an object of class "analyze.coherency". The vertical axis shows the Fourier periods. The horizontal axis shows the averages. User-defined minimum and maximum levels can be applied to cross-wavelet power averages, minimum levels can be applied to coherence averages. Also, an option is given to individualize the period axis and/or axis of averages by specifying tick marks and labels.

There is an option to label periods according to significance of averages (if p-values are provided) at given levels of significance. Labels are point symbols along the line of averages which can be assigned individually.

The idea to show significance levels by colors of plotting characters and its implementation has been adopted from Huidong Tian and Bernard Cazelles (archived R package WaveletCo).

Usage

wc.avg(WC, which.avg = "wp", exponent = 1,
   show.siglvl = TRUE, 
   siglvl = c(0.05, 0.1), 
   sigcol = c("red", "blue"), sigpch = 20, sigcex = 1,
   minimum.level = NULL, maximum.level = NULL,
   label.avg.axis = TRUE, 
   averagelab = NULL, averagetck = 0.02, averagetcl = 0.5,
   spec.avg.axis = list(at = NULL, labels = TRUE, 
                        las = 1, hadj = NA, padj = NA),
   label.period.axis = TRUE, 
   periodlab = NULL, periodtck = 0.02, periodtcl = 0.5,
   spec.period.axis = list(at = NULL, labels = TRUE, 
                           las = 1, hadj = NA, padj = NA),
   show.legend = TRUE, legend.coords = "topright", 
   main = NULL, 
   lwd = 1, col = 1, 
   lwd.axis = 1,
   verbose = FALSE)

Arguments

Author(s)

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

References

Aguiar-Conraria L., and Soares M.J., 2011. Business cycle synchronization and the Euro: A wavelet analysis. Journal of Macroeconomics 33 (3), 477–489.

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

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 P.C., 1994. Wavelet spectrum analysis and ocean wind waves. In: Foufoula-Georgiou E., and Kumar P., (eds.), Wavelets in Geophysics, Academic Press, San Diego, 151–166.

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.

Veleda D., Montagne R., and Araujo M., 2012. Cross-Wavelet Bias Corrected by Normalizing Scales. Journal of Atmospheric and Oceanic Technology 29, 1401–1408.

See Also

analyze.coherency, wc.image, wc.sel.phases, wc.phasediff.image, wt.image, wt.avg,
wt.sel.phases, wt.phase.image, reconstruct

Examples

## Not run: 
## The following example is modified from Veleda et al., 2012:

series.length <- 3*128*24
x1 <- periodic.series(start.period = 1*24, length = series.length)
x2 <- periodic.series(start.period = 2*24, length = series.length)
x3 <- periodic.series(start.period = 4*24, length = series.length)
x4 <- periodic.series(start.period = 8*24, length = series.length)
x5 <- periodic.series(start.period = 16*24, length = series.length)
x6 <- periodic.series(start.period = 32*24, length = series.length)
x7 <- periodic.series(start.period = 64*24, length = series.length)
x8 <- periodic.series(start.period = 128*24, length = series.length)

x <- x1 + x2 + x3 + x4 + 3*x5 + x6 + x7 + x8 + rnorm(series.length)
y <- x1 + x2 + x3 + x4 - 3*x5 + x6 + 3*x7 + x8 + rnorm(series.length)

matplot(data.frame(x, y), type = "l", lty = 1, xaxs = "i", col = 1:2, 
 xlab = "index", ylab = "",
 main = "hourly series with periods of 1, 2, 4, 8, 16, 32, 64, 128 days", 
 sub = "(out of phase at period 16, different amplitudes at period 64)")
legend("topright", legend = c("x","y"), col = 1:2, lty = 1)

## 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, y = y)

## Computation of cross-wavelet power and wavelet coherence, x over y:
## 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.
my.wc <- analyze.coherency(my.data, c("x","y"), 
                           loess.span = 0, 
                           dt = 1/24, dj = 1/20, 
                           window.size.t = 1, window.size.s = 1/2, 
                           lowerPeriod = 1/4,
                           make.pval = TRUE, n.sim = 10)

## Plot of cross-wavelet power, 
## with color breakpoints according to quantiles:
wc.image(my.wc, main = "cross-wavelet power spectrum, x over y",
   legend.params = list(lab = "cross-wavelet power levels (quantiles)"),
   periodlab = "period (days)")
## Note:
## The default time axis shows an index of given points in time, 
## which is the count of hours in our example.      
## By default, arrows are plotted which show the phase differences 
## of x over y at respective significant periods. 
## (Please see our guide booklet for further explanation.)
         
## With time elapsed in days 
## (starting from 0 and proceeding in steps of 50 days)
## instead of the (default) time index:
index.ticks  <- seq(1, series.length, by = 50*24)
index.labels <- (index.ticks-1)/24

## Insert your specification of the time axis: 
wc.image(my.wc, color.key = "i", 
   main = "cross-wavelet power spectrum, x over y", 
   legend.params = list(lab = "cross-wavelet power levels (quantiles)"),
   periodlab = "period (days)", timelab = "time elapsed (days)",
   spec.time.axis = list(at = index.ticks, labels = index.labels))
       
## Plot of average cross-wavelet power:
wc.avg(my.wc, siglvl = 0.05, sigcol = "red", periodlab = "period (days)")

## The same plot, but with enhanced symbol size, user-defined axes, 
## minimum and a maximum plot level of averages:
wc.avg(my.wc, siglvl = 0.05, sigcol = "red", sigcex = 1.5, 
   minimum.level = 0, maximum.level = 17, 
   periodlab = "period (days)",
   spec.period.axis = list(at = c(1,2,4,8,16,32,64,128)), 
   spec.avg.axis = list(at = seq(0,16,2)), 
   lwd = 1.5)
       
## Another style of the plot:
## 'cex.axis' in 'par' controls for the size of axis tick labels, 
## while 'cex.lab' controls for the size of axis and legend labels. 
## Note that scaling by 'cex' would also affect 'sigcex'. 
op <- par(no.readonly = TRUE)
par(cex.lab = 1.3, cex.axis = 1.1)
wc.avg(my.wc, siglvl = 0.05, sigcol = "red", sigcex = 1.5, 
   minimum.level = 0, maximum.level = 17, 
   periodlab = "period (days)",
   spec.period.axis = list(at = c(1,2,4,8,16,32,64,128)), 
   spec.avg.axis = list(at = seq(0,16,2)),
   lwd = 1.5)
par(op)      

## Plot of wavelet coherence 
## (with color breakpoints according to quantiles):
wc.image(my.wc, which.image = "wc", main = "wavelet coherence, x over y",
   legend.params = list(label.digits = 3),
   periodlab = "period (days)")
         
## Plot of average wavelet coherence:
wc.avg(my.wc, which.avg = "wc", 
   siglvl = 0.05, sigcol = "red", legend.coords = "topleft", 
   periodlab = "period (days)", 
   lwd = 1.5)

## The same plot, setting the minimum plot level of averages to 0: 
wc.avg(my.wc, which.avg = "wc", 
   siglvl = 0.05, sigcol = "red", legend.coords = "topleft", 
   minimum.level = 0, 
   periodlab = "period (days)", 
   lwd = 1.5)
 
## Please see also the examples in our guide booklet,
## URL http://www.hs-stat.com/projects/WaveletComp/WaveletComp_guided_tour.pdf.


## End(Not run)

Image plot of the cross-wavelet power spectrum and wavelet coherence spectrum of two time series

Description

This function plots the cross-wavelet power image, or alternatively the wavelet coherence image, of two time series, which are provided by an object of class "analyze.coherency". The vertical axis shows the Fourier periods. The horizontal axis shows time step counts, but can be easily transformed into a calendar axis if dates are provided in either row names or a variable named "date" in the data frame at hand. Both axes can be relabeled. In particular, an option is given to individualize the period and/or time axis by specifying tick marks and labels.

An option is given to raise cross-wavelet power (or wavelet coherence) values to any (positive) exponent before plotting in order to accentuate the contrast of the image.

The color levels can be defined according to quantiles of values or according to equidistant breakpoints (covering the interval from 0 to maximum level), with the number of levels as a further parameter. A user-defined maximum level can be applied to cross-wavelet power images. In addition, there is an option to adopt an individual color palette.

Further plot design options concern: plot of the cone of influence, plot of contour lines to border areas of significance, plot of the ridge, and plot of arrows (optional: "smoothed" arrows computed from smoothing filters as defined in analyze.coherency) to reflect phase differences.

For that matter, the significance level of contour lines can be defined separately. The plot of the ridge can be restricted to a high-level region ("high" according to a given level of plotted values). In particular, the area to be filled with arrows can be determined in several ways: to reflect significance (at a given level) with respect to cross-wavelet power, wavelet coherence, or individual wavelet power, and/or to flag a high-value region. Furthermore, there is an option to clear out the area where the p-values of cross-wavelet power (coherence, respectively) exceed a given level.

Finally, there is an option to format and insert a color legend (a right-hand vertical color bar) and to set the plot title. For further processing of the plot, graphical parameters of plot regions are provided as output.

The name and parts of the layout were inspired by a similar function developed by Huidong Tian and Bernard Cazelles (archived R package WaveletCo). The code for the arrow design to reflect phase differences has been adopted from Huidong Tian.

Usage

wc.image(WC, 
     which.image = "wp", exponent = 1, 
     plot.coi = TRUE, 
     plot.contour = TRUE, siglvl.contour = 0.1, col.contour = "white",
     plot.ridge = FALSE, lvl = 0, col.ridge = "black",
     plot.arrow = TRUE, use.sAngle = FALSE, 
     p = 1, 
     which.arrow.sig = which.image, 
     siglvl.arrow = 0.05, col.arrow = "black",
     clear.area = FALSE, 
     which.area.sig = which.image, siglvl.area = 0.2, 
     color.key = "quantile", 
     n.levels = 100, 
     color.palette = "rainbow(n.levels, start = 0, end = .7)", 
     maximum.level = NULL,
     useRaster = TRUE, max.contour.segments = 250000,
     plot.legend = TRUE,
     legend.params = list(width=1.2, shrink = 0.9, mar = 5.1, 
                          n.ticks = 6, 
                          label.digits = 1, label.format = "f", 
                          lab = NULL, lab.line = 2.5),
     label.time.axis = TRUE, 
     show.date = FALSE, date.format = NULL, date.tz = NULL,
     timelab = NULL, timetck = 0.02, timetcl = 0.5,
     spec.time.axis = list(at = NULL, labels = TRUE, 
                           las = 1, hadj = NA, padj = NA),
     label.period.axis = TRUE, 
     periodlab = NULL, periodtck = 0.02, periodtcl = 0.5,
     spec.period.axis = list(at = NULL, labels = TRUE, 
                             las = 1, hadj = NA, padj = NA),
     main = NULL,
     lwd = 2, lwd.axis = 1, 
     graphics.reset = TRUE,
     verbose = FALSE)

Arguments

Value

A list of class graphical parameters with the following elements:

Author(s)

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

References

Aguiar-Conraria L., and Soares M.J., 2011. Business cycle synchronization and the Euro: A wavelet analysis. Journal of Macroeconomics 33 (3), 477–489.

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

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 P.C., 1994. Wavelet spectrum analysis and ocean wind waves. In: Foufoula-Georgiou E., and Kumar P., (eds.), Wavelets in Geophysics, Academic Press, San Diego, 151–166.

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.

Veleda D., Montagne R., and Araujo M., 2012. Cross-Wavelet Bias Corrected by Normalizing Scales. Journal of Atmospheric and Oceanic Technology 29, 1401–1408.

See Also

analyze.coherency, wc.avg, wc.sel.phases, wc.phasediff.image, wt.image, wt.avg,
wt.sel.phases, wt.phase.image, reconstruct

Examples

## Not run: 
## The following example is modified from Veleda et al., 2012:

series.length <- 3*128*24
x1 <- periodic.series(start.period = 1*24, length = series.length)
x2 <- periodic.series(start.period = 2*24, length = series.length)
x3 <- periodic.series(start.period = 4*24, length = series.length)
x4 <- periodic.series(start.period = 8*24, length = series.length)
x5 <- periodic.series(start.period = 16*24, length = series.length)
x6 <- periodic.series(start.period = 32*24, length = series.length)
x7 <- periodic.series(start.period = 64*24, length = series.length)
x8 <- periodic.series(start.period = 128*24, length = series.length)

x <- x1 + x2 + x3 + x4 + 3*x5 + x6 + x7 + x8 + rnorm(series.length)
y <- x1 + x2 + x3 + x4 - 3*x5 + x6 + 3*x7 + x8 + rnorm(series.length)

matplot(data.frame(x, y), type = "l", lty = 1, xaxs = "i", col = 1:2, 
 xlab = "index", ylab = "",
 main = "hourly series with periods of 1, 2, 4, 8, 16, 32, 64, 128 days", 
 sub = "(out of phase at period 16, different amplitudes at period 64)")
legend("topright", legend = c("x","y"), col = 1:2, lty = 1)

## 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, y = y)

## Computation of cross-wavelet power and wavelet coherence, x over y:
## 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.
my.wc <- analyze.coherency(my.data, c("x","y"), 
                           loess.span = 0, 
                           dt = 1/24, dj = 1/20, 
                           window.size.t = 1, window.size.s = 1/2, 
                           lowerPeriod = 1/4,
                           make.pval = TRUE, n.sim = 10)

## Plot of cross-wavelet power spectrum, 
## with color breakpoints according to quantiles:
wc.image(my.wc, 
   main = "cross-wavelet power spectrum, x over y",
   legend.params = list(lab = "cross-wavelet power levels (quantiles)"),
   periodlab = "period (days)")
## Note:
## The default time axis shows an index of given points in time, 
## which is the count of hours in our example.      
## By default, arrows are plotted which show the phase differences 
## of x over y at respective significant periods. 
## (Please see our guide booklet for further explanation.)

## The same plot, but with equidistant color breakpoints: 
wc.image(my.wc, color.key = "i", 
   main = "cross-wavelet power spectrum, x over y",
   legend.params = list(lab = "cross-wavelet power levels (equidistant)"),
   periodlab = "period (days)")

## The same plot, but adopting a palette of gray colors, 
## omitting the arrows:
wc.image(my.wc, color.key = "i", 
   main = "cross-wavelet power spectrum, x over y", 
   legend.params = list(lab = "cross-wavelet power levels (equidistant)"),
   color.palette = "gray( (1:n.levels)/n.levels )", 
   plot.arrow = FALSE,
   periodlab = "period (days)")
         
## The same plot, now with ridge of power:
wc.image(my.wc, color.key = "i", 
   main = "cross-wavelet power spectrum, x over y", 
   legend.params = list(lab = "cross-wavelet power levels (equidistant)"),
   color.palette = "gray( (1:n.levels)/n.levels )", 
   plot.arrow = FALSE,
   plot.ridge = TRUE, col.ridge = "red",
   periodlab = "period (days)")             
         
## The plot, turning back to arrows, now in yellow color:
wc.image(my.wc, color.key = "i", 
   main = "cross-wavelet power spectrum, x over y", 
   legend.params = list(lab = "cross-wavelet power levels (equidistant)"),
   color.palette = "gray( (1:n.levels)/n.levels )", 
   col.arrow = "yellow",
   periodlab = "period (days)")   
         
## Alternate styles of the time axis:          
    
## The plot with time elapsed in days, starting from 0 and proceeding 
## in steps of 50 days (50*24 hours), instead of the (default) time index:
index.ticks  <- seq(1, series.length, by = 50*24)
index.labels <- (index.ticks-1)/24    

## Insert your specification of time axis: 
wc.image(my.wc, color.key = "i", 
   main = "cross-wavelet power spectrum, x over y", 
   legend.params = list(lab = "cross-wavelet power levels (equidistant)"),
   color.palette = "gray( (1:n.levels)/n.levels )", 
   col.arrow = "yellow",
   periodlab = "period (days)", timelab = "time elapsed (days)",
   spec.time.axis = list(at = index.ticks, labels = index.labels))

## The plot with (automatically produced) calendar axis:
wc.image(my.wc, color.key = "i", 
   main = "cross-wavelet power spectrum, x over y", 
   legend.params = list(lab = "cross-wavelet power levels (equidistant)"),
   color.palette = "gray( (1:n.levels)/n.levels )", 
   col.arrow = "yellow",
   periodlab = "period (days)",
   show.date = TRUE, date.format = "%F %T")

## Individualizing your calendar axis (works with show.date = TRUE)...
## How to obtain, for example, monthly date ticks and labels:

## The sequence of tick positions:
monthly.ticks <- seq(as.POSIXct("2014-11-01 00:00:00", format = "%F %T"), 
                     as.POSIXct("2015-11-01 00:00:00", format = "%F %T"), 
                     by = "month")
## Observe that the following specification may produce an error:
## 'seq(as.Date("2014-11-01"), as.Date("2015-11-01"), by = "month")'
## Time of the day is missing here!

## The sequence of labels (e.g. information on month and year only):
monthly.labels <- strftime(monthly.ticks, format = "%b %Y")

## Insert your specification of time axis as parameter to wc.image: 
wc.image(my.wc, color.key = "i", 
   main = "cross-wavelet power spectrum, x over y", 
   legend.params = list(lab = "cross-wavelet power levels (equidistant)"),
   color.palette = "gray( (1:n.levels)/n.levels )", 
   col.arrow = "yellow", 
   periodlab = "period (days)",
   show.date = TRUE, date.format = "%F %T", 
   spec.time.axis = list(at = monthly.ticks, labels = monthly.labels, 
                         las = 2))
## Note: 
## The monthly ticks specify the midpoints of the colored cells and 
## match the location of corresponding (default) time index ticks.

## A cross-wavelet power plot with individualized period axis and exponent 
## to accentuate contrast in the image:
wc.image(my.wc, exponent = 0.5, color.key = "i", 
   main = "cross-wavelet power spectrum, x over y",
   legend.params = list(lab = "cross-wavelet power levels 
   (raised by exponent 0.5, equidistant levels)"),
   color.palette = "gray( (1:n.levels)/n.levels )", 
   col.arrow = "yellow", 
   periodlab = "period (days)",
   spec.period.axis = list(at = c(1,2,4,8,16,32,64,128))) 
         
## An option to switch to the corresponding frequency axis:
my.periods <- c(1,2,4,8,16,32,64,128)
my.frequencies <- paste("1/",my.periods, sep = "")
wc.image(my.wc, exponent = 0.5, color.key = "i", 
   main = "cross-wavelet power spectrum, x over y",
   legend.params = list(lab = "cross-wavelet power levels
   (raised by exponent 0.5, equidistant levels)"),
   color.palette = "gray( (1:n.levels)/n.levels )", 
   col.arrow = "yellow", 
   periodlab = "frequency (per day)", 
   spec.period.axis = list(at = my.periods, labels = my.frequencies))            

## Adding, for example, horizontal lines at period ticks...

## There is an option to add further objects to the image plot region, 
## by setting 'graphics.reset = FALSE' 
## (but recall previous par settings after plotting):

op <- par(no.readonly = TRUE)
wc.image(my.wc, exponent = 0.5, color.key = "i", 
   main = "cross-wavelet power spectrum, x over y",
   legend.params = list(lab = "cross-wavelet power levels
   (raised by exponent 0.5, equidistant levels)"),
   color.palette="gray( (1:n.levels)/n.levels )", 
   col.arrow = "yellow", 
   periodlab = "frequency (per day)", 
   spec.period.axis = list(at = my.periods, labels = my.frequencies),
   timelab = "",
   show.date = TRUE, date.format = "%F %T",
   graphics.reset = FALSE)            
abline(h = log2(my.periods))
year2015 <- as.POSIXct("2015-01-01 00:00:00", format = "%F %T")
abline(v = year2015)
axis(1, at = year2015, labels = 2015, padj = 1)
par(op)

## For further axis plotting options: 
## Please see the examples in our guide booklet,
## URL http://www.hs-stat.com/projects/WaveletComp/WaveletComp_guided_tour.pdf.    
            
## Plot of wavelet coherence of x over y, 
## with color breakpoints according to quantiles:
wc.image(my.wc, which.image = "wc", 
   main = "wavelet coherence, x over y",
   legend.params = list(lab = "wavelet coherence levels (quantiles)", 
                        lab.line = 3.5, label.digits = 3),
   periodlab = "period (days)")

## Plot of wavelet coherence, but with equidistant color breakpoints:
wc.image(my.wc, which.image = "wc", color.key = "i",
   main = "wavelet coherence, x over y",
   legend.params = list(lab = "wavelet coherence levels (equidistant)"),
   periodlab = "period (days)")  
   

## End(Not run)

Image plot of phase differences of periodic components for two time series

Description

This function plots the phase difference image of two time series, which is provided by an object of class "analyze.coherency". The vertical axis shows the Fourier periods. The horizontal axis shows time step counts, but can be easily transformed into a calendar axis if dates are provided in either row names or a variable named "date" in the data frame at hand. Both axes can be relabeled. In particular, an option is given to individualize the period and/or time axis by specifying tick marks and labels.

The color levels are defined according to equidistant breakpoints (covering the interval from -pi to +pi), with the number of levels as a further parameter. In addition, there is an option to adopt an individual color palette.

If the default palette is retained, colors indicate the following. Green: phase differences close to zero, which means that the two time series are in phase at the respective period. Yellowgreen: in phase, series 1 leading. Turquoise: in phase, series 2 leading. Red: phase differences are close to +pi, out of phase, series 2 leading. Blue: phase differences are close to -pi, out of phase, series 1 leading.

Further plot design options concern: plot of the cone of influence, plot of contour lines to border areas of significance with respect to cross-wavelet power or wavelet coherence at a given significance level.

Finally, there is an option to insert and format a color legend (a right-hand vertical color bar) and to set the plot title. For further processing of the plot, graphical parameters of plot regions are provided as output.

Usage

wc.phasediff.image(WC, use.sAngle = FALSE, 
   plot.coi = TRUE, 
   plot.contour = TRUE, which.contour = "wp", 
   siglvl = 0.1, col.contour = "white",
   n.levels = 100, 
   color.palette = "rainbow(n.levels, start = 0, end = .7)", 
   useRaster = TRUE, max.contour.segments = 250000,
   plot.legend = TRUE,
   legend.params = list(width = 1.2, shrink = 0.9, mar = 5.1, 
                        n.ticks = 6, 
                        pi.style = TRUE, 
                        label.digits = 1, label.format = "f", 
                        lab = NULL, lab.line = 3),
   label.time.axis = TRUE, 
   show.date = FALSE, date.format = NULL, date.tz = NULL,
   timelab = NULL, timetck = 0.02, timetcl = 0.5,
   spec.time.axis = list(at = NULL, labels = TRUE, 
                         las = 1, hadj = NA, padj = NA),
   label.period.axis = TRUE, 
   periodlab = NULL, periodtck = 0.02, periodtcl = 0.5,
   spec.period.axis = list(at = NULL, labels = TRUE, 
                           las = 1, hadj = NA, padj = NA),
   main = NULL,
   lwd = 2, lwd.axis = 1,
   graphics.reset = TRUE,
   verbose = FALSE)

Arguments

Value

A list of class graphical parameters with the following elements:

Author(s)

Angi Roesch and Harald Schmidbauer

References

Aguiar-Conraria L., and Soares M.J., 2011. Business cycle synchronization and the Euro: A wavelet analysis. Journal of Macroeconomics 33 (3), 477–489.

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

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 P.C., 1994. Wavelet spectrum analysis and ocean wind waves. In: Foufoula-Georgiou E., and Kumar P., (eds.), Wavelets in Geophysics, Academic Press, San Diego, 151–166.

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

Veleda D., Montagne R., and Araujo M., 2012. Cross-Wavelet Bias Corrected by Normalizing Scales. Journal of Atmospheric and Oceanic Technology 29, 1401–1408.

See Also

analyze.coherency, wc.image, wc.avg, wc.sel.phases, wt.image, wt.avg, wt.sel.phases, wt.phase.image, reconstruct

Examples

## Not run: 
## The following example is modified from Veleda et al., 2012:

series.length <- 3*128*24
x1 <- periodic.series(start.period = 1*24, length = series.length)
x2 <- periodic.series(start.period = 2*24, length = series.length)
x3a <- periodic.series(start.period = 4*24, length = series.length)
x3b <- periodic.series(start.period = 4*24, length = series.length,
                       phase = 24)
x4 <- periodic.series(start.period = 8*24, length = series.length)
x5 <- periodic.series(start.period = 16*24, length = series.length)
x6 <- periodic.series(start.period = 32*24, length = series.length)
x7 <- periodic.series(start.period = 64*24, length = series.length)
x8 <- periodic.series(start.period = 128*24, length = series.length)

x <- x1 + x2 + x3a + x4 + 3*x5 + x6 + x7 + x8 + rnorm(series.length)
y <- x1 + x2 + x3b + x4 - 3*x5 + x6 + 3*x7 + x8 + rnorm(series.length)

matplot(data.frame(x, y), type = "l", lty = 1, xaxs = "i", col = 1:2, 
 xlab = "index", ylab = "",
 main = "hourly series with periods of 1, 2, 4, 8, 16, 32, 64, 128 days", 
 sub = "(different phases at periods 4 and 16)")
legend("topright", legend = c("x","y"), col = 1:2, lty = 1)

## 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, y = y)

## Computation of cross-wavelet power and wavelet coherence, x over y:
## 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.
my.wc <- analyze.coherency(my.data, c("x","y"), 
                           loess.span = 0, 
                           dt = 1/24, dj = 1/20, 
                           window.size.t = 1, window.size.s = 1/2, 
                           lowerPeriod = 1/4,
                           make.pval = TRUE, n.sim = 10)

## Plot of cross-wavelet power spectrum, 
## with color breakpoints according to quantiles:
wc.image(my.wc, 
   main = "cross-wavelet power spectrum, x over y",
   legend.params = list(lab = "cross-wavelet power levels (quantiles)"),
   periodlab = "period (days)")

## Default plot of phase differences 
## (with contour lines referring to cross-wavelet power)
wc.phasediff.image(my.wc, which.contour = "wp", 
   main = "image of phase differences, x over y",
   periodlab = "period (days)")
      
## With time elapsed in days 
## (starting from 0 and proceeding in steps of 50 days) 
## instead of the (default) time index:
index.ticks  <- seq(1, series.length, by = 50*24)
index.labels <- (index.ticks-1)/24      
wc.phasediff.image(my.wc, which.contour = "wp", 
   main = "image of phase differences, x over y",
   periodlab = "period (days)",
   timelab = "time elapsed (days)",
   spec.time.axis = list(at = index.ticks, labels = index.labels))

## The same plot, but with (automatically produced) calendar axis:  
wc.phasediff.image(my.wc, which.contour = "wp", 
   main = "image of phase differences, x over y",
   periodlab = "period (days)",
   show.date = TRUE, date.format = "%F %T")
     
## For further axis plotting options:
## Please see the examples in our guide booklet,
## URL http://www.hs-stat.com/projects/WaveletComp/WaveletComp_guided_tour.pdf. 

## Plot of phase difference with numerals as labels of the color legend bar: 
wc.phasediff.image(my.wc, 
   legend.params=list(pi.style = FALSE, label.digits = 2))
   

## End(Not run)

Comparison plot of phases for selected periodic components of two time series

Description

This function plots the phases for periodic components of two time series, which are provided by an object of class "analyze.coherency".

Periodic components can be selected by specification of a single Fourier period or of a period band. In the latter case, and in the default case (no specification at all), phases are averaged across periods for each time series. Other options: restriction to the cone of influence, restriction to an area of significance (with respect to cross-wavelet power, wavelet coherence or individual wavelet power). Phase differences (i.e. angles, smoothed or not smoothed) can be added to the plot.

(The time axis can be altered to display dates, see e.g. wt.image. In particular, an option is given to individualize the phase and/or time axis by specifying tick marks and labels.)

Usage

wc.sel.phases(WC, sel.period = NULL, sel.lower = NULL, sel.upper = NULL, 
      only.coi = FALSE, 
      only.sig = TRUE, which.sig = "wp", siglvl = 0.05, 
      phase.cols = c("red", "blue"), 
      show.Angle = TRUE, use.sAngle = FALSE, Angle.col = "black", 
      show.legend = TRUE, legend.coords = "topleft", legend.horiz = TRUE,
      label.time.axis = TRUE, 
      show.date = FALSE, date.format = NULL, date.tz = NULL, 
      timelab = NULL, timetck = 0.02, timetcl = 0.5,
      spec.time.axis = list(at = NULL, labels = TRUE, 
                            las = 1, hadj = NA, padj = NA),
      label.phase.axis = TRUE, 
      phaselab = NULL, phasetck = 0.02, phasetcl = 0.5,
      spec.phase.axis = list(at = NULL, labels = TRUE, 
                             las = 1, hadj = NA, padj = NA),
      phaselim = c(-pi,pi+show.legend*ifelse(legend.horiz,0.8,2)),      
      main = NULL, sub = NULL,
      lwd = 1, lwd.Angle = 2, lwd.axis = 1, 
      verbose = FALSE)

Arguments

Value

A list of class "sel.phases" with the following elements:

Author(s)

Angi Roesch and Harald Schmidbauer

References

Aguiar-Conraria L., and Soares M.J., 2011. Business cycle synchronization and the Euro: A wavelet analysis. Journal of Macroeconomics 33 (3), 477–489.

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

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 P.C., 1994. Wavelet spectrum analysis and ocean wind waves. In: Foufoula-Georgiou E., and Kumar P., (eds.), Wavelets in Geophysics, Academic Press, San Diego, 151–166.

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

Veleda D., Montagne R., and Araujo M., 2012. Cross-Wavelet Bias Corrected by Normalizing Scales. Journal of Atmospheric and Oceanic Technology 29, 1401–1408.

See Also

analyze.coherency, wc.image, wc.avg, wc.phasediff.image, wt.image, wt.avg,
wt.sel.phases, wt.phase.image, reconstruct

Examples

## Not run: 
## The following example is modified from Veleda et al., 2012:

series.length <- 3*128*24
x1 <- periodic.series(start.period = 1*24, length = series.length)
x2 <- periodic.series(start.period = 2*24, length = series.length)
x3a <- periodic.series(start.period = 4*24, length = series.length)
x3b <- periodic.series(start.period = 4*24, length = series.length,
                       phase = 24)
x4 <- periodic.series(start.period = 8*24, length = series.length)
x5 <- periodic.series(start.period = 16*24, length = series.length)
x6 <- periodic.series(start.period = 32*24, length = series.length)
x7 <- periodic.series(start.period = 64*24, length = series.length)
x8 <- periodic.series(start.period = 128*24, length = series.length)

x <- x1 + x2 + x3a + x4 + 3*x5 + x6 + x7 + x8 + rnorm(series.length)
y <- x1 + x2 + x3b + x4 - 3*x5 + x6 + 3*x7 + x8 + rnorm(series.length)

matplot(data.frame(x, y), type = "l", lty = 1, xaxs = "i", col = 1:2, 
 xlab = "index", ylab = "",
 main = "hourly series with periods of 1, 2, 4, 8, 16, 32, 64, 128 days", 
 sub = "(different phases at periods 4 and 16)")
legend("topright", legend = c("x","y"), col = 1:2, lty = 1)

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, y = y)

## Computation of cross-wavelet power and wavelet coherency of x over y:
## 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.
my.wc <- analyze.coherency(my.data, c("x","y"), loess.span = 0, 
                          dt = 1/24, dj = 1/20,
                          window.size.t = 1, window.size.s = 1/2, 
                          lowerPeriod = 1/4, 
                          make.pval = TRUE, n.sim = 10)

## Plot of cross-wavelet power spectrum, 
## with color breakpoints according to quantiles:
wc.image(my.wc, main = "cross-wavelet power spectrum, x over y",
   legend.params = list(lab = "cross-wavelet power levels (quantiles)"),
   periodlab = "period (days)")

## Select period 64 and compare plots of corresponding phases, including 
## the phase differences (angles) in their non-smoothed (default) version:
wc.sel.phases(my.wc, sel.period = 64, show.Angle = TRUE)

## With time elapsed in days
## (starting from 0 and proceeding in steps of 50 days) 
## instead of the (default) time index:
index.ticks  <- seq(1, series.length, by = 50*24)
index.labels <- (index.ticks-1)/24
wc.sel.phases(my.wc, sel.period = 64, show.Angle = TRUE,  
   timelab = "time elapsed (days)",
   spec.time.axis = list(at = index.ticks, labels = index.labels))

## The same plot, but with (automatically produced) calendar axis:
wc.sel.phases(my.wc, sel.period = 64, show.Angle = TRUE,   
   show.date = TRUE, date.format = "%F %T")

## For further axis plotting options:
## Please see the examples in our guide booklet,
## URL http://www.hs-stat.com/projects/WaveletComp/WaveletComp_guided_tour.pdf.    

## Now, select period 16...
## and observe that corresponding components are out of phase:
wc.sel.phases(my.wc, sel.period = 16, show.Angle = TRUE,
   show.date = TRUE, date.format = "%F %T")
## ... compare to period 4...
wc.sel.phases(my.wc, sel.period = 4, show.Angle = TRUE,
   show.date = TRUE, date.format = "%F %T")   

## In the following, no periods are selected. 
## In this case, instead of individual phases, the plot shows 
## average phases for each series: 
wc.sel.phases(my.wc)


## End(Not run)

Series of weather data and ambient gamma dose rate readings

Description

Ten years (from January 2005 through December 2014; 3652 values) of daily mean temperature, relative humidity, and ambient gamma dose rate (German abbreviation: ODL) readings from Mannheim-Rheinau (Germany).

ODL data retrieved in 2015 from ODL-INFO of the German Federal Office for Radiation Protection (in German: Bundesamt fuer Strahlenschutz, BfS), weather data from the German Meteorological Office (in German: Deutscher Wetterdienst, DWD).
— We owe this example to our former student Nadiya Appelhans.

Usage

data("weather.radiation.Mannheim")

Format

A data frame of four columns:

Details

The ambient gamma dose rate is an equivalent dose representing the stochastic health effects of low levels of ionizing radiation on the human body. According to BfS, radioactivity is to be found everywhere in the environment. It may be of natural or artificial origin. — We drew Mannheim randomly from a set of places for which weather and radiation data were available; there is no conspicuous radiation in Mannheim.

Source

The Federal Office for Radiation Protection (BfS) Germany, https://www.bfs.de/EN/home/home_node.html, and in particular https://www.bfs.de/EN/topics/ion/environment/environment_node.html

ODL-INFO, https://odlinfo.bfs.de/ODL/DE/home/home_node.html

DWD (Deutscher Wetterdienst), https://www.dwd.de/EN/Home/home_node.html

Examples

data(weather.radiation.Mannheim)

plot(as.Date(weather.radiation.Mannheim$date, tz = "UTC"),
     weather.radiation.Mannheim$radiation, type = "l",
     xlab = "day", 
     ylab = "mean daily ambient gamma dose rate (ODL) in microsieverts per hour")

Wavelet transformation, and a simulation algorithm

Description

This function provides Morlet wavelet transformation results of the given time series, performed within the lower-order function WaveletTransform, subject to criteria concerning the time and frequency resolution, and an (optional) lower and/or upper Fourier period. It performs a simulation algorithm to test against a specified alternative hypothesis, which can be chosen from a variety of options, and provides p-values. The selected model will be fitted to the data and simulated according to estimated parameters in order to provide surrogate time series. This function is called by function analyze.wavelet.

The name and parts of the layout were inspired by a similar function developed by Huidong Tian and Bernard Cazelles (archived R package WaveletCo). 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.

Usage

wt(x, start = 1, dt = 1, dj = 1/20, 
   lowerPeriod = 2*dt, upperPeriod = floor(length(x)*dt/3), 
   make.pval = TRUE, method = "white.noise", params = NULL, 
   n.sim = 100, save.sim = FALSE)

Arguments

Value

A list with the following elements:

Author(s)

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

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 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.

See Also

WaveletTransform, analyze.wavelet

Plot of wavelet power averages across time of a single time series

Description

This function plots wavelet power averages across time of a single time series, which are provided by an object of class "analyze.wavelet", or alternatively of class "analyze.coherency". (In the latter case, the series number or name must be specified.) The vertical axis shows the Fourier periods. The horizontal axis shows the averages. User-defined minimum and maximum average levels can be applied. Also, an option is given to individualize the period axis and/or axis of averages by specifying tick marks and labels.

There is an option to label periods according to significance of averages (if p-values are provided) at given levels of significance. Labels are point symbols along the line of averages which can be assigned individually.

The idea to show significance levels by colors of plotting characters and its implementation has been adopted from Huidong Tian and Bernard Cazelles (archived R package WaveletCo).

Usage

wt.avg(WT, my.series = 1, exponent = 1,
  show.siglvl = TRUE, siglvl = c(0.05, 0.1), 
  sigcol = c("red", "blue"), sigpch = 20, sigcex = 1,
  minimum.level = NULL, maximum.level = NULL,
  label.avg.axis = TRUE, 
  averagelab = NULL, averagetck = 0.02, averagetcl = 0.5,
  spec.avg.axis = list(at = NULL, labels = TRUE, 
                       las = 1, hadj = NA, padj = NA),
  label.period.axis = TRUE, 
  periodlab = NULL, periodtck = 0.02, periodtcl = 0.5,
  spec.period.axis = list(at = NULL, labels = TRUE, 
                          las = 1, hadj = NA, padj = NA),
  show.legend = TRUE, legend.coords = "topright", 
  main = NULL, 
  lwd = 1, col = 1, 
  lwd.axis = 1,
  verbose = FALSE)

Arguments

Author(s)

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

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 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.

See Also

analyze.wavelet, wt.image, wt.sel.phases, wt.phase.image, reconstruct

Examples

## 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")
        
my.data <- data.frame(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. 
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)

## Plot of wavelet power spectrum (with equidistant color breakpoints): 
wt.image(my.wt, color.key = "i", main = "wavelet power spectrum",
   legend.params = list(lab = "wavelet power levels (equidistant levels)"),
   periodlab = "period (days)")
## Note:
## The default time axis shows an index of given points in time, 
## which is the count of hours in our example.
         
## With time elapsed in days 
## (starting from 0 and proceeding in steps of 50 days) 
## instead of the (default) time index:
index.ticks  <- seq(1, series.length, by = 50*24)
index.labels <- (index.ticks-1)/24

## Insert your specification of time axis: 
wt.image(my.wt, color.key = "i", main = "wavelet power spectrum",
   legend.params = list(lab = "wavelet power levels (equidistant levels)"),
   periodlab = "period (days)", timelab = "time elapsed (days)",
   spec.time.axis = list(at = index.ticks, labels = index.labels))
         
## Plot of average wavelet power:
wt.avg(my.wt, siglvl = 0.05, sigcol = "red", periodlab = "period (days)")

## The same plot, but with enhanced symbol size, user-defined period axis, 
## and horizontal grid:
wt.avg(my.wt, siglvl = 0.05, sigcol = "red", sigcex = 1.3, 
   periodlab = "period (days)",
   spec.period.axis = list(at = c(1,8,32,128)), 
   periodtck = 1, periodtcl = NULL, 
   lwd = 1.5, lwd.axis = 0.25)
       
## Another style of the plot: 
## With user-defined period axis and axis of averages, 
## minimum and maximum plot levels of averages:

op <- par(no.readonly = TRUE)
par(cex.lab = 1.3, cex.axis = 1.1)
wt.avg(my.wt, siglvl = 0.05, sigcol = "red", sigcex = 1.3, 
   minimum.level = 0, maximum.level = 11,
   periodlab = "period (days)", 
   spec.period.axis = list(at = c(1,8,32,128)), 
   spec.avg.axis = list(at = 0:10),
   lwd = 1.5)
par(op)   
## Note:
## 'cex.axis' in 'par' controls for the size of axis tick labels, 
## while 'cex.lab' controls for the size of axis and legend labels. 
## Scaling by 'cex' would also affect 'sigcex'. 


## Please see also the examples in our guide booklet,
## URL http://www.hs-stat.com/projects/WaveletComp/WaveletComp_guided_tour.pdf.


## End(Not run)

Image plot of the wavelet power spectrum of a single time series

Description

This function plots the wavelet power spectrum of a single time series, which is provided by an object of class "analyze.wavelet", or alternatively of class "analyze.coherency". (In the latter case, the series number or name must be specified.) The vertical axis shows the Fourier periods. The horizontal axis shows time step counts, but can be easily transformed into a calendar axis if dates are provided in either row names or as a variable named "date" in the data frame at hand. Both axes can be relabeled. In particular, an option is given to individualize the period and/or time axis by specifying tick marks and labels.

An option is given to raise wavelet power values to any (positive) exponent before plotting in order to accentuate the contrast of the image.

The color levels can be defined according to quantiles of values or according to equidistant breakpoints (covering the interval from 0 to maximum level), with the number of levels as a further parameter. A user-defined maximum level can be applied. In addition, there is an option to adopt an individual color palette.

Further plot design options concern: plot of the cone of influence, plot of wavelet power contour lines at a specified level of significance, plot of power ridges.

Finally, there is an option to insert and format a color legend (a right-hand vertical color bar) and to set the plot title. For further processing of the plot, graphical parameters of plot regions are provided as output.

The name and parts of the layout were inspired by a similar function developed by Huidong Tian and Bernard Cazelles (archived R package WaveletCo).

Usage

wt.image(WT, my.series = 1, exponent = 1,
     plot.coi = TRUE, 
     plot.contour = TRUE, siglvl = 0.1, col.contour = "white",
     plot.ridge = TRUE, lvl = 0, col.ridge = "black", 
     color.key = "quantile", 
     n.levels = 100, 
     color.palette = "rainbow(n.levels, start = 0, end = .7)",
     maximum.level = NULL,
     useRaster = TRUE, max.contour.segments = 250000,
     plot.legend = TRUE,
     legend.params = list(width = 1.2, shrink = 0.9, mar = 5.1, 
                          n.ticks = 6, 
                          label.digits = 1, label.format = "f", 
                          lab = NULL, lab.line = 2.5),
     label.time.axis = TRUE, 
     show.date = FALSE, date.format = NULL, date.tz = NULL, 
     timelab = NULL, timetck = 0.02, timetcl = 0.5,
     spec.time.axis = list(at = NULL, labels = TRUE, 
                           las = 1, hadj = NA, padj = NA),
     label.period.axis = TRUE, 
     periodlab = NULL, periodtck = 0.02, periodtcl = 0.5,
     spec.period.axis = list(at = NULL, labels = TRUE, 
                             las = 1, hadj = NA, padj = NA),
     main = NULL,
     lwd = 2, lwd.axis = 1, 
     graphics.reset = TRUE,
     verbose = FALSE)

Arguments

Value

A list of class graphical parameters with the following elements:

Author(s)

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

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 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.

See Also

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

Examples

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

## Plot of wavelet power spectrum 
## with color breakpoints referring to quantiles:  
wt.image(my.wt, main = "wavelet power spectrum",
   legend.params = list(lab = "wavelet power levels (quantiles)", 
                        lab.line = 3.5, 
                        label.digits = 2),
   periodlab = "period (days)")
## Note:
## The default time axis shows an index of given points in time, 
## which is the count of hours in our example.

## The same plot, but with equidistant color breakpoints: 
wt.image(my.wt, color.key = "i", main = "wavelet power spectrum",
   legend.params = list(lab = "wavelet power levels (equidistant)"),
   periodlab = "period (days)")
         
## Alternative styles of the time axis:          
         
## The plot with time elapsed in days, starting from 0 and proceeding 
## in steps of 50 days (50*24 hours),
## instead of the (default) time index:
index.ticks  <- seq(1, series.length, by = 50*24)
index.labels <- (index.ticks-1)/24
## Insert your specification of the time axis: 
wt.image(my.wt, color.key = "i", main = "wavelet power spectrum",
   legend.params = list(lab = "wavelet power levels (equidistant)"),
   periodlab = "period (days)", timelab = "time elapsed (days)",
   spec.time.axis = list(at = index.ticks, labels = index.labels))

## The plot with (automatically produced) calendar axis:
wt.image(my.wt, color.key = "i", main = "wavelet power spectrum",
   legend.params = list(lab = "wavelet power levels (equidistant)"), 
   periodlab = "period (days)",
   show.date = TRUE, date.format = "%F %T")

## Individualizing your calendar axis (works with 'show.date = TRUE')...
## How to obtain, for example, monthly date ticks and labels:

## The sequence of tick positions:
monthly.ticks <- seq(as.POSIXct("2014-11-01 00:00:00", format = "%F %T"), 
                     as.POSIXct("2016-11-01 00:00:00", format = "%F %T"), 
                     by = "month")
## Observe that the following specification may produce an error:
## 'seq(as.Date("2014-11-01"), as.Date("2016-11-01"), by = "month")'
## Time of the day is missing here!

## The sequence of labels (e.g. information on month and year only):
monthly.labels <- strftime(monthly.ticks, format = "%b %Y")

## Insert your specification of the time axis: 
wt.image(my.wt, color.key = "i", main = "wavelet power spectrum",
   legend.params = list(lab = "wavelet power levels (equidistant)"), 
   periodlab = "period (days)", 
   show.date = TRUE, date.format = "%F %T", 
   spec.time.axis = list(at = monthly.ticks, labels = monthly.labels, 
                         las = 2))
## Note: 
## The monthly ticks specify the midpoints of the colored cells and match 
## the location of corresponding (default) time index ticks.
         
## Furthermore, the plot with an individualized period axis:
wt.image(my.wt, color.key = "i", main = "wavelet power spectrum",
   legend.params = list(lab = "wavelet power levels (equidistant)"), 
   periodlab = "period (days)", 
   show.date = TRUE, date.format = "%F %T",
   spec.time.axis = list(at = monthly.ticks, labels = monthly.labels, 
                         las = 2),
   spec.period.axis = list(at = c(1,8,32,128))) 
         
## Switching the time axis from index to time elapsed in hours 
## (starting from 0, and proceeding in steps of 500 hours), 
## and the period axis from days to hours:
index.ticks  <- seq(1, series.length, by = 500)
index.labels <- index.ticks - 1
wt.image(my.wt, color.key = "i", main = "wavelet power spectrum",
   legend.params = list(lab = "wavelet power levels (equidistant)"), 
   timelab = "time elapsed (hours)", periodlab = "period (hours)",
   spec.time.axis = list(at = index.ticks, labels = index.labels),
   spec.period.axis = list(at = c(1,8,32,128), labels = c(1,8,32,128)*24))            

## A plot with different colors:
wt.image(my.wt, main = "wavelet power spectrum",
   legend.params = list(lab = "wavelet power levels (quantiles)", 
                        lab.line = 3.5, label.digits = 2), 
   color.palette = "gray((1:n.levels)/n.levels)", col.ridge = "yellow",
   periodlab = "period (days)")

## In the case of monthly (or quarterly) data, the time axis should be 
## labeled at equally spaced time points. An example:

monthyear <- seq(as.Date("2014-01-01"), as.Date("2018-01-01"),
                 by = "month")
monthyear <- strftime(monthyear, format = "%b %Y")

xx <- periodic.series(start.period = 6, length = length(monthyear))
xx <- xx + 0.2*rnorm(length(monthyear))

plot(xx, type = "l", xlab = "index", ylab = "", xaxs = "i",
     main = "monthly series with period of 6 months")

monthly.data <- data.frame(date = monthyear, xx = xx)
 
my.wt <- analyze.wavelet(monthly.data, "xx", loess.span = 0, 
                         dt = 1, dj = 1/250, 
                         make.pval = TRUE, n.sim = 250)
## Note: 
## The natural choice of 'dt' in this example is 'dt = 1', 
## resulting in periods measured in months. 
## (Setting 'dt = 1/12' would result in periods measured in years.)

## The default wavelet power plot then shows the monthly:
wt.image(my.wt, main = "wavelet power spectrum", 
   periodlab = "period (months)")

## The following plot shows the elapsed time, measured in months: 
wt.image(my.wt, main = "wavelet power spectrum", 
   periodlab = "period (months)", timelab = "time elapsed (months)",
   spec.time.axis = list(at = 1:length(monthyear), 
                         labels = (1:length(monthyear))-1))
                                     
## In case you prefer the monthyear labels themselves: 
wt.image(my.wt,  main = "wavelet power spectrum", 
   periodlab = "period (months)", timelab = "month and year",
   spec.time.axis = list(at = 1:length(monthyear), labels = monthyear)) 

## You may sometimes wish to enhance your plot with additional information. 
## There is an option to add further objects to the image plot region, 
## by setting 'graphics.reset = FALSE' 
## (but recall previous par settings after plotting):

op <- par(no.readonly = TRUE)
wt.image(my.wt, main = "wavelet power spectrum", 
   periodlab = "period (months)", 
   spec.period.axis = list(at = c(2,4,6,8,12)), 
   spec.time.axis = list(at = 1:length(monthyear),
                         labels = substr(monthyear,1,3)),
   graphics.reset = FALSE)
abline(h = log2(6), lty = 3)         
abline(v = seq(1, length(monthyear), by = 12), lty = 3)
mtext(2014:2018, side = 1, 
      at = seq(1, length(monthyear), by = 12), line = 2)
par(op)

## For further axis plotting options:
## Please see the examples in our guide booklet,
## URL http://www.hs-stat.com/projects/WaveletComp/WaveletComp_guided_tour.pdf.


## End(Not run)

Image plot of the phases of periodic components for a single time series

Description

This function plots the wavelet phase image for a time series, which is provided by an object of class "analyze.wavelet", or alternatively of class "analyze.coherency". (In the latter case, the series number or name must be specified.) The vertical axis shows the Fourier periods. The horizontal axis shows time step counts, but can be easily transformed into a calendar axis if dates are provided in either row names or as a variable named "date" in the data frame at hand. Both axes can be relabeled. In particular, an option is given to individualize the period and/or time axis by specifying tick marks and labels.

The color levels are defined according to equidistant breakpoints (covering the interval from -pi to +pi), with the number of levels as a further parameter. In addition, there is an option to adopt an individual color palette.

If the default palette is retained, colors indicate the following. Green: Phases close to zero. Red: phases close to +pi. Blue: phases close to -pi.

Further plot design options concern: plot of the cone of influence, plot of contour lines to border areas of significance with respect to cross-wavelet power or wavelet coherency at a given significance level, plot of power ridges.

Finally, there is an option to insert and format a color legend (a right-hand vertical color bar) and to set the plot title. For further processing of the plot, graphical parameters of plot regions are provided as output.

Usage

wt.phase.image(WT, my.series = 1,
       plot.coi = TRUE, 
       plot.contour = TRUE, 
       siglvl = 0.1, col.contour = "white",
       plot.ridge = TRUE, col.ridge = "black", 
       n.levels = 100, 
       color.palette = "rainbow(n.levels, start = 0, end = .7)", 
       useRaster = TRUE, max.contour.segments = 250000,
       plot.legend = TRUE,
       legend.params = list(width = 1.2, shrink = 0.9, mar = 5.1, 
                            n.ticks = 6, 
                            pi.style = TRUE,
                            label.digits = 1, label.format = "f", 
                            lab = NULL, lab.line = 3),
       label.time.axis = TRUE, 
       show.date = FALSE, date.format = NULL, date.tz = NULL, 
       timelab = NULL, timetck = 0.02, timetcl = 0.5,
       spec.time.axis = list(at = NULL, labels = TRUE, 
                             las = 1, hadj = NA, padj = NA),
       label.period.axis = TRUE, 
       periodlab = NULL, periodtck = 0.02, periodtcl = 0.5,
       spec.period.axis = list(at = NULL, labels = TRUE, 
                               las = 1, hadj = NA, padj = NA),
       main = NULL,
       lwd = 2, lwd.axis = 1, 
       graphics.reset = TRUE,
       verbose = FALSE)

Arguments

Value

A list of class graphical parameters with the following elements:

Author(s)

Angi Roesch and Harald Schmidbauer

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.

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

analyze.wavelet, wt.image, wt.avg, wt.sel.phases, reconstruct

Examples

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

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

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

## Default image of phases:
wt.phase.image(my.wt, 
   main = "image of phases", 
   periodlab = "period (days)")

## With time elapsed in days 
## (starting from 0 and proceeding in steps of 50 days) 
## instead of the (default) time index:
index.ticks  <- seq(1, series.length, by = 50*24)
index.labels <- (index.ticks-1)/24
wt.phase.image(my.wt, 
   main = "image of phases", 
   periodlab = "period (days)",
   timelab = "time elapsed (days)",
   spec.time.axis = list(at = index.ticks, labels = index.labels))

## The same plot, but with (automatically produced) calendar axis:
wt.phase.image(my.wt, 
   main = "image of phases", periodlab = "period (days)", 
   show.date = TRUE, date.format = "%F %T")

## For further axis plotting options:
## Please see the examples in our guide booklet,
## URL http://www.hs-stat.com/projects/WaveletComp/WaveletComp_guided_tour.pdf.    

## Image plot of phases with numerals as labels of the color legend bar: 
wt.phase.image(my.wt, 
   legend.params=list(pi.style = FALSE, label.digits = 2))
   

## End(Not run)

Plot phases for selected periodic components of a single time series

Description

This function plots the phases for selected periodic components of a time series, which are provided by an object of class "analyze.wavelet", or alternatively of class "analyze.coherency". (In the latter case, the series number or name must be specified.)

Periodic components can be selected by specification of a single Fourier period or of a period band. In the latter case, and in the default case (no specification at all), there is an option to average the phases across periods. Other options: restriction to the cone of influence, restriction to an area of significance (with respect to wavelet power).

(The time axis can be altered to give dates, see e.g. wt.image. In particular, an option is given to individualize the phase and/or time axis by specifying tick marks and labels.)

Usage

wt.sel.phases(WT, my.series = 1,
      sel.period = NULL, sel.lower = NULL, sel.upper = NULL, 
      only.coi = FALSE,
      only.sig = TRUE, siglvl = 0.05, 
      show.avg.phase = FALSE, phase.avg.col = "black", 
      label.time.axis = TRUE, 
      show.date = FALSE, date.format = NULL, date.tz = NULL, 
      timelab = NULL, timetck = 0.02, timetcl = 0.5,
      spec.time.axis = list(at = NULL, labels = TRUE, 
                            las = 1, hadj = NA, padj = NA),
      label.phase.axis = TRUE, 
      phaselab = NULL, phasetck = 0.02, phasetcl = 0.5,
      spec.phase.axis = list(at = NULL, labels = TRUE, 
                             las = 1, hadj = NA, padj = NA),
      main = NULL, sub = NULL, 
      lwd = 1, lwd.axis = 1,
      verbose = FALSE)

Arguments

Value

A list of class "sel.phases" with the following elements:

Author(s)

Angi Roesch and Harald Schmidbauer

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.

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

analyze.wavelet, wt.image, wt.avg, wt.phase.image, reconstruct

Examples

## 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")
     
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.
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)

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

## Select period 8 and plot corresponding phases across time:
wt.sel.phases(my.wt, sel.period = 8)

## With time elapsed in days 
## (starting from 0 and proceeding in steps of 50 days) 
## instead of the (default) time index:
index.ticks  <- seq(1, series.length, by = 50*24)
index.labels <- (index.ticks-1)/24
wt.sel.phases(my.wt, sel.period = 8,  
   timelab = "time elapsed (days)",
   spec.time.axis = list(at = index.ticks, labels = index.labels))

## The same plot, but with (automatically produced) calendar axis:
wt.sel.phases(my.wt, sel.period = 8,  
   show.date = TRUE, date.format = "%F %T")

## For further axis plotting options:
## Please see the examples in our guide booklet,
## URL http://www.hs-stat.com/projects/WaveletComp/WaveletComp_guided_tour.pdf.        

## In the following, no period is selected.
## By setting 'show.avg.phase = TRUE', the plot shows average phases 
## instead of individual phases: 
wt.sel.phases(my.wt, show.avg.phase = TRUE)


## End(Not run)