Tests of Stationarity

<<< locally stationary time series page

Why Should You Test?

Okay. So, your time series is loaded into your favourite statistical package. What now? Probably the first thing you need to do is produce a plot of your time series. The plot will give you an idea of the overall levels and variability of the series. The plot will give you an idea of any trends or seasonality in the series. This kind of evaluation is part of an initial data analysis and an excellent description can be found in Chapter 2 of Chatfield, listed in the references. After trend and seasonality are assessed they are often removed and the residuals are then further analyzed for stochastic structure.

Often, the next step commonly advocated is to compute autocorrelations or autocovariance (again, see the brief introduction to stationary series for more details on this). However, these statistics rely on the assumption that the series is stationary, i.e. has statistical properties that do not change with time.

We were a little vague in the introduction about the definition of a stationary series. In fact, there are different, related (and more technical) definitions of stationarity. We will attempt to describe them here in informal terms.

Strict Stationarity

Strict stationarity is the strongest form of stationarity. It means that the joint statistical distribution of any collection of the time series variates never depends on time. So, the mean, variance and any moment of any variate is the same whichever variate you choose. The formal mathematical definition of strictly stationary series can be found on the Wiki page. However, for day to day use strict stationarity is too strict. Hence, the following weaker definition is often used instead.

Stationarity of order 2

For everyday use we often consider time series that have:

Such time series are known as second-order stationary or stationary of order 2.

From now on, whenever we mention stationarity, we mean second-order stationarity.

Note: it is possible to consider a weaker form of stationarity still: a series that is first-order stationary which means that the mean is a constant function of time. Economists are keen on this kind of stationarity, particularly in how to combine time series with time-varying means to obtain one which is first-order stationary (for example). This latter concept is known as cointegration.

The Role of Stationarity

So far we have explained that stationarity (second order or strict) is about imposing constancy of certain time series quantities. Why is this a useful concept? Certainly, much data seem to obey this rule in that future statistical behaviour is identical to past behaviour.

On the other hand, much data is not stationary or at least only approximately stationary. A real problem is that although there are tests for stationarity we submit that they are not used much in practice. Why is this?

Why do analysts persist with stationary models that are not appropriate and potentially risky? We offer four reasons. 1. Fear of diversity. There is a single mathematical model (the Fourier-Cramer model) for stationary time series. For nonstationary series the situation can be complex and the diversity of potential models can be daunting. 2. Education. Many undergraduate time series courses only have time or the ambition to consider stationary models, and 3. Mathematical expediency. Stationary models are mathematically easier to study and develop asymptotic theories for (that is, mathematically we understand how our modelling works for larger and larger samples). 4. Maturity. Stationary theory is mature, widely known and widely applicable.

However, it is the case that many real time series are just not stationary. Series often display trends (invalidating first-order stationarity), or seasonalities, or changes in variance (invalidating second-order stationarity). Hence, statistics and related fields have a second armoury of techniques that can manipulate time series to become stationary (differencing, variable transformations such as taking logs or square roots). After manipulation the series can be treated as stationary and standard methods used.

Tests of Second-order Stationarity

This section assume that the trend and seasonality has been modeled and removed from your time series and you wish to test whether it is second-order stationary. Here, by TEST we mean a statistically rigorous hypothesis test. We will examine two tests here based on methods described in Priestley and Subba Rao (1969) and Nason (2013). The reason for focussing on these tests is that there are freely available implementations in the freeware programming language R. However, it is important to note that there are several other possible tests that could be used and some of them are listed in Nason (2013).

The Priestley-Subba Rao (PSR) Test

Description. The PSR test centres around the time-varying Fourier spectrum ft(w) where t is time and w is frequency. For a stationary time series the time-varying spectrum is, not surprisingly, a constant function of time. The PSR test investigates how "non-constant" ft(w) is a function of time. It does this by looking at the logarithm of an estimator of ft(w), i.e. obtaining

Y(t, w) = log{Ft(w)},

where Ft(w) is an estimator of f. Then, approximately,

E {Y(t, w)} = log{ft(w)}.

and the variance of Y(t, w) is approximately constant. Here the logarithm acts as a variance stabilizer permitting us to focus on changes in the mean structure of Y. These actions permit us to write Y(t, w) as a linear model with constant variance and test the constancy of f using a standard one-way analysis of variance (ANOVA).

Implementation. The PSR test is implemented in the fractal package in R available from the CRAN package repository (at the time of writing the package is currently in the Archive). The function that actually implements the test is called stationarity.

Example. Let's look at how to use the stationarity function to carry out a test of stationarity. First, we will get a time series to use as an example and test for stationarity. We will use the Earthquake/Explosion data set that is described in Shumway and Stoffer (2011). This can be acquired through the astsa package.

First, install the fractal, astsa packages into R in your normal way. Then load the packages and make the earthquake/explosion data available:

library("fractal")
library("astsa")
data("eqexp")

The eqexp object contains 17 columns corresponding to recording of 8 earthquake signals, 8 explosion signals and a final signal from "the Novaya Zemlya event of unknown origin". Each column consists of two signals: the P-wave which occupies the first 1024 entries, and the Q-wave which occupies the second 1024 entries. We will extract signal 14's P-wave (which corresponds to the explosion P-wave depicted in Nason (2013)) and then plot it.

exP <- eqexp[1:1024, 14]
ts.plot(exP)

Explosion P-wave for signal 14

Let us now apply the Priestley-Subba Rao test of stationarity:

stationarity(exP)

and this results in numerical output ending with:

Priestley-Subba Rao stationarity Test for exP
---------------------------------------------
Samples used : 1024
Samples available : 1020
Sampling interval : 1
SDF estimator : Multitaper
Number of (sine) tapers : 5
Centered : TRUE
Recentered : FALSE
Number of blocks : 10
Block size : 102
Number of blocks : 10
p-value for T : 0
p-value for I+R : 0.0003388925
p-value for T+I+R : 0


The key line to examine in this output is the p-value for the T component (which is the test of stationarity p-value for variation over time). In this example the p-value is essentially zero which indicates that there is very strong evidence to reject the null hypothesis of stationarity.

Wavelet Spectrum Test

Description. This test of stationarity looks at a quantity called βj(t) which is closely related to a wavelet-based time-varying spectrum of the time series (it is a linear transform of the evolutionary wavelet spectrum of the locally stationary wavelet processes of Nason, von Sachs and Kroisandt, 2000). Again, we need to see whether the βj(t) function varies over time or is constant. Naturally, though, since we have data we need to examine an estimate of βj(t).

To verify constancy of βj(t) we use the technique due to von Sachs and Neumann (2000) which examines the Haar wavelet coefficients of the estimate of βj(t). A function f(t) is constant if and only if all of its Haar wavelet coefficients are zero. So, von Sachs and Neumann (2000) perform a multiple hypothesis test on all the Haar wavelet coefficients (of a potentially time-varying quantity) and is they are all deemed to be zero then the function is deemed to be constant. Von Sachs and Neumann (2000) develop powerful theory that establishes the asymptotic normality of the Haar wavelet coefficients under mild assumptions and for some time series with heavy tails also. Von Sachs and Neumann (2000) introduce this idea for local Fourier spectra and Nason (2013) applies it to wavelet spectra.

Implementation and Example. This test of stationarity can be found in the locits package (shortly to appear on CRAN) as the hwtos2 function. We will continue our example from above and apply the test to the exP data set.
First, load the locits package, then apply the hwtos2 test and print out the results:

library("locits")
ans <- hwtos2(exP)
ans

The results are:

Class 'tos' : Stationarity Object :
       ~~~~ : List with 9 components with names
nreject rejpval spvals sTS AllTS AllPVal alpha x xSD


summary(.):
----------
There are 441 hypothesis tests altogether
There were 11 FDR rejects
The rejection p-value was 0.0002681456
Using Bonferroni rejection p-value is 0.0001133787
And there would be 9 rejections.
Listing FDR rejects... (thanks Y&Y!)
P: 7 HWTlev: 0 indices on next line...[1] 1
P: 7 HWTlev: 1 indices on next line...[1] 1
P: 7 HWTlev: 2 indices on next line...[1] 1
P: 7 HWTlev: 4 indices on next line...[1] 2
P: 7 HWTlev: 5 indices on next line...[1] 3
P: 8 HWTlev: 0 indices on next line...[1] 1
P: 8 HWTlev: 1 indices on next line...[1] 1
P: 8 HWTlev: 4 indices on next line...[1] 2
P: 9 HWTlev: 0 indices on next line...[1] 1
P: 9 HWTlev: 1 indices on next line...[1] 1
P: 9 HWTlev: 4 indices on next line...[1] 2

As mentioned above the test performs a multiple hypothesis test. There are various ways to assess the significance of multiple hypothesis tests and the output above shows assessment via Bonferroni correction and false discovery rate (FDR). It indicates that 11 hypothesis were rejected against FDR assessment and 9 according to Bonferroni. In either case the composite null hypothesis of stationarity is rejected and multiple null hypotheses are rejected.

This style of stationarity test can also indicate where the nonstationarities are located in the series. This is because the test knows the scale and location of the Haar wavelet coefficients whose null hypotheses were rejected. The estimated nonstationarities can be plotted by locits by simple application of plot by:

plot(ans)

which results in the following plot.
Explosion P-wave for signal 14 nonstationarities

This plot shows the exP time series (as plotted above) in grey. Each red double-headed arrow corresponds to one of the 11 FDR nonstationarities identified by the test. The length of the arrow corresponds to the scale of the Haar wavelet coefficient whose null hypothesis was rejected and the location of that Haar wavelet coefficient is fixed by the mid-point of the arrow. The numbers 6, 7, 8 on the right-hand side of the plot correspond to the scale of the wavelet periodogram, j where the nonstationarity was found. It is particularly noticeable that most of the nonstationarities seem centred around the t=100 point where there appears to be a significant burst of power.

If not Stationary, then what?

The explosion data set is an extreme example of a time series that is almost certainly not stationary. It would be clearly inappropriate to use methods designed for stationary series on the explosion P wave. For example, the autocovariance structure of the series is clearly different at the beginning and end of the series. Also, it would not be appropriate to apply the regular spectrum (periodogram) estimator to the whole series as it would erroneously average out the differences in the series.

In this case, and others of nonstationarity, it would be better to use local methods of autocovariance or local spectral quantities to estimate the second-order structure of the series. Information on how to do this can be found in the sections on local autocovariance and local autocorrelation or local spectral analysis.

References

Nason, G.P. (2013) A test for second-order stationarity and approximate confidence intervals for localized autocovariances for locally stationary time series, Journal of the Royal Statistical Society, Series B, 75, (to appear).

Nason, G.P., von Sachs, R. and Kroisandt, G. (2000) Wavelet processes and adaptive estimation of the evolutionary wavelet spectrum., Journal of the Royal Statistical Society, Series B, 62, 271-292.

Priestley, M.B. and Subba Rao, T. (1969) A Test for Non-Stationarity of Time-Series. Journal of the Royal Statistical Society, Series B, 31, 140-149.

Shumway, R.H. and Stoffer, D.S. (2011) Time Series Analysis and Its Applications (With R Examples). Springer: New York.

von Sachs, R. and Neumann, M.H. (2000) A wavelet-based test for stationarity. Journal of Time Series Analysis, 21, 597-613.


<<< locally stationary time series page