The details of our testing methodology are set out in the
Cardinali and Nason paper referenced below. Essentially, the testing process works as follows. First, one
has to define a test statistic. Given a time series this has
return a statistic that measures `degree of nonstationarity'.
For example, estimating the EWS, and then computing the
sum of the sample variances of each scale is such as
measure (and known as the $T_{vS}$ statistic).
This statistic is zero for a constant spectrum and positive
for non-constant spectrum (and generally larger for larger variations
of the spectrum).
Once a test statistic T is selected then a parametric Monte Carlo
test can be used. First, T is computed on the series itself.
Then, for statistical assessment of the `significance' of the test
statistic the following procedure is carried out. Assuming, for
a moment that the time series is stationary, we estimate its
evolutionary wavelet spectrum (EWS) and then average this over
time ($\bar{S}_j$). Then we use the function
LSWsim to simulate a time series whose EWS is the
constant, stationary, spectral estimate. Then we compute our
test statistic, $$T_b$$, on this simulated series.
Then we calculate $$T_b$$ for Bsim-1 simulations. The function
then returns BSim numbers. The first is the test statistic
computed on the actual data. The remaining ones are the test
statistic computed on the simulated stationary series.
The idea being that if the time series is really stationary then
the first value will be comparable to the ones obtained by simulation.
If the time series is not stationary then the first test statistic
will be much larger than the ones obtained by simulation (since
the actual data T will have been computed on a time series with
varying spectrum, whereas the simulated ones are all computed on
constant spectra, and their variation is only due to sampling
variation).
The test statistic supplied to this function (as argument
theTS) should take an EWS object as an argument.
For example, the WaveThresh function ewspec
produces a suitable spectral estimate in its $S
argument (both objects are actually examples of a
non-decimated wavelet transform object, class wd).
The function plotBS can be used the present the results
of this function in an interpretable form and calculate the
p-value of the test.