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

For everyday use we often consider time series that have:

- a constant mean
- a constant variance
- an autocovariance that does not depend on time.

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.

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.

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

*Description.* The PSR test centres around the time-varying Fourier
spectrum *f _{t}(w)* where

*Y(t, w) = log{F _{t}(w)}*,

where *F _{t}(w)* is an estimator of

*E {Y(t, w)} = log{f _{t}(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)

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.

*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

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

*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")
```

The results are:

ans <- hwtos2(exP)

ans

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

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.

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.

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