Introduction to Equity Pair Trading
Overview
This article will cover the Equity Pair trading use case and go from a single pair example to an industrialised search process with commonly used statistical tests, combinatoric generators and computed in a brute force fashion using vectorization and parallelization where appropriate. Finally, we will use some filtering and visualization tools to eyeball the results and suggest potential next steps. Link to github source
Sections
Introducing More Statistical Rigour
Industrializing Our Search Process
Calculation Details — Vectorization or Parallisation
Final Visual Check — Seaborn FacetGrid
Pre-requisites:
Thomson Reuters Eikon / Refinitiv Workspace with access to Eikon Data APIs (Free Trial Available)
Python 2.x/3.x
Required Python Packages: eikon, pandas, numpy, matplotlib, statsmodels, swifter, seaborn
Introduction to Pair Trading
Pair trading is a time-honoured market-neutral type of statistical arbitrage trading strategy. It is market-neutral because once you have identified a suitable pair (two stocks for example) — you would be long one stock and short the other. This means that whatever the broader market does ie goes up or down — you are theoretically somewhat insulated. Another interesting characteristic of this strategy is its low capital outlay as the short side of the trade funds the long side of the trade — so this type of approach is very efficient from a capital as well as a risk point of view.
So how can we find these pairs? Well there are many approaches — some concentrate on finding companies in similar sectors or business areas or geographies or markets. Some focus on identifying groups of fundamental criteria.
Lets start by taking two well-known automakers — say Tesla and General Motors in the United States and see what we can find out about the price movement of these.
Single Pair Example
First lets import some of the packages we will be using.
%matplotlib inline
import eikon as ek
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.tsa.stattools as ts
from statsmodels.tsa.stattools import adfuller
import itertools
import swifter
import seaborn as sns
ek.set_app_key('YOUR APP KEY HERE')/Users/jasonram/anaconda3/lib/python3.7/site-packages/dask/dataframe/utils.py:13: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
import pandas.util.testing as tm
First, lets load some timeseries data for both stocks using the ek.get_timeseries function.
rics = ['TSLA.O','GM']
s = '2013-12-23'
e = '2019-12-08'
i = 'daily'
df = ek.get_timeseries(rics, # RICs
fields=['CLOSE'], # fields to be retrieved
start_date=s, # start time
end_date=e, # end time
interval=i)
df
Next we want to get rid of any non-complete (ie N/A) rows as thiis will not play well with some our routines downstream.
df.dropna(inplace=True)
dfLets add both Spread and Ratio columns as we will want to chart these — particularly the Ratio — as this is the series we will be hoping to trade on and we also print a Pearson Correlation Coefficient.
df['Spread'] = df['TSLA.O'] - df['GM']
df['Ratio'] = df['TSLA.O'] / df['GM']
mean = df['Ratio'].mean()
corr = df['TSLA.O'].corr(df['GM'])
print(corr)
df0.44026737012835365
The code below charts our ratio between TSLA.O and GM. We add a mean value of the ratio — this is the mean to which we hope the ratio will revert to. Adding a mean and standard deviation thresholds onto the chart can give us interesting entry and exit signals over time.
fig = plt.figure()
ax = plt.axes()
x = df.index.values
plt.plot(x, df['Ratio'], label='Ratio')
plt.plot(x, (df['Ratio'].values.mean().repeat(len(x))), label='Mean')
stdplus = df['Ratio'].values.mean() + (1.5*df['Ratio'].values.std())
stdminus = df['Ratio'].values.mean() - (1.5*df['Ratio'].values.std())
#plt.plot(x, (df['Ratio'].values.mean().repeat(len(x))), label='SD1')
plt.plot(x, (stdplus.repeat(len(x))), label='SD+1.5')
plt.plot(x, (stdminus.repeat(len(x))), label='SD-1.5');
From the above chart we can see that adding the mean and standard deviation thresholds onto the chart might give us interesting entry and exit signals over time. The lower bound of 1.5 times standard deviation seems to be a pretty reliable entry point for this pair — and has been tested on 5 occasions over the last decade — all with minimal drawdown and large reversion to mean in each case.
The upper bound of 1.5 times standard deviation seems to be a reasonable threshold till 2017 — but since then it has been breached on a number of occasions — so we may wish to expand the upper bound to 2 standard deviations or more.
Whilst interesting — this approach suffers from 2 drawbacks — firstly, its still a bit imprecise and secondly it relies on me knowing what the appropriate pair should be. Lets deal with the first of these drawbacks.
Introduce More Statistical Rigour
Lets consider some more rigorous statistical tests.
For this particular strategy we want to see a relatively high level of correlation between the two stock price series. Correlation is a measure of directional relatedness — ie if GM rises does Tesla rise. The range of correlation is from -1 to +1. Correlation in non-stationary timeseries may just be reflecting autoregression (AR processes) — so we need conduct some additional tests.
It is important that the ratio series is mean reverting. It is this property that we wish to identify & exploit. To help us in this task we can test for stationarity of the series. We will use the Augmented Dickey Fuller (ADF) test to see if we can reject the null hypothesis that the Ratio timeseries has a unit root (ie the series is non-stationary).
Another statistical property we wish to be present is that the two stocks should be cointegrated. This is the co-movement of both time series (which could be non-stationary). We will use the Engle-Granger two-step cointegration test to see if we can reject the null hypothesis that the pair are NOT cointegrated. This is important as both stocks in the pair should display some quality of cointegratedness for this type of strategy to be exploitable. Thankfully the excellent statsmodels package provides this functionality for us.
print(adfuller(df['Ratio'])[0:5])(-3.5915050452251873, 0.00592552532947228, 0, 1499, {'1%': -3.4347199356122493, '5%': -2.86347004827819, '10%': -2.567797534300163})result = ts.coint(df['TSLA.O'], df['GM'])
print(result)(-3.6980090738968627, 0.018433668450996164, array([-3.90376106, -3.34020915, -3.04728056]))
From the ADF output we can see that the critical value is -3.59 with a p-value of 0.005. The critical value is less than the t-values at the 1%,5% & 10% confidence levels and the p-value is less than 0.05 (95% confidence level). So we can reject the null hypothesis.
From the cointegration output we can see that the critical value is -3.69 with a p-value of 0.018. The critical value is less than the t-values at the 5% & 10% confidence levels, but not at the 1% level and the p-value of 0.018 is less than 0.05 (95% confidence interval). So we can reject the null hypothesis.
Industrialising Our Search Process
So we have seen how it works for one pair — lets try and deal with the second drawback and try to industrialise this search for potential pairs to a whole universe of stocks — say the FTSE 100 index — but it could of course be any universe or list of instruments. So lets get a list of constituents for our chosen index using the ek.get_data function. The get_data function can be passed a chain RIC or a list of RICs — here we pass a chain RIC for the FTSE 100.
rics, err = ek.get_data('0#.FTSE','CF_Name')
rics
Here we take our list of RIC codes for FTSE 100 index and then we add a start and end date parameter, select the daily interval and the ‘CLOSE’ field and get the underlying time series for each RIC in the list of instruments. We then build a larger dataframe by concatenating each Price series to the base dataframe, whilst renaming the column to the relevant RIC.
instruments = rics['Instrument'].astype(str).values.tolist()
s = '2010-01-06'
e = '2019-12-08'
inv = 'daily'
data1 = pd.DataFrame()
for i in instruments:
df1 = ek.get_timeseries(i, # RICs
fields=['CLOSE'], # fields to be retrieved
start_date=s, # start time
end_date=e, # end time
interval=inv)
df1.rename(columns = {'CLOSE': i}, inplace = True)
if len(data1):
data1 = pd.concat([data1, df1], axis=1)
else:
data1 = df1So lets have a look at our larger base dataframe.
data1
We need to check our dataframe to see if the are any NAN values present as this will not play well with the statsmodels package we wish to use downstream.
data2=data1.copy()Drop any rows that are entirely blank — I noticed a few of these being returned
data2.dropna(how='all', inplace=True)Drop columns which have a nan in them
nalist = data2.columns[data2.isna().any()].tolist()
nalist['ICAG.L',
'EVRE.L',
'MNG.L',
'AUTOA.L',
'GVC.L',
'JETJ.L',
'CCH.L',
'AVST.L',
'POLYP.L',
'GLEN.L',
'AVV.L',
'OCDO.L']data2.drop(nalist, axis=1, inplace=True)
data2
Calculation details — Vectorization or Parallelization — decisions, decisions? Who cares with swifter!
We now wish to derive a list of all combinations of RICs — we can do this using the excellent itertools.combinations combinatoric generator function — with no repeated elements.
Then for each of the RIC pairs/ratios in this list we want to calculate and store the output from the ADF test, Engle-Granger cointegration test and Pearsons correlation coefficient. As this is a somewhat brute force approach I tried to follow some of the performance advice I had seen — and finally opted for the syntactically compact apply method using a lambda function. This is supposed to be more efficient than simple looping, iterrows or itertuples.
Still not satisfied — I searched for some parallelisation or speedup potential and ran across the Switfer package. It seemed to offer a great deal — supercharging any vanilla df.apply function by either:
“a) trying to vectorize the function where possible or, where vectorization is not possible:
b) automatically decide on the basis of speed whether to use Dask Parallel Processing or the vanilla Pandas apply”
Well I must say — what fantastic promise — and just one package to import and a single modifier to add — instead of df.apply I use df.swifter.apply. What is even better is we have a progress bar while we wait — super awesome. I really ❤ python sometimes. I trust you are digging the slackness here :)
After some power has been consumed and heat generated — we hopefully end up with some light in the shape of complete combinatoric calculations for over 4000 pairs.
all_pairs = list(itertools.combinations(data2.columns,2)) #combinatoric generator no repeated elements, no reverses
coint_df = pd.DataFrame(all_pairs,columns =['RIC1','RIC2'])
coint_df['ADF'] = coint_df.swifter.apply(lambda row: adfuller((data2[row.RIC1]/data2[row.RIC2])), axis=1)
coint_df['EG'] = coint_df.swifter.apply(lambda row: ts.coint(data2[row.RIC1], data2[row.RIC2]), axis=1)
coint_df['PCorr'] = coint_df.swifter.apply(lambda row: np.corrcoef(data2[row.RIC1],data2[row.RIC2]), axis =1)
coint_dfHBox(children=(FloatProgress(value=0.0, description='Pandas Apply', max=3916.0, style=ProgressStyle(descriptio…
HBox(children=(FloatProgress(value=0.0, description='Pandas Apply', max=3916.0, style=ProgressStyle(descriptio…
HBox(children=(FloatProgress(value=0.0, description='Pandas Apply', max=3916.0, style=ProgressStyle(descriptio…
We now wish to narrow this list down by selecting a subset that might display promising characteristics but first we need to unpack selected values from the returns of the statistical tests we ran previously.
Here again — we use the syntatically compact apply method with a lambda function to unpack our desired values — heck why not speed this up as well by adding some swifter super-awesomeness!
coint_df['ADF_Crit'] = coint_df.swifter.apply(lambda row: row.ADF[0], axis=1)
coint_df['ADF_t5'] = coint_df.swifter.apply(lambda row: list(row.ADF[4].values())[1],axis=1)
coint_df['ADF_pVal'] = coint_df.swifter.apply(lambda row: row.ADF[1], axis=1)
coint_df['EG_Crit'] = coint_df.swifter.apply(lambda row: row.EG[0], axis=1)
coint_df['EG_t5'] = coint_df.swifter.apply(lambda row: row.EG[2][1], axis=1)
coint_df['EG_pVal'] = coint_df.swifter.apply(lambda row: row.EG[1], axis=1)
coint_df['PCorrVal'] = coint_df.swifter.apply(lambda row: row.PCorr[0][1], axis=1)HBox(children=(FloatProgress(value=0.0, description='Pandas Apply', max=3916.0, style=ProgressStyle(descriptio…
HBox(children=(FloatProgress(value=0.0, description='Pandas Apply', max=3916.0, style=ProgressStyle(descriptio…
HBox(children=(FloatProgress(value=0.0, description='Pandas Apply', max=3916.0, style=ProgressStyle(descriptio…
HBox(children=(FloatProgress(value=0.0, description='Pandas Apply', max=3916.0, style=ProgressStyle(descriptio…
HBox(children=(FloatProgress(value=0.0, description='Pandas Apply', max=3916.0, style=ProgressStyle(descriptio…
HBox(children=(FloatProgress(value=0.0, description='Pandas Apply', max=3916.0, style=ProgressStyle(descriptio…
HBox(children=(FloatProgress(value=0.0, description='Pandas Apply', max=3916.0, style=ProgressStyle(descriptio…coint_df
Filtering our universe down to identify potentially attractive candidates
Now we have unpacked all the relevant values we can simply filter the larger dataframe for pairs that satisfy our specific — in this case — additive, conditions below:
a) ADF Critical Value being less than t-test value at the 95% confidence level AND
b) ADF pVal being less than 0.02 AND
c) EG Critical Value being less than t-test value at the 95% confidence level AND
d) EG pVal being less than 0.02 AND
e) PCorrVal being greater than 0.4
This gives us a more manageable list of potential candidates. A point to note I have selected the threshold values intuitively — but you could widen or tighten these based on your requirements.
low_coint = coint_df[(coint_df['ADF_Crit'] < coint_df['ADF_t5']) & (coint_df['ADF_pVal'] < 0.03) &
(coint_df['EG_Crit'] < coint_df['EG_t5']) & (coint_df['EG_pVal'] < 0.03) & (coint_df['PCorrVal'] > 0.4)]
low_coint.reset_index(drop=True)
Final visualisation check — Seaborn FacetGrid
Now we have our promising list of pairs — lets move onto the final stage which for me is always visualisation. We will retrieve the ratio for each pair as a series and store to new dataframe then hand off this data to the excellent Seaborn facetgrid chart.
pairlist = pd.DataFrame()
for row in low_coint.itertuples():
pairlist[str(row.RIC1) + '/' + str(row.RIC2)] = data2[row.RIC1] / data2[row.RIC2]
pairlist
fig = plt.figure(figsize=(20,20))
for i,col in enumerate(pairlist.columns):
ax=plt.subplot(10,5,i+1)
pairlist[[col]].plot(ax=ax)
plt.show()
Conclusion
In this article — we have covered the basic idea behind pair-trading — which is a time-honoured market-neutral type of statistical arbitrage (stat-arb) strategy. We explained some of the salient points using a single pair as an example. Given the 2-fold weakness in this manual strategy — we needed to know what an appropriate pair might be — and the lack of a robust statistical selection process — we then embarked on a journey to industrialise our search process using brute force calculations of all relevant combinations of index constituents — in our case the FTSE 100.
Along the way we implemented iterative API calls to build a larger dataframe than a single call would allow for. We also implemented various statistical testing routines to try to add some robustness to our search and then used some simple filters to narrow down our universe to the most promising candidates. Along the way we met a new friend — Swifter — which offers great potential speedups via vectorisation or Dask parallelisation with basically no additional knowledge, code or effort. To be honest I already thought multiprocessing/multiprocess was fairly concise and straightforward — but Swifter really blew me away. I strongly recommend adding it to your toolkit for pandas intensive calculation workflows.
Finally — we graphed the filtered pair ratios en-masse so we could see the rough shapes of the ratios using the easy to use and efficient FacetGrid charting.
From the charts above we can see some very promising candidates which appear to have stable mean-reverting shapes. Annoyingly — also present are a number of shapes which although seemingly satisfying our statistical criterial, do not appear to be of an appropriate character. We would need to investigate this further which is beyond the scope of this article — suffice to say there are statistical exceptions to our basic rule filters which we may need to take account for or modify.
From here we could select appropriate pairs and then go to the next level of selecting appropriate +/- standard-deviation thresholds — by eye — but better still automating it all via algorithm. Again that is beyond the scope of this article.
For my eye — the most promising pairs from the list are INF.L/ULVR.L, UU.L/Land.L, WTB.L/LGEN.L, BP.L/RSA.L, PHNX.L/BAES.L, SPX.L/ICP.L. Of course — these are not in any way to be construed as recommendations or solicitations — just candidates worthy of further research and testing — again beyond the scope of this article.
I hope to have introduced the basic concepts of this type of strategy and provided you with a practical codebase to investigate further.
Further Resources for Eikon Data API
For Content Navigation in Eikon — please use the Data Item Browser Application: Type ‘DIB’ into Eikon Search Bar.
