PAIRS TRADING WITH A KALMAN FILTER¶
In this notebook we are going look at the concept of building a trading strategy backtest based on mean reverting, co-integrated pairs of assets (Stock and ETFs). So to restate the theory in in terms of US equities, assets that are statistically co-integrated move in a way that means when their prices start to diverge by a certain amount (i.e. the spread between the 2 assets prices increases), we would expect that divergence to
eventually revert back to the mean. In this instance we would look to sell the outperforming stock,and buy the under performing stock under the notion that the under performing stock would eventually recover with the over performing stock and rise in price, or vice versa the over performing stock would in time suffer from the same downward pressure of the under performing stock and fall in relative value.
Hence, pairs trading is a market neutral trading strategy enabling traders to profit from virtually any market conditions: Bull Markets, Bear Markets, or Sideways Markets.
So in our search for co-integrated assets, economic theory would suggest that we are more likely to find pairs of that are driven by the same factors, or similar business practices. After all, it is logical to consider 2 assets in the same industry to be similar products, to be at the mercy of the same general ups and downs of the volatile market.
Kalman Filter¶
So what is aKalman Filter? Well this site (http://www.bzarg.com/p/how-a-kalman-filter-works-in-pictures/) explains and states the following:
You can use a Kalman filter in any place where you have uncertain information about some dynamic system, and you can make an educated guess about what the system is going to do next. Even if messy reality comes along and interferes with the clean motion you guessed about, the Kalman filter will often do a very good job of figuring out what actually happened. And it can take advantage of correlations between crazy phenomena that you maybe wouldn’t have thought to exploit!
Kalman filters are ideal for systems which are continuously changing. They have the advantage that they are light on memory (they don’t need to keep any history other than the previous state), and they are very fast, making them well suited for real time problems and embedded systems.
So lets start to import the relevant modules we will need for our strategy backtest:
Imports¶
from time import time
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib.cm as cm
from scipy import stats
import datetime as dt
import pandas as pd
import math
import os.path
import time
import json
import requests
import pandas_market_calendars as mcal
from datetime import timedelta, datetime
from dateutil import parser
import seaborn as sns
import matplotlib as mpl
import quantstats as qs
import statsmodels.api as sm
from pykalman import KalmanFilter
from math import sqrt
import warnings
import ffn
import pyfolio as pf
from pandas_datareader import data as web
pd.set_option('display.max_columns', None)
warnings.filterwarnings('ignore')
Download ticker data¶
symbols = ['GDX','GDXJ','GLD', 'AAPL','GOOGL', 'FB','TWTR','AMD',
'NVDA','CSCO', 'ORCL', 'ATVI', 'TTWO', 'EA', 'HYG',
'LQD', 'JNK', 'SLV', 'USLV', 'SIVR', 'USO', 'UWT',
'QQQ', 'SPY', 'VOO', 'VDE', 'VTI', 'EMLP', 'VDC',
'FSTA', 'KXI', 'IBB', 'VHT','VNQ', 'IYR', 'MSFT',
'PG', 'TMF', 'UPRO', 'WFC', 'JPM', 'GS', 'CVX',
'XOM', 'INTC', 'COST', 'WMT', 'T', 'VZ', 'CMCSA', 'AMZN']
def get_symbols(symbols,data_source,ohlc,begin_date=None,end_date=None):
out = []
new_symbols = []
for symbol in symbols:
df = web.DataReader(symbol, data_source,begin_date, end_date)
df = df[ohlc]
new_symbols.append(symbol)
out.append(df.astype('float'))
data = pd.concat(out, axis = 1)
data.columns = new_symbols
data = data.dropna(axis=1)
return data.dropna(axis=1)
start = pd.Timestamp('2014-01-01')
end = pd.Timestamp('2020-03-05')
prices = get_symbols(symbols,data_source='yahoo',ohlc='Close',\
begin_date=start,end_date=end)
Plot the resulting DataFrame of price data just to make sure we have what we need and as a quick sanity check:
combo = prices.copy()
combo.index = pd.DatetimeIndex(combo.index)
combo.head()
combo.info()
num_stocks = len(combo.columns)
print('Number of Stocks =', num_stocks)
Plot price series¶
n_secs = len(combo.columns)
colors = cm.rainbow(np.linspace(0, 1, n_secs))
combo.div(combo.iloc[0,:]).plot(color=colors, figsize=(12, 6))# Normalize Prices
plt.title('All Stocks Normalized Price Series')
plt.xlabel('Date')
plt.ylabel('Price (USD$)')
plt.grid(b=None, which=u'major', axis=u'both')
plt.legend(bbox_to_anchor=(1.01, 1.1), loc='upper left', ncol=1)
plt.show();
What Is Cointegration?¶
The most common test for Pairs Trading is the cointegration test. Cointegration is a statistical property of two or more time-series variables which indicates if a linear combination of the variables is stationary.
Stationary process: parameters such as mean and variance also do not change over time.
Ok so it looks from the chart as if we have downloaded price data for around 50 assets; this should be more than enough to find at least a couple of co-integrated pairs to run our backtest over.
We will now define a quick function that will run our assets, combining them into pairs one by one and running co-integration tests on each pair. That result will then be stored in a matrix that we initialise,
and then we will be able to plot that matrix as a heatmap. Also, if the co-integration test meets our threshold statistical significance (in our case 5%), then that pair of tickers will be stored in a list for later retrieval.
# NOTE CRITICAL LEVEL HAS BEEN SET TO 5% FOR COINTEGRATION TEST
def find_cointegrated_pairs(dataframe, critial_level = 0.05):
n = dataframe.shape[1] # the length of dateframe
pvalue_matrix = np.ones((n, n)) # initialize the matrix of p
keys = dataframe.columns # get the column names
pairs = [] # initilize the list for cointegration
for i in range(n):
for j in range(i+1, n): # for j bigger than i
stock1 = dataframe[keys[i]] # obtain the price of "stock1"
stock2 = dataframe[keys[j]]# obtain the price of "stock2"
result = sm.tsa.stattools.coint(stock1, stock2) # get conintegration
pvalue = result[1] # get the pvalue
pvalue_matrix[i, j] = pvalue
if pvalue < critial_level: # if p-value less than the critical level
pairs.append((keys[i], keys[j], pvalue)) # record the contract with that p-value
return pvalue_matrix, pairs
Let’s now run our data through our function, save the results and plot the heatmap:
df = combo
binance_symbols = df.columns
# Set up the split point for our "training data" on which to perform the co-integration test (the remaining dat awill be fed to our backtest function)
split = int(len(df) * 0.3)
# Run our dataframe (up to the split point) of ticker price data through our co-integration function and store results
pvalue_matrix, pairs = find_cointegrated_pairs(df[:split])
# Convert our matrix of stored results into a DataFrame
pvalue_matrix_df = pd.DataFrame(pvalue_matrix)
# Use Seaborn to plot a heatmap of our results matrix
sns.clustermap(pvalue_matrix_df, xticklabels=binance_symbols,yticklabels=binance_symbols, figsize=(12, 12))
plt.title('Stock P-value Matrix')
plt.tight_layout()
plt.show();
So we can see from the very dark squares that it looks as though there are indeed a few pairs of assets who’s co-integration score is below the 5% threshold
hardcoded into the function we defined. To see more explicitly which pairs these are, let’s print out our list of stored pairs that was part of the fucntion results we stored:
for pair in pairs:
print("Asset {} and Asset {} has a co-integration score of {}".format(pair[0],pair[1],round(pair[2],4)))
We will now use the “pykalman” module to set up a couple of functions that will allow us to generate Kalman filters which we will apply to our data and in turn our regression that is fed the said data.
def KalmanFilterAverage(x):
# Construct a Kalman filter
kf = KalmanFilter(transition_matrices = [1],
observation_matrices = [1],
initial_state_mean = 0,
initial_state_covariance = 1,
observation_covariance=1,
transition_covariance=.01)
# Use the observed values of the price to get a rolling mean
state_means, _ = kf.filter(x.values)
state_means = pd.Series(state_means.flatten(), index=x.index)
return state_means
# Kalman filter regression
def KalmanFilterRegression(x,y):
delta = 1e-3
trans_cov = delta / (1 - delta) * np.eye(2) # How much random walk wiggles
obs_mat = np.expand_dims(np.vstack([[x], [np.ones(len(x))]]).T, axis=1)
kf = KalmanFilter(n_dim_obs=1, n_dim_state=2, # y is 1-dimensional, (alpha, beta) is 2-dimensional
initial_state_mean=[0,0],
initial_state_covariance=np.ones((2, 2)),
transition_matrices=np.eye(2),
observation_matrices=obs_mat,
observation_covariance=2,
transition_covariance=trans_cov)
# Use the observations y to get running estimates and errors for the state parameters
state_means, state_covs = kf.filter(y.values)
return state_means
def half_life(spread):
spread_lag = spread.shift(1)
spread_lag.iloc[0] = spread_lag.iloc[1]
spread_ret = spread - spread_lag
spread_ret.iloc[0] = spread_ret.iloc[1]
spread_lag2 = sm.add_constant(spread_lag)
model = sm.OLS(spread_ret,spread_lag2)
res = model.fit()
halflife = int(round(-np.log(2) / res.params[1],0))
if halflife <= 0:
halflife = 1
return halflife
Now let us define our main “Backtest” function that we will run our data through. The fucntion takes one pair of tickers at a time, and then returns several outputs, namely the DataFrame of cumulative returns,
the Sharpe Ratio and the Compound Annual Growth Rate (CAGR). Once we have defined our function, we can iterate over our list of pairs and feed the relevant data, pair by pair, into the function, storing the outputs for each pair for later use and retrieval.
def backtest(df,s1, s2):
#############################################################
# INPUT:
# DataFrame of prices (df)
# s1: the symbol of asset one
# s2: the symbol of asset two
# x: the price series of asset one
# y: the price series of asset two
# OUTPUT:
# df1['cum rets']: cumulative returns in pandas data frame
# sharpe: Sharpe ratio
# CAGR: Compound Annual Growth Rate
x = df[s1]
y = df[s2]
# Run regression (including Kalman Filter) to find hedge ratio and then create spread series
df1 = pd.DataFrame({'y':y,'x':x})
df1.index = pd.to_datetime(df1.index)
state_means = KalmanFilterRegression(KalmanFilterAverage(x),KalmanFilterAverage(y))
df1['hr'] = - state_means[:,0]
df1['spread'] = df1.y + (df1.x * df1.hr)
# calculate half life
halflife = half_life(df1['spread'])
# calculate z-score with window = half life period
meanSpread = df1.spread.rolling(window=halflife).mean()
stdSpread = df1.spread.rolling(window=halflife).std()
df1['zScore'] = (df1.spread-meanSpread)/stdSpread
##############################################################
# trading logic
entryZscore = 1.25
exitZscore = -0.08
#set up num units long
df1['long entry'] = ((df1.zScore < - entryZscore) & ( df1.zScore.shift(1) > - entryZscore))
df1['long exit'] = ((df1.zScore > - exitZscore) & (df1.zScore.shift(1) < - exitZscore))
df1['num units long'] = np.nan
df1.loc[df1['long entry'],'num units long'] = 1
df1.loc[df1['long exit'],'num units long'] = 0
df1['num units long'][0] = 0
df1['num units long'] = df1['num units long'].fillna(method='pad')
#set up num units short
df1['short entry'] = ((df1.zScore > entryZscore) & ( df1.zScore.shift(1) < entryZscore))
df1['short exit'] = ((df1.zScore < exitZscore) & (df1.zScore.shift(1) > exitZscore))
df1.loc[df1['short entry'],'num units short'] = -1
df1.loc[df1['short exit'],'num units short'] = 0
df1['num units short'][0] = 0
df1['num units short'] = df1['num units short'].fillna(method='pad')
#set up totals: num units and returns
df1['numUnits'] = df1['num units long'] + df1['num units short']
df1['spread pct ch'] = (df1['spread'] - df1['spread'].shift(1)) / ((df1['x'] * abs(df1['hr'])) + df1['y'])
df1['port rets'] = df1['spread pct ch'] * df1['numUnits'].shift(1)
df1['cum rets'] = df1['port rets'].cumsum()
df1['cum rets'] = df1['cum rets'] + 1
##############################################################
try:
sharpe = ((df1['port rets'].mean() / df1['port rets'].std()) * sqrt(252))
except ZeroDivisionError:
sharpe = 0.0
##############################################################
start_val = 1
end_val = df1['cum rets'].iat[-1]
start_date = df1.iloc[0].name
end_date = df1.iloc[-1].name
days = (end_date - start_date).days
CAGR = (end_val / start_val) ** (252.0/days) - 1
df1[s1+ " "+s2+'_cum_rets'] = df1['cum rets']
return df1[s1+ " "+s2+'_cum_rets'], sharpe, CAGR
So now let’s run our full list of pairs through our Backtest function, and print out some results along the way, and finally after storing the equity curve for each pair, produce a chart that plots out each curve.
results = []
for pair in pairs:
rets, sharpe, CAGR = backtest(df[split:],pair[0],pair[1])
results.append(rets)
print("The pair {} and {} produced a Sharpe Ratio of {} and a CAGR of {}".format(pair[0],pair[1],
round(sharpe,2),
round(CAGR,4)))
rets0 = pd.concat(results, axis=1)
rets0.plot(figsize=(12,6),legend=True)
plt.legend(bbox_to_anchor=(1.01, 1.1), loc='upper left', ncol=1)
plt.grid(b=None, which=u'major', axis=u'both')
plt.title('Pairs Returns')
plt.xlabel('Date')
plt.ylabel('Returns');
Now we run a few extra lines of code to combine, equally weight, and print our our final equity curve:
filename = 'pairs_rets.csv'
rets0.to_csv(filename)
#concatenate together the individual equity curves into a single DataFrame
results_df = pd.concat(results,axis=1).dropna()
#equally weight each equity curve by dividing each by the number of pairs held in the DataFrame
results_df /= len(results_df.columns)
#sum up the equally weighted equity curves to get our final equity curve
final_res = results_df.sum(axis=1)
# square root of sample size for correct number of bins for returns distribution
print('Bin Count =', np.sqrt(len(final_res)))
Pair_Rets = ffn.to_returns(final_res)
Pair_Rets = pd.DataFrame(Pair_Rets)
Pair_Rets = Pair_Rets.fillna(0)
Pair_Rets.columns = ['Pairs_Returns']
fig, (ax1, ax2) = plt.subplots(1,2,figsize=(12, 5))
sns.distplot(Pair_Rets, hist = True, kde = True, bins=35,
hist_kws = {'linewidth': 1, 'alpha':.5},
label='Pairs Returns', color='#4b91fa', ax=ax1)
ax1.axvline(x=0.000709, color='#ff0000', linewidth=1.25, linestyle='dashed',label = 'Returns Mean')
ax1.set_title('Pairs Returns Distribution')
ax1.margins(0.001)
ax1.set_xlabel('Returns (%)')
ax1.set_ylabel('Density')
stats.probplot(Pair_Rets.Pairs_Returns, plot=ax2)
plt.tight_layout()
plt.show();
Pair_Rets.Pairs_Returns.describe()
perf = final_res.calc_stats()
num_pairs = len(results_df.columns)
print('Number of Pairs =', num_pairs)
# set SPY as benchmark
bench = df.loc[str(Pair_Rets.index[0]):str(Pair_Rets.index[-1])].SPY.pct_change().dropna()
Pair_Rets0 = Pair_Rets.loc[str(bench.index[0]):str(bench.index[-1])]
Total Returns Performance¶
fig = pf.create_returns_tear_sheet(Pair_Rets.Pairs_Returns, benchmark_rets=bench)
plt.figure(figsize=(14, 7))
pf.plot_perf_stats(factor_returns=Pair_Rets.Pairs_Returns, returns=bench)
plt.show();
The Sharpe Ratio is very good. Also take into consideration trading cost was not factored in these performance stats. The testing period was only from August 2015 until March 2020 only looking at specific pairs returns when traded in the training period using 1.25 - (-0.08) zscore entry/exit logic. Then creating an equal weight portfolio of those high pairs to trade in the test period. Cummulative returns were almost 58% for about 4yrs with limited drawdown.
qs.extend_pandas()
stock = Pair_Rets.Pairs_Returns
#qs.reports.html(stock, "SPY", title='Kalman Filter Pairs Strategy')