この記事は主に Auquan のペアトレーディング戦略を転載および整理したもので（非オリジナル）、主に 2 つの点を修正しました：
- 米国株を A 株（主に銀行株）に置き換えました。（A 株は実際には空売りできませんが、研究には問題ありません。）
- 共整合の ADF 検定を追加し、読者が共整合の定義をより深く理解できるようにしました。
これにより、比較的完全なペアトレーディングの入門資料が構成され、皆さんに読んでいただけます。この記事のすべてのコードは JoinQuant で実行済みです。
ペアトレーディングは、数学的分析に基づく戦略の良い例です。
原則は次のとおりです：
X と Y という 2 つの証券のペアがあり、何らかの基礎的な経済的リンクがあります。例としては、同じ製品を製造する 2 つの企業、たとえばペプシとコカ・コーラが挙げられます。これら 2 つの間のスプレッド（価格の比率または差）が時間とともに一定であると期待します。しかし、時折、これら 2 つのペア間のスプレッドに乖離が生じることがあります。ペア内の乖離は、一時的な供給 / 需要の変化、1 つの証券に対する大きな買い / 売り注文、企業の重要なニュースへの反応などによって引き起こされる可能性があります。一時的に 2 つの証券間に乖離がある場合、つまり 1 つの株が上昇し、もう 1 つが下落する場合、ペアトレードは、パフォーマンスの良い株をショートし、パフォーマンスの悪い株をロングすることで、2 つの間の「スプレッド」が最終的に収束することに賭けることになります。
ペアトレーディングは、市場中立の取引戦略であり、トレーダーが実質的にどんな市場条件でも利益を得ることを可能にします：上昇トレンド、下降トレンド、または横ばいの動き。
まず、人工的な例を構築します。
1. コンセプトの説明 —2 つの偽の証券を生成することから始める#
X のデイリーリターンを正規分布から抽出してモデル化します。次に、累積和を計算して、各日の X の値を取得します。
import numpy as np
import pandas as pd
import statsmodels
from statsmodels.tsa.stattools import coint
import matplotlib.pyplot as plt
# ランダム数生成器のシードを設定
np.random.seed(107)
# デイリーリターンを生成
Xreturns = np.random.normal(0, 1, 100)
# 合計してすべての価格を上にシフト
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. 共整合#
共整合は、非常に大まかに言えば、「異なる」形の相関です。2 つの系列が共整合している場合、それらの比率は平均値の周りで変動します。ペアトレーディングが 2 つの時系列間で機能するためには、時間の経過とともに比率の期待値が平均に収束する必要があります。つまり、共整合している必要があります。上記で構築した時系列は共整合しています。
共整合性は、相関性とは異なり、2 つの条件を満たす必要があります：
1）株式 X と株式 Y の価格対数系列は一次単位根系列である。
もし株式 X の対数価格が非定常時系列であり（これは常識です）、株式 X の単純単期収益率（対数価格の差を単純単期収益率と見なす）系列が定常時系列であるならば、株式 X の対数価格は一次単位根系列である。
2）株式 X と株式 Y で構成された回帰方程式の残差系列は定常である。
上記の 2 つの条件を満たす場合、株式 X と株式 Y は共整合性を持つと言えます。
高い相関性が必ずしも共整合性を意味するわけではありません。2 つの株の相関がほぼ同じであっても、共整合関係の確率には大きな差があります。
(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）：単位根が存在する（非定常）
- 対立仮説（Alternative Hypothesis）：単位根が存在しない（定常）
ADF（y, lags, trend, max_lags, method）関数の説明：
（参考文献『量的投資を Python で』P338）
- パラメータ lags の値は検定モデルを決定します。lags の値は一次差分の遅延次数を示します。
- 検定結果の Test Statistic 値と Critical Values 値を比較し、どの区間に落ちるかを確認します。有意水準（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値が得られるはずです。なぜなら、物理的に可能な限り共整合した2つの系列を人工的に作成したからです。
```python
# 共整合性テストのp値を計算
# 2つの時系列間の比率が平均値の周りで定常であるかどうかを知らせます
score, pvalue, _ = coint(X, Y)
p 値を確認するだけで、p 値が十分に小さい場合、株式 X と株式 Y は共整合関係にあると見なされます。（方法 2 は方法 1 よりも簡単だと思いませんか～～ただし、方法 1 は共整合の概念を理解するのに役立ちます）
4. ペアトレードを行う方法は？#
2 つの共整合した時系列（上記の 600015.XSHG と 601169.XSHG など）は、お互いに近づいたり離れたりするため、スプレッドが高い時と低い時があります。ペアトレードは、1 つの証券を買い、もう 1 つを売ることによって行います。このようにして、両方の証券が一緒に下落したり上昇したりしても、私たちは利益も損失も出さない — 市場中立です。
上記の X と Y に戻ると、Y = ⍺ X + e とし、比率（Y/X）がその平均値⍺の周りで動くと、2 つの比率が平均に戻ることで利益を得ます。これを行うために、X と Y が遠く離れているときを監視します。つまり、⍺が高すぎるか低すぎるときです：
- 比率をロングする：これは、比率⍺が通常より小さく、増加することを期待する場合です。上記の例では、Y を買い、X を売ることで賭けを行います。
- 比率をショートする：これは、比率⍺が大きく、減少することを期待する場合です。上記の例では、Y を売り、X を買うことで賭けを行います。
注意：常に「ヘッジポジション」を持っています：ショートポジションは、売却した証券が価値を失うと利益を得ますし、ロングポジションは、証券が価値を得ると利益を得ますので、全体の市場の動きには影響されません。私たちは、証券 X と Y が相対的に動く場合にのみ利益または損失を得ます。
5. このように振る舞う証券を見つけるためにデータを使用する#
これを行う最良の方法は、共整合している可能性がある証券から始めて、統計的テストを実行することです。すべてのペアに対して統計的テストを実行すると、多重比較バイアスに陥る可能性があります。
多重比較バイアスとは、多くのテストを実行するために、多くのテストを実行する際に有意な p 値を誤って生成する可能性が高まることを指します。ランダムデータに対して 100 のテストを実行すると、5 つの p 値が 0.05 未満になることが予想されます。n の証券を共整合性のために比較する場合、n (n-1)/2 の比較を行い、多くの誤って有意な p 値が見られることが期待されます。これを避けるために、共整合している可能性がある少数のペアを選び、それぞれを個別にテストします。これにより、多重比較バイアスへの曝露が減少します。
それでは、共整合性を示す証券を見つけてみましょう。中国の銀行株のバスケットで作業しましょう。「600000.XSHG」、「600015.XSHG」、「600016.XSHG」、「600036.XSHG」、「601009.XSHG」などです。これらの株は同じセグメントで運営されており、共整合した価格を持つ可能性があります。証券のリストをスキャンし、すべてのペア間の共整合性をテストします。共整合性テストスコア行列、p 値行列、および p 値が 0.05 未満のペアを返します。この方法は多重比較バイアスに陥りやすく、実際には証券は 2 回目の検証ステップを受けるべきです。この例のためにこれを無視しましょう。
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"]
# 各株の共整合性テストのp値を示すヒートマップ
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 スコアとして扱う方がより有用です。z スコアは次のように定義されます：
Z スコア（値）=（値 - 平均）/ 標準偏差
警告
実際には、これはデータにスケールを与えるために通常行われますが、これは基礎となる分布を仮定します。通常は正規分布です。しかし、金融データの多くは正規分布ではなく、統計を生成する際に単に正規性や特定の分布を仮定しないように非常に注意する必要があります。比率の真の分布は非常に太い尾を持ち、極端な値がモデルを混乱させ、大きな損失を引き起こす可能性があります。
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 スコアが - 1.0 未満のときは比率を「ロング」する
- z スコアが 1.0 を超えるときは比率を「ショート」する
- z スコアがゼロに近づくときはポジションを終了する
これは氷山の一角に過ぎず、概念を説明するための非常に単純な例です。
- 実際には、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:]
一般的に、全体のサンプルサイズに対して統計を取ることは悪い場合があります。たとえば、市場が上昇しており、両方の証券も上昇している場合、過去 3 年間の平均価格は今日の代表的なものではないかもしれません。このため、トレーダーは通常、最近のデータのローリングウィンドウに依存する統計を使用します。
比率の値を使用する代わりに、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 スコアを見てみましょう。
# 60日間の標準偏差を取得
std_60 = train.rolling(window=60, center=False).std()
std_60.name = 'std 60d'
# 各日のzスコアを計算
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 スコアは文脈なしではあまり意味がありません。価格と並べてプロットして、どのように見えるかを確認しましょう。
# 比率とzスコアからの売買信号をプロット
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()
これは、実際に取引している株にとって何を意味するのでしょうか？見てみましょう。
# 価格とzスコアからの売買信号をプロット
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()
# 比率を買うときはS1を買い、S2を売る
buyR[buy != 0] = S1[buy != 0]
sellR[buy != 0] = S2[buy != 0]
# 比率を売るときはS1を売り、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 比率（600015.XSHG の株を 1 つ買い、601169.XSHG の株を売る）を購入し、高いときに 1 比率（600015.XSHG の株を売り、601169.XSHG の株を 1 つ買う）を売却し、これらの取引の PnL を計算するシンプルなバックテスターを作成します。
def trade(S1, S2, window1, window2):
# ウィンドウの長さが0の場合、アルゴリズムは意味をなさないので終了
if (window1 == 0) or (window2 == 0):
return 0
# ローリング平均とローリング標準偏差を計算
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
# 取引をシミュレート
# お金もポジションも持たない状態から開始
money = 0
countS1 = 0
countS2 = 0
for i in range(len(ratios)):
# zスコアが1を超えたらショート
if zscore[i] > 1:
money += S1[i] - S2[i] * ratios[i]
countS1 -= 1
countS2 += ratios[i]
# zスコアが-1未満ならロング
elif zscore[i] < -1:
money -= S1[i] - S2[i] * ratios[i]
countS1 += 1
countS2 -= ratios[i]
# zスコアが-0.5と0.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)
戦略は利益が出ているようです！（結果 = 12.716）これで、移動平均のウィンドウを変更したり、売買およびポジション終了の閾値を変更したりして、パフォーマンスの改善を確認できます。また、ロジスティック回帰や SVM などのより洗練されたモデルを試して、1/-1 の予測を行うこともできます。
テストデータでのパフォーマンスを確認してみましょう（結果 = 1.487）。
trade(data['600015.XSHG'].iloc[le:], data['601169.XSHG'].iloc[le:], 5, 60)
7. 過剰適合を避ける#
過剰適合は、取引戦略の最も危険な落とし穴です。私たちのモデルでは、ローリングパラメータ推定を使用し、ウィンドウの長さを最適化したいと考えています。すべての可能な合理的なウィンドウ長を反復し、モデルが最も良く機能する長さを選択できます。以下に、トレーニングデータの PnL に基づいてウィンドウ長をスコアリングし、最適なものを見つけるシンプルなループを書きます。結果は 58 です。
# この戦略を使用して最高のリターンを得るウィンドウ長0-254を見つける
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)
次に、テストデータでのモデルのパフォーマンスを確認し、このウィンドウ長が最適ではないことを見つけます！これは、元の選択が明らかにサンプルデータに過剰適合していたためです。
# このデータセットに基づいて最適なウィンドウ長を見つけ、
# このウィンドウ長を使用したリターンを計算
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])
# このデータセットに基づいて最適なウィンドウ長を見つけ、
# このウィンドウ長を使用したリターンを計算
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 です。
トレーニングデータとテストデータの PnL をウィンドウ長ごとに別々にプロットすると、これが確認できます。
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 のペアトレーディング戦略
- 蔡立耑『量的投資 —Python をツールとして』