本文主要基於 Auquan 的 Pairs Trading Strategy 轉載和整理（非原創），主要修改了兩個地方：
- 用 A 股（主要是銀行股）替換了美股。（雖然 A 股事實上無法做空，但是用於研究是沒有問題的。）
- 增加了協整的 ADF 檢驗，旨在幫助讀者更深刻理解協整的定義。
如此，構成比較完整的配對交易的入門資料，供大家閱讀。本文的所有代碼已在 JoinQuant 上運行過。
配對交易是一個基於數學分析的策略的良好示例。
原理如下：
假設你有一對證券 X 和 Y，它們之間有某種潛在的經濟聯繫。舉例來說，可能是兩家公司生產相同的產品，例如百事可樂和可口可樂。你預期這兩者之間的價差（價格的比率或差異）會隨時間保持不變。然而，時不時地，這對的價差可能會出現偏離。這對中的偏離可能是由於暫時的供需變化、對某一證券的大額買入 / 賣出訂單、對某一公司的重要新聞的反應等等。當兩個證券之間出現暫時的偏離，即一隻股票上漲而另一隻下跌時，配對交易的策略是做空表現較好的股票，並做多表現較差的股票，押注於這兩者之間的 “價差” 最終會收斂。
配對交易是一種市場中立的交易策略，使交易者能夠從幾乎任何市場條件中獲利：上漲趨勢、下跌趨勢或橫盤運動。
我們將從構建一個人造示例開始。
1. 解釋概念 —— 首先生成兩個虛擬證券#
我們通過從正態分佈中抽樣來建模 X 的每日回報。然後我們進行累積和以獲得 X 在每一天的值。
import numpy as np
import pandas as pd
import statsmodels
from statsmodels.tsa.stattools import coint
import matplotlib.pyplot as plt
# just set the seed for the random number generator
np.random.seed(107)
# Generate the daily returns
Xreturns = np.random.normal(0, 1, 100)
# sum them and shift all the prices up
X = pd.Series(np.cumsum(Xreturns), name='X') + 50
X.plot(figsize=(15,7))
plt.show()
現在我們生成 Y。Y 應該與 X 有深厚的經濟聯繫，因此 Y 的價格應該變化得相當相似。我們通過取 X，將其上移並添加一些從正態分佈中抽取的隨機噪聲來建模。
noise = np.random.normal(0, 1, 100)
Y = X + 5 + noise
Y.name = 'Y'
pd.concat([X, Y], axis=1).plot(figsize=(15,7))
plt.show()
2. 協整#
協整，簡單來說，是一種 “不同” 的相關性。如果兩個序列是協整的，則它們之間的比率將圍繞一個均值變化。為了使配對交易在兩個時間序列之間有效，隨時間的比率的期望值必須收斂到均值，即它們應該是協整的。我們上面構建的時間序列是協整的。
協整性，區別於相關性，需要滿足兩個條件：
1）股票 X 和股票 Y 的價格對數序列是一階單整序列。
如果股票 X 的對數價格是非平穩時間序列（毫無疑問，這是常識）& 股票 X 的簡單單期收益率（對數價格相減結果可看做簡單單期收益率）序列是平穩時間序列，則股票 X 的對數價格是一階單整序列。
2）股票 X 和股票 Y 構造的回歸方程，其殘差序列是平穩的。
滿足以上兩個條件，股票 X 和股票 Y 可以說具有協整性。
高相關性並不一定意味著協整性。即使兩個股票的相關性差不多，但協整關係的概率差別比較大。
(Y/X).plot(figsize=(15,7))
plt.axhline((Y/X).mean(), color='red', linestyle='--')
plt.xlabel('Time')
plt.legend(['Price Ratio', 'Mean'])
plt.show()
3. 測試協整性#
方法 1：按照定義去判斷是否具有協整性。#
首先，驗證股票 X 和股票 Y 的對數價格是否非平穩序列。
平穩性檢驗，我們用 ADF 檢驗：
- 原假設（Null Hypothesis）：有單位根（unit root)，即非平穩
- 备择假设（Alternative Hypothesis）：無單位根，即平穩
ADF（y,lags,trend,max_lags,method）函數說明：
（參考《量化投資以 Python 為工具》P338）
- 參數 lags 的值確定檢驗模型。lags 的值表示一階差分的滯後階數。
- 檢驗結果的 Test Statistic 值和 CriticalValues 值進行比較，落在哪個區間。小於顯著性水平（1% 或 5%）的話，則拒絕原假設。（或者只看 p-value 即可，p-value 足夠小，就拒絕原假設。）
stock_list = ['600015.XSHG', '601169.XSHG']
sh = get_price(stock_list, start_date="2019-01-01", end_date="2020-04-30", frequency="daily", fields=['close'])['close']
PAf=sh['600015.XSHG']
PBf=sh['601169.XSHG']
PAflog=np.log(PAf)
PBflog=np.log(PBf)
from arch.unitroot import ADF
adfA=ADF(PAflog)
print(adfA.summary().as_text())
Augmented Dickey-Fuller Results#
Test Statistic -0.972
P-value 0.763
Lags 0#
Trend: Constant
Critical Values: -3.45 (1%), -2.87 (5%), -2.57 (10%)
Null Hypothesis: The process contains a unit root.
Alternative Hypothesis: The process is weakly stationary.
從上述結果來看，Test Statistic=-0.972，大於原假設分佈下的 1、5、10 分位數，從而不能拒絕原假設，進而說明 600015（華夏銀行）的對數價格序列是非平穩的。同理，601169 (北京銀行) 的對數價格序列也是非平穩的（驗證略）。
第二步，驗證股票 X 和股票 Y 的簡單單期收益率序列是否平穩。
retA=PAflog.diff()[1:]
adfretA=ADF(retA)
print(adfretA.summary().as_text())
Augmented Dickey-Fuller Results#
Test Statistic -17.403
P-value 0.000
Lags 0#
Trend: Constant
Critical Values: -3.45 (1%), -2.87 (5%), -2.57 (10%)
Null Hypothesis: The process contains a unit root.
Alternative Hypothesis: The process is weakly stationary.
從上述結果來看，600015（華夏銀行）的對數價格的差分（即簡單單期收益率）序列是平穩的。同理，601169 (北京銀行) 的簡單單期收益率序列也是平穩的（驗證略）。
最後，驗證股票 X 和股票 Y 構造的回歸方程，其殘差序列是平穩的。
import statsmodels.api as sm
model = sm.OLS(PBflog, sm.add_constant(PAflog))
results = model.fit()
print(results.summary())
#回歸截距項
alpha=results.params[0]
#提取回歸系數
beta=results.params[1]
#殘差
spread=PBflog-beta*PAflog-alpha
#繪製殘差序列時序圖
spread.plot()
plt.title (' 價差序列 ')
#對殘差做平穩性檢驗
adfSpread = ADF(spread,trend='nc') #殘差均值為0，所以trend設為nc
print(adfSpread.summary().as_text())
Augmented Dickey-Fuller Results
=====================================
Test Statistic -3.555
P-value 0.000
Lags 3
-------------------------------------
Trend: No Trend
Critical Values: -2.57 (1%), -1.94 (5%), -1.62 (10%)
Null Hypothesis: The process contains a unit root.
Alternative Hypothesis: The process is weakly stationary.
從上述結果可知，殘差是平穩的。
至此為止，我們可以說股票X和股票Y具有協整關係。
##### 方法2：有一個方便的測試在statsmodels.tsa.stattools中。我們應該看到一個非常低的p-value，因為我們人為地創建了兩個系列，使其協整的可能性最大。
```python
# compute the p-value of the cointegration test
# will inform us as to whether the ratio between the 2 timeseries is stationary
# around its mean
score, pvalue, _ = coint(X,Y)
看 p-value 即可，p-value 足夠小，即可認為股票 X 和股票 Y 具有協整關係。（是不是覺得方法 2 比方法 1 簡單多了～～不過，方法 1 更幫助理解協整的概念）
4. 如何進行配對交易？#
因為兩個協整的時間序列（例如上面的 600015.XSHG 和 601169.XSHG）會向彼此漂移並分開，因此會有時候價差很高，有時候價差很低。我們通過買入一種證券並賣出另一種證券來進行配對交易。這樣，如果兩個證券一起下跌或一起上漲，我們既不賺也不賠 —— 我們是市場中立的。
回到上面的 X 和 Y，遵循 Y = ⍺ X + e，使得比率（Y/X）圍繞其均值⍺波動，我們在兩者回歸均值時賺取利潤。為了做到這一點，我們將觀察 X 和 Y 何時相距較遠，即⍺過高或過低：
- 做多比率 當比率⍺低於通常水平時，我們預期它會上升。在上面的例子中，我們通過買入 Y 並賣出 X 來押注。
- 做空比率 當比率⍺較大時，我們預期它會變小。在上面的例子中，我們通過賣出 Y 並買入 X 來押注。
注意 我們始終保持 “對沖頭寸”：做空頭寸如果賣出的證券失去價值則獲利，做多頭寸如果證券增值則獲利，因此我們對整體市場運動免疫。我們只有在證券 X 和 Y 相對於彼此移動時才會賺或賠。
5. 使用數據尋找類似行為的證券#
最好的方法是從你懷疑可能協整的證券開始，並進行統計檢驗。如果你只是對所有配對進行統計檢驗，你將會受到多重比較偏差的影響。
多重比較偏差 只是指在進行多次檢驗時，錯誤生成顯著 p-value 的機會增加，因為我們進行了很多檢驗。如果在隨機數據上進行 100 次檢驗，我們應該預期看到 5 個 p-value 低於 0.05。如果你在比較 n 個證券的協整性，你將進行 n (n-1)/2 次比較，並且應該預期看到許多錯誤顯著的 p-value，這將隨著你增加而增加。為了避免這種情況，選擇少量你有理由懷疑可能協整的配對，並單獨測試每一對。這將減少多重比較偏差的暴露。
所以讓我們試著找到一些顯示協整的證券。讓我們使用一籃中國的銀行股，例如 “600000.XSHG”、“600015.XSHG”、“600016.XSHG”、“600036.XSHG”、“601009.XSHG” 等等。這些股票在相似的細分市場運作，可能具有協整的價格。我們掃描一個證券列表，並測試所有配對之間的協整性。它返回一個協整檢驗分數矩陣、一個 p-value 矩陣，以及任何 p-value 小於 0.05 的配對。這種方法容易受到多重比較偏差的影響，實際上這些證券應該經過第二次驗證步驟。為了這個例子，我們忽略這一點。
def find_cointegrated_pairs(data):
# data的長度
n = data.shape[1]
# 初始化
score_matrix = np.zeros((n, n))
pvalue_matrix = np.ones((n, n))
# 抽取列的名稱
keys = data.keys()
# 初始化強協整組
pairs = []
for i in range(n):
for j in range(i+1, n):
S1 = data[keys[i]]
S2 = data[keys[j]]
result = coint(S1, S2)
score = result[0]
pvalue = result[1]
score_matrix[i, j] = score
pvalue_matrix[i, j] = pvalue
if pvalue < 0.05:
pairs.append((keys[i], keys[j]))
return score_matrix, pvalue_matrix, pairs
stock_list = ["600000.XSHG", "600015.XSHG", "600016.XSHG", "600036.XSHG", "601009.XSHG","601166.XSHG", "601169.XSHG", "601328.XSHG", "601398.XSHG", "601988.XSHG", "601998.XSHG"]
data = get_price(stock_list, start_date="2019-01-01", end_date="2020-04-30", frequency="daily", fields=["close"])["close"]
# Heatmap to show the p-values of the cointegration test
# between each pair of stocks
scores, pvalues, pairs = find_cointegrated_pairs(data)
import seaborn
m = [0,0.2,0.4,0.6,0.8,1]
seaborn.heatmap(pvalues,
xticklabels=stock_list,
yticklabels=stock_list,
cmap='RdYlGn_r' ,
mask = (pvalues >= 0.98))
plt.show()
print(pairs)
[('600015.XSHG', '601169.XSHG'), ('600036.XSHG', '601328.XSHG'), ('601166.XSHG', '601398.XSHG')]
看起來 '600015.XSHG' 和 '601169.XSHG' 是協整的。讓我們看看價格，以確保沒有異常情況發生。
S1 = data['600015.XSHG']
S2 = data['601169.XSHG']
score, pvalue, _ = coint(S1, S2)
print(pvalue)
ratios = S1 / S2
ratios.plot(figsize=(15,7))
plt.axhline(ratios.mean())
plt.legend(['Price Ratio'])
plt.show()
這個比率看起來確實圍繞著穩定的均值波動。絕對比率在統計上並不是很有用。通過將其視為 z-score 來標準化我們的信號會更有幫助。Z 分數定義為：
Z Score (Value) = (Value — Mean) / Standard Deviation
警告
在實踐中，這通常是為了給數據提供某種尺度，但這假設了潛在的分佈。通常是正態的。然而，許多金融數據並不是正態分佈的，我們必須非常小心，不要僅僅假設正態性或任何特定的分佈來生成統計數據。比率的真實分佈可能是非常厚尾的，容易受到極端值的影響，從而破壞我們的模型並導致巨大的損失。
def zscore(series):
return (series - series.mean()) / np.std(series)
zscore(ratios).plot(figsize=(15,7))
plt.axhline(zscore(ratios).mean(), color='black')
plt.axhline(1.0, color='red', linestyle='--')
plt.axhline(-1.0, color='green', linestyle='--')
plt.legend(['Ratio z-score', 'Mean', '+1', '-1'])
plt.show()
6. 簡單策略：#
- 當 z-score 低於 - 1.0 時，進行 “做多” 比率
- 當 z-score 高於 1.0 時，進行 “做空” 比率
- 當 z-score 接近零時退出頭寸
這只是冰山一角，僅僅是一個非常簡單的示例來說明這些概念。
- 在實踐中，你會想計算一個更優化的權重，以決定持有多少 S1 和 S2 的股份
- 你還會想使用不斷更新的統計數據進行交易。
首先將數據分為 70% 的訓練集和 30% 的測試集。
ratios = data['600015.XSHG'] / data['601169.XSHG']
print(len(ratios))
le = int(len(ratios)*7/10)
train = ratios[:le]
test = ratios[le:]
一般來說，對整個樣本大小進行統計可能是有害的。例如，如果市場上漲，且兩個證券也隨之上漲，那麼你過去三年的平均價格可能無法代表今天的情況。因此，交易者通常使用依賴於最近數據的滾動窗口的統計數據。
而不是使用比率值，讓我們使用 5 日移動平均來計算 z 分數，並使用 60 日移動平均和 60 日標準差作為均值和標準差。
ratios_mavg5 = train.rolling(window=5,center=False).mean()
ratios_mavg60 = train.rolling(window=60,center=False).mean()
std_60 = train.rolling(window=60,center=False).std()
zscore_60_5 = (ratios_mavg5 - ratios_mavg60)/std_60
plt.figure(figsize=(15,7))
plt.plot(train.index, train.values)
plt.plot(ratios_mavg5.index, ratios_mavg5.values)
plt.plot(ratios_mavg60.index, ratios_mavg60.values)
plt.legend(['Ratio','5d Ratio MA', '60d Ratio MA'])
plt.ylabel('Ratio')
plt.show()
我們可以使用移動平均來計算每個給定時間的比率的 z 分數。這將告訴我們比率有多極端，以及此時進入頭寸是否明智。讓我們看看現在的 z 分數。
# Take a rolling 60 day standard deviation
std_60 = train.rolling(window=60,center=False).std()
std_60.name = 'std 60d'
# Compute the z score for each day
zscore_60_5 = (ratios_mavg5 - ratios_mavg60)/std_60
zscore_60_5.name = 'z-score'
plt.figure(figsize=(15,7))
zscore_60_5.plot()
plt.axhline(0, color='black')
plt.axhline(1.0, color='red', linestyle='--')
plt.axhline(-1.0, color='green', linestyle='--')
plt.legend(['Rolling Ratio z-Score', 'Mean', '+1', '-1'])
plt.show()
z 分數在沒有上下文的情況下意義不大，讓我們將其與價格一起繪製，以了解它的樣子。
# Plot the ratios and buy and sell signals from z score
plt.figure(figsize=(15,7))
train[60:].plot()
buy = train.copy()
sell = train.copy()
buy[zscore_60_5>-1] = 0
sell[zscore_60_5<1] = 0
buy[60:].plot(color='g', linestyle='None', marker='^')
sell[60:].plot(color='r', linestyle='None', marker='^')
x1,x2,y1,y2 = plt.axis()
plt.axis((x1,x2,ratios.min(),ratios.max()))
plt.legend(['Ratio', 'Buy Signal', 'Sell Signal'])
plt.show()
這對我們正在交易的實際股票意味著什麼？讓我們看看。
# Plot the prices and buy and sell signals from z score
plt.figure(figsize=(18,9))
S1 = data['600015.XSHG'].iloc[:le]
S2 = data['601169.XSHG'].iloc[:le]
S1[60:].plot(color='b')
S2[60:].plot(color='c')
buyR = 0*S1.copy()
sellR = 0*S1.copy()
# When buying the ratio, buy S1 and sell S2
buyR[buy!=0] = S1[buy!=0]
sellR[buy!=0] = S2[buy!=0]
# When selling the ratio, sell S1 and buy S2
buyR[sell!=0] = S2[sell!=0]
sellR[sell!=0] = S1[sell!=0]
buyR[60:].plot(color='g', linestyle='None', marker='^')
sellR[60:].plot(color='r', linestyle='None', marker='^')
x1,x2,y1,y2 = plt.axis()
plt.axis((x1,x2,min(S1.min(),S2.min())*0.95,max(S1.max(),S2.max())*1.05))
plt.legend(['華夏銀行','北京銀行', 'Buy Signal', 'Sell Signal'])
plt.show()
注意我們有時在短腿上賺錢，有時在長腿上賺錢，有時兩者都賺。
讓我們看看這個信號能產生什麼樣的利潤。我們寫一個簡單的回測器，當比率低時買入 1 個比率（買入 1 股 600015.XSHG 並賣出比率 x 601169.XSHG 的股票），當比率高時賣出 1 個比率（賣出 1 股 600015.XSHG 並買入比率 x 601169.XSHG 的股票），並計算這些交易的盈虧。
def trade(S1, S2, window1, window2):
# If window length is 0, algorithm doesn't make sense, so exit
if (window1 == 0) or (window2 == 0):
return 0
# Compute rolling mean and rolling standard deviation
ratios = S1/S2
ma1 = ratios.rolling(window=window1,center=False).mean()
ma2 = ratios.rolling(window=window2,center=False).mean()
std = ratios.rolling(window=window2,center=False).std()
zscore = (ma1 - ma2)/std
# Simulate trading
# Start with no money and no positions
money = 0
countS1 = 0
countS2 = 0
for i in range(len(ratios)):
# Sell short if the z-score is > 1
if zscore[i] > 1:
money += S1[i] - S2[i] * ratios[i]
countS1 -= 1
countS2 += ratios[i]
# Buy long if the z-score is < 1
elif zscore[i] < -1:
money -= S1[i] - S2[i] * ratios[i]
countS1 += 1
countS2 -= ratios[i]
# Clear positions if the z-score between -.5 and .5
elif abs(zscore[i]) < 0.5:
money += countS1*S1[i] + S2[i] * countS2
countS1 = 0
countS2 = 0
# print('Z-score: '+ str(zscore[i]), countS1, countS2, S1[i] , S2[i])
return money
trade(data['600015.XSHG'].iloc[:le], data['601169.XSHG'].iloc[:le], 5, 60)
這個策略似乎是有利可圖的！（result=12.716）現在我們可以通過改變移動平均窗口、改變買入 / 賣出和退出頭寸的閾值等進一步優化，並檢查驗證數據上的性能改進。我們還可以嘗試更複雜的模型，如邏輯回歸、SVM 等來進行 1/-1 的預測。
讓我們看看它在測試數據上的表現（result=1.487）
trade(data['600015.XSHG'].iloc[le:], data['601169.XSHG'].iloc[le:], 5, 60)
7. 避免過擬合#
過擬合是交易策略最危險的陷阱。在我們的模型中，我們使用了滾動參數估計，並可能希望優化窗口長度。我們可以簡單地遍歷所有可能的合理窗口長度，並根據我們的模型表現最佳的長度進行選擇。下面我們寫一個簡單的循環，根據訓練數據的盈虧來評分窗口長度，並找到最佳的長度，結果是 58。
# Find the window length 0-254
# that gives the highest returns using this strategy
length_scores = [trade(data['600015.XSHG'].iloc[:le],
data['601169.XSHG'].iloc[:le], 5, l)
for l in range(255)]
best_length = np.argmax(length_scores)
print ('Best window length:', best_length)
現在我們檢查我們的模型在測試數據上的表現，發現這個窗口長度與最佳的不同！這是因為我們的原始選擇顯然是過擬合到樣本數據。
# Find the returns for test data
# using what we think is the best window length
length_scores2 = [trade(data['600015.XSHG'].iloc[le:],
data['601169.XSHG'].iloc[le:],5, l)
for l in range(255)]
# 用上述計算得出的best_length來計算length_scores2
print (best_length, 'day window:', length_scores2[best_length])
# Find the best window length based on this dataset,
# and the returns using this window length
best_length2 = np.argmax(length_scores2)
print (best_length2, 'day window:', length_scores2[best_length2])
我們可以看到上述的第一個結果是基於之前計算的 58，得出 length_scores2=1.487，而實際最佳的結果是 best_length=52，對應的 length_scores2=2.09。
我們可以通過單獨繪製訓練和測試數據的盈虧隨窗口長度的變化來看到這一點。
plt.figure(figsize=(15,7))
plt.plot(length_scores)
plt.plot(length_scores2)
plt.xlabel('Window length')
plt.ylabel('Score')
plt.legend(['Training', 'Test'])
plt.show()
參考文獻：
- Auquan 的Pairs Trading Strategy
- 蔡立耑《量化投資 — 以 Python 為工具》