Using wavelet-based methods, this function estimates the trend and evolutionary wavelet spectrum (EWS) of a nonstationary time series.
Two methods are implemented (see references), the direct estimator (T.est.type="linear" and
S.do.diff=FALSE), and the difference estimator (T.est.type="nonlinear") and S.do.diff=TRUE)
The defaults give the direct estimator.
All the defaults are set carefully. Key times to change defaults are
if the data contains "cusps", then the difference estimator is preferred.
to assess stability of the estimate to the wavelet, change the wavelet number T.filter.number and
S.filter.number and/or the wavelet type T.family and S.family, see details.
The arguments affecting trend are preceded by T. and those affecting spectral estimation are preceded
by S..
TLSW(
x,
do.trend.est = TRUE,
do.spec.est = TRUE,
T.est.type = "linear",
T.filter.number = 4,
T.family = "DaubExPhase",
T.transform = "nondec",
T.boundary.handle = TRUE,
T.max.scale = floor(log2(length(x)) * 0.7),
T.thresh.type = "hard",
T.thresh.normal = TRUE,
T.CI = FALSE,
T.sig.lvl = 0.05,
T.reps = 200,
T.CI.type = "normal",
T.lacf.max.lag = floor(10 * (log10(length(x)))),
S.filter.number = T.filter.number,
S.family = T.family,
S.smooth = TRUE,
S.smooth.type = "mean",
S.binwidth = floor(6 * sqrt(length(x))),
S.max.scale = floor(log2(length(x)) * 0.7),
S.boundary.handle = TRUE,
S.inv.mat = NULL,
S.do.diff = FALSE,
S.lag = 1,
S.diff.number = 1,
gen.filter.number = S.filter.number,
gen.family = S.family
)An object of class "TLSW", a list that contains the following components:
Input data
Input parameter, logical variable specifying if spectral estimation was performed.
A list object, returned if do.spec.est = TRUE. Contains relevant input parameters
and the following fields related to the spectrum estimate:
S: The evolutionary wavelet spectral (smoothed and corrected) estimate of the input data. This object is of
class wd and so can be plotted and printed in the usual way using wavethresh
functionality.
WavPer: The raw wavelet periodogram of the input
data. The EWS estimate (S, above) is the smoothed corrected version of this raw
wavelet periodogram.
SmoothWavPer: The smoothed, but uncorrected raw
wavelet periodogram of the input data.
Input parameter, logical variable specifying if trend estimation was performed.
A list object, returned if do.trend.est = TRUE. Contains relevant input parameters
and the following fields related to the trend estimate:
T: A vector of length length(x) containing the trend estimate.
lower.CI: Returned if T.CI = TRUE. The lower limit of the pointwise confidence interval.
upper.CI: Returned if T.CI = TRUE. The upper limit of the pointwise confidence interval.
The time series you wish to analyse.
Logical variable, indicating whether trend estimation is to be performed on the time series.
Logical variable, indicating whether spectral estimation is to be performed on the time series.
String indicating type of wavelet thresholding used. Can be "linear" (default), which means
that all non-boundary wavelet coefficients are set to zero, or "nonlinear", where
each wavelet coefficient is thresholded using a time-varying, noise-dependent threshold.
The index number for the wavelet used for trend estimation.
The family of the wavelet used for trend estimation.
String giving the type of wavelet transform used for trend estimation.
Can be "dec", in which case a standard (decimated) wavelet transform is used, or "nondec" (default),
in which case a nondecimated transform is used.
Logical variable, if TRUE, the time series is
boundary corrected when estimating the trend.
Integer variable, selects the number of scales of the wavelet transform to apply thresholding to for trend estimation.
String variable, used only if T.est.type = "nonlinear"; the type of
thresholding function used in the trend estimation. Can be
"soft" or "hard" (default).
Logical variable, used only if T.est.type = "nonlinear";
if TRUE, uses a threshold assuming the data are normally
distributed. If FALSE, uses a larger threshold to reflect non-normality.
Logical variable. If TRUE, a (1-T.sig.lvl) pointwise confidence interval is
computed for the trend estimate. When T.transform = "dec" and T.est.type = "linear", this is
computed using the asymptotic distribution of the trend estimator.
Otherwise, it is computed via bootstrapping.
Used only if T.CI = TRUE; a numeric value
(0 <= T.sig.lvl <= 1) with which a (1-T.sig.lvl) pointwise
confidence interval for the trend estimate is generated.
Used only if T.transform = "nondec" and T.CI = TRUE; the number of bootstrap
replications used to calculate the confidence interval.
Used only if T.transform = "nondec" and T.CI = TRUE; the type of confidence
interval computed. Can be "percentile", in which case empirical percentiles are used, or
"normal" (default), in which case the (symmetric) normal approximation is used.
Used only if T.est.type = "linear" and T.CI = TRUE;
the maximum lag of the autocovariance to compute needed for calculating the asymptotic confidence interval.
The index number for the wavelet used for spectrum estimation.
The family of the wavelet used for spectrum estimation.
A logical variable to indicate whether smoothing is performed on the wavelet periodogram.
String indicating which type of smoothing to use on wavelet periodogram. Can be one of
"mean": (default) running mean smoother.
"median": running median smoother.
"epan": Epanechnikov kernel smoother.
The bin width of the smoother used to smooth the raw wavelet periodogram.
The coarsest wavelet scale used to estimate the spectrum. Should be a positive integer less than \(J\), where \(n=2^J\) is the length of the time series.
Logical variable, if TRUE, the time series is boundary corrected, to get a more accurate spectrum estimate at the boundaries of the times series. If FALSE, no boundary correction is applied. Recommended to use TRUE.
The user can pre-calculate and supply the appropriate correction matrix used to correct the raw wavelet periodogram. If left blank, then the correction matrix is calculated when performing spectral estimation.
Logical variable, indicating if the time series is to be differenced before spectral estimation is performed.
The lag of differencing used, only applicable if S.do.dif = TRUE.
The number of differencing operations performed,
only applicable if S.do.diff = TRUE. A first difference is strongly recommended as default.
The index number for the wavelet that generates the
stochastic component of the time series. For the "DaubExPhase" family, the filter number can be between
1 to 10. For the "DaubLeAsymm" family, the filter number can be between 4 to 10.
Recommended to leave as the default, set to the same as S.filter.number.
The family of the generating wavelet. It is recommended to
use either the Daubechies Extremal Phase family, or the Daubechies Least
Asymmetric family, corresponding to the "DaubExPhase" or the "DaubLeAsymm"
options. Recommended to leave as the default, set to the same as S.family.
The fitted trend LSW process \(X_{t,n} \), \(t = 0, \ldots , n-1\), and \(n = 2^J\) is a doubly-indexed stochastic process with the following representation in the mean square sense: $$X_{t} = T_t + \varepsilon_t = T_t + \sum_{j = 1}^{\infty} \sum_{k \in \mathbb{Z}} w_{j,k;n} \psi_{j,k} (t) \xi_{j,k} ,$$
where \(\{\xi_{j,k} \}\) is a random, uncorrelated, zero-mean orthonormal increment sequence, \(\{w_{j,k;n} \}\) is a set of amplitudes, and \(\{\psi_{j, k} \}_{j,k}\) is a set of discrete non-decimated wavelets. The trend component \(T_t := T(t/n)\) is assumed to be a general smooth (Holder) continuous function. See the referenced papers for full details of the model. The key considerations for users are:
The model assumes smooth trend and spectral components. The larger the T.filter.number
the smoother the assumption on the underlying trend and similarly for S.filter.number and the
spectral estimate.
The choice of wavelet (smoothness assumption) does affect the estimation so one should check the
robustness of their conclusions to the choice of wavelet (T.filter.number and S.filter.number).
This is akin to selecting the kernel in nonparametric modelling.
The underlying methods are designed for signals of length \(n=2^J\) and so modifications are made to signals which are not of this form. A natural approach is to extend the data (at both ends) and the default approach does this by reflection with a trend correction to avoid discontinuities.
McGonigle, E. T., Killick, R., and Nunes, M. (2022a). Trend locally stationary wavelet processes. Journal of Time Series Analysis, 43(6), 895-917.
McGonigle, E. T., Killick, R., and Nunes, M. (2022b). Modelling time-varying first and second-order structure of time series via wavelets and differencing. Electronic Journal of Statistics, 6(2), 4398-4448.
McGonigle, E. T., Killick, R., and Nunes, M. (2025). TrendLSW: Trend and Spectral Estimation of Nonstationary Time Series in R. Journal of Statistical Software, 115(10), 1-30, doi:10.18637/jss.v115.i10.
plot.TLSW, summary.TLSW, print.TLSW, wd, ewspec3
# simulates an example time series and estimates its trend and evolutionary wavelet spectrum
spec <- matrix(0, nrow = 9, ncol = 2^9)
spec[1, ] <- seq(from = 1, to = 9, length = 512)
trend <- sin(pi * (seq(from = 0, to = 4, length = 512)))
set.seed(1)
x <- TLSWsim(trend = trend, spec = spec)
plot.ts(x)
x.TLSW <- TLSW(x)
summary(x.TLSW)
plot(x.TLSW) # by default plots both the trend and spectrum estimates
Run the code above in your browser using DataLab