Strategy Library

Pairs Trading-Copula vs Cointegration

Abstract

In this research, We investigate two pairs trading methods and compare the result. Pairs trading involves in investigating the dependence structure between two highly correlated assets. With the assumption that mean reversion will occur, long or short positions are entered in the opposite direction when there is a price divergence. Typically the asset price distribution is modeled by a  Gaussian distribution of return series but the joint normal distribution may fail to catch some key features of the dependence of stock pairs' price like tail dependence. We investigate using copula theory to identify these trading opportunities.

We will discuss the basic framework of copula from the mathematical perspective and explain how to apply the approach in pairs trading. The implementation of the algorithm is based on based on the paper Trading strategies with copulas from Stander Y, Marais D, Botha I(2013). We compare the performance of the copula pairs trading strategy with the co-integration pairs trading method based on the paper Statistical arbitrage trading strategies and high-frequency trading from Hanson T A, Hall J R. (2012). The co-integration technique assumes a co-integration relationship between paired equities to identify profitable trading opportunities. The empirical results suggest that the copula-based strategy is more profitable than the traditional pairs trading techniques.

Framework of Copula

1. Definition

Given a random vector \(X_1,X_2,...,X_p\), its marginal cumulative distribution functions (CDFs) are \(F_i(x) = P[X_i \leq x]\). By applying the probability integral transform to each component, the marginal distributions of \((U_1,U_2,...,U_p) = (F_1(X_1),F_2(X_2),...,F_p(X_p))\) are uniform (from Wikipedia). Then the copula of \(X_1,X_2,...,X_p\) is defined as the joint cumulative distribution function of \(U_1,U_2,...,U_p\), for which the marginal distribution of each variable U is uniform as  \(U(0,1)\).

\[C(u_1,u_2,...,u_p) = P[U_1\leq u_1,U_2\leq u_2,..., U_1\leq u_1]\]

Copulas function contains all the dependency characteristics of the marginal distributions and will better describe the linear and non-linear relationship between variables, using probability. They allow the marginal distributions to be modeled independently from each other, and no assumption on the joint behavior of the marginals is required.

2. Bivariate Copulas

Since this research focuses on bivariate copulas (for pairs trading we have 2 random variables) some probabilistic properties are specified. Let X and Y be two random variables with cumulative probability function \(F_1(X)\) and \(F_2(Y)\). \(U=F_1(X), V=F_2(Y)\) which are uniformly distributed.  Then the copula function is \(C(u,v)=P(U\leq u,V\leq v)\). Taking the partial derivative of the copula function over U and V would give the conditional distribution function as follows:

\[P(U\leq u\mid V= v)=\frac{\partial C(u,v)}{\partial v}\] \[P(V\leq v\mid U= u)=\frac{\partial C(u,v)}{\partial u}\]

3. Archimedean Copulas

There are many copula functions that enable us to describe dependence structures between variables, other than the Gaussian assumption. Here we will focus three of these; the Clayton, Gumbel and Frank copula formulas from the Archimedean class. Archimedean copulas are based on the Laplace transforms φ of univariate distribution functions. They are constructed by a particular generator function \(\phi\).

\[C(u,v)=\phi^{-1}( \phi(u),\phi(v) )\]

The probability density function is:

\[c(u,v)=\phi_{(2)}^{-1}(\phi(u)+\phi(v))\phi^{'}(u)\phi^{'}(v)\]

Where \(\phi_{(2)}^{-1}\) is the inverse of the second derivative of the generator function.

Copula Copula function C(u,v;θ)
Clayton Copula \[(u^{-\theta}+v^{-\theta}-1)^{-1/\theta}\]
Gumbel Copula \[exp(-[(-\ln u)^{\theta}+(-\ln v)^{\theta}]^{1/\theta})\]
Frank Copula\[-\theta^{-1}\ln\left[1+\frac{(exp(-\theta u)-1)(exp(-\theta v)-1)}{exp(-\theta)-1}\right]\]

Genest and MacKay proved that the relation between the copula generator function and Kendall rank correlation tau in the bivariate case can be given by:

\[\tau=1+4\int_{0}^{1} \frac{\partial \phi (v)}{\partial \phi^{'}(v)}dv\]

So we can easily estimate the parameter in Archimedean copulas if we know Kendall’s tau rank measure and the generator function. Please refer to step 3 to see the formulas.

Part I - Copula Method

ETFs have many different stock sectors and asset classes which provide us a wide range of pairs trading candidates. Our data set consists of daily data of the ETFs traded on the NASDAQ or the NYSE.

We use the first 3 years of data to choose the best fitting copula and asset pair ("training formation period"). Next, we use a period of 5 years from 2011 to 2017 ("the trading period"), to execute the strategy. During the trading period we use a rolling 12 month window of data to get the copula parameters ("rolling formation period").

Step 1: Selecting the Paired Stocks

The general method of pair selection is based on both fundamental and statistical analysis.

1) Assemble a list of potentially related pairs

Any random pairs could be correlated. It is possible that those variables are not causally related to each other, but because of a spurious relationship due to either coincidence or the presence of a certain third, unseen factor. Thus, it is important for us to start with a list of securities that have something in common. For this demonstration, we choose some of the most liquid ETFs traded on the Nasdaq or the NYSE.  The relationship for those potentially related pairs could be due to an index, sector or asset class overlap. e.g. QQQ and XLK are two ETFs which track the market leading indices.

2) Filter the trading pair with statistical correlation

To determine which stock pairs to include in the analysis, correlations between the pre-selected ETF pairs are analyzed. Below are three types of correlation measures we usually use in statistics:

Correlation Measurement Techniques
Pearson correlation\[r = \frac{\sum (x_i- \bar{x})(y_i- \bar{y})}{\sqrt{\sum (x_i- \bar{x})^2)\sum (y_i- \bar{y})^2)} }\]
Kendall rank correlation\[\tau=\frac{n_c-n_d}{\frac{1}{2}n(n-1)}\]
Spearman rank correlation\[\rho=1-\frac{6\sum d_i^2}{n(n^2-1)}\]
 \(n\) = number of value in each data set \(n_c\) = number of concordant \(n_d\) = number of discordant \(d_i\) = the difference between the ranks of corresponding values \(x_i\) and \(y_i\)

We can get these coefficients in Python using functions from the stats library in SciPy. The correlations have been calculated using daily log stock price returns during the training formation period. We found the 3 correlation techniques give the paired ETFs the same correlation coefficient ranking. The Pearson correlation assumes that both variables should be normally distributed. Thus here we use Kendall rank as the correlation measure and choose the pairs with the highest Kendall rank correlation to implement the pairs trading. We get the daily historical closing price of our ETFs pair by using the History function and converting the prices to a log return series. Let \(P_x\) and \(P_y\) denote the historical stock price series for stock x and stock y. The log returns for the ETFs pair are given by:

\[R_x = ln(\frac{P_{x,t}}{P_{x,t-1}}),   R_y = ln(\frac{P_{y,t}}{P_{y,t-1}})\]   t = 1,2,...,n where n is the number of price data
def _pair_selection(self):

    tick_syl =  [["QQQ","XME","TNA","FAS","XLF","EWC","QLD"],
                 ["XLK","EWG","TLT","FAZ","XLU","EWA","QID"]]
    logreturn={}
    for i in range(2):
        syl = [self.AddSecurity(SecurityType.Equity, x, Resolution.Daily).Symbol.Value for x in tick_syl[i]]
        history = self.History(syl, self.lookbackdays,Resolution.Daily)
        # generate the log return series of paired stocks
        close = history['close'].unstack(level=0)
        df_logreturn = (np.log(close) - np.log(close.shift(1))).dropna()
        for j in tick_syl[i]:
            logreturn[j] = df_logreturn[j]
    # estimate coefficients of different correlation measures
    tau_coef,pr_coef,sr_coef= [],[],[]
    for i in range(len(tick_syl[i])):
        tik_x, tik_y= logreturn[tick_syl[0][i]], logreturn[tick_syl[1][i]]
        tau_coef.append(kendalltau(tik_x, tik_y)[0])
        pr_coef.append(pearsonr(tik_x, tik_y)[0])
        sr_coef.append(spearmanr(tik_x, tik_y)[0])
    index_max = tau_coef.index(max(tau_coef))
    self.ticker = [tick_syl[0][index_max],tick_syl[1][index_max]]

Step 2: Estimating Marginal Distributions of log-return

In order to construct the copula, we need to transform the log-return series \(R_x\) and \(R_y\) to two uniformly distributed values u and v. This can be done by estimating the marginal distribution functions of \(R_x\) and \(R_y\) and plugging the return values into a distribution function. As we make no assumptions about the distribution of the two log-return series, here we use the empirical distribution function to approach the marginal distribution \(F_1(R_x)\) and \(F_2(R_y)\). The Python ECDF function from the statsmodel library gives us the Empirical CDF as a step function.

Step 3: Estimating Copula Parameters

As discussed above, we estimate the copula parameter theta by the relationship between the copula and the dependence measure Kendall’s tau, for each of the Archimedean copulas.

Copula Kendall's tau parameter θ
Clayton Copula\[\frac{\theta}{\theta +2}\]\[\theta=2\tau(1-\tau)^{-1}\]
Gumbel Copula\[1-\theta^{-1}\]\[\theta=(1-\tau)^{-1}\]
Frank Copula\[1+4[D_1(\theta)-1]/\theta\]\[arg min\left(\frac{\tau-1}{4}-\frac{D_1(\theta)-1}{\theta}\right)^2\]
\[D_1(\theta)=\frac{1}{\theta}\int_{0}^{\theta}\frac{t}{exp(t)-1}dt \]
def _parameter(self, family, tau):
    if  family == 'clayton':
        return 2*tau/(1-tau)
    elif family == 'frank':
        # debye = quad(integrand, sys.float_info.epsilon, theta)[0]/theta  is first order Debye function
    	# frank_fun is the squared difference
    	# Minimize the frank_fun would give the parameter theta for the frank copula
      integrand = lambda t: t/(np.exp(t)-1)
    	frank_fun = lambda theta: ((tau - 1)/4.0  - (quad(integrand, sys.float_info.epsilon, theta)[0]/theta - 1)/theta)**2
      return minimize(frank_fun, 4, method='BFGS', tol=1e-5).x
    elif family == 'gumbel':
        return 1/(1-tau)

Step 4: Selecting the Best Fitting Copula

Once we get the parameter estimation for the copula functions, we use the AIC criteria to select the copula that provides the best fit in algorithm initialization.

\[AIC=-2L(\theta)+2k\]

where \(L(\theta)=\sum_{t=1}^T\log c(u_t,v_t;\theta)\) is the log-likelihood function and k is the number of parameters, here k=1.

The density functions of each copula function are as follows:

Copula Density function c(u,v;θ)
Clayton Copula \[(\theta+1)(u^{-\theta}+v^{-\theta}-1)^{-2-1/\theta}u^{-\theta-1}v^{-\theta-1}\]
Gumbel Copula\[C(u,v;\theta)(uv)^{-1}A^{-2+2/\theta}[(\ln u)(\ln v)]^{\theta -1}[1+(\theta-1)A^{-1/\theta}]\]
Frank Copula \[\frac{-\theta(exp(-\theta)-1)(exp(-\theta(u+v)))}{((exp(-\theta u)-1)(exp(-\theta v)-1)+(exp(-\theta)-1))^2}\]
\[A=(-\ln u)^{\theta}+(-\ln v)^{\theta}\]
def _lpdf_copula(self, family, theta, u, v):
    ''' estimate the log probability density function of three kinds of Archimedean copulas '''
    if  family == 'clayton':
        pdf = (theta+1) * ((u**(-theta)+v**(-theta)-1)**(-2-1/theta)) * (u**(-theta-1)*v**(-theta-1))
    elif family == 'frank':
        num = -theta * (np.exp(-theta)-1) * (np.exp(-theta*(u+v)))
        denom = ((np.exp(-theta*u)-1) * (np.exp(-theta*v)-1) + (np.exp(-theta)-1))**2
        pdf = num/denom
    elif family == 'gumbel':
        A = (-np.log(u))**theta + (-np.log(v))**theta
        c = np.exp(-A**(1/theta))
        pdf = c * (u*v)**(-1) * (A**(-2+2/theta)) * ((np.log(u)*np.log(v))**(theta-1)) * (1+(theta-1)*A**(-1/theta))
    return np.log(pdf)

The copula that provides the best fit is the one that corresponds to the lowest value of AIC criterion. The chosen pair is "QQQ" & "XLK".

self.family = ['clayton', 'frank', 'gumbel']
tau = kendalltau(x, y)[0]  # estimate Kendall'rank correlation
AIC ={}  # generate a dict with key being the copula family, value = [theta, AIC]
for i in self.family:
    lpdf = [self._lpdf_copula(i, self._parameter(i,tau), x, y) for (x, y) in zip(u, v)]
    # Replace nan with zero and inf with finite numbers in lpdf list
    lpdf = np.nan_to_num(lpdf)
    loglikelihood = sum(lpdf)
    AIC[i] = [self._parameter(i,tau), -2*loglikelihood + 2]
    # choose the copula with the minimum AIC
    self.copula = min(AIC.items(), key = lambda x: x[1][1])[0]

Step 5: Generating the Trading Signals

The copula functions include all the information about the dependence structures of two return series. According to Stander Y, Marais D, Botha I. in Trading strategies with copulas, the fitted copula is used to derive the confidence bands for the conditional marginal distribution function of \(C(v\mid u)\) and \(C(u\mid v)\), that is the mispricing indexes. When the market observations fall outside the confidence band, it is an indication that pairs trading opportunity is available. Here we choose 95%  as the upper confidence band, 5% as the lower confidence band as indicated in the paper. The confidence level was selected based on a back-test analysis in the paper that shows using 95% seems to lead to appropriate trading opportunities to be identified.

Given current returns \(R_x, R_y\) of stock X and stock Y, we define the "mis-pricing indexes" are:

\[MI_{X|Y}=P(U\leq u\mid V\leq v)=\frac{\partial C(u,v)}{\partial v}\] \[MI_{Y|X}=P(V\leq v\mid U\leq u)=\frac{\partial C(u,v)}{\partial u}\]

For further mathematical proof, please refer to Xie W, Wu Y. Copula-based pairs trading strategy. The conditional probability formulas of bivariate copulas can be derived by taking partial derivatives of copula functions shown in Table 1. The results are as follows:

Gumbel Copula

\[C(v\mid u)=C(u,v;\theta)[(-\ln u)^\theta+(-\ln v)^\theta]^{\frac{1-\theta}{\theta}}(-\ln u)^{\theta-1}\frac{1}{u}\] \[C(u\mid v)=C(u,v;\theta)[(-\ln u)^\theta+(-\ln v)^\theta]^{\frac{1-\theta}{\theta}}(-\ln v)^{\theta-1}\frac{1}{v}\]

Clayton Copula

\[C(v\mid u)=u^{-\theta-1}(u^{-\theta}+v^{-\theta}-1)^{-\frac{1}{\theta}-1}\] \[C(u\mid v)=v^{-\theta-1}(u^{-\theta}+v^{-\theta}-1)^{-\frac{1}{\theta}-1}\]

Frank Copula

\[C(v\mid u)=\frac{(exp(-\theta u)-1)(exp(-\theta v)-1)+(exp(-\theta v)-1)}{(exp(-\theta u)-1)(exp(-\theta v)-1)+(exp(-\theta)-1)}  \] \[C(u\mid v)=\frac{(exp(-\theta u)-1)(exp(-\theta v)-1)+(exp(-\theta u)-1)}{(exp(-\theta u)-1)(exp(-\theta v)-1)+(exp(-\theta)-1)} \]

After selection of trading pairs and the best-fitted copulas, we take the following steps for trading. Please note we implement the Steps 1, 2, 3 and 4 on the first day of each month using the daily data for the last 12 months, which means our empirical distribution functions and copula parameters theta estimation are updated once a month. In summary each month:

  • During the 12 months' rolling formation period, daily close prices are used to calculate the daily log returns for the pair of ETFs and then compute Kendall's rank correlation.
  • Estimate the marginal distribution functions of log returns of X and Y, which are ecdf_x and ecdf_y separately.
  • Plug Kendall's tau into copula parameter estimation functions to get the value of theta.
  • Run linear regression over the two price series. The coefficient is used to determine how many shares of stock X and Y to buy and sell. For example, if the coefficient is 2, for every X share that is bought or sold, 2 units of Y are sold or bought.
def _set_signal(self):
    history = self.History(self.ticker, self.lookbackdays,Resolution.Daily)
    # generate the log return series of paired stocks
    close = history['close'].unstack(level=0)
    logreturn = (np.log(close) - np.log(close.shift(1))).dropna()
    x, y = logreturn[self.ticker[0]], logreturn[self.ticker[1]]
    # estimate Kendall'rank correlation each trading day
    tau = kendalltau(x, y)[0]
    # etstimate the copula parameter: theta
    self.theta = self._parameter(self.copula, tau)
    # simulate the empirical distribution function for returns of two paired stocks
    self.ecdf_x, self.ecdf_y  = ECDF(x), ECDF(y)
    # run linear regression over the two history return series
    self.coef = stats.linregress(x,y).slope

Finally during the trading period, each day we convert today's returns to u and v by using empirical distribution functions ecdf_x and ecdf_y. After that, two mispricing indexes are calculated every trading day by using the estimated copula C.  The algorithm constructs short positions in X and long positions in Y on the days that \(MI_{Y|X}<0.05\) and \(MI_{X|Y}>0.95\). It constructs short position in Y and long positions in X on the days that \(MI_{Y|X}>0.95\) and \(MI_{X|Y}<0.05\).

def OnData(self,data):
    for i in self.syl:
        self.price_list[i].append(self.Portfolio[i].Price)
    # compute today's log return of 2 stocks
    if len(self.price_list[self.syl[0]]) < 2 or len(self.price_list[self.syl[1]]) < 2: return
    else:
        return_x = np.log(float(self.price_list[self.syl[0]][-1]/self.price_list[self.syl[0]][-2]))
        return_y = np.log(float(self.price_list[self.syl[1]][-1]/self.price_list[self.syl[1]][-2]))
    # Convert the two returns to uniform values u and v using the empirical distribution functions
    u_value = self.ecdf_x(return_x)
    v_value = self.ecdf_y(return_y)
    # Compute the mispricing indices for u and v by using estimated copula
    self._misprice_index(self.copula, self.theta, u_value, v_value)
    quantity = self.CalculateOrderQuantity(self.syl[1],0.4)
    if self.MI_u_v < self.floor_CL and self.MI_v_u > self.cap_CL:
        if self.Portfolio[self.syl[0]].Quantity < 0 and self.Portfolio[self.syl[1]].Quantity > 0:
            self.Liquidate(self.syl[0])
            self.Liquidate(self.syl[1])
            quantity = self.CalculateOrderQuantity(self.syl[1],0.4)
            self.Sell(self.syl[1], 1 * quantity )
            self.Buy(self.syl[0], self.coef * quantity)
        else:
            self.Sell(self.syl[1], 1 * quantity )
            self.Buy(self.syl[0], self.coef * quantity)
    elif self.MI_u_v > self.cap_CL and self.MI_v_u < self.floor_CL:
        if self.Portfolio[self.syl[0]].Quantity > 0 and self.Portfolio[self.syl[1]].Quantity < 0:
            self.Liquidate(self.syl[0])
            self.Liquidate(self.syl[1])
            quantity = self.CalculateOrderQuantity(self.syl[1],0.4)
            self.Buy(self.syl[1], 1 * quantity )
            self.Sell(self.syl[0], self.coef * quantity)
        else:
            self.Buy(self.syl[1], 1 * quantity )
            self.Sell(self.syl[0], self.coef * quantity)

Part II - Cointegration Method

For the cointegration pairs trading method, we choose the same ETF pair "GLD" and "DGL".  There is no need to choose a copula function so there is only a 12 month rolling formation period. The trading period is 5 years from January 2011 to  May 2017.

Step 1: Generate the Spread Series

At the start of each month, we generate the log price series of two ETFs with the daily close. Then the spread series is estimated using regression analysis based on log price series data. For equities X and Y, we run linear regression over the log price series and get the coefficient β.

\[spread_t=\log(price_t^y)-\beta \log(price_t^x)\]

Step 2: Compute the Threshold

Using the standard deviation of spread during the rolling formation period, a threshold of two standard deviations is set up for the trading strategy as indicated in the paper.

price_x = pd.Series([float(i.Close) for i in self.symbols[0].hist_window],
                     index = [i.Time for i in self.symbols[0].hist_window])

price_y = pd.Series([float(i.Close) for i in self.symbols[1].hist_window],
                     index = [i.Time for i in self.symbols[1].hist_window])
if len(price_x) < 250: return
spread = self.regr(np.log(price_x), np.log(price_y))
mean = np.mean(spread)
std = np.std(spread)
ratio = floor(self.Portfolio[self.symbols[1]].Price / self.Portfolio[self.symbols[0]].Price)
if spread[-1] > mean + self.threshold * std:
    if not self.Portfolio[self.symbols[0]].Quantity > 0 and not self.Portfolio[self.symbols[0]].Quantity < 0:
        self.Sell(self.symbols[1], 100)
        self.Buy(self.symbols[0],  ratio * 100)

elif spread[-1] < mean - self.threshold * std:
    if not self.Portfolio[self.symbols[0]].Quantity < 0 and not self.Portfolio[self.symbols[0]].Quantity > 0:
        self.Sell(self.symbols[0], 100)
        self.Buy(self.symbols[1], ratio * 100)
else:
    self.Liquidate()

Step 3: Set up the Trading Signals

On each trading day, we enter a trade whenever the spread moves more than two standard deviations away from its mean. In other words, we construct short positions in X and long positions in Y on the day that spread mean+2*std. We construct short positions in Y and long positions in X on the day that spread<mean-2*std. The trade is exited if the spread reverts to its equilibrium (defined as less than half a standard deviation from zero spread). The value of mean and standard deviation are calculated from the rolling formation period and will be updated once a month.

Summary

Ultimately pairs trading intends to capture the price divergence of two correlated assets through mean reversion. Our results demonstrate that the copula approach for pairs trading is superior to the conventional cointegration method because it is based on the probability of the dependence structure, vs cointegration which relies on simple linear regression variance from normal pricing. We found through testing the performance of the copula method less sensitive to the starting parameters. Because the cointegration method relies on standard distribution and the ETF pairs had low volatility there were few trading opportunities.

method Transactions Profit Sharpe Ratio Drawdown
Copula 346 274.293% 1.022 19.4%
Cointegration 91 26.358% 0.298 23.7%

Generally, ETFs are not very volatile and so mean-reversion did not provide many trading opportunities. There are only 39 trades during 5 years for cointegration method. It is observed that the use of copula in pairs trading provides more trading opportunities as it does not require any rigid assumptions according to Liew R Q, Wu Y. - Pairs trading A copula approach.

Algorithm

Backtest for copula method

Backtest for cointegration method

References

  1. Stander Y, Marais D, Botha I. Trading strategies with copulas[J]. Journal of Economic and Financial Sciences, 2013, 6(1): 83-107. Online Copy
  2. Hanson T A, Hall J R. Statistical arbitrage trading strategies and high-frequency trading[J]. 2012. Online Copy
  3. Mahfoud M, Michael M. Bivariate Archimedean copulas: an application to two stock market indices[J]. BMI Paper, 2012.  Online Copy
  4. Rad H, Low R K Y, Faff R. The profitability of pairs trading strategies: distance, cointegration and copula methods[J]. Quantitative Finance, 2016, 16(10): 1541-1558.online copy
  5. Mahfoud M, Michael M. Bivariate Archimedean copulas: an application to two stock market indices[J]. BMI Paper, 2012.  Online Copy
  6. LANDGRAF N, SCHOLTUS K, DIRIS D R B. High-Frequency copula-based pairs trading on US Goldmine Stocks[J]. 2016.
  7. Genest, C. and MacKay, J., 1986, The Joy of Copulas: Bivariate Distributions with Uniform Marginals, The American Statistician, 40, 280-283
  8. Jean Folger, Pairs Trading Example Online Copy
  9. Xie W, Wu Y. Copula-based pairs trading strategy[C] Asian Finance Association (AsFA) 2013 Conference. doi. 2013, 10.
  10. Liew R Q, Wu Y. Pairs trading: A copula approach[J]. Journal of Derivatives & Hedge Funds, 2013, 19(1): 12-30.

You can also see our Documentation and Videos. You can also get in touch with us via Chat.

Did you find this page Helpful ?