Currently, there exists several good tests for stationarity. Probably, the first test for stationarity was proposed by Priestley and Subba Rao (JRSSB, 1969). Since then, several other good tests have come on the market. Examples, include the test proposed in Pararoditis (Bernoulli, 2009), Dahlhaus (JoE, 2009, Handbook of Statistics, 2012), Dette, Preuss and Vetter (2011) and Nason (JRSSB, 2013) . Please note this is not an exhaustive list and if I have omitted a test send me an email.

Here we describe the test for stationarity proposed in Dwivedi and Subba Rao (JTSA, 2011) and for multivariate time series in Jentsch and Subba Rao (2013). This is the very basic code for the Gaussian univariate case. Below I briefly how it works.

These tests are based on the idea that for a second order stationary processes the correlation between the DFTs is close to zero (this can be understood by considering the Spectral representation of a second order time series, the DFTs can be treated as an approximation of the orthogonal increment process). Whereas for nonstationary time series there does exist correlation between the DFTs (if one could construct a spectral representation of a nonstationary time series, the increment process would not be orthogonal).

Both the test statistic proposed in Dwivedi and Subba Rao (2011) and Jentsch and Subba Rao are based on looking for correlations in DFTs. The test statistic proposed in Dwivedi and Subba Rao test is for univariate time series and is really designed for Gaussian time series. In the case of non-Gaussianity the test statistic contains an (fourth order cumulant) additional term, that needs to be estimated to ensure that the test does not lead to too many false positives. This issue is addressed in the Jentsch and Subba Rao (2013), where a bootstrap procedure is proposed to estimate this variance (this means the test can be applied to nonlinear time series). In addition, in Jentsch and Subba Rao (2013) we show that the correlations in the DFTs characterise the nonstationarity, therefore we propose a test statistic that can capture a wider class of alternatives (by inclusion of a lag parameter l).

Here is the basic code (assuming Gaussianity, but allowing for the lag parameter l) can be found at (it include both the DFT covariance, the test, though this can be modified and some basic options for plotting).

Personally, rather than just doing a test I like to `see' the covariances. Therefore, I would recommend plotting the DFT covariance for different l (see Jentsch and Subba Rao (2013) for the precise definition) over r (it's a bit like an ACF plot). Several large spikes (over the 1.96 lines) suggest nonstationarity in the time series. As a rule of thumb, if the time series is locally stationary the DFT covariance should decrease with the lag r. On the other hand, if the time series is periodically stationary one should see large spikes at evenly spaced distances. The number of spikes over the length of the time series correspond to the number of periods in the periodically stationary time series (in other words it it was periodically correlated with period 2, you should see two spikes between one and T). To understand how correlations in the DFT are related to periodically stationary processes one should look at the literature on harmonizable time series Goodman (1965) is a good reference.

Please note, that the above is really designed for Gaussian time series. In addition, a very crude estimator of the spectral density is used (rectangle kernel based on an average of 13 frequencies). For the general case, you should contact my collaborator, Carsten Jentsch, for the code.