Same algo just start 9 month earlier.

Hey Justin,

Thanks for the share. I've noticed that there is a coint function in statsmodels.tsa.stattools. Is there a significant difference between the coint function and the ADF test? Any sense in using both?

I've attached a backtest below that attempts to find the pvalue of both tests for each pair, every day. Disclaimer: what I often think is happening in python is actually not.

@Jamie,
I haven't tried the coint function in stattools yet, though I imagine it's very similar. I just took a quick glimpse at the code, and it's effectively running a regression of the lagged version of the input timeseries versus the unlagged version which is quite similar to ADF. The difference may lie in how the critical values are computed.

The Engle-Granger test is also sometimes used to test for co-integration, but I haven't looked at that implementation yet.

Thumbs up (y)

Great Algo. Its amazing. Very Helpful

Hello Justin / All
Could you suggest how I can run this algo on multiple pairs, rather than just one pair?

Thanks!
Preston

Try making a pairs trading class that keeps track of all the bookkeeping for each given pair. See David's generalized Kalman filters pair trading algo for a great example of class based pairs trading

Thanks for sharing info

Hi all,
I cloned Justin's algo, however when I run a backtest, the performance remains at 0% for the entirety of the backtest window.
I made no changes to the original source code.

Any ideas why this would be occurring?

Thanks!

Adam,

You probably run algo in daily mode and it only work in minute mode.
Here is my latest backtest of original Justin's Lent algo started just 9 month earlier.

All,
It's worth noting that when I post backtests, code, and research notebooks, the intent is to illustrate a methodology, and provide some code templates to spur the creative thought process of the community and save folks some time by providing cut-and-paste code fragments that can be integrated into their own code. By no means am I posting something that has been fully vetted, and immediately investable in it's exact form, by any stretch of the imagination. I often bias for simpler, rather than overly complex, examples as well, so as to benefit a broader spectrum of readers.

Vladimir,
I see you've recognized that the backtest I posted above seems to fail pretty badly over a different timeframe. We see this a lot with strategies we look at, many of which are overfit to just the 2 year period in the contests we run. We try to work with the algo owner and provide advice as to why it may have broken down over the different timeframes. Perhaps you can extend your analysis to provide me some advice as to improving this strategy? Maybe you have some recommendations as to how to incorporate a regime switching model which is very likely to help a strategy such as this given the time frame it seems to fail (the financial/commodity futures crises that occurred in late 2008). Perhaps a stochastic volatility regime switching model might help significantly. If you have experience in this area I'm sure the community would find it a solid addition to incorporate into strategies such as these to make them more robust. I know I would.

Justin,
Why did you choose the pair USO and GLD? I guess a broader question is can you suggest a process for scanning through a basket of stocks and determining if there are tradable pairs? I'm assuming tests for cointegration would be one method such as ADF as you used. It would be nice if there could be an algo to run through a basket of stocks and auto determine which would make "good" pairs.

Just a thought.
Adam

Hey Adam,
I just chose USO/GLD in order to replicate this example that uses those same tickers, from this book: http://www.amazon.com/Quantitative-Trading-Build-Algorithmic-Business/dp/0470284889/

That book is a really good intro to stat arb pair trading (as well as his other books). All the code in the book is in Matlab, so my algo was an attempt to implement it in Python, in our backtester, and incorporate some of the other statistical techniques described throughout the book.

You are correct, that screening a bunch of potential pairs is a reasonable research idea, but you should be cognizant of simply datamining. You first want to determine a sensible economic basis for which the pairs of stocks should be tied (e.g. pairs of stocks in the same sector would be reasonable pairs of stocks to search across). Writing an algo in our backtester to accomplish this would be fairly straightforward: First you can use our Morningstar fundamentals database to grab all stocks in the Energy sector, perhaps even filtering down to stocks of companies of a certain band of marketcap (e.g. only mid-cap energy stocks), then in before_trading_starts(), you loop over each stock pair computing the ADF p-value (or other cointegration stat), keep all the stock pairs that meet your criteria, and then in handle_data() you just run the ones that meet the criteria through an algo similar to the one I shared to enter/exit the trades.

Myself or someone on our team here at Q can try to develop a template for this and share it.

As well you can look at this forum post that shows how to develop a single algo that trades a portfolio of multiple pairs:
It's the algo backtest in the first comment from David Edwards, here:
https://www.quantopian.com/posts/quantopian-lecture-series-this-time-youre-more-wrong?c=1

Justin,
I noticed in the blog section you have a notebook on using a Bayesian optimizer...would you know how i can pull it into Q? its currently on github..thanks!

@Adam, At present it's not possible to use the Bayesian optimizer from the blog post in the Q environment. It was more of a proof of concept implementation idea. As you mentioned, the code I used for the blog post is on github and you can sign up for a trial with SigOpt to get a username/API key to work with it in your own python/zipline environment locally. Offering some of these alternative methods of optimization as a service is an interesting concept which we will have to think about as we develop our Q platform in the future. Thanks for the feedback!

Thanks Justin! Would be neat to be able to do that type of optimization and / or a particle swarm technique in Q. :)

Hello Justin,
I believe I found a gap in the trading logic. In the statistic filtering section (lines ~155-176) the algorithm immediately exits if a test fails. That prevents new trades from being opened but does nothing to handle existing trades. Open trades stay open until all of the statistical tests pass again and the algorithm reaches its standard exit logic.
By design we should also have a high likelihood of being in a trade when this happens so the impact could be quite high. The problem in detecting this is that if the relationship re-establishes quickly the performance won't suffer. But if we include a time period in which the relationship doesn't return quickly, as Vladimir did, the results are noticeable.
I added a few lines to close any positions that are open when the statistical tests break down. There are probably better ways of handling the exit logic, but this simple change shows the benefit of having it there. The algorithm doesn't do as well during the original test period but the performance improves over the extended period.

(I also made minor change on lines 20 and 21 to use sid() function to set x and y assets rather than symbol(). The rest of the algorithm is unchanged.)

Pair trading using Copula methods instead of cointegration is the new rage. Anyone tried it?

Pair trading using Copula methods instead of cointegration is the new
rage. Anyone tried it?

This paper offers a systematic comparison of copula methods and cointegration methods when applied to US goldmine stocks. Also, the paper contrasts the pair selection criteria based on the ADF statistic, Kendall's tau, Spearman's rho and distance metric... One note: I'm not the author.

Thanks Julian. I had a go at it and results are looking very good.

Anything you can share Aqua, to play with?

I put in lot of resources (time and money) to get copulas working on Q. But you can use this for a start:

https://www.quantconnect.com/forum/discussion/2011/pairs-trading---copula-method-vs-cointegration

Thanks, I was more looking for Q code to play with ....the spirit of sharing ;)

I put a lot of time (and money) in zipline-live and the algos's I developed and I still share ... One day karma will come and pay graciously

Peter, I will try to put up something and post it without disclosing my secret sauce :)

Hi Peter,

Here is the Gaussian Copula class. You can use it to trade a pair or a basket. Let me know if you have any questions:

def cov2corr( A ):  
    d = np.sqrt(A.diagonal())  
    A = ((A.T/d).T)/d  
    return A

from sklearn.svm import SVR  
from scipy.stats import mvn  
#from statsmodels.sandbox.distributions.extras import mvnormcdf

# avioding boundary trouble  
loc_eps = 0.0000001


# ================================================================  
# copula classes  
# ================================================================


class GaussianCopulaConditional:  
    """ """

    def __init__(self, corr_matrix=np.eye(2), conditionals=None):  
        """  
        Gaussian covariance matrix and variables variables list;  
        target variables will be kept in 'targets' list  
        """  
        self.cr = pd.DataFrame(corr_matrix)  
        cls = self.cr.columns.values

        # by default conditionals contains everything except the last component  
        if conditionals is None:  
            conditionals = cls[:-1]

        # the rest of components are targets  
        targets = [s for s in cls if s not in conditionals]

        # keep the variables lists  
        self.targets = targets  
        self.conditionals = conditionals

        # covariance submatrices  
        sII = np.array(self.cr.loc[conditionals, conditionals])  
        sIJ = np.array(self.cr.loc[conditionals, targets])  
        sJI = np.array(self.cr.loc[targets, conditionals])  
        sJJ = np.array(self.cr.loc[targets, targets])

        # conditional mean and covariance  
        pr = np.linalg.solve(sII, sIJ)  
        self.cond_cov = sJJ - sJI.dot(pr)  
        self.mupr = np.linalg.solve(sII, sIJ).T

    def make_input(self, u, u_cond):  
        """ """  
        if u is None:  
            u = 0.5 * np.ones((1, len(self.targets)))  
        if u_cond is None:  
            u_cond = 0.5 * np.ones((1, len(self.conditionals)))  
        u[u == 0] = loc_eps  
        u[u == 1] = 1 - loc_eps  
        return u, u_cond

    def pdf(self, u=None, u_cond=None):  
        """  
        u is 2-dim array n_obs x n_targets  
        u_cond is 2-dim array 1 x n_conditionals  
        """  
        u, u_cond = self.make_input(u, u_cond)  
        x_cond = ss.norm.ppf(u_cond)  
        x = ss.norm.ppf(u)  
        mn = self.mupr.dot(x_cond.T).T[0]  
        if len(self.targets) <= 1:  
            res = ss.norm.pdf(x, mn[0], self.cond_cov[0, 0])  
            fx = ss.norm.pdf(x)  
            res = pd.Series((res / fx)[:, 0], index=u[:, 0])  
        else:  
            res = ss.multivariate_normal.pdf(x, mean=mn, cov=self.cond_cov)  
            fx = ss.norm.pdf(x).prod(axis=1)  
            res = pd.Series(res / fx, index=range(x.shape[0]))  
        res.name = 'Cond PDF of ' + ', '.join(self.targets)  
        return res

    def cdf(self, u=None, u_cond=None):  
        """  
        u is 2-dim array n_obs x n_targets  
        u_cond is 2-dim array 1 x n_conditionals  
        """  
        u, u_cond = self.make_input(u, u_cond)  
        x_cond = ss.norm.ppf(u_cond)  
        x = ss.norm.ppf(u)  
        mn = self.mupr.dot(x_cond.T).T[0]  
        if len(self.targets) <= 1:  
            # univariate normal  
            res = ss.norm.cdf(x, mn[0], self.cond_cov[0, 0])[0,0]  
        else:  
            # mvnormcdf does not accept multiple input points  
            res = np.zeros(x.shape[0])  
            for i in range(x.shape[0]):  
                res[i] = mvnormcdf(x[i, :], mn, self.cond_cov)  
        return res


class GaussianCopula:  
    """ """

    def __init__(self, corr_matr=pd.DataFrame(np.eye(2))):  
        """  
        initialize copula (as two dimensional indipendent one by default)  
        """  
        self.set_corr(corr_matr)

    def set_corr(self, cr):  
        """ set correlation matrix parameter and compute its inverse """  
        self.cr = pd.DataFrame(cr)  
        self.dim = self.cr.shape[0]  
        cls = self.cr.columns  
        self.cri = pd.DataFrame(np.linalg.inv(self.cr), columns=cls, index=cls)

    def fit(self, data):  
        """  
        fit the copula from data  
        the data should be 'm times n' Numpy array or Pandas DataFrame,  
        m = number of samples (observations), n = number of components  
        """  
        cr = pd.DataFrame(data).corr(method='spearman')  
        self.set_corr(pearson_from_spearman(cr))

    def rvs(self, size=100000):  
        """  
        simulate sample for Monte Carlo method  
        size = sample volume, 100,000 by default  
        """  
        sam = ss.multivariate_normal.rvs(cov=self.cr, size=size)  
        return pd.DataFrame(ss.norm.cdf(sam), columns=self.cr.columns)

    def pdf(self, u):  
        """ probability density function of the copula """  
        x = ss.norm.ppf(u)  
        mat = self.cri - np.eye(self.dim)  
        pok = -0.5 * x.dot(mat).dot(x)  
        mult = 1 / np.sqrt(np.linalg.det(self.cri))  
        return mult * np.exp(pok)

    def probability(self, event_func, n_obs=100000):  
        """  
        calculates probability of an event A, which is defined via its  
        indicator function "event_func"  
        """  
        u = self.rvs(size=n_obs)  
        return event_func(u).mean()

    def conditional_probability(self, event_func, cond_func, n_obs=100000):  
        """  
        calculates conditional probability of an event A, given event B,  
        which are defined via there indicator functions "event_func" and  
        "cond_func", respectively  
        """  
        u = self.rvs(size=n_obs)  
        prob_of_cond = cond_func(u).mean()  
        prob_of_event_and_cond = (event_func(u) & cond_func(u)).mean()  
        if prob_of_cond > 0:  
            res = prob_of_event_and_cond / prob_of_cond  
        else:  
            res = 0  
        return res

    def conditional_mean_variance_event(self, y_func, cond_func, n_obs=100000):  
        """  
        calculates conditional mean and variance of a random variable Y = f(U),  
        given condition A, with function f defined by "y_func", and  
        condition A defined by its indicator function "cond_func"  
        """  
        u = self.rvs(size=n_obs)  
        cond_ser = cond_func(u).astype('float')  
        prob_of_cond = cond_ser.mean()  
        event_and_cond = (y_func(u) * cond_ser).mean()  
        event2_and_cond = (y_func(u) ** 2 * cond_ser).mean()  
        if prob_of_cond > 0:  
            me = event_and_cond / prob_of_cond  
            va = event2_and_cond / prob_of_cond - me ** 2  
        else:  
            me, va = 0, 0  
        return me, va

    def conditional_mean_variance_vector(self, y_func, x_func, n_obs=5000):  
        """  
        creates an SVR model for computing expected mean of Y given X,  
        where X, Y are functions of U;  
        returns the model 'mod';  
        to predict value of Y in a new point X, call mod.predict(X), as usual  
        """  
        u = self.rvs(size=n_obs)

        # expected value portion  
        y = y_func(u)  
        x = x_func(u)  
        mod_ev = SVR()  
        mod_ev.fit(x, y)  
        ym = mod_ev.predict(x)

        # variance portion  
        yv = (y - ym) ** 2  
        mod_va = SVR()  
        mod_va.fit(x, yv)

        return mod_ev, mod_va

# =================================================================  
# conversion from Pearson correlations to Spearman correlations  
# and vice versa  
# =================================================================


def spearman_from_pearson(r):  
    """ convertion for normal multivariate """  
    s = 6 / np.pi * np.arcsin(0.5 * r)  
    return s


def pearson_from_spearman(s):  
    """ conversion for normal multivariate """  
    r = 2 * np.sin(s * np.pi / 6)  
    return r


# =====================================================================  
# Gaussian multivariate normal cdf  
# =====================================================================

informcode = {0: 'normal completion with ERROR < EPS',  
              1: '''completion with ERROR > EPS and MAXPTS function values used;  
                    increase MAXPTS to decrease ERROR;''',  
              2: 'N > 500 or N < 1'}

def mvstdnormcdf(lower, upper, corrcoef, **kwds):  
    n = len(lower)  
    lower = np.array(lower)  
    upper = np.array(upper)  
    corrcoef = np.array(corrcoef)

    correl = np.zeros(int(n*(n-1)/2.0))  #dtype necessary?

    if (lower.ndim != 1) or (upper.ndim != 1):  
        raise ValueError('can handle only 1D bounds')  
    if len(upper) != n:  
        raise ValueError('bounds have different lengths')  
    if n==2 and corrcoef.size==1:  
        correl = corrcoef  
    elif corrcoef.ndim == 1 and len(corrcoef) == n*(n-1)/2.0:  
        correl = corrcoef  
    elif corrcoef.shape == (n,n):  
        correl = corrcoef[np.tril_indices(n, -1)]  
    else:  
        raise ValueError('corrcoef has incorrect dimension')

    if not 'maxpts' in kwds:  
        if n >2:  
            kwds['maxpts'] = 10000*n

    lowinf = np.isneginf(lower)  
    uppinf = np.isposinf(upper)  
    infin = 2.0*np.ones(n)

    np.putmask(infin,lowinf,0)# infin.putmask(0,lowinf)  
    np.putmask(infin,uppinf,1) #infin.putmask(1,uppinf)  
    #this has to be last  
    np.putmask(infin,lowinf*uppinf,-1)

    error, cdfvalue, inform = sp.stats.mvn.mvndst(lower,upper,infin,correl,**kwds)  
    if inform:  
        print('something wrong', informcode[inform], error)  
    return cdfvalue

def mvnormcdf(upper, mu, cov, lower=None,  **kwds):  
    upper = np.array(upper)  
    if lower is None:  
        lower = -np.ones(upper.shape) * np.inf  
    else:  
        lower = np.array(lower)  
    cov = np.array(cov)  
    stdev = np.sqrt(np.diag(cov)) # standard deviation vector  
    #do I need to make sure stdev is float and not int?  
    #is this correct to normalize to corr?  
    lower = (lower - mu)/stdev  
    upper = (upper - mu)/stdev  
    divrow = np.atleast_2d(stdev)  
    corr = cov/divrow/divrow.T  
    #v/np.sqrt(np.atleast_2d(np.diag(covv)))/np.sqrt(np.atleast_2d(np.diag(covv))).T

    return mvstdnormcdf(lower, upper, corr, **kwds)  

Attached notebook on usage.

Thanks Aqua Rooster, really interesting information, I came to it while reading "Algoritmic trading" by Ernie Chan. How much improvement do you see in copula compared with the Cointegration CADF and Johansen referred by Mr Chan?

@Jose, Not much improvement because the underlying principles of pair trading are the same. Only advantage is that you can trade baskets instead of pairs.

So how are the baskets identified? Or is that what the Gaussian copula thingy does?

@Grant, the copula gives you the conditional expectation of a stock's return given the return of other stocks in the basket. You can base on this decision the weights to use in the basket.