How to Check Time Series Stationarity in Python

By Leo Smigel

Updated on

You can use visual inspection, global vs. local analysis, and statistics to analyze stationarity. The Augmented Dickey-Fuller (ADF) test is the most commonly used parametric test, and the Zivot-Andrews test is better than the ADF at detecting stationarity through structural breaks.

This post assumes you understand what stationarity is and why it’s important.

Note: I use matplotlib plots for analysis and Plotly charts for presentation in this post. The complete code for this post is on the Analyzing Alpha GitHub Repository.

Visual Inspection

It’s easy to see if a process is creating a stationary time series. If we see the mean or the distribution changing, it’s non-stationary.

Let’s first create a random walk model using Python, and then we’ll plot it in Plotly.

# Grab imports
import numpy as np
import pandas as pd

# Create index
periods = 2000
vol = .002
dti = pd.date_range("2021-01-01",
                    periods=periods,
                    freq="H")

# Generate normal distribution
random_walk = pd.DataFrame(index=dti,
                   data=np.random.normal(size=periods)  * vol)

# Plot it using a custom function
make_plot(random_walk,
          "Stationary "Random Walk" Time Series").show()

We can identify a stochastic process by analyzing a time series plot and gain an understanding of if there’s a shift in the mean or distribution through time. This is called global vs. local analysis, which we’ll discuss further below.

However, in practice, it’s trivial to use more advanced plot functions to identify a stationary process.

Decomposition

Time series decomposition separates the signal, trend, seasonality, and error automatically. Let’s decompose Apple’s revenue. We can see the clear seasonality and trend in this non-stationary data.

from statsmodels.tsa.seasonal import seasonal_decompose
sd = seasonal_decompose(apple_revenue_history)
sd.plot()

Autocorrelation Function

The autocorrelation and partial-autocorrelation functions analyze a data set for statistical significance between the first data point and prior data points.

A non-stationary time series data will show significance between itself and its lagged values, and that significance will decay to zero slowly as in the first plot.

The second plot, a stationary time series, will quickly drop to zero. The red lines are confidence intervals where values above or below the lines are more significant than two standard deviations.

from statsmodels.graphics.tsaplots import plot_acf
plot_acf(df2, lags=40, alpha=0.05)
plot_acf(df2.diff().dropna(), lags=40, alpha=0.05)

Sometimes it’s evident after visual inspection that we have stationary or non-stationary data. But even in these circumstances, we want to test our hypothesis using statistical tests to understand our series better.

Global vs. Local Analysis

We can break our time series into multiple segments and analyze the summary statistics of each against the time series or another partition to see if our time series data is changing through time.

def get_quads(df):
    quadlen = int(len(df1) * 0.25)
    ss = df[:quadlen].describe()
    ss[1] = df[quadlen:quadlen*2].describe()
    ss[2] = df[quadlen*2:quadlen*3].describe()
    ss[3] = df[quadlen*3:].describe()
    return ss

Random Walk

Notice how the random walk time series has a constant mean, variance, and the distribution between the min, percentages, and max are relatively consistent (no seasonality)?

random_walk = pd.DataFrame(index=dti,
                   data=np.random.random(size=periods) * vol)
get_quads(random_walk)
   0           1           2           3
count  500.000000  500.000000  500.000000  500.000000
mean     0.000136   -0.000054    0.000070   -0.000216
std      0.002067    0.001909    0.002008    0.002020
min     -0.005851   -0.005614   -0.006003   -0.006325
25%     -0.001263   -0.001344   -0.001262   -0.001492
50%      0.000095   -0.000094    0.000045   -0.000229
75%      0.001503    0.001225    0.001331    0.001073
max      0.007002    0.006702    0.006338    0.005969

Trend Stationary

You can see the mean is increasing in the trending time series.

trending = pd.DataFrame(index=dti,
                   data=np.random.random(size=periods) * vol).cumsum()
get_quads(trending)
                0           1           2           3
count  500.000000  500.000000  500.000000  500.000000
mean     0.253949    0.754164    1.245338    1.752112
std      0.148433    0.138083    0.153263    0.144376
min      0.001630    0.509963    0.989912    1.509835
25%      0.122746    0.637347    1.111457    1.629752
50%      0.261204    0.749095    1.243160    1.748031
75%      0.376252    0.878295    1.376870    1.879781
max      0.508350    0.989398    1.509197    2.000878

Volatile Time Series

The standard deviation is changing in the volatile time series.

varying = pd.DataFrame(index=dti,
                   data=np.random.normal(size=periods) * vol \
              * np.logspace(1,5,num=periods, dtype=int))
              0           1           2           3
count  500.000000  500.000000  500.000000  500.000000
mean     0.000510    0.082075   -0.206946   -1.990338
std      0.094116    0.935861   10.263379   99.674838
min     -0.664262   -3.887168  -46.277549 -472.434713
25%     -0.039151   -0.280147   -4.025722  -45.672561
50%     -0.001081    0.067228   -0.068825   -4.059768
75%      0.037203    0.476481    3.295115   35.828933
max      0.424776    5.000719   61.889090  455.791059

Seasonal Time Series

You can identify seasonality by analyzing the distribution through min, max, and the percentages in between. Notice how the 50% area is changing.

def simulate_seasonal_term(periodicity, total_cycles, noise_std=1.,
                           harmonics=None):
    duration = periodicity * total_cycles
    assert duration == int(duration)
    duration = int(duration)
    harmonics = harmonics if harmonics else int(np.floor(periodicity / 2))

    lambda_p = 2 * np.pi / float(periodicity)

    gamma_jt = noise_std * np.random.randn((harmonics))
    gamma_star_jt = noise_std * np.random.randn((harmonics))

    total_timesteps = 100 * duration # Pad for burn in
    series = np.zeros(total_timesteps)
    for t in range(total_timesteps):
        gamma_jtp1 = np.zeros_like(gamma_jt)
        gamma_star_jtp1 = np.zeros_like(gamma_star_jt)
        for j in range(1, harmonics + 1):
            cos_j = np.cos(lambda_p * j)
            sin_j = np.sin(lambda_p * j)
            gamma_jtp1[j - 1] = (gamma_jt[j - 1] * cos_j
                                 + gamma_star_jt[j - 1] * sin_j
                                 + noise_std * np.random.randn())
            gamma_star_jtp1[j - 1] = (- gamma_jt[j - 1] * sin_j
                                      + gamma_star_jt[j - 1] * cos_j
                                      + noise_std * np.random.randn())
        series[t] = np.sum(gamma_jtp1)
        gamma_jt = gamma_jtp1
        gamma_star_jt = gamma_star_jtp1
    wanted_series = series[-duration:] # Discard burn in

    return wanted_series

duration = 100 * 3
periodicities = [10, 100]
num_harmonics = [3, 2]
std = np.array([2, 3])
np.random.seed(8678309)

terms = []
for ix, _ in enumerate(periodicities):
    s = simulate_seasonal_term(
        periodicities[ix],
        duration / periodicities[ix],
        harmonics=num_harmonics[ix],
        noise_std=std[ix])
    terms.append(s)
terms.append(np.ones_like(terms[0]) * 10.)
seasonal = pd.DataFrame(index=dti[:duration],data=np.sum(terms, axis=0))

                0            1            2            3
count    75.000000    75.000000    75.000000    75.000000
mean    328.055459  -280.450561  -149.840302   120.579444
std     733.650363   750.299351   958.572453  1007.596395
min    -863.464617 -1694.424423 -1815.786373 -1733.291362
25%    -239.600008  -770.199418  -828.906000  -694.550888
50%      72.425299  -346.127307  -405.201792    59.358360
75%     923.255524    59.875525   485.926284   945.035545
max    1813.015082  1642.090923  1916.794045  1832.084148

Parametric Tests

Statistical tests identify specific types of stationarity. Most parametric tests analyze if the time series has a unit root, which is a systematic pattern that is random in the presence of serial correlation — which is a fancy way of saying a lagged version of itself. The reason it’s called a unit root test is due to the math of the process.

Testing if the time series has a unit root and is therefore not stationary is the null hypothesis, which is a fancy way of saying the commonly accepted belief. We want to reject the null hypothesis with a level of certainty to state the time series is stationary — or more accurately, the process generating the time series is stationary.

Let’s look at the most common unit root tests.

Dickey-Fuller

The Dickey-Fuller test is the first statistical test that analyzes if a unit root exists in an autoregressive model of a time series.

It runs into issues with serial correlation, which is why there’s an Augmented Dickey-Fuller test.

Augmented Dickey-Fuller (ADF)

The Augmented Dickey-Fuller tests for a unit root in a univariate process in the presence of serial correlation. The ADF test handles more complex models and is the typical go-to for most analysts.

Non-Stationary Process

Let’s use the ADF test on Apple’s revenue, which is non-stationary.

from statsmodels.tsa.stattools import adfuller
t_stat, p_value, _, _, critical_values, _  = adfuller(df6['observed'].values, autolag='AIC')
print(f'ADF Statistic: {t_stat:.2f}')
for key, value in critical_values.items():
     print('Critial Values:')
     print(f'   {key}, {value:.2f}')
ADF Test Statistic: -2.12
Critial Values:
   1%, -3.69
Critial Values:
   5%, -2.97
Critial Values:
   10%, -2.63

You can compare the ADF test statistic of -2.12 against the critical values. We see that the test statistic is less than all of the critical values, so we cannot reject the null hypothesis — in other words, as we already knew, Apple’s revenue is non-stationary; we see trends, and it’s mean and variance are changing.

Stationary Process

Now let’s do the same for the random walk data, which exhibits stationarity.

from statsmodels.tsa.stattools import adfuller
result = adfuller(random_walk[0].values, autolag='AIC')
t_stat, p_value, _, _, critical_values, _  = adfuller(random_walk[0].values, autolag='AIC')
print(f'ADF Statistic: {t_stat:.2f}')
for key, value in critical_values.items():
     print('Critial Values:')
     print(f'   {key}, {value:.2f}')

print(f'\np-value: {p_value:.2f}')
print("Non-Stationary") if p_value > 0.05 else print("Stationary")
ADF Statistic: -43.83
Critial Values:
   1%, -3.43
Critial Values:
   5%, -2.86
Critial Values:
   10%, -2.57

p-value: 0.00
Stationary

Notice that this time I’ve included the p-value. When a p-value is greater than 5%, we accept the null hypothesis. We also see an extreme ADF test statistic from this random walk model.

I like the critical value approach compared to the p-value as I can see to what degree I can reject that the time series is not stationary.

Kwiakowski-Phillips-Schmidt-Shin (KPSS)

The Kwiatkowski-Phillips-Schmidt-Shin (KPSS) tests for the null hypothesis that the time series process is stationary or trend stationarity. In other words, the KPSS test null hypothesis is the opposite of the ADF test.

Exhibits Non-Stationarity

from statsmodels.tsa.stattools import kpss
t_stat, p_value, _, critical_values = kpss(df6['observed'].values, nlags='auto')

print(f'ADF Statistic: {t_stat:.2f}')
for key, value in critical_values.items():
     print('Critial Values:')
     print(f'   {key}, {value:.2f}')

print(f'\np-value: {p_value:.2f}')
print("Stationary") if p_value > 0.05 else print("Non-Stationary")

ADF Statistic: 0.91
Critial Values:
   10%, 0.35
Critial Values:
   5%, 0.46
Critial Values:
   2.5%, 0.57
Critial Values:
   1%, 0.74

p-value: 0.01
Non-Stationary

Exhibits Stationarity

from statsmodels.tsa.stattools import kpss
t_stat, p_value, _, critical_values = kpss(random_walk[0].values, nlags='auto')

print(f'ADF Statistic: {t_stat:.2f}')
for key, value in critical_values.items():
     print('Critial Values:')
     print(f'   {key}, {value:.2f}')

print(f'\np-value: {p_value:.2f}')
print("Stationary") if p_value > 0.05 else print("Non-Stationary")
from statsmodels.tsa.stattools import kpss
t_stat, p_value, _, critical_values = kpss(random_walk[0].values, nlags='auto')

print(f'ADF Statistic: {t_stat:.2f}')
for key, value in critical_values.items():
     print('Critial Values:')
     print(f'   {key}, {value:.2f}')

print(f'\np-value: {p_value:.2f}')
print("Non-Stationary") if p_value > 0.05 else print("Stationary")

Zivot and Andrews

Zivot-Andrews tests for the same thing as the ADF and KPSS tests for the presence of a structural break.

Let’s add a break to our random walk model. Notice how the variance and cyclicality properties are relatively constant, but there’s a significant shift in the mean.

Now let’s run both an ADF test and a Zivot-Andrews to test for non-stationarity.

t_stat, p_value, _, _, critical_values, _  = adfuller(stationary_with_break[0].values, autolag='AIC')
print(f'ADF Statistic: {t_stat:.2f}')
for key, value in critical_values.items():
     print('Critial Values:')
     print(f'   {key}, {value:.2f}')

print(f'\np-value: {p_value:.2f}')
print("Non-Stationary") if p_value > 0.05 else print("Stationary")
ADF Statistic: -1.05
Critial Values:
   1%, -3.43
Critial Values:
   5%, -2.86
Critial Values:
   10%, -2.57

p-value: 0.73
Non-Stationary
from statsmodels.tsa.stattools import zivot_andrews
t_stat, p_value, critical_values, _, _ = zivot_andrews(stationary_with_break[0].values)
print(f'Zivot-Andrews Statistic: {t_stat:.2f}')
for key, value in critical_values.items():
     print('Critial Values:')
     print(f'   {key}, {value:.2f}')

print(f'\np-value: {p_value:.2f}')
print("Non-Stationary") if p_value > 0.05 else print("Stationary")
Zivot-Andrews Statistic: -60.35
Critial Values:
   1%, -5.28
Critial Values:
   5%, -4.81
Critial Values:
   10%, -4.57

p-value: 0.00
Stationary

We see different results between the two stationarity tests. The ADF test fails to incorporate the structural break while the Zivot-Andrew identifies the series as stationary.

The Bottom Line

The process for testing for stationarity in time series data is relatively straightforward. You analyze the information visually, perform a decomposition, review the summary statistics, and then select a parametric test to gain confidence in your assumptions.

The Augmented Dickey-Fuller or ADF test is the most commonly used test for stationarity; however, it’s not always the best. To analyze stationarity in time series with structural breaks, you’ll want to use the Zivot-Andrews test. As always, it’s essential to have an intuitive understanding of the problem set and to have a deep knowledge of the tools at your disposal.

Leave a Comment

Get Exclusive Tips

Sign up for the newsletter to get tips and strategies I don't share anywhere else.

Analyzing Alpha
2009 Mackenzie Way
Suite 100
Cranberry Twp, PA 16066

P: 412-212-3240
E: info [ at ] analyzingalpha.com

© Analyzing Alpha

Contact Privacy DCMA Terms Disclosure Disclaimer Sitemap

Analyzing Alpha is reader-supported. We may earn an affiliate commission when you buy through links on our site. As an Amazon Associate, we earn from qualifying purchases.