*I'd like to thank Dr. Tom Starke for providing the inspiration for this article series. The code below is a modification of that which used to be found on his website leinenbock.com, which later became drtomstarke.com.*

A while back we began discussing statistical mean reversion testing. In that article we looked at a couple of techniques that helped us determine whether a time series was mean reverting or not. In particular we looked at the Augmented Dickey-Fuller Test and the Hurst Exponent. In this article we will consider another test for mean reversion, namely the **Cointegrated Augmented Dickey Fuller** (CADF) test.

Firsty, it should be noted that it is actually very difficult to find a directly tradable asset that possesses mean-reverting behaviour. For instance, equities broadly behave like GBMs and hence render the mean-reverting trade strategies relatively useless. However, there is nothing stopping us from creating a *portfolio* of price series that is stationary. Hence we can apply mean-reverting trading strategies to the portfolio.

The simplest form of mean-reverting trade strategies is the classic "pairs trade", which usually involves a dollar-neutral long-short pair of equities. The theory goes that two companies in the same sector are likely to be exposed to similar market factors, which affect their businesses. Occasionally their relative stock prices will diverge due to certain events, but will revert to the long-running mean.

Let's consider two energy sector equities Approach Resources Inc given by the ticker AREX and Whiting Petroleum Corp given by the ticker WLL. Both are exposed to similar market conditions and thus will likely have a stationary pairs relationship. We are now going to create some plots, using pandas and the Matplotlib libraries to demonstrate the cointegrating nature of AREX and WLL. The first plot (Figure 1) displays their respective price histories for the period Jan 1st 2012 to Jan 1st 2013.

**Fig 1 - Time series plots of AREX and WLL**

If we create a scatter plot of their prices, we see that the relationship is broadly linear (see Figure 2) for this period.

**Fig 2 - Scatter plot of AREX and WLL prices**

The pairs trade essentially works by using a linear model for a relationship between the two stock prices:

\begin{eqnarray} y(t) = \beta x(t) + \epsilon(t) \end{eqnarray}Where $y(t)$ is the price of AREX stock and $x(t)$ is the price of WLL stock, both on day $t$.

If we plot the residuals $\epsilon(t) = y(t) - \beta x(t)$ (for a particular value of $\beta$ that we will determine below) we create a new time series that, at first glance, looks relatively stationary. This is given in Figure 3:

**Fig 3 - Residual plot of AREX and WLL linear combination**

We will describe the code for each of these plots below.

## Cointegrated Augmented Dickey-Fuller Test

In order to statistically confirm whether this series is mean-reverting we could use one of the tests that we considered in the previous article, namely the Augmented Dickey-Fuller Test or the Hurst Exponent. However, neither of these tests will actually help us determine $\beta$, the hedging ratio needed to form the linear combination, they will only tell us whether, for a particular $\beta$, the linear combination is stationary.

This is where the Cointegrated Augmented Dickey-Fuller (CADF) test comes in. It determines the optimal hedge ratio by performing a linear regression against the two time series and then tests for stationarity under the linear combination.

### Python Implementation

We will now use Python libraries to test for a cointegrating relationship between AREX and WLL for the period of Jan 1st 2012 to Jan 1st 2013. We will use Yahoo Finance for the data source and Statsmodels to carry out the ADF test, as above.

The first task is to create a new file, `cadf.py`

, and import the necessary libraries. The code makes use of NumPy, Matplotlib, Pandas and Statsmodels. In order to correctly label the axes and download data from Yahoo Finance via pandas, we import the matplotlib.dates module and the pandas.io.data module. We also make use of the Ordinary Least Squares (OLS) function from pandas:

```
#!/usr/bin/python
# -*- coding: utf-8 -*-
# cadf.py
import datetime
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import pandas as pd
import pandas.io.data as web
import pprint
import statsmodels.tsa.stattools as ts
from pandas.stats.api import ols
```

The first function, `plot_price_series`

, takes a pandas DataFrame as input, with two columns given by the placeholder strings "ts1" and "ts2". These will be our pairs equities. The function simply plots the two price series on the same chart. This allows us to visually inspect whether any cointegration may be likely.

We use the Matplotlib dates module to obtain the months from the datetime objects. Then we create a figure and a set of axes on which to apply the labelling/plotting. Finally, we plot the figure:

```
# cadf.py
def plot_price_series(df, ts1, ts2):
months = mdates.MonthLocator() # every month
fig, ax = plt.subplots()
ax.plot(df.index, df[ts1], label=ts1)
ax.plot(df.index, df[ts2], label=ts2)
ax.xaxis.set_major_locator(months)
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b %Y'))
ax.set_xlim(datetime.datetime(2012, 1, 1), datetime.datetime(2013, 1, 1))
ax.grid(True)
fig.autofmt_xdate()
plt.xlabel('Month/Year')
plt.ylabel('Price ($)')
plt.title('%s and %s Daily Prices' % (ts1, ts2))
plt.legend()
plt.show()
```

The second function, `plot_scatter_series`

, plots a scatter plot of the two prices. This allows us to visually inspect whether a linear relationship exists between the two series and thus whether it is a good candidate for the OLS procedure and subsequent ADF test:

```
# cadf.py
def plot_scatter_series(df, ts1, ts2):
plt.xlabel('%s Price ($)' % ts1)
plt.ylabel('%s Price ($)' % ts2)
plt.title('%s and %s Price Scatterplot' % (ts1, ts2))
plt.scatter(df[ts1], df[ts2])
plt.show()
```

The third function, `plot_residuals`

, is designed to plot the residual values from the fitted linear model of the two price series. This function requires that the pandas DataFrame has a "res" column, representing the residual prices:

```
# cadf.py
def plot_residuals(df):
months = mdates.MonthLocator() # every month
fig, ax = plt.subplots()
ax.plot(df.index, df["res"], label="Residuals")
ax.xaxis.set_major_locator(months)
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b %Y'))
ax.set_xlim(datetime.datetime(2012, 1, 1), datetime.datetime(2013, 1, 1))
ax.grid(True)
fig.autofmt_xdate()
plt.xlabel('Month/Year')
plt.ylabel('Price ($)')
plt.title('Residual Plot')
plt.legend()
plt.plot(df["res"])
plt.show()
```

Finally, the procedure is wrapped up in a `__main__`

function. The first task is to download the OHLCV data for both AREX and WLL from Yahoo Finance. Then we create a separate DataFrame, df, using the same index as the AREX frame to store both of the adjusted closing price values. We then plot the price series and the scatter plot.

After the plots are complete the residuals are calculated by calling the pandas ols function on the WLL and AREX series. This allows us to calculate the $\beta$ hedge ratio. The hedge ratio is then used to create a "res" column via the formation of the linear combination of both WLL and AREX.

Finally the residuals are plotted and the ADF test is carried out on the calculated residuals. We then print the results of the ADF test:

```
# cadf.py
if __name__ == "__main__":
start = datetime.datetime(2012, 1, 1)
end = datetime.datetime(2013, 1, 1)
arex = web.DataReader("AREX", "yahoo", start, end)
wll = web.DataReader("WLL", "yahoo", start, end)
df = pd.DataFrame(index=arex.index)
df["AREX"] = arex["Adj Close"]
df["WLL"] = wll["Adj Close"]
# Plot the two time series
plot_price_series(df, "AREX", "WLL")
# Display a scatter plot of the two time series
plot_scatter_series(df, "AREX", "WLL")
# Calculate optimal hedge ratio "beta"
res = ols(y=df['WLL'], x=df["AREX"])
beta_hr = res.beta.x
# Calculate the residuals of the linear combination
df["res"] = df["WLL"] - beta_hr*df["AREX"]
# Plot the residuals
plot_residuals(df)
# Calculate and output the CADF test on the residuals
cadf = ts.adfuller(df["res"])
pprint.pprint(cadf)
```

The output of the code (along with the Matplotlib plots) is as follows:

```
(-2.9607012342275936,
0.038730981052330332,
0,
249,
{'1%': -3.4568881317725864,
'10%': -2.5729936189738876,
'5%': -2.8732185133016057},
601.96849256295991)
```

It can be seen that the calculated test statistic of -2.96 is smaller than the 5% critical value of -2.87, which means that we can reject the null hypothesis that there isn't a cointegrating relationship at the 5% level. Hence we can conclude, with a reasonable degree of certainty, that AREX and WLL possess a cointegrating relationship, at least for the time period sample considered.

### Why Statistical Testing?

Fundamentally, as far as algorithmic trading is concerned, the statistical tests outlined above are only as useful as the profits they generate when applied to trading strategies. Thus, surely it makes sense to simply evaluate performance at the strategy level, as opposed to the price/time series level? Why go to the trouble of calculating all of the above metrics when we can simply use trade level analysis, risk/reward measures and drawdown evaluations?

Firstly, any implemented trading strategy based on a time series statistical measure will have a far larger sample to work with. This is simply because when calculating these statistical tests, we are making use of each *bar* of information, rather than each *trade*. There will be far less round-trip trades than bars and hence the statistical significance of any trade-level metrics will be far smaller.

Secondly, any strategy we implement will depend upon certain parameters, such as look-back periods for rolling measures or z-score measures for entering/exiting a trade in a mean-reversion setting. Hence strategy level metrics are only appropriate *for these parameters*, while the statistical tests are valid for the underlying time series sample.

In practice we want to calculate both sets of statistics. Python, via the statsmodels and pandas libraries, make this extremely straightforward. The additional effort is actually rather minimal!